José Ferreira, Alan Oliveira, André Souto, José Cecílio LASIGE, Departamento
de Informática, Faculdade de Ciências da Universidade Lisboa, Campo Grande
016, 1749-016 Lisboa
fc53311<EMAIL_ADDRESS>
1 1 March 2023

# Software-based security approach for networked embedded devices

###### Abstract

As the Internet of Things (IoT) continues to expand, data security has become
increasingly important for ensuring privacy and safety, especially given the
sensitive and, sometimes, critical nature of the data handled by IoT devices.
There exist hardware-based trusted execution environments used to protect
data, but they are not compatible with low-cost devices that lack hardware-
assisted security features. The research in this paper presents software-based
protection and encryption mechanisms explicitly designed for embedded devices.
The proposed architecture is designed to work with low-cost, low-end devices
without requiring the usual changes on the underlying hardware. It protects
against hardware attacks and supports runtime updates, enabling devices to
write data in protected memory. The proposed solution is an alternative data
security approach for low-cost IoT devices without compromising performance or
functionality. Our work underscores the importance of developing secure and
cost-effective solutions for protecting data in the context of IoT.

Keywords: IoT Security, Trusted Execution Environment, Code Protection, Memory
Integrity.

##  1 Introduction

As a result of the Internet of Things (IoT) popularization, millions of
embedded devices are being deployed and connected to the worldwide network
[1]. These devices allow IoT to extend the boundaries of the Internet. At the
same time, they connect digital processes to the physical world [2]. Although
the resulting systems create added-value services for the respective
applications, they inherently also bring vulnerabilities [3]. The threats must
be mitigated since these systems might collect and process critical and
sensitive data.

Among the techniques proposed to protect data are the Trusted Execution
Environments (TEEs). TEEs are designed to provide mechanisms to protect
applications (code and critical data), ensuring confidentiality and integrity
[3]. However, most existing TEEs proposals rely on hardware [3], such as
Trusted Platform Modules (TPMs) [4], Intel SGX [5] and ARM TrustZone [6]. Due
to cost and size constraints, those hardware features mainly exist on high-end
platforms and are unavailable on cost-effective and low-end embedded devices
[3].

Alternatively, software-based approaches are being proposed, such as PISTIS
[3], Security MicroVisor ($S\mu V$) [2], SofTEE [7], and Virtual Ghost [8].
These solutions have some advantages when compared with hardware-assisted
ones. One advantage is related to update costs. It is easier and cheaper to
update a software-based TEE. Another advantage is hardware portability, given
that those TEEs do not require specific hardware features (such as ARM
TrustZone or Intel SGX) [7]. Despite these advantages, most of those software-
based secure architectures do not consider hardware-based attacks in their
scheme, i.e., they do not consider scenarios in which the attacker has access
to the device where the application is running and can change its code.

In the literature, it is possible to find solutions hybrid architectures,
i.e., architectures that combine hardware and software modules to achieve the
best of both solutions [9, 10, 11, 12]. Although those solutions offer strong
security guarantees for applications running on low-end devices, to be
implemented in those devices, they require hardware modification [3]. Hardware
modification is impracticable in real-world scenarios, considering that every
device needs the addition of customized hardware [3], and it is difficult for
legacy devices to take advantage of them [7]. Moreover, the hybrid solution
approaches usually consider that hardware-based attacks are out of scope.

Considering the lack of hardware security features in low-cost embedded
systems and the benefits offered by software-based TEEs, this work aims to
develop a software-based security approach for networked embedded devices
(SbS4NED) that provides a set of lightweight mechanisms to protect software
and data integrity (continuously verifying the integrity of memory) and offer
correction in case of unexpected changes. The application code will also be
protected using encryption. This way, it becomes tempered resistant and offers
more reliability to the verification process. Moreover, it will be supported
by lightweight cryptography algorithms presented at the National Institute of
Standards and Technology (NIST) competition [13]. In particular, Xoodyak [14],
one of NIST’s finalists, will be used to encrypt data.

##  2 Thread Model and Assumptions

Defining the assumptions about the attacker’s capabilities and goals, the
system’s components and interactions, and the security goals that need to be
achieved are essential. In this work, we consider adversaries with the
following capabilities:

  * •

The adversaries have access to the device. They may modify the application
code running on the device to read or change the data the application handles.

  * •

The adversaries can sniff the network, modify messages exchanged between
devices, and perform man-in-the-middle attacks.

  * •

Software-based adversaries may be present on the device where the architecture
will be deployed. Their goal may be to change the data available in memory and
consequently control the entities that rely on data accuracy.

We assume that the SbS4NED is correctly installed on the devices by a trusted
party. We also assume it is bug-free, encrypted, and working as expected.
Therefore, the adversary can not surpass the code, and the verification
process carried out by its components. The final assumption is that each
device has mechanisms to compute the encryption key used to protect the local
files where SbS4NED keys are stored.

##  3 Software-based security Architecture

As mentioned in the Introduction, this work aims to build a software-based
tempered-resistant solution that protects the software and data in networked
embedded devices. Driven by this goal, we design the SbS4NED proposed
architecture.

Figure 1: SbS4NED architecture.

Figure 1 shows a high-level description of the proposed architecture, where
the SbS4NED Computing Module (SbS4NED_CM) runs inside the Gateway. It is
responsible for monitoring applications running on the nodes connected to the
Gateway. Each node will have an agent (SbS4NED_Agent) generating the signature
of the application code running on the node and sending it to the SbS4NED_CM
for code integrity check. The SbS4NED_CM and the node code’s application are
encrypted to increase security and to offer more protection to SbS4NED_CM
internal processes. Moreover, the messages exchanged between SbS4NED_CM and
its agent are also encrypted. Next, we describe the SbS4NED components:

  * •

App Manager – It interacts with the applications deployed in the node and aims
to perform application updates and send and receive data from the nodes.

  * •

Key Manager – This component is responsible for managing (i.e., generating and
renewing) the keys used internally by SbS4NED_CM and for external
communication (with a SbS4NED_Agent running on the node). It uses Diffie-
Hellman (DH) key-exchange protocol for external communication to generate or
renew the key.

  * •

Crypto – Provides cryptography services inside the SbS4NED_CM and the
SbS4NED_Agent. It can encrypt, decrypt, and compute the message authentication
code (MAC). In the SbS4NED_CM side, the App Manager can also use this
component to encrypt the app-compiled code before sending it to the node. This
way, secure code update is ensured.

  * •

Integrity Checker – Designed for memory integrity checking. It writes the data
from App Manager in the memory and holds the (randomized) position where it is
written. The Integrity Checker is also responsible for remotely checking the
integrity of the nodes’ code.

  * •

Logger – It is responsible for keeping the log files updated regarding the
memory integrity state, which app the data came from, which nodes are
connected, and any network activity that must be logged to easily detect if an
attacker is trying to join the network or injecting any data on the network.

  * •

App Thread – It is used for executing the application code developed by the
user. It offers an API to interact with the node underlayer software and
hardware layer. All the interaction must be done using the API to ensure, for
instance, that the exchanged data is encrypted.

Figure 2: SbS4NED modules.

The system architecture of SbS4NED_CM and its agents is illustrated in Figure
2. The figure provides an overview of the individual components and their
interconnections within the system.

The Key Manager and Crypto components are the essential modules used in the
SbS4NED_CM and also in its agents. These components are deployed on both sides
of the system, providing the necessary encryption services and ensuring that
data transmitted between the agents and SbS4NED_CM remain secure and
confidential.

The App Manager, Integrity Checker, and Logger are part of SbS4NED_CM. These
must be deployed in the Gateway. The App Manager plays a crucial role in
dealing with the application code deployed in the agents, providing the
necessary services to manage, update, and configure them. On the other hand,
the Integrity Checker ensures that the code and messages transmitted by the
agent are authentic and have not been tampered with.

###  3.1 Data Protection

Memory integrity is nowadays a crucial security concern. The integrity of the
data stored in memory is essential to ensure the system’s proper functioning
and to prevent unauthorized access or manipulation of sensitive information.
When data is written in the memory, two pieces of information are stored: the
value ($v$) that can be accessed by any other external entity, such as an
actuator, and the data integrity ($I$) needed to check the integrity of the
data. Data integrity $I$ is computed using the MAC, and its purpose is to
ensure that the data in memory has not been tampered with or modified. The
formula to compute it is $I=$MAC($v\oplus t$), where $t$ can be a timestamp or
a random number known only by the integrity checker. The position of the value
$I$ in memory is arbitrary, and only the Integrity Checker knows where it is
placed. The choice of MAC over hash functions is because MAC uses a secret key
to generate the authentication code. Assuming that SbS4NED_CM has exclusive
access to this private key, only authorized accesses to SbS4NED_CM can
generate the correct MAC result. In contrast, anyone can compute the hash
value by identifying the function used. A secret key prevents anyone from
computing the hash value and faking the integrity data. This way, SbS4NED_CM
provides a strong level of security against malicious attacks to the memory
and an easy and efficient verification of the memory’s integrity.

###  3.2 Application Protection

Besides confidentiality protection provided by encryption, the application
code is locally stored in the node and decrypted only during execution.
Application code also has integrity protection to ensure it remains unmodified
even when other parties can access the node. To check the code’s integrity, a
challenge-response protocol is used, in which the SbS4NED_CM will send
$enc(t)$, an encrypted challenge to the node, where $t$ can be a timestamp or
a value randomly generated by the SbS4NED_CM. The node has to reply with the
hash of the whole application code ($ac$) XOR-ed with $t$, i.e.,
$enc(hash(ac\oplus t))$. Then, the SbS4NED_CM checks the validity of the
result. If the node fails the validation, SbS4NED_CM could force an update to
restore the node application. If the node gives no response, SbS4NED_CM
assumes that the node is lost or compromised and its data is dropped.

###  3.3 Keys Renewal

A renewal cycle mechanism can extend the system’s lifetime. Therefore, before
any application update, the key used to encrypt the application code and for
communication between the SbS4NED_CM and node is renewed using a DH protocol.
However, the application could take a long time without needing an update.
Therefore, the node is provided with an encryption application code, and the
SbS4NED_CM can initialize the DH key exchange with the node to generate the
new key even when there is no call of the functionally. The generated key will
be used to renew the encrypt of the application code and in further
communication with SbS4NED_CM.

###  3.4 Encryption Algorithms

Since the proposed architecture uses an encryption algorithm and targets low-
end devices, the algorithms that can be used must be lightweight, including
the cryptographic ones. Thus, this architecture will use NIST lightweight
encryption algorithms to protect code and data during message exchange [13].
Although by the time of writing this manuscript, the final stage of the NIST
competition comprises ten finalists, we are only interested in algorithms that
can deal with stream encryption. The main reason is that we want the encrypted
code and messages to have the same size as the original ones. Among the ten
finalists, few support stream encryption natively. For the SbS4NED
architecture, Xoodyak and ISAP schemes with keyed mode association are
considered. Since the SbS4NED architecture is modular, other algorithms may
also be used.

##  4 Proof of concept

To verify and characterize the proposed architecture, we are currently
implementing it. We plan to deploy the architecture in a prototype, enabling
system testing in a distributed environment. The prototype will consist of a
Gateway on which the SbS4NED_CM will be running, with connections to multiple
nodes where the SbS4NED agents and applications will be deployed. By
conducting tests in a distributed environment, we can ensure that the
architecture can effectively handle the communication and data exchange
requirements between the Gateway and the nodes. Additionally, we will be able
to identify potential issues or limitations during deployment, which will help
us refine the architecture further.

During the early stages of designing the proposed architecture, we conducted
experiments using two NIST lightweight encryption algorithms, ISAP and
Xoodyak, to determine which would be more suitable for our research. Our
experiment used a Raspberry Pi Model 3+ platform with a Quad-core @1.4GHz and
1GB LPDDR2 SRAM. We tested both algorithms for their execution time and memory
usage using the same key length of 16 bytes, a nounce of 16 bytes, and file
sizes from 1 to 65 kilobytes (kB) (Figure 3). The tests also included both the
encryption (Figure3(a)) and decryption (Figure 3(b)) processes.

(a) [ISAP] and [Xoodyak] Encryption on Raspberry Pi 3+.

(b) [ISAP] and [Xoodyak] Decryption on Raspberry Pi 3+.

Figure 3: [ISAP] and [Xoodyak] on Raspberry Pi 3+.

Our experiment shows that Xoodyak was more suitable for our research than ISAP
in terms of both average execution time and memory usage. The average
execution time for Xoodyak was consistently lower, especially for larger file
sizes (generally, Xoodyak is 30 to 60 times faster than ISAP with a maximum
standard deviation of 0.004 ms). Moreover, both algorithms required about 370
kb of RAM to perform encryption and decryption tasks. It was noticed that file
size does not affect the algorithm’s memory usage. In summary, our experiments
with ISAP and Xoodyak algorithms, conducted on a Raspberry Pi Model 3+
platform, helped us to determine which algorithm is more suitable for our
research.

##  5 Conclusion

This work proposes a software-based secure execution environment architecture
that is lightweight, requires no hardware modifications and is resistant to
the most common hardware attacks. Currently, this architecture is being
implemented as a proof of concept. After the implementation, tests will be
conducted to analyze its performance and characterize its efficiency.

## Acknowledgments

This work was supported by the LASIGE Research Unit (ref. UIDB/00408/2020 and
ref. UIDP/00408/2020), and by the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 871259 (ADMORPH project).

## References

  * [1] M. Ammar and B. Crispo, “Verify&revive: Secure detection and recovery of compromised low-end embedded devices,” in Annual Computer Security Applications Conference, (New York, NY), p. 717–732, ACM, 2020.
  * [2] M. Ammar, B. Crispo, B. Jacobs, D. Hughes, and W. Daniels, “S$\mu$v—the security microvisor: A formally-verified software-based security architecture for the internet of things,” IEEE Transactions on Dependable and Secure Computing, vol. 16, no. 5, pp. 885–901, 2019.
  * [3] M. Grisafi, M. Ammar, M. Roveri, and B. Crispo, “PISTIS: Trusted computing architecture for low-end embedded systems,” in 31st USENIX Security Symposium, (Boston, MA), pp. 3843–3860, Aug. 2022.
  * [4] T. C. Group, “TCG specification architecture overview.” https://trustedcomputinggroup.org/wp-content/uploads/TCG\\_1\\_4\\_Architecture\\_Overview.pdf, 2007\.  [Online - Accessed on 11-11-2022].
  * [5] V. Costan and S. Devadas, “Intel sgx explained.” Cryptology ePrint Archive, Paper 2016/086, 2016.  https://eprint.iacr.org/2016/086.
  * [6] B. Ngabonziza, D. Martin, A. Bailey, H. Cho, and S. Martin, “Trustzone explained: Architectural features and use cases,” in 2nd IEEE Int. Conf. on Collaboration and Internet Computing, pp. 445–451, 2016.
  * [7] U. Lee and C. Park, “Softee: Software-based trusted execution environment for user applications,” IEEE Access, vol. 8, pp. 121874–121888, 2020.
  * [8] J. Criswell, N. Dautenhahn, and V. Adve, “Virtual ghost: Protecting applications from hostile operating systems,” SIGARCH Comput. Archit. News, vol. 42, p. 81–96, feb 2014.
  * [9] I. Nunes, K. Eldefrawy, N. Rattanavipanon, M. Steiner, and G. Tsudik, “VRASED: A verified Hardware/Software Co-Design for remote attestation,” in 28th USENIX Security Symposium, (Santa Clara, CA), pp. 1429–1446, 2019.
  * [10] K. Eldefrawy, A. Francillon, D. Perito, and G. Tsudik, “Smart: Secure and minimal architecture for (establishing a dynamic) root of trust,” in 19th Annual Network and Distributed System Security Symposium, February 5-8, San Diego, USA (ISOC, ed.), (San Diego), 2012.
  * [11] P. Koeberl, S. Schulz, A.-R. Sadeghi, and V. Varadharajan, “Trustlite: A security architecture for tiny embedded devices,” in Proceedings of the Ninth European Conference on Computer Systems, EuroSys ’14, (New York, NY, USA), ACM, 2014.
  * [12] F. Brasser, B. El Mahjoub, A.-R. Sadeghi, C. Wachsmann, and P. Koeberl, “Tytan: Tiny trust anchor for tiny devices,” in 2015 52nd ACM/EDAC/IEEE Design Automation Conference (DAC), pp. 1–6, 2015.
  * [13] NIST, “Lightweight cryptography.” https://csrc.nist.gov/projects/lightweight-cryptography.  Accessed: 2022-12-04.
  * [14] J. Daemen, S. Hoffert, S. Mella, M. Peeters, G. Van Assche, and R. Van Keer, “Xoodyak, a lightweight cryptographic scheme,” NIST, 05 2021.


# Dynamical fractal and anomalous noise in a clean magnetic crystal

Jonathan N. Hallén  TCM Group, Cavendish Laboratory, University of Cambridge,
Cambridge CB3 0HE, UK  Max Planck Institute for the Physics of Complex
Systems, 01187 Dresden, Germany  Santiago A. Grigera  Instituto de Física de
Líquidos y Sistemas Biológicos, UNLP-CONICET, La Plata, Argentina  D. Alan
Tennant  Department of Physics and Astronomy, University of Tennessee,
Knoxville, TN 37996, USA  Department of Materials Science and Engineering,
University of Tennessee, Knoxville, TN 37996, USA  Claudio Castelnovo  TCM
Group, Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, UK
Roderich Moessner  Max Planck Institute for the Physics of Complex Systems,
01187 Dresden, Germany

(November 1, 2022)

###### Abstract

Fractals – objects with non-integer dimensions – occur in manifold settings
and length scales in nature, ranging from snowflakes and lightning strikes to
natural coastlines. Much effort has been expended to generate fractals for use
in many-body physics. Here, we identify an emergent dynamical fractal in a
disorder-free, stoichiometric three-dimensional magnetic crystal in
thermodynamic equilibrium. The phenomenon is born from constraints on the
dynamics of the magnetic monopole excitations in spin ice, which restrict them
to move on the fractal. This observation explains the anomalous exponent found
in magnetic noise experiments in the spin ice compound Dy2Ti2O7, and it
resolves a long standing puzzle about its rapidly diverging relaxation time.
The capacity of spin ice to exhibit such striking phenomena holds promise of
further surprising discoveries in the cooperative dynamics of even simple
topological many-body systems.

##  I Introduction

The current intense research efforts on the behaviour of topological matter,
besides unearthing many exciting phenomena, attempt to yield an understanding
of these systems on the same level of both generality and detail as is
available for conventional systems [1, 2]. One particular frontier concerns
the dynamical properties, especially of topological systems that host exotic
‘fractionalised’ excitations such as Laughlin quasiparticles with anyonic
statistics in the quantum Hall effect [3], or emergent magnetic monopoles in
the topological spin liquid known as spin ice [4] (see Fig. 1). Indeed, the
dynamical behaviour of the latter has been an enigma since its discovery [5,
6, 7]. Most recently, as well as strikingly, ultrasensitive low-temperature
SQUID experiments have thrown up a new puzzle: the magnetic noise spectral
density exhibits an anomalous power law as a function of frequency, with the
low temperature exponent $\alpha\approx 1.5$ deviating strongly from the well-
known $\alpha=2$ of a paramagnet [8, 9] (see also Ref. [10, 11, 12]). Within
the – generally successful – framework of what we call the ‘standard model’ of
spin ice dynamics [13, 14], this behaviour cannot be accounted for using
broadly accepted model Hamiltonian parameters [9].

Figure 1:  In spin ice, classical Ising-like moments reside on a pyrochlore
lattice of corner-sharing tetrahedra. In each of the exponentially numerous
ground states, two spins point into a tetrahedron, and two out. Violations of
this ice rule take the form of magnetic monopoles (marked by a red sphere).
These monopoles move on the diamond lattice formed by the centres of the
tetrahedra. From a given position, monopoles can hop by flipping any one of
the three majority spins (bonds highlighted in green and grey), but not the
minority spin (unmarked bonds). Here we show parts of a spin ice
configuration, including all tetrahedra that this monopole could reach within
3 hops. Some spins experience a vanishing transverse field [15], and their
dynamics is substantially suppressed. If moves through such spins are
forbidden, the monopole is restricted to move on the paths highlighted in
green – these form the emergent dynamical fractal.

In this work, we provide the missing ingredient: besides the constraints
imposed by the emergent gauge field in spin ice, there is a further dynamical
bottleneck on account of the local (transverse) field distribution, which
suppresses the dynamics of another quarter of the spins [15]. Both
restrictions reflect the random yet correlated orientation of the spins in the
spin ice ground states, and force the monopoles to move on an effectively
disordered, slowly-time-evolving cluster in real space even in the absence of
quenched disorder.

Crucially, this cluster is close to a percolation transition. It therefore
exhibits a well-developed fractal structure on short and intermediate length
scales, which we characterise in detail. We show that it is through hosting
monopole motion that the fractal structure bequeathes anomalous exponents to
the magnetic noise. Our numerical modelling of this process allows us to
quantitatively reproduce the experimental noise curves, with only a single
global fitting parameter for a microscopic timescale. In the process, we also
shed light onto a further, long-standing puzzle in spin ice: the steeper than
expected rise of the macroscopic relaxation time upon cooling. Our theory
explains this phenomenon naturally, in a clean (i.e., stoichiometric and
uniform) system, as a reflection of the sparseness and structure of the
dynamical fractal.

Although fractals have been found or artificially constructed in a range of
systems such as porous materials [16], polymers [17], synthetic atom lattices
[18], molecular systems [19], and the reciprocal space of quasi-crystals [20],
we are not aware of any previous examples of fractals in clean stoichiometric
bulk crystals. Notably, the fractal geometry in spin ice influences the
dynamics in a qualitative and experimentally observable way while leaving no
signatures in the thermodynamics – which is presumably why it has eluded
discovery for so long.

##  II Model

Spin ice is a topological magnet [21, 4] with fractionalised quasiparticles in
the form of mobile magnetic monopoles [22], as illustrated in Fig. 1. The low-
temperature behaviour of spin ice in and out of equilibrium can largely be
recast in terms of the dynamics of a dilute gas of monopoles [13], and how
they interact with the background spin configuration. Specifically, the widely
used ‘standard model’ (SM) of incoherent spin ice dynamics forbids spin flips
which create – rather than hop – monopoles (Fig. 1). For a monopole in a
tetrahedron, this constraint systematically blocks one direction out of four
[14].

This model has successfully described important features of the dynamics –
especially the exponentially divergent relaxation time at low temperatures.
However, it has failed to account for the large energy scale in the leading
exponential growth of this time scale [23, 24], and it has in particular been
challenged by susceptibility [25, 26, 27] and anomalous magnetic noise
experiments [8, 9]. These puzzles and their resolution are discussed below.

Our new ‘beyond the standard’ model (bSM) incorporates the observation that
the internal field distribution on spins across which monopoles hop is
peculiarly bimodal [15]. In particular, one-third of the flippable spins
experience a near-vanishing transverse field (see App. A). We model these
spins as flipping at a lower rate $1/\tau_{\mathrm{slow}}$. By contrast, the
other spins experience a finite transverse field and flip at some reference
rate $1/\tau_{\mathrm{fast}}$. (In the SM, all spins attempt to flip at a
single rate $1/\tau_{0}$.) We take all $\tau$’s to be independent of
temperature. The scale $\tau_{\mathrm{fast}}$ defines our unit of time, and –
as we will see – it is in fact the only fitting parameter in our analysis.

Regarding the interaction parameters, we use a model Hamiltonian, ${\cal
H}_{\mathrm{OP}}$, an extension of the conventional dipolar spin ice
Hamiltonian that was previously obtained from a combined fit to neutron
scattering, magnetic susceptibility, and specific heat measurements [28]. It
comprises long-range dipolar interactions and first, second, and third
nearest-neighbour exchange terms. More details are given in App. A.1 (where a
comparison is drawn to the maximally simple nearest-neighbour spin ice model).

##  III Results

The left panel of Fig. 2 shows the magnetic noise measured using SQUID
magnetometry on a single crystal of DTO [9] in the temperature range
$0.64~{}{\rm K}\leq T\leq 1.04~{}{\rm K}$, high enough to be above the
’freezing’ of spin ice [7], and low enough for monopoles to be sparse, weakly-
interacting quasiparticles. The noise is expressed in terms of the power
spectral density (PSD), $S(\nu)$, defined as the temporal Fourier transform of
the magnetisation, $M(t)$, autocorrelation function: $S(\nu)={\cal
F}\left[\langle M(0)M(t)\rangle\right]$.

Comparison with our simulations using ${\cal H}_{\mathrm{OP}}$ shows that bSM
dynamics reproduces experiments over four orders of magnitude in frequency and
six orders of magnitude in noise power (Fig. 2). The only fitting parameter is
$\tau_{\mathrm{fast}}=85$ $\mu$s (we estimate
$\tau_{\mathrm{slow}}/\tau_{\mathrm{fast}}\gtrsim 10^{3}$ (see App. A). This
is indistinguishable from $\tau_{\mathrm{slow}}=\infty$ in these plots). By
comparison, SM dynamics with fitting parameter $\tau_{0}=200$ $\mu$s is unable
to describe the experimental data at low temperature.

Figure 2:  Left panel: Power spectral density (PSD) of magnetisation
fluctuations extracted from SQUID measurements on Dy2Ti2O7 (filled circles)
[9]. Monte Carlo results for ${\cal H}_{\mathrm{OP}}$ with SM dynamics (grey)
and bSM dynamics (black) are shown with solid lines. The curves have been
shifted vertically (see App. A). An overall time scale factor,
$\tau_{\mathrm{fast}}=85\ \mu$s and $\tau_{0}=200\ \mu$s, was applied for the
bSM and SM results respectively to visually match the experimental data, with
$\tau_{\mathrm{slow}}=\infty$ for simplicity. We note that the deviation at
high frequencies is a known feature of Monte Carlo dynamics, and that other
minor discrepancies are to be expected in a model with a necessarily sharper
parameter distribution than the experimental system (see App. A). Right panel:
Relaxation times $\tau$ were extracted from fits of the PSD curves to the
function $A\left[1+(2\pi f\tau)^{\alpha}\right]^{-1}$ for both experimental
and numerical data.

The extracted relaxation time from the bSM (right panel of Fig. 2) also agrees
remarkably well with experiment [9], especially when contrasted with the SM,
which yields much too short a relaxation time at low temperature. This has
been a puzzle in the community for many years [6, 7, 13, 29, 23, 26, 30, 27,
31]: in a gas of freely moving monopoles, the relaxation time of the
magnetisation scales with their inverse density [13],
$1/\rho\sim\exp(\Delta_{m}/T)$, set by the energy cost $\Delta_{m}$ of an
isolated monopole. Increasing the energy in the Arrhenius law to
$\Delta_{\tau}>\Delta_{m}$ is largely precluded by basic statistical
mechanics, whereas estimates suggest that a $\Delta_{\tau}$ in excess of twice
$\Delta_{m}$ is actually required to fit the experimental growth of the
relaxation time [29]. Previous theories of the steep rise of the relaxation
time upon cooling invoked extrinsic contributions due to open boundary
effects, disorder, and an autonomously temperature-dependent microscopic time
scale [27, 26, 32, 33, 34, 35]; the identification of an intrinsic mechanism
leading to a parametrically faster growth of the relaxation time than $1/\rho$
has been lacking.

To explain the anomalous noise and susceptibility, and the corresponding
strongly diverging relaxation time, we turn to what is perhaps the most
remarkable discovery of this work. An isolated monopole in spin ice with bSM
dynamics has between zero and three choices of sites to move to, with the
statistical average being two. Linking the sites reachable by successive
monopole hops, as in Fig. 3, yields a fractal cluster! Its fractal exponents
are then picked up in the experimental anomalous noise signal, thus altering
the relaxation properties of the system.

In more detail, to understand the fluctuations of the magnetisation, we need
to analyse the statistical properties of the monopole motion, as (i) their
motion proceeds via flipping of spins and (ii) they are the natural (sparse,
weakly interacting) quasiparticles in the regime $T\lesssim 1$ K. The motion
of the monopoles takes place on a ‘dynamical cluster’ in real space, defined
by excluding the spins which are not flippable energetically, or because of a
small local transverse field.

Figure 3:  Top panel: Number of sites that a monopole can visit within a given
number of steps $n$, the so-called chemical distance, with SM (orange circles)
and bSM (blue circles) dynamics. The corresponding results for a walker on the
diamond lattice at various bond filling fractions are shown by green, red, and
cyan crosses. In the bSM case, there is a crossover between fractal $\sim
n^{1.85}$ and standard $n^{3}$ scaling occurring at $n_{cr}\approx 14$, when
the monopole has travelled a distance of the order of the correlation length.
Within the fractal regime, i.e. up to $n_{cr}$, the monopole can access
approximately 130 sites. The error bars, shown only for $n=6$ for simplicity,
reflect the finite width of the distribution and are not due to sampling (the
statistical uncertainty is smaller than the size of the symbols). Bottom
panel: Sites that a specific monopole, starting at the site marked in orange,
can reach with SM dynamics (upper row) and bSM dynamics (lower row) within $n$
steps. The colour of the sites shifts from blue to green with the number of
steps away from the starting site. The fractal cluster is much sparser,
underpinning the rapidly diverging relaxation time (Fig. 2).

For $\tau_{\mathrm{slow}}=\infty$, this defines a percolation problem close to
the critical point, so that the fractal structure is visible on small and
intermediate scales (Fig. 3). In this regime, the number of sites the
structure contains grows anomalously slowly with the chemical distance $n$ –
i.e., the minimum number of steps on the lattice needed to join two sites. We
find that the growth agrees very well with $n^{1.85}$, the predicted exponent
from critical percolation in $d=3$ [36, 37] (see Fig. 3), in contrast to the
conventional $n^{d}$. Note the sparseness of the fractal: up to $n_{cr}=14$,
it contains only about 130 of the 2071 sites available on the diamond lattice.
Around $n_{cr}$, the system crosses over to conventional three-dimensional
behaviour at longer lengthscales. Although our dynamical rules yield a non-
standard correlated percolation problem, this does not impact its critical
behaviour, and the properties of the cluster are in fact accurately reproduced
by a random walker on a random percolation cluster on the diamond lattice at
filling fraction $p=0.43$, only slightly above its critical value
$p_{c}\approx 0.39$ [38] (see App. B).

Our discussion so far has largely been based on the behaviour of a single
monopole. Beyond this, the motion of other monopoles can change the local spin
and transverse field configurations, thereby endowing the percolation cluster
with a slow dynamics of its own (see App. B). A higher monopole density
(alongside a finite $\tau_{\mathrm{slow}}$) therefore diminishes the crossover
value $n_{cr}$.

It is the crossover in real space that terminates the anomalous regime in the
PSD towards low frequencies in the time domain. At yet lower frequencies, a
plateau appears on account of the trajectories of different monopoles
overlapping. From the theory of random walks on percolation clusters above the
percolation threshold [39, 40, 41], the PSD of a monopole should show
anomalous decay $S\sim\nu^{-(1+\sigma)}$ at high frequency and conventional
decay $S\sim\nu^{-2}$ at low frequencies. Exact enumeration in the random
percolation problem was previously used to obtain $\sigma=0.50\pm 0.01$ in
$d=3$ [42]. This largely explains the experimentally observed anomalous power
law which hovers near $\nu^{-1.5}$ in the regime under consideration [9].

##  IV  Discussion

As is so often the case with discoveries in condensed matter physics, their
beauty derives not only from the notability of the new phenomena in itself,
but also from the often unexpected ingredients of the explanation, and the way
they snugly fit together.

From the point of view of percolation, we have discovered a purely
dynamically-generated fractal object in a uniform, stoichiometric, disorder-
free bulk crystal. Its existence is predicated on the emergence of point-like
mobile objects – the magnetic monopoles – in a three-dimensional topological
spin liquid which remains fluctuating down to low temperatures [22]. The
monopole motion is subject in turn to a twofold set of constraints. One,
imposed by the emergent gauge field represented by the spin background,
reflects the large-scale topological nature of the spin ice state. The other
is purely microscopic in origin, resulting from an interplay of the geometric
spin arrangement and the short-range statistical spin correlations imposed by
the ice rules. It is the combination of these constraints that makes the
monopoles move on a fractal structure. This in turn is the origin of the
previously puzzling anomalous magnetic noise and rapidly diverging relaxation
time. Conversely, this highly non-trivial phenomenology provides a rare
validation of a concrete model of the microscopic dynamics, in particular the
bimodal distribution of local transverse fields.

Most remarkably, this is all accessible by probing magnetisation response and
fluctuations, specifically using AC-susceptibility or noise measurements,
whose frequency dependence reflects spatial information, linked by the
(sub-)diffusive motion of the monopoles. The emergence of a dynamical fractal
in a disorder-free crystal presents a novel mechanism for the existence of
anomalous noise, a subject which has been studied in many other contexts,
within materials science [43, 44, 45, 46] and elsewhere [47, 48, 49, 50]. Our
work provides the first bulk magnetic measurement of a fractal.

Some of the above ingredients can be manipulated in a controlled manner, all
of which provide promising avenues for future experiments. Trivially,
different spin ice compounds have different interaction parameters, resulting
in broadly different transverse field distributions. Moreover, the microscopic
dynamics can also be altered by applying uniaxial strain [51] or via the
crystal field scheme, by switching from (Kramers) Dy- to, say, (non-Kramers)
Ho-based spin ice. For example, HTO is predicted to have a much smaller ratio
$\tau_{\mathrm{slow}}/\tau_{\mathrm{fast}}$ [15].

Note that our effective percolation problem has turned out to be at a sweet
spot near criticality – reducing constraints would eliminate the anomalous
nature of the signal, whereas increasing them would likely eliminate
equilibration. Indeed, it will be interesting to see if doing so (e.g., by the
controlled introduction of some form of quenched or dynamical disorder) may
shed light on the ‘glassy’ physics below $650$ mK [52].

Needless to say, the advent of highly-tunable NISQ platforms opens up an
entirely new set of directions, especially in two dimensions [53, 54], where
anomalous noise has already been seen in the classical nanomagnetic artificial
spin ice [55]. This includes questions related to quantum diffusion of
monopoles and the role of increasingly coherent many-body quantum dynamics.
The latter could alter the noise in a characteristic way (see App. B) and
therefore be used as evidence for the presence of so-called ring exchange
processes in candidate quantum spin ice compounds.

In closing, we note that it was thanks to the wealth of information
accumulated over 25 years of work by the magnetism community [52, 4] that we
were able to analyse the cooperative dynamics of spin ice on the level of
detail required for discovering the dynamical fractal. This underlines the
value of well-characterised model systems for driving the discovery of
striking yet subtle phenomena.

## Acknowledgments

We express our gratitude to Cristiano Nisoli, Peter Schiffer, Shivaji L.
Sondhi and Shu Zhang for useful discussions. This work was supported in part
by the Engineering and Physical Sciences (EPSRC) grants No. EP/P034616/1 and
EP/T028580/1 (CC and JNH), the Deutsche Forschungsgemeinschaft under grants
SFB 1143 (project-id 247310070) (RM and JNH), and the cluster of excellence
ct.qmat (EXC 2147, project-id 390858490) (RM and JNH). The work of DAT was
supported by the U.S. Department of Energy, Office of Science, Office of Basic
Energy Sciences, under Award Number DE-SC0022311. Part of this work was
carried out within the framework of a Max Planck independent research group on
strongly correlated systems.

## References

  * Chaikin and Lubensky [1995] P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Moessner and Moore [2021] R. Moessner and J. E. Moore,  _Topological Phases of Matter_ (Cambridge University Press, 2021).
  * Stern [2008] A. Stern, Annals of Physics 323, 204 (2008).
  * Udagawa and Jaubert [2021] M. Udagawa and L. Jaubert,  _Spin Ice_ (Springer, 2021).
  * Snyder  _et al._ [2001] J. Snyder, J. Slusky, R. Cava, and P. Schiffer, Nature 413, 48 (2001).
  * Matsuhira  _et al._ [2001] K. Matsuhira, Y. Hinatsu, and T. Sakakibara, Journal of Physics: Condensed Matter 13, L737 (2001).
  * Snyder  _et al._ [2004] J. Snyder, B. Ueland, J. Slusky, H. Karunadasa, R. Cava, and P. Schiffer, Physical Review B 69, 064414 (2004).
  * Dusad  _et al._ [2019] R. Dusad, F. K. Kirschner, J. C. Hoke, B. R. Roberts, A. Eyal, F. Flicker, G. M. Luke, S. J. Blundell, and J. S. Davis, Nature 571, 234 (2019).
  * Samarakoon  _et al._ [2022] A. M. Samarakoon, S. A. Grigera, D. A. Tennant, A. Kirste, B. Klemke, P. Strehlow, M. Meissner, J. N. Hallén, L. Jaubert, C. Castelnovo, and R. Moessner, Proceedings of the National Academy of Sciences 119, e2117453119 (2022), https://www.pnas.org/doi/pdf/10.1073/pnas.2117453119 .
  * Kirschner  _et al._ [2018] F. K. K. Kirschner, F. Flicker, A. Yacoby, N. Y. Yao, and S. J. Blundell, Phys. Rev. B 97, 140402 (2018).
  * Watson  _et al._ [2019] C. A. Watson, I. Sochnikov, J. R. Kirtley, R. J. Cava, and K. A. Moler, arXiv preprint arXiv:1903.11465 (2019).
  * Nisoli [2021] C. Nisoli, Europhysics Letters 135, 57002 (2021).
  * Ryzhkin [2005] I. Ryzhkin, Journal of Experimental and Theoretical Physics 101, 481 (2005).
  * Jaubert and Holdsworth [2009] L. D. Jaubert and P. C. Holdsworth, Nature Physics 5, 258 (2009).
  * Tomasello  _et al._ [2019] B. Tomasello, C. Castelnovo, R. Moessner, and J. Quintanilla, Physical review letters 123, 067204 (2019).
  * Yu [2008] B. Yu, Applied Mechanics Reviews 61 (2008).
  * Newkome  _et al._ [2006] G. R. Newkome, P. Wang, C. N. Moorefield, T. J. Cho, P. P. Mohapatra, S. Li, S.-H. Hwang, O. Lukoyanova, L. Echegoyen, J. A. Palagallo, V. Iancu, and S.-W. Hla, Science 312, 1782 (2006), https://www.science.org/doi/pdf/10.1126/science.1125894 .
  * Kempkes  _et al._ [2019] S. N. Kempkes, M. R. Slot, S. E. Freeney, S. J. Zevenhuizen, D. Vanmaekelbergh, I. Swart, and C. M. Smith, Nature physics 15, 127 (2019).
  * Rothemund  _et al._ [2004] P. W. K. Rothemund, N. Papadakis, E. Winfree, and A. Condon, PLoS biology 2, e424 (2004).
  * Viebahn  _et al._ [2019] K. Viebahn, M. Sbroscia, E. Carter, J.-C. Yu, and U. Schneider, Physical Review Letters 122, 110404 (2019).
  * Castelnovo  _et al._ [2012] C. Castelnovo, R. Moessner, and S. Sondhi, Annual Review of Condensed Matter Physics 3, 35 (2012).
  * Castelnovo  _et al._ [2008] C. Castelnovo, R. Moessner, and S. L. Sondhi, Nature 451, 42 (2008).
  * Matsuhira  _et al._ [2011] K. Matsuhira, C. Paulsen, E. Lhotel, C. Sekine, Z. Hiroi, and S. Takagi, Journal of the Physical Society of Japan 80, 123711 (2011).
  * Pomaranski  _et al._ [2013] D. Pomaranski, L. Yaraskavitch, S. Meng, K. Ross, H. Noad, H. Dabkowska, B. Gaulin, and J. Kycia, Nature Physics 9, 353 (2013).
  * Quilliam  _et al._ [2011] J. Quilliam, L. Yaraskavitch, H. Dabkowska, B. Gaulin, and J. Kycia, Physical Review B 83, 094424 (2011).
  * Yaraskavitch  _et al._ [2012] L. Yaraskavitch, H. Revell, S. Meng, K. Ross, H. Noad, H. Dabkowska, B. Gaulin, and J. Kycia, Physical Review B 85, 020410 (2012).
  * Revell  _et al._ [2013] H. Revell, L. Yaraskavitch, J. Mason, K. Ross, H. Noad, H. Dabkowska, B. Gaulin, P. Henelius, and J. Kycia, Nature Physics 9, 34 (2013).
  * Samarakoon  _et al._ [2020] A. M. Samarakoon, K. Barros, Y. W. Li, M. Eisenbach, Q. Zhang, F. Ye, V. Sharma, Z. Dun, H. Zhou, S. A. Grigera, _et al._ , Nature communications 11, 1 (2020).
  * Jaubert and Holdsworth [2011] L. D. C. Jaubert and P. C. W. Holdsworth, Journal of Physics: Condensed Matter 23, 164222 (2011).
  * Bramwell [2012] S. T. Bramwell, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, 5738 (2012).
  * Raban  _et al._ [2022] V. Raban, L. Berthier, and P. C. Holdsworth, Physical Review B 105, 134431 (2022).
  * Sala  _et al._ [2014] G. Sala, M. Gutmann, D. Prabhakaran, D. Pomaranski, C. Mitchelitis, J. Kycia, D. Porter, C. Castelnovo, and J. Goff, Nature Materials 13, 488 (2014).
  * Paulsen  _et al._ [2014] C. Paulsen, M. J. Jackson, E. Lhotel, B. Canals, D. Prabhakaran, K. Matsuhira, S. Giblin, and S. Bramwell, Nature Physics 10, 135 (2014).
  * Paulsen  _et al._ [2016] C. Paulsen, S. Giblin, E. Lhotel, D. Prabhakaran, G. Balakrishnan, K. Matsuhira, and S. Bramwell, Nature Physics 12, 661 (2016).
  * Edberg  _et al._ [2022] R. Edberg, I. M. Fjellvåg, L. Hoffman, P. P. Deen, K. Lefmann, P. Henelius, and A. Rydh, arXiv preprint arXiv:2202.02801 (2022).
  * Havlin and Nossal [1984] S. Havlin and R. Nossal, Journal of Physics A: Mathematical and General 17, L427 (1984).
  * Stauffer and Aharony [2018] D. Stauffer and A. Aharony,  _Introduction to percolation theory_ (Taylor & Francis, 2018).
  * Xu  _et al._ [2014] X. Xu, J. Wang, J.-P. Lv, and Y. Deng, Frontiers of Physics 9, 113 (2014).
  * de Gennes [1976] P. G. de Gennes, La recherche 7, 919 (1976).
  * Havlin and Ben-Avraham [1987] S. Havlin and D. Ben-Avraham, Advances in Physics 36, 695 (1987), https://doi.org/10.1080/00018738700101072 .
  * Hughes [2021] B. D. Hughes, Complex Media and Percolation Theory , 191 (2021).
  * Eduardo Roman [1990] H. Eduardo Roman, Journal of statistical physics 58, 375 (1990).
  * Aliev  _et al._ [2007] F. Aliev, R. Guerrero, D. Herranz, R. Villar, F. Greullet, C. Tiusan, and M. Hehn, Applied physics letters 91, 232504 (2007).
  * Balandin [2013] A. A. Balandin, Nature nanotechnology 8, 549 (2013).
  * Najafi  _et al._ [2021] M. Najafi, S. Tizdast, Z. Moghaddam, and M. Samadpour, Physica Scripta 96, 125259 (2021).
  * Costanzi and Dahlberg [2017] B. N. Costanzi and E. D. Dahlberg, Physical review letters 119, 097201 (2017).
  * Voss and Clarke [1975] R. F. Voss and J. Clarke, Nature 258, 317 (1975).
  * Press [1978] W. H. Press, Comments on Astrophysics 7, 103 (1978).
  * Varotsos  _et al._ [2013] C. A. Varotsos, I. Melnikova, M. N. Efstathiou, and C. Tzanis, Theoretical and applied climatology 111, 641 (2013).
  * Vernotte and Lantz [2015] F. Vernotte and E. Lantz, Metrologia 52, 222 (2015).
  * Pili  _et al._ [2021] L. Pili, A. Steppke, M. Barber, F. Jerzembeck, C. Hicks, P. Guruciaga, D. Prabhakaran, R. Moessner, A. Mackenzie, S. Grigera, _et al._ , arXiv preprint arXiv:2105.07977 (2021).
  * Bramwell and Gingras [2001] S. T. Bramwell and M. J. Gingras, Science 294, 1495 (2001).
  * King  _et al._ [2021] A. D. King, C. Nisoli, E. D. Dahl, G. Poulin-Lamarre, and A. Lopez-Bezanilla, Science 373, 576 (2021), https://www.science.org/doi/pdf/10.1126/science.abe2824 .
  * Stern  _et al._ [2021] M. Stern, C. Castelnovo, R. Moessner, V. Oganesyan, and S. Gopalakrishnan, Phys. Rev. B 104, 115114 (2021).
  * Goryca  _et al._ [2021] M. Goryca, X. Zhang, J. Li, A. L. Balk, J. D. Watts, C. Leighton, C. Nisoli, P. Schiffer, and S. A. Crooker, Phys. Rev. X 11, 011042 (2021).
  * Even [2011] S. Even,  _Graph algorithms_ (Cambridge University Press, 2011).
  * Prakash  _et al._ [1992] S. Prakash, S. Havlin, M. Schwartz, and H. E. Stanley, Physical Review A 46, R1724 (1992).

##  Appendix A Methods

###  A.1 Model Hamiltonians

We use two different Hamiltonians for the spin-spin interactions. The first
one was obtained as a combined fit to neutron scattering, magnetic
susceptibility, and specific heat measurements on Dy2Ti2O7 [28]. It consists
of long-range dipolar interactions and short-range exchange between first,
second, and third-neighbour spins:

$\displaystyle{\cal H}_{\rm{OP}}\ \\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\
Da^{3}\sum_{i<j}\left[\frac{\vec{S}_{i}\cdot\vec{S}_{j}}{r_{ij}^{3}}-\frac{3\left(\vec{S}_{i}\cdot\vec{r}_{ij}\right)\left(\vec{S}_{j}\cdot\vec{r}_{ij}\right)}{r_{ij}^{5}}\right]$
(1) $\displaystyle+$ $\displaystyle\\!\\!\\!\
J_{1}\sum_{\left<i,j\right>}\vec{S}_{i}\cdot\vec{S}_{j}+J_{2}\sum_{\left<i,j\right>_{2}}\vec{S}_{i}\cdot\vec{S}_{j}$
$\displaystyle+$
$\displaystyle\\!\\!\\!J_{3}\sum_{\left<i,j\right>_{3}}\vec{S}_{i}\cdot\vec{S}_{j}+J_{3}^{\prime}\sum_{\left<i,j\right>_{3^{\prime}}}\vec{S}_{i}\cdot\vec{S}_{j}\,.$

The dipolar interaction strength is $D=1.3224\ \mathrm{K}/a^{3}$, and the
exchange strengths are $J_{1}=3.41$ K, $J_{2}=0.0$ K, $J_{3}=-0.00466$ K, and
$J_{3}^{\prime}=0.0439$ K [9].

For reference, we contrast this with the simplest possible model: the nearest-
neighbour Hamiltonian,

${\cal
H}_{\rm{NN}}=-J_{\rm{eff}}\sum_{\left<i,j\right>}\vec{S}_{i}\cdot\vec{S}_{j}\,$
(2)

that has an effective exchange parameter $J_{\mathrm{eff}}$. In our work
$J_{\mathrm{eff}}$ was set to $5.7$ K, chosen so that the monopole number
densities of ${\cal H}_{\mathrm{OP}}$ and ${\cal H}_{\mathrm{NN}}$
approximately match at $T=0.75$ K.

In regards to the phenomena of interest in our work – namely the magnetic
noise, the relaxation time scale, and the emergent dynamical fractal– we find
that ${\cal H}_{\rm{OP}}$ and ${\cal H}_{\rm{NN}}$ give qualitatively similar
results. For this reason, part of the analysis presented in the appendices
uses ${\cal H}_{\rm{NN}}$ for computational convenience and to access larger
system sizes.

###  A.2 Transverse fields and spin flipping rates

In Ref. [15] it was observed that spin ice Hamiltonians generally lead to a
bimodal distribution of effective transverse fields at the site of a spin
whose flip leads to the hop of a monopole (see Fig. 4). As explained in detail
there, this arises because the six spins neighbouring a given site provide
environments with two different symmetries. In the high-symmetry case (right
panel of Fig. 4), which occurs in a third of the instances, the site is an
inversion centre, leading to a vanishing net field due to the six neighbours.
In the low-symmetry case, by contrast, it is only the longitudinal field
component which vanishes, while the transverse one remains finite. Farther
range interactions, of dipolar or exchange origin, are found to only introduce
minor corrections. (Specifically to ${\cal H}_{\mathrm{OP}}$, the second
neighbour exchange terms vanish, and the third neighbour ones are not expected
to have a transverse component.)

Figure 4: Two tetrahedra in spin ice with a single monopole. The easy axes of
spins $a_{1}$ and $a_{2}$, $b_{1}$ and $b_{2}$, and $c_{1}$ and $c_{2}$ are
parallel to one another. The longitudinal field on the central spin from its
six nearest-neighbours vanishes in both cases and the central spin can
therefore flip without energy penalty, enabling the monopole to hop. In the
left scenario, there is a non-zero effective transverse field from the
nearest-neighbour spins, corresponding to $\tau_{\rm fast}$. In the
configuration on the right, the effective transverse field vanishes and there
are no matrix elements that enable the spin to flip [15], corresponding to
$\tau_{\rm slow}$.

In the standard model (SM) of spin ice dynamics, this difference is not
included and all spins are assumed to flip stochastically with an average
characteristic rate $1/\tau_{0}$. In this work we consider instead ‘beyond the
standard’ model (bSM) dynamics, which attempts to flip spins at a slow rate
$1/\tau_{\mathrm{slow}}$ if the local field has a vanishing transverse
component, and at a fast rate $1/\tau_{\mathrm{fast}}$ otherwise. In all
cases, the time scales are assumed to be independent of temperature.

###  A.3 Methods – spin ice

All numerical results were obtained using Monte Carlo simulations with the
Metropolis algorithm, or variations thereof as explained below, with periodic
boundary conditions.

We denote the total number of spins in the system by $N$ (note that the number
of spins on the pyrochlore lattice equals the number of bonds on the diamond
lattice). A $16$ site cubic unit cell was used, and we considered systems of
linear size $10$ ($N=16000$) for the thermal noise and diffusion measurements,
and of linear size $16$ ($N=65536$) for the characterisation of the
percolation cluster and measurements that involve a single mobile monopole
(see below).

Spin ice configurations at different temperatures were prepared by simulated
annealing from a high-temperature random configuration. The Ewald summation
technique was used to include long-ranged dipolar interactions. Worm updates
were used to ensure equilibration at the lowest temperatures, and when
preparing the system for single mobile monopole simulations (more on this
later). Worm updates consist of finding a chain of spins aligned head to tail,
forming either a closed loop or terminating at a monopole at one end, and
attempting to flip the corresponding spins simultaneously. The update is
accepted with the customary Metropolis probability
$P=\min{\left[1,\exp\left(\Delta E/T\right)\right]}$, where $\Delta E$ is the
overall change in energy. If the update is accepted, it neither creates nor
annihilates any monopoles.

The dynamics of thermal systems was simulated using single spin flip updates.
A random spin is selected and the energy cost $\Delta E$ associated with
flipping this spin is calculated. In the case of bSM dynamics, the move is
accepted with probability $P=G\,\min{\left[1,\exp\left(\Delta
E/T\right)\right]}$, where $G=1$ if the spin flips with time scale
$\tau_{\mathrm{fast}}$, and $G=\tau_{\mathrm{fast}}/\tau_{\mathrm{slow}}$ if
the spin flips with time scale $\tau_{\mathrm{slow}}$. Slow and fast time
scales were determined based on the configuration of the six nearest-
neighbouring spins, according to the discussion above – see also Ref. [15].

A Monte Carlo sweep consists of $N$ attempted updates, and one sweep
corresponds to the time scale $\tau_{0}$ in SM dynamics, and to the time scale
$\tau_{\mathrm{fast}}$ in bSM dynamics. In the monopole picture this allows
for monopole creation, annihilation, and motion. Double monopoles are allowed,
although they quickly become energetically excluded at low temperatures.

In the thermal simulations, we are also able to follow specific monopoles from
creation to annihilation. When doing this, we first equilibrated the system
and then tracked the monopoles created afterwards, up to the point when they
are annihilated. From the measured trajectories, it is possible to extract
properties such as the monopole lifetime or radius of gyration distributions.
One can also compute the PSD of each individual monopole trajectory; however,
as this requires more intensive simulations than to compute the PSD of the
magnetisation of the entire system, we limited ourselves to consider only the
higher portion of the frequency spectrum, $\nu>10^{-3}/\tau_{\rm fast}$.

####  A.3.1 Single monopole simulations

The dynamical magnetic properties of spin ice, in the temperature regime under
consideration, are best understood in terms of the statistical properties of
monopole motion. Monopole moves occur through spin flips, and magnetisation
fluctuations are thus directly linked to monopole displacement.

Single mobile monopole simulations are a good approximation of ${\cal
H}_{\mathrm{NN}}$ at low temperatures (low monopole density), when one
neglects monopole creation and annihilation events. The simulations were
initiated from a random spin ice state with zero monopoles. A pair of
monopoles was then generated, one of which was kept at a fixed location
throughout the simulation (annihilation events between the pair were
thereafter forbidden). The dynamics is generated from updates consisting of
randomly selecting one of the four spins neighbouring the mobile monopole, and
flipping this if the spin (a) is a majority spin and (b) has a non-zero
transverse field (for bSM dynamics only). Four such attempted updates
correspond to one Monte Carlo sweep in the thermal simulations, and we
therefore equate this to $\tau_{0}$ (SM) or $\tau_{\mathrm{fast}}$ (bSM). This
assumes that the monopole can move at no energy cost, as is the case in the
nearest-neighbour model. Before any measurements were taken, the mobile
monopole was allowed to perform $6.6\times 10^{7}$ attempted moves to ensure
that the positions of the mobile and static monopoles had become sufficiently
uncorrelated.

To control that the single monopole simulations do indeed reproduce the
magnetic noise of the thermal system at low temperatures, we compared the PSDs
obtained from the thermal simulations, single mobile monopole simulations, and
individual monopole trajectories in the thermal simulations (Fig. 5). After
rescaling them to match in magnitude, the three cases agree well, indicating
that the low frequency noise can indeed be understood from the motion of
individual monopoles.

Figure 5:  Comparison of the PSD computed using Monte Carlo simulations of
${\cal H}_{\mathrm{NN}}$ at $T=0.6$ K (blue), using single mobile monopole
simulations (red), and using the trajectories of individual monopoles in the
Monte Carlo simulations (black, obtained only for frequencies $\nu\gtrsim
10^{-3}/\tau_{\rm fast}$). The red and black lines have been rescaled so that
the three cases match in magnitude.

###  A.4 Methods – bond percolation

Uncorrelated bond percolation configurations at different filling fractions
were generated by randomly placing bonds on the diamond lattice. As a way to
further mimic spin ice, a specific type of correlation was also considered,
where the bonds are placed subject to the constraint that no vertex is
connected to more than three bonds. This mimics the constraint that every site
containing a single monopole in spin ice has one minority spin that the
monopole cannot move through. The constraint was implemented by first creating
an ordered bond configuration with fixed connectivity, ensuring that every
vertex has exactly three bonds. The configuration was then scrambled by
randomly applying dimer loop updates, and further bonds were subsequently
removed at random to reach the required filling fraction. We find that this
constraint increases the percolation threshold to $p\approx 0.41$, above which
a structured percolation cluster is created. (Here the word _structured_ is
used to refer to the extra correlations introduced by the constraint.)

The random walk on the uncorrelated percolation cluster (PC) or the structured
percolation cluster (SPC) works similarly to the single mobile monopole
simulations. A walker is placed at a random vertex on the percolation cluster.
The updates consist of randomly choosing one of the four bonds at that vertex
and moving the walker along it, if the bond exists. This is the diffusion
model commonly referred to as “the blind ant in the labyrinth” [39, 40]. To
enable comparison with bSM spin ice, we consider four such attempted updates
to be equivalent to $\tau_{\mathrm{fast}}$.

###  A.5 Connected clusters

An adapted depth-first search [56] was used to identify clusters of connected
sites. Following the standard algorithm, a monopole or walker is moved from an
initial reference site to all sites it can reach within a specific chemical
distance (the minimum number of steps on the lattice) from the reference site.
In the case of monopoles, the spin configuration is updated as the search is
performed, and backtracking from dead-ends or sites at the largest chemical
distance must be performed explicitly to maintain the correct spin
configuration.

###  A.6 Power spectral density

PSDs were calculated using Welch’s method. The high frequency aliasing tail
(namely, beyond the Nyquist frequency) has been cut off in the examples shown.
Relaxation times $\tau$ were extracted by fitting the low frequency knee with
the function

$f=A\frac{\tau}{1+\left(\tau\nu\right)^{\alpha}}\,.$ (3)

###  A.7 Experimental and numerical units

The experimental magnetic noise data shown in Fig. 2 of the main text were
recorded in units of $\Phi_{0}^{2}/{\rm Hz}$ [9] ($\Phi_{0}=h/2e\simeq
2.068\cdot 10^{-15}$ Vs being the magnetic flux quantum).

In our simulations we measure the magnetisation in units of the single spin
magnetic moment strength $\mu$, which for Dy2Ti2O7 is approximately
$10\mu_{B}$ (Bohr magnetons) [4]. The natural unit of time is
$\tau_{\mathrm{fast}}$ in simulations with bSM dynamics and the $\tau_{0}$ in
simulations with the SM dynamics. When comparing the models below we use
$\tau_{0}=\tau_{\mathrm{fast}}$. The PSD computed from simulations is measured
in units of $\mu^{2}\tau_{\mathrm{fast}}$.

The numerical results in Fig. 2A in the main text have been scaled to match
the experimental ones; unfortunately the latter depend on details such as the
size of the squid and the portion of sample contributing to the magnetic flux
signal, which carry too much uncertainty to be able to derive the appropriate
rescaling factor from first principles. Based on visual inspection, factors of
850 and 795 respectively were used to rescale the SM and bSM numerical PSD
curves. The same factor was used for all displayed temperatures. For visual
clarity the PSD curves were shifted vertically by multiplicative factor
$2^{0}$, $2^{4}$ and $2^{8}$ for temperature 0.65, 0.84, and 1.04 K
respectively.

##  Appendix B Supplementary Text

###  B.1 Comparison between simulations and experiments

The agreement between bSM simulations and experimental results (displayed in
Fig. 2 of the main text) is quite remarkable, considering our model only has
one global fitting parameter. There are however some discrepancies which we
comment on here.

Deviations at high frequencies are a known feature of simple Monte Carlo
dynamics (see for example [9]). Monte Carlo dynamics relies on the
approximations that a single spin attempts to flip at uncorrelated times and
that it on average attempts one flip in time $\tau_{0}$ (or
$\tau_{\mathrm{fast}}$ in the bSM case); it models a Poissonian process. As a
result of this, the PSD for frequencies $\nu>1/\tau_{0}$ always decays as
$\nu^{-2}$. At the same time, the experimental signal at the highest
frequencies is very small (for $T=0.65$ K, it is about eight orders of
magnitude below the low-frequency signal), and the relative upturn in the
experimental curve there is likely also a signal-to-noise issue.

The theoretical model naturally has a sharper parameter distribution than the
experimental system, and this leads to a signal with sharper crossovers
between different regimes, which would account for why the numerical curve
overshoots the experimental one somewhat at low frequencies, and then drops
slightly below it at intermediate frequencies. Future work may help fine-tune
the description and reduce these discrepancies; however, this will require a
separate microscopic determination of any additional parameters (e.g., of the
value of spin flip time scales like $\tau_{\mathrm{fast}}$), lest it be a mere
fit improvement by increasing the number of adjustable parameters.

###  B.2  Constrained spin ice as a correlated percolation problem

The fractal nature of the clusters the monopoles move on arises from the
proximity to the percolation threshold when one maps blocked spins onto
missing bonds in a bond percolation problem. In the standard language of
percolation theory [37], the percolation cluster first appears at the
percolation threshold $p_{c}$. For bond percolation on the diamond lattice,
this has been numerically estimated at $p_{c}\approx 0.39$ [38]. At this
point, the percolation cluster is self-similar on all length scales. If one
increases the filling fraction beyond the threshold value, the percolation
cluster remains fractal only on length scales shorter than the correlation
length $\xi$, which decreases as the filling fraction increases. A random
walker moving on the percolation cluster will move sub-diffusively on short
length scales, with a mean-squared displacement $\langle R^{2}(t)\rangle\sim
t^{\sigma}$ and $\sigma<1$. For three-dimensional lattices, $\sigma=0.50\pm
0.01$ has been computed using exact enumeration [42]. (Note that our $\sigma$
is directly related to $d_{w}$ in the notation of Ref. [42], with
$\sigma=2/d_{w}$ and $d_{w}$ defined by $\sqrt{\langle R^{2}(t)\rangle}\sim
t^{1/d_{w}}$.). On scales greater than $\xi$ the motion of the walker becomes
diffusive, $\langle R^{2}(t)\rangle\sim t$. This behaviour is observed for a
monopole in bSM spin ice (Fig. 6). After a time of approximately $10^{3}\
\tau_{\mathrm{fast}}$, the crossover from sub-diffusive to diffusive behaviour
occurs. The corresponding length scale of the crossover is approximately $8\
a_{d}$, where $a_{d}$ is the nearest-neighbour distance on the diamond
lattice.

Figure 6:  The root-mean-squared displacement of a monopole in bSM spin ice,
measured in units of the diamond lattice nearest-neighbour distance $a_{d}$.
From the crossover behaviour we can extract a typical time and length scale
for the fractal behaviour, and we find that they are approximately
$1200\tau_{\mathrm{fast}}$ and $8a_{d}$ respectively.

In Fig. 7, the PSD of a single mobile monopole (see Methods) with bSM dynamics
in ${\cal H}_{\mathrm{NN}}$ is shown together with the PSD of a random walker
moving on the percolation cluster on the diamond lattice at various filling
fractions. At high frequencies, the PSD of the monopole matches that of a
random walker on the percolation cluster at $p_{c}$, and decays as
$\nu^{-1.5}$. On longer time scales, the monopole moves on length scales
greater than $\xi$, at which point the lattice is no longer fractal and normal
diffusive behaviour is recovered. This leads to the $\nu^{-2}$ decay observed
at low frequencies. From the cross-over between the two decay regimes, we
again conclude that the fractal motion of the monopole occurs on time scales
smaller than approximately $10^{3}\ \tau_{\mathrm{fast}}$.

Figure 7:  PSDs computed from Monte Carlo simulations of a random walk on the
uncorrelated percolation cluster (PC) and on the structured percolation
cluster (SPC) on the diamond lattice. The cyan line shows that the PSD indeed
goes as $\nu^{-1.5}$ when $p=p_{c}$. The PSD computed for the position of a
single mobile monopole with bSM nearest-neighbour spin ice (blue line) is
reproduced by a random walker at filling fractions $p=0.43$ and $p=0.46$ for
the PC and SPC respectively. Results for filling fraction $p=0.50$,
corresponding to the number of blocked spins in bSM spin ice, are shown in the
right panel. A monopole moving with SM dynamics has a PSD that decays
approximately as $\nu^{-2}$ at all frequencies (left panel; orange line).

In the bSM version of spin ice, the ice rule constraint blocks $1/4$ of the
spins (if double monopoles are energetically ruled out) and the transverse
field constraint blocks $1/3$ of the remaining spins, leaving the system at an
effective bond filling $p=1/2$. While this value is larger than $p_{c}$ and
one would thus expect a three-dimensional cluster, the blocked spins are not
uncorrelated; most prominently, there is always at least one blocked spin at
every site. This type of correlation has been shown to increase the
percolation threshold in $d=2$ [57], and it is reasonable to expect an
effective percolation threshold for bSM spin ice larger than that of
uncorrelated bond correlation.

Indeed, a walker on the random-bond percolation cluster at $p=1/2$ shows a
smaller regime of fractal behaviour as compared to a monopole in bSM spin ice
(see Fig. 3 in the main text, and Fig. 8). The bSM results are instead well-
reproduced around $p=0.43$.

Figure 8:  Number of accessible sites within chemical distance $n$ of a random
reference site. Results are included for a monopole on spin ice, SM/bSM
(circles); for uncorrelated percolation, PC (crosses); and for percolation
with at most three bonds at every vertex, SPC (up and down triangles; note
that the down triangles are almost perfectly hidden by the red crosses). The
grey triangles show the behaviour when every vertex has exactly $3$ bonds,
distributed at random. Results for filling fraction $p=0.50$, corresponding to
the number of blocked spins in bSM spin ice, are shown in brown crosses and
magenta triangles. The filling fractions that give the best correspondence
with bSM dynamics, $p=0.43$ for the uncorrelated percolation cluster and
$p=0.46$ for the structured percolation cluster, are shown in red crosses and
yellow triangles respectively.

Following the method described above, one can ensure that no vertex in the
percolation model is connected by more than three bonds (generating structured
percolation clusters, SPC). Removing $1/4$ of the bonds in this way reproduces
the cluster growth of SM spin ice (see Fig. 8). Reducing the filling fraction
further one finds that the bSM results are now reproduced around $p=0.46$,
i.e. closer to $0.50$. It is reasonable to assume that the transverse field
distributions contribute some further correlations, and that these are
responsible for the remaining discrepancy. We have not studied these
correlations and their effects further, and they remain a potential topic of
future research.

It is perhaps worthy of notice how good an agreement can be found between
magnetic noise from single monopole motion in bSM nearest-neighbour spin ice
and an appropriately correlated bond percolation problem. Recall that the
kinematic constraints in the former depend on the local spin correlations,
which evolve in time; on the contrary, the bond percolation problem is
entirely static (quenched disorder). However, the difference appears to be
minimal within the time scales and system sizes accessed in our study.

###  B.3 Monopoles and percolation clusters

Here we also take the opportunity to expand on some of the less obvious
properties of monopoles with bSM dynamics, and how they compare to the
behaviour of a random walker on a percolation cluster. Firstly, one should
keep in mind that the percolation structure in spin ice is only meaningful in
relation to the monopoles. Furthermore, there is not a single percolation
structure; every monopole has its individual set of spins it can move through
and sites it can reach. As there is never a single monopole in the system,
these sets are also continuously changed by the motion of other monopoles, and
they may block and unblock paths for each other. In fact, this can lead to
monopoles being constrained to a small number of sites, where they are either
forced to remain until another monopole moves by to open up a new path for
them, or until they annihilate. For ${\cal H}_{\mathrm{NN}}$ at $T=0.6$ K, we
have found that at any given moment in time approximately $17$% of the
monopoles are constrained to a small number of sites ($\lesssim 60$), whereas
the remaining monopoles are able to access approximately $95$% of the sites in
the system. Monopoles constrained to a small number of sites generate a small
contribution to the magnetic noise, and only at relatively high frequencies.
As these are also less numerous, we have chosen to focus our analysis on
monopoles moving on the “percolation cluster”, i.e., the ones that can access
a majority of the system sites.

How fast the spin configuration around a monopole changes, i.e., the time
scale on which the monopole’s fractal cluster evolves, is determined by the
monopole density. A lower estimate of the time scale on which the fractals
evolve is given by measuring the time $\Delta_{t}$ between consecutive changes
of the transverse field across a spin. $\Delta_{t}$ is distributed according
to a temperature dependent exponential distribution
$\sim\exp{\left(\Delta_{t}/\tau_{B}(T)\right)}$, from which we can extract a
temperature dependent characteristic time $\tau_{B}(T)$ (see Fig. 9).
$\tau_{B}(T)$ is the typical time scale on which the local transverse field
evolves, and we find that $\log\left(\tau_{B}\right)\propto A/T$, with
$A\approx 6$ K ($A$ is similar to the monopole generation cost).

Figure 9: The characteristic time $\tau_{B}(T)$ on which the local transverse
field configuration changes. Extracted from fitting the distribution of times
$\Delta_{t}$ between consecutive changes of the transverse field across a spin
with a function $C\exp{\left[\Delta_{t}/\tau_{B}(T)\right]}$.

The evolution of the fractal cluster destroys the memory in the system, and
will set a time scale above which the monopole motion becomes diffusive.
However, as this time scale grows exponentially with inverse temperature, it
becomes asymptotically unimportant when compared to the dominant effect that
limits anomalous scaling: the fact that the system sits above the percolation
threshold. How memory is destroyed is explained further below, where we
discuss the effect of loop updates as an example of this.

###  B.4 Local one-dimensional motion

It is interesting to note that, between blocked spins in the absence of double
monopoles and spins that are prevented from flipping due to vanishing internal
transverse fields, a monopole in bSM spin ice moves on average along a random
one-dimensional path embedded on a three-dimensional diamond lattice.
Remarkably, random walk motion $n(t)$ along such path $\vec{x}_{n}$ produces
indeed anomalous diffusion: $\langle x^{2}_{n(t)}\rangle\sim\sqrt{\langle
n^{2}(t)\rangle}\sim\sqrt{t}$ and a PSD that decays as $\nu^{-1.5}$. (We thank
Shu Zhang at the Max Planck Institute for the Physics of Complex Systems for
pointing out a possible analogy to random walks on random one-dimensional
paths.) This is suggestive of a potential connection. However, the relevance
of correlations (e.g., self-avoidance) in one-dimensional-like paths in
monopole motion in bSM spin ice and/or on a critical percolation cluster
remains to be assessed, as well as the importance of a finite density of dead
ends and branching points which ought to be considered in these systems. It is
an intriguing direction for future work.

###  B.5  Effects of finite $\tau_{\mathrm{slow}}$ and loop updates

In Fig. 10, the PSD computed with bSM dynamics is shown for different ratios
$\tau_{\mathrm{slow}}/\tau_{\mathrm{fast}}$ between $1$ (equivalent to SM
dynamics) and $10^{4}$ (approximately the ratio predicted for Dy2Ti2O7 in Ref.
[15], using ${\cal H}_{\mathrm{NN}}$ at temperature $T=0.7$ K.

The same figure also shows the PSDs obtained when updates of the slow spins
are altogether forbidden, demonstrating that it is a good approximation within
the time scales of interest to our study for
$\tau_{\mathrm{slow}}/\tau_{\mathrm{fast}}\gtrsim 10^{3}$. For other spin ice
compounds, the difference between the two time scales may be significantly
smaller. For example, the ratio predicted for HTO is
$\tau_{\mathrm{slow}}/\tau_{\mathrm{fast}}\approx 10^{2}$ [15], and both time
scales ought to be accounted for in the simulations.

Figure 10:  PSD computed for ${\cal H}_{\mathrm{NN}}$ in Monte Carlo
simulations for different values of $\tau_{\mathrm{slow}}$, at $T=0.7$ K,
where $\tau_{\mathrm{slow}}/\tau_{\mathrm{fast}}=1$ is equivalent to SM
dynamics. For $\tau_{\mathrm{slow}}\gtrsim 10^{3}\tau_{\mathrm{fast}}$ the
spins with no transverse field can be well approximated as static, as
demonstrated by comparison to the dashed orange curve, which shows the PSD for
the bSM where updates of the slow spins are not allowed.

Evidence that the blocked spins retain memory, and that this is key to
generating the fractal structures the monopoles move on, becomes apparent if
we alternate bSM single spin flip dynamics with closed loop worm updates that
do not change the magnetisation nor move monopoles. In this case, we notice
that the PSD becomes less and less anomalous (see Fig. 11). In the results
presented here, we only considered hexagonal loops consisting of six spins,
but the general result remains unchanged if one includes larger closed loops.
To approximately recover the noise of SM dynamics for a system consisting of
$16\times 8^{3}=8192$ spins, it is sufficient to apply $4096$ hexagonal loop
updates per Monte Carlo sweep. However, the noise already becomes less
anomalous when one applies even a single hexagonal loop update per sweep. This
suggests that the magnetic noise in quantum spin ice, where a significant rate
of loop updates is expected, may be less anomalous – an effect that could be
taken as evidence of quantum coherent ring exchange processes in candidate
quantum spin ice systems.

Figure 11:  Applying closed loop updates in addition to single spin flip
attempts in bSM spin ice removes ‘spin memory’ from the system and
approximately recovers the noise curve of SM dynamics (black curve). The
coloured lines here show the PSD of bSM dynamics with different numbers of
hexagonal loop updates applied per Monte Carlo sweep, varying from $0$ to
$8^{5}$. These results were obtained in thermal equilibrium at $T=0.7$ K for
${\cal H}_{\mathrm{NN}}$, and we used a system consisting of $8192$ spins.