# The role of correlated hopping in many-body physics of flat-band systems:
Nagaoka ferromagnetism

Tom Westerhout<EMAIL_ADDRESS>Mikhail I. Katsnelson
<EMAIL_ADDRESS>Institute for Molecules and Materials, Radboud
University, Heijendaalseweg 135, 6525 AJ Nijmegen, The Netherlands

###### Abstract

In narrow-band systems, correlated contribution to effective hopping becomes
important, and one needs to carefully consider three types of hopping
processes: between doubly and singly occupied sites, between singly occupied
and empty sites, and between doubly occupied and empty or two singly occupied
sites. All three hopping parameters cannot vanish simultaneously, and one
should specify for which of these processes the band becomes flat. By means of
exact diagonalization of finite systems we demonstrate that three hoppings
play qualitatively different roles in many-body effects, in particular, in the
formation of half-metallic ferromagnetic (Nagaoka) states.

In solids, the narrower the electronic band the more important the correlation
effects are. A prototypical example of an electron system with a flat band is
the quantum Hall regime for a two-dimensional electron gas in a magnetic field
when the kinetic energy is suppressed by Landau quantization [1]. In this
situation many-body effects become dominant leading, for example, to
fractional quantum Hall effect, and all relevant terms in the Hamiltonian
originate from the Coulomb interaction [2]. The conventional perturbative
techniques, such as Feynman diagrammatic expansion, are then not applicable
and more advanced and subtle theoretical approaches are required [3].

Flat bands can also be observed in solids even without a quantized magnetic
field. The most well-known example is the magic-angle twisted bilayer graphene
[4, 5, 6] where the formation of a flat band is accompanied by a broad variety
of strongly correlated effects and electronic phase transitions including
unconventional superconductivity [7, 8, 9, 10]. Another example is a single
layer of InSe where an almost flat hole band arises without any twisting. The
appearance of this band is accompanied, among other things, by an anomalously
strong electron-phonon interaction and closeness to ferromagnetic and charge-
density-wave instabilities [11, 12]. Strong correlation effects are thus
practically unavoidable in narrow-band systems.

Recently, it was emphasized [13] that in the flat-band regime, kinetic energy
contribution from the Coulomb interaction, a so called _correlated hopping_ ,
becomes relevant. This concept is very old, originating from the polar model
of Shubin and Vonsovsky in 1934 [14]. The polar model, including the
correlated hopping effects, has since then been reconsidered with contemporary
theoretical tools such as Hubbard X-operators in Refs. [15, 16].
Independently, the concept was introduced in Ref. [17]. A possible role of
correlated hopping in the physics of high-temperature cuprate
superconductivity was discussed by Hirsch [18, 19]. In one dimension, the
Hubbard-like model with correlated hopping can be solved exactly in some cases
[20, 21]. The two-dimensional case has been considered numerically by means of
exact diagonalization [22, 23] and Quantum Monte Carlo [24] methods.

The discovery of real two-dimensional materials with flat (or almost flat)
bands requires more careful and detailed studies of the correlation effects,
one by one. In this Letter we focus on the appearance of a saturated, half-
metallic ferromagnetic state typical for an almost half-filled narrow-band
case. We call this phenomenon Nagaoka ferromagnetism due to the seminal work
[25]. Nagaoka has proven rigorously that in the infinite-$U$ Hubbard model
near half-filling, with just one extra or missing electron, the ground state
is a saturated ferromagnet with maximum possible total spin. The result holds
for various two- and three-dimensional crystal structures, but does not work
in one dimension. This initial formulation is, strictly speaking, insufficient
since the gap between this ground state and excited states vanishes in
thermodynamic limit. However, it was proven by Irkhin and Katsnelson via
virial expansion for X-operator Green’s functions [26] that Nagaoka’s result
remains valid for finite but small concentrations of extra (or missing)
electrons. Further increase in electron (or hole) concentration results in a
transition to a non saturated ferromagnetic state [27, 28, 29], but rigorous
results about this transition are still missing. The case of a one-dimensional
system with next- and next-nearest-neighbor hoppings was considered in Ref.
30. Rigorous results for ferromagnetism in the Hubbard model on a special
class of decorated lattices with degenerate ground states were obtained in
Ref. 31. Ferromagnetism in the Hubbard model was also studied using dynamical
mean-field theory which is formally associated to the limit of infinite
dimension [32, 33, 34]. It is impossible to give a full list of references on
the topic, but important works [35, 36, 37, 38, 39, 40, 41, 42] should be
mentioned to illustrate different approaches to the problem.

We start our analysis with the standard Hubbard model with nearest-neighbor
hopping:

$H_{\text{Hubbard}}=t\sum_{\begin{subarray}{c}\langle i,j\rangle,\;\sigma\\\
i<j\end{subarray}}(c^{\dagger}_{i\sigma}c^{\vphantom{\dagger}}_{j\sigma}+\text{H.c.})+U\sum_{i}n_{i\uparrow}n_{i\downarrow}\,.$
(1)

Here, $c^{\dagger}_{i\sigma}$ and $c^{\vphantom{\dagger}}_{i\sigma}$ are
fermionic creation and annihilation operators respectively, creating or
annihilating an electron at site $i$ with spin
$\sigma=\\{\uparrow,\downarrow\\}$;
$n_{i\sigma}=c^{\dagger}_{i\sigma}c^{\vphantom{\dagger}}_{i\sigma}$; $\langle
i,j\rangle$ denote the nearest neighbors.

We focus on the almost half-filled case when the number of electrons is one
fewer than the number of sites. For $U\to\infty$ the ground state of the
system is ferromagnetic, at least, for some types of crystal lattices [25]. In
the following, we will consider the effect of narrow or even flat bands on the
stability of this ferromagnetic phase.

Figure 1: Phase diagrams of (a) $2\times 2$ plaquette with 3 electrons and (b)
$3\times 3$ plaquette with 8 electrons as functions of $F$ and $U$. Points are
numerical data and lines are a guide to eye. Different colors correspond to
different values of the total spin $S$ defined in Eq. (4).

To this end, instead of the Hamiltonian (1) with a single hopping parameter
$t$, we consider a more complete polar model [14, 15, 16] which also contains
hopping like terms originating from the electron-electron interaction. These
terms contain three fermionic operators belonging to one site $i$ and one
fermionic operator belonging to another site $j$. In the nearest-neighbor
approximation they are described by one parameter,

$F\\!=\left\langle ii\left|v\right|ij\right\rangle\\!=\\!\\!\int d{\bf r}d{\bf
r}^{\prime}\psi_{i}^{\ast}({\bf r})\psi_{i}^{\ast}({\bf
r}^{\prime})v\left({\bf r-r}^{\prime}\right)\psi_{i}({\bf r})\psi_{j}({\bf
r}^{\prime})\,,$ (2)

where $v\left({\bf r-r}^{\prime}\right)$ is the potential of inter-electron
interaction, $\psi_{i}({\bf r})$ is the state of the electron at site $i$ and
the sites $i$, $j$ are nearest neighbors.

After adding such terms to the Hamiltonian (1), we obtain

$H=U\sum_{i}n_{i\uparrow}n_{i\downarrow}+t\\!\\!\sum_{\begin{subarray}{c}\langle
i,j\rangle,\;\sigma\\\
i<j\end{subarray}}\\!(c^{\dagger}_{i\sigma}c^{\vphantom{\dagger}}_{j\sigma}+\text{H.c.})\\\
+F\\!\\!\sum_{\begin{subarray}{c}\langle i,j\rangle,\;\sigma\\\
i<j\end{subarray}}\\!(n_{i\bar{\sigma}}+n_{j\bar{\sigma}})(c^{\dagger}_{i\sigma}c^{\vphantom{\dagger}}_{j\sigma}+\text{H.c.})\,.$

Here, $\bar{\sigma}$ denotes the inverse of $\sigma$ such that if
$\sigma=\,\uparrow$ then $\bar{\sigma}=\,\downarrow$, and if
$\sigma=\,\downarrow$ then $\bar{\sigma}=\,\uparrow$. By considering different
occupations of the sites connected by the hopping process, we can rewrite the
Hamiltonian as

$H\\!=U\sum_{i}n_{i\uparrow}n_{i\downarrow}+\\!\\!\sum_{\begin{subarray}{c}\langle
i,j\rangle,\;\sigma\\\
i<j\end{subarray}}\\!\bigg{[}\beta^{\prime\prime}(1-n_{i\bar{\sigma}})(1-n_{j\bar{\sigma}})+\beta^{\prime}n_{i\bar{\sigma}}n_{j\bar{\sigma}}\\\
+\gamma\big{(}n_{i\bar{\sigma}}(1-n_{j\bar{\sigma}})+(1-n_{i\bar{\sigma}})n_{j\bar{\sigma}}\big{)}\bigg{]}(c^{\dagger}_{i\sigma}c^{\vphantom{\dagger}}_{j\sigma}+\text{H.c.})\,,$
(3)

where

$\left\\{\begin{aligned} \beta^{\prime\prime}&=t\\\ \beta^{\prime}&=t+2F\\\
\gamma&=t+F\end{aligned}\right.\,.$

In the above, we have used the Wannier basis which is orthonormal by
construction. The corresponding expressions in a more general case of non-
orthogonal basis functions (for example, atomic orbitals) are derived in Ref.
[43].

Generally speaking, the parameter $F$ is by no means small. Recently, it was
estimated from the first-principle constrained random phase approximation
(cRPA) method for benzene molecule [44]. The results are: $F=0.53\;$eV,
$t=-2.75\;$eV. For comparison, Hubbard $U$ for this system is $U=9.38\;$eV,
nearest-neighbor Coulomb interaction is $V_{\text{NN}}=6.12\;$eV, and nearest-
neighbor direct exchange interaction is $J=0.18\;$eV. Note that the parameters
for benzene molecule are similar to those of graphene [45].

The Hamiltonian (3) is part of the complete polar model [15] which originally
also contained the intersite Coulomb interaction. Here $\beta^{\prime\prime}$,
$\beta^{\prime}$, and $\gamma$ are parameters describing hopping between an
empty and a singly occupied site, between a doubly and singly occupied site,
and between two singly occupied sites (or between doubly occupied and an empty
site), respectively. To avoid confusion, note that our definition of the
parameter $\beta^{\prime\prime}$ differs by a minus sign from the original one
in Ref. [15]. Importantly, at $F\neq 0$ the three hopping parameters cannot
all be zero, but in the following we analyze the effect of zeroing one of them
at a time. We will see that these three hopping parameters play different
roles in the formation of magnetic ground state of the system under
consideration. First, we present the numerical results and then give their
qualitative explanation.

We focus on the square lattice and perform exact diagonalization for small
clusters to study how the total spin in the ground state depends on $t$, $F$,
and $U$. The QuSpin package [46] is used for exact diagonalization. The total
spin is defined by

$S(S+1)=\sum_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}=\sum_{ij}(S^{x}_{i}S^{x}_{j}+S^{y}_{i}S^{y}_{j}+S^{z}_{i}S^{z}_{j})\,,$
(4)

where

$S^{x}_{i}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}c^{\dagger}_{i\sigma^{\prime}}\sigma^{x}_{\sigma^{\prime}\sigma}c^{\vphantom{\dagger}}_{i\sigma}$

and $S^{y}$ and $S^{z}$ are defined similarly. $\sigma^{x}$, $\sigma^{y}$,
$\sigma^{z}$ are Pauli matrices.

Let us first consider a $2\times 2$ cluster. Even though the system is very
small, it still reproduces an experimentally realizable situation of a
quadruple quantum dot where Nagaoka ferromagnetic ground state was recently
observed [47]. In Fig. 1(a) we show the phase diagram of the system. Since for
three electrons there are only two possible values of the total spin, $1/2$
and $3/2$, we find two phases. The horizontal line at $F=0$ corresponds to the
standard Hubbard model, and, as expected, the system is antiferromagnetic at
$U=0$ and becomes ferromagnetic for $U/t\gtrsim 18.4$.

The main difference between the Hubbard model with $U=0$ and $U\to\infty$, is
the “Gutzwiller projector” forbidding the appearance of doubly occupied sites
[48]. Putting $\gamma=0$ has the same consequences as $U\to\infty$ since it
makes impossible the creation (as well as annihilation) of doubly occupied and
empty sites from a pair of singly occupied sites. This means that with
$\gamma=0$ the system behaves as an infinitely strongly correlated one even
for small $U$, and one could expect that Nagaoka ferromagnetism arises even if
the Coulomb interaction is not large. We indeed observe this behavior in Fig.
1(a) where for $\gamma=0$ (in other words, at $F=-1$) ferromagnetism is
favored even at very small $U$. At the same time, setting $F$ close enough to
$0$, that is considering small $\beta^{\prime\prime}$, makes ferromagnetism
impossible.

The simple example of a $2\times 2$ cluster already clearly shows
qualitatively different roles of the hoppings parameters in Eq. (3): whereas
small $\gamma$ favors ferromagnetism, small $\beta^{\prime\prime}$ destroys
it. In Fig. 1(b) we show the phase diagram for a $3\times 3$ cluster with
eight electrons. Here, the situation is very similar: for $\gamma=0$ and
$U/t\gtrsim 4.2$, Nagaoka ferromagnetic ground state is realized, and for
$\beta^{\prime\prime}=0$ the system is always antiferromagnetic. One
additional feature which we observe for a larger system is the stabilization
of a non-saturated ferromagnetic phase for small $U$ and $\gamma$. In the
Supplemental Material [49] we also provide data for clusters of $5$, $6$, $8$,
$10$, $12$, and $13$ sites—they all demonstrate qualitatively similar
behavior.

Let us now discuss qualitatively the physical origin of Nagaoka
ferromagnetism. First, note that $\beta^{\prime}$ and $\beta^{\prime\prime}$
are dual with respect to particle-hole inversion. Hence in the following,
without loss of generality, we discuss the role of $\beta^{\prime\prime}$
since in the numerical simulations we had one extra hole.

Nagaoka ferromagnetism can be understood in terms of effective narrowing of
the electronic band when an electron moves in a non-ferromagnetic
environment—the mechanism known as double exchange [50]. The band edges remain
the same for any spin configuration, but the density of states drops
exponentially in case of random or antiferromagnetic spin configurations,
forming so called Lifshitz tails [51, 52, 53]. Thus, for any small but finite
charge carrier occupation, the lowest average band energy is obtained for a
ferromagnet. An explicit spin Hamiltonian describing this tendency can be
derived in the static approximation for the narrow-band Hubbard model on Bethe
lattice [54]; it confirms the previously known non-Heisenberg form of the
corresponding exchange interaction which is proportional to
$J\cos(\theta_{ij}/2)$, where $\theta_{ij}$ is the angle between magnetic
moments at sites $i$ and $j$ [50].

Importantly, the coupling constant $J$ is actually $\beta^{\prime\prime}$
since this band narrowing effect is related to the hopping between singly
occupied and empty sites. This explains the tendency to ferromagnetism at
finite $\beta^{\prime\prime}$ and vanishing $\gamma$; the latter makes the
system strongly correlated and the effective $U$ infinitely large, forbidding
the appearance of doubly occupied sites.

The role of finite $\gamma$ is the exact opposite. Creation and annihilation
of doubly occupied and empty site pairs decrease the energy of the
antiferromagnetic state resulting in indirect (“kinetic”) antiferromagnetic
exchange [50]. Its value in the large-$U$ Hubbard model is proportional to
$-t^{2}/U$ and in the full polar model it becomes $-\gamma^{2}/U$. Therefore,
finite $\gamma$ favors antiferromagnetism.

One charge carrier (either an electron or a hole) put in an antiferromagnetic
environment forms a magnetic polaron (also known as fluctuon or ferron)
creating around itself a local ferromagnetic environment [55, 56, 57, 58, 59,
60, 52, 53, 61]. In the simplest model of an infinite potential well in two
dimensions [61], the energy of formation of such a ferromagnetic droplet can
be estimated as

$\mathcal{E}(R)=-\Delta+\frac{\hbar^{2}z_{0}^{2}}{2m^{*}R^{2}}+J\frac{\pi
R^{2}}{S_{0}}\,,$

where $\Delta$ is the energy difference between the cases of antiferromagnetic
and ferromagnetic environments, $R$ is the radius of the well, $z_{0}=2.40483$
is the first zero of the Bessel function, $J_{0}(z)$, $m^{*}$ is the effective
mass, $J$ is the antiferromagnetic indirect exchange, and $S_{0}$ is the area
per magnetic atom. After optimization with respect to $R$, we estimate the
magnetic polaron energy as

$\mathcal{E}=-\Delta+C\sqrt{WJ},$

where $W$ is the total bandwidth and $C$ is some numerical factor (its value
can be found in Ref. [61]). To be stable, the magnetic polaron should have
negative energy. Assuming that both $\Delta$ and $W$ are proportional to
$\beta^{\prime\prime}$ and that $J\propto\gamma^{2}/U$, one obtains the
stability criterion

$\gamma^{2}<D|\beta^{\prime\prime}|U,$ (5)

where $D$ is another numerical factor. Eq. (5) is too rough to be
quantitatively correct, especially for small clusters, but qualitatively it
does explain the opposite roles of $\beta^{\prime\prime}$ and $\gamma$—the
first one stabilizing and the latter destabilizing Nagaoka ferromagnetism.

To conclude, in this work we have analyzed the role of correlated hopping in
narrow-band systems focusing on one particular many-body phenomenon—Nagaoka
ferromagnetism. By performing exact diagonalization of small clusters we have
shown that for one extra hole (electron) $\beta^{\prime\prime}$
($\beta^{\prime}$) favors the ferromagnetic ground state whereas $\gamma$
favors the antiferromagnetic phase. Smallness of $\gamma$ acts as the
Gutzwiller projector making the system strongly correlated even at small
values of the Hubbard $U$. We also provide qualitative explanation of our
numerical results.

The code to carry out the analysis is publicly available at
http://github.com/twesterhout/correlated-hoppings.

This work is supported by European Research Council via Synergy Grant No.
854843 - FASTCORR. Numerical simulations in this work were carried out on the
Dutch national e-infrastructure with the support of SURF Cooperative.

## References

  * Girvin and Prange [1987] S. M. Girvin and R. Prange,  _The quantum Hall effect_ (1987).
  * Murthy and Shankar [2003] G. Murthy and R. Shankar, Hamiltonian theories of the fractional quantum Hall effect, Rev. Mod. Phys. 75, 1101 (2003).
  * Hansson  _et al._ [2017] T. H. Hansson, M. Hermanns, S. H. Simon, and S. F. Viefers, Quantum Hall physics: Hierarchies and conformal field theory techniques, Rev. Mod. Phys. 89, 025005 (2017).
  * Suárez Morell  _et al._ [2010] E. Suárez Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Flat bands in slightly twisted bilayer graphene: Tight-binding calculations, Phys. Rev. B 82, 121407 (2010).
  * Bistritzer and MacDonald [2011] R. Bistritzer and A. H. MacDonald, Moiré bands in twisted double-layer graphene, Proc. Natl. Acad. Sci. 108, 12233 (2011).
  * Li  _et al._ [2009] G. Li, A. Luican, J. M. B. Lopes dos Santos, A. H. Castro Neto, A. Reina, J. Kong, and E. Y. Andrei, Observation of van Hove singularities in twisted graphene layers, Nat. Phys. 6, 109 (2009).
  * Cao  _et al._ [2018a] Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori, and P. Jarillo-Herrero, Correlated insulator behaviour at half-filling in magic-angle graphene superlattices, Nature 556, 80 (2018a).
  * Cao  _et al._ [2018b] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, Nature 556, 43 (2018b).
  * Yankowitz  _et al._ [2019] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, Science 363, 1059 (2019).
  * Sharpe  _et al._ [2019] A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Emergent ferromagnetism near three-quarters filling in twisted bilayer graphene, Science 365, 605 (2019).
  * Lugovskoi  _et al._ [2019] A. V. Lugovskoi, M. I. Katsnelson, and A. N. Rudenko, Strong electron-phonon coupling and its influence on the transport and optical properties of hole-doped single-layer InSe, Phys. Rev. Lett. 123, 176401 (2019).
  * Stepanov  _et al._ [2022] E. A. Stepanov, V. Harkov, M. Rösner, A. I. Lichtenstein, M. I. Katsnelson, and A. N. Rudenko, Coexisting charge density wave and ferromagnetic instabilities in monolayer InSe, npj Comput. Mater. 8 (2022).
  * Cea  _et al._ [2022] T. Cea, P. A. Pantaleón, N. R. Walet, and F. Guinea, Electrostatic interactions in twisted bilayer graphene, Nano Materials Science 4, 27 (2022).
  * Schubin and Vonsovsky [1934] S. Schubin and S. V. Vonsovsky, On the electron theory of metals, Proc. R. Soc. London A. 145, 159 (1934).
  * Vonsovsky and Katsnelson [1979a] S. V. Vonsovsky and M. I. Katsnelson, Some types of instabilities in the electron energy spectrum of the polar model of the crystal. I. The maximum-polarity state, J. Phys. C: Solid State Phys. 12, 2043 (1979a).
  * Vonsovsky and Katsnelson [1979b] S. V. Vonsovsky and M. I. Katsnelson, Some types of instabilities in the electron energy spectrum of the polar model of the crystal. II. The criterion of stability of a metallic state, J. Phys. C: Solid State Phys. 12, 2055 (1979b).
  * Foglio and Falicov [1979] M. E. Foglio and L. M. Falicov, New approach to the theory of intermediate valence. I. General formulation, Phys. Rev. B 20, 4554 (1979).
  * Hirsch [1989a] J. Hirsch, Hole superconductivity, Phys. Lett. A 134, 451 (1989a).
  * Hirsch [1989b] J. E. Hirsch, Bond-charge repulsion and hole superconductivity, Phys. C Supercond. Its Appl. 158, 326 (1989b).
  * Arrachea  _et al._ [1996] L. Arrachea, A. A. Aligia, and E. Gagliano, Anomalous flux quantization in a Hubbard ring with correlated hopping, Phys. Rev. Lett. 76, 4396 (1996).
  * Aligia  _et al._ [2000] A. A. Aligia, K. Hallberg, C. D. Batista, and G. Ortiz, Phase diagrams from topological transitions: The Hubbard chain with correlated hopping, Phys. Rev. B 61, 7883 (2000).
  * Gagliano  _et al._ [1995] E. R. Gagliano, A. A. Aligia, L. Arrachea, and M. Avignon, Single-particle spectral function of a generalized Hubbard model: Metal-insulator transition, Phys. Rev. B 51, 14012 (1995).
  * Arrachea and Aligia [2000] L. Arrachea and A. A. Aligia, Pairing correlations in a generalized Hubbard model for the cuprates, Phys. Rev. B 61, 9686 (2000).
  * Liberto  _et al._ [2014] M. D. Liberto, C. E. Creffield, G. I. Japaridze, and C. M. Smith, Quantum simulation of correlated-hopping models with fermions in optical lattices, Phys. Rev. A 89, 013624 (2014).
  * Nagaoka [1966] Y. Nagaoka, Ferromagnetism in a narrow, almost half-filled $s$ band, Phys. Rev. 147, 392 (1966).
  * Irkhin and Katsnelson [1985] V. Y. Irkhin and M. I. Katsnelson, Spin waves in narrow band ferromagnet, J. Phys. C: Solid State Phys. 18, 4173 (1985).
  * Linden and Edwards [1991] W. v. d. Linden and D. M. Edwards, Ferromagnetism in the Hubbard model, J. Phys. Condens. Matter 3, 4917 (1991).
  * Irkhin and Katsnelson [1990] V. Y. Irkhin and M. I. Katsnelson, Ground state and electron-magnon interaction in an itinerant ferromagnet: Half-metallic ferromagnets, J. Phys. Condens. Matter 2, 7151 (1990).
  * Irkhin and Zarubin [2004] V. Y. Irkhin and A. V. Zarubin, Density-of-states picture and stability of ferromagnetism in the highly correlated Hubbard model, Phys. Rev. B 70, 035116 (2004).
  * Müller-Hartmann [1995] E. Müller-Hartmann, Ferromagnetism in Hubbard models: Low density route, J. Low Temp. Phys. 99, 349 (1995).
  * Mielke and Tasaki [1993] A. Mielke and H. Tasaki, Ferromagnetism in the Hubbard model, Comm. Math. Phys. 158, 341 (1993).
  * Uhrig [1996] G. S. Uhrig, Exact single spin flip for the Hubbard model in $d=\infty$, Phys. Rev. Lett. 77, 3629 (1996).
  * Held  _et al._ [1996] K. Held, M. Ulmke, and D. Vollhardt, Correlated-electron theory of metamagnetism in strongly anisotropic antiferromagnets, Mod. Phys. Lett. B 10, 203 (1996).
  * Vollhardt  _et al._ [1996] D. Vollhardt, N. Blümer, K. Held, M. Kollar, J. Schlipf, and M. Ulmke, Non-perturbative approaches to magnetism in strongly correlated electron systems, Z. Phys. B Condens. Matter 103, 283 (1996).
  * Tasaki [1998] H. Tasaki, From Nagaoka’s ferromagnetism to flat-band ferromagnetism and beyond: An introduction to ferromagnetism in the Hubbard model, Prog. Theor. Phys. 99, 489 (1998).
  * Maksymenko  _et al._ [2012] M. Maksymenko, A. Honecker, R. Moessner, J. Richter, and O. Derzhko, Flat-band ferromagnetism as a Pauli-correlated percolation problem, Phys. Rev. Lett. 109, 096404 (2012).
  * Mielke [1999] A. Mielke, Ferromagnetism in single-band Hubbard models with a partially flat band, Phys. Rev. Lett. 82, 4312 (1999).
  * Tamura  _et al._ [2002] H. Tamura, K. Shiraishi, T. Kimura, and H. Takayanagi, Flat-band ferromagnetism in quantum dot superlattices, Phys. Rev. B 65, 085324 (2002).
  * Roth [1969] L. M. Roth, Electron correlation in narrow energy bands. II. One reversed spin in an otherwise fully aligned narrow $s$ band, Phys. Rev. 186, 428 (1969).
  * Jarrett  _et al._ [1968] H. Jarrett, W. Cloud, R. Bouchard, S. Butler, C. Frederick, and J. Gillson, Evidence for itinerant d-electron ferromagnetism, Phys. Rev. Lett. 21, 617 (1968).
  * Plischke [1974] M. Plischke, Ferromagnetism in the single-band Hubbard model: An exact high-temperature expansion, J. Stat. Phys. 11, 159 (1974).
  * Hanisch  _et al._ [1997] T. Hanisch, G. S. Uhrig, and E. Müller-Hartmann, Lattice dependence of saturated ferromagnetism in the Hubbard model, Phys. Rev. B 56, 13960 (1997).
  * Irkhin [1998] V. Y. Irkhin, Derivation of the $t$–$J$ model: Electron spectrum and exchange interactions in narrow energy bands, Phys. Rev. B 57, 13375 (1998).
  * van Loon  _et al._ [2021] E. G. C. P. van Loon, M. Rösner, M. I. Katsnelson, and T. O. Wehling, Random phase approximation for gapped systems: Role of vertex corrections and applicability of the constrained random phase approximation, Phys. Rev. B 104, 045134 (2021).
  * Schüler  _et al._ [2013] M. Schüler, M. Rösner, T. O. Wehling, A. I. Lichtenstein, and M. I. Katsnelson, Optimal Hubbard models for materials with nonlocal Coulomb interactions: Graphene, silicene, and benzene, Phys. Rev. Lett. 111, 036601 (2013).
  * Weinberg and Bukov [2019] P. Weinberg and M. Bukov, QuSpin: A python package for dynamics and exact diagonalisation of quantum many body systems. Part II: Bosons, fermions and higher spins, SciPost Physics 7, 020 (2019).
  * Dehollain  _et al._ [2020] J. P. Dehollain, U. Mukhopadhyay, V. P. Michal, Y. Wang, B. Wunsch, C. Reichl, W. Wegscheider, M. S. Rudner, E. Demler, and L. M. K. Vandersypen, Nagaoka ferromagnetism observed in a quantum dot plaquette, Nature 579, 528 (2020).
  * Vollhardt [1984] D. Vollhardt, Normal 3He: An almost localized Fermi liquid, Rev. Mod. Phys. 56, 99 (1984).
  * [49] See Supplemental Material for data on other cluster sizes.
  * Anderson [1963] P. W. Anderson, Theory of magnetic exchange interactions: Exchange in insulators and semiconductors, in  _Solid State Physics_, Vol. 14 (Elsevier, 1963) pp. 99–214.
  * Brinkman and Rice [1970] W. F. Brinkman and T. M. Rice, Single-particle excitations in magnetic insulators, Phys. Rev. B 2, 1324 (1970).
  * Auslender and Katsnelson [2005] M. I. Auslender and M. I. Katsnelson, Electron self-trapping and the fluctuation density-of-states tail at the critical point, Phys. Rev. B 72, 113107 (2005).
  * Auslender and Katsnelson [2006] M. I. Auslender and M. I. Katsnelson, Electron self-trapping at quantum and classical critical points, Ann. Phys-new. York. 321, 1762 (2006).
  * Auslender and Katsnelson [1982] M. I. Auslender and M. I. Katsnelson, The effective spin Hamiltonian and phase separation instability of the almost half-filled Hubbard model and the narrow-band s-f model, Solid State Commun. 44, 387 (1982).
  * Krivoglaz [1973] M. A. Krivoglaz, Fluctuon states of electrons, Uspekhi Fizicheskih Nauk 111, 617 (1973).
  * Nagaev [1983] E. L. Nagaev, Physics of magnetic semiconductors MIR publishers (1983).
  * Nagaev [2001] E. L. Nagaev, Colossal-magnetoresistance materials: Manganites and conventional ferromagnetic semiconductors, Phys. Rep. 346, 387 (2001).
  * Visscher [1974] P. B. Visscher, Phase separation instability in the Hubbard model, Phys. Rev. B 10, 943 (1974).
  * Mott [1973] N. F. Mott, Metal-insulator transitions, Contemp. Phys. 14, 401 (1973).
  * Auslender and Katsnelson [1981] M. I. Auslender and M. I. Katsnelson, Magnetic susceptibility of the spin polaron states in the s-f exchange model above curie temperature, J. Magn. Magn. Mater. 24, 117 (1981).
  * Soriano and Katsnelson [2020] D. Soriano and M. I. Katsnelson, Magnetic polaron and antiferromagnetic-ferromagnetic transition in doped bilayer CrI3, Phys. Rev. B 101, 041402 (2020).


Current address: ]Department of Physics, University of Cagliari, S. P.
Monserrato-Sestu km 0,700, 09042 Monserrato (CA), Italy.

Current address: ]Department of Materials Science and Engineering, Uppsala
University, Box 534, SE-751 21 Uppsala, Sweden.

Current address: ]Fachbereich Physik and Forschungszentrum OPTIMAS, Technische
Universität Kaiserslautern, 67663 Kaiserslautern, Germany.

# Collective magnetic dynamics in artificial spin ice probed by AC
susceptibility

Merlin Pohlit<EMAIL_ADDRESS>Department of Physics and
Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden  Giuseppe
Muscas  [  Department of Physics and Astronomy, Uppsala University, Box 516,
SE-75120 Uppsala, Sweden  Ioan-Augustin Chioar  Department of Physics and
Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden  Henry
Stopfel  [  Department of Physics and Astronomy, Uppsala University, Box 516,
SE-75120 Uppsala, Sweden  Agne Ciuciulkaite  Department of Physics and
Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden  Erik Östman
Department of Physics and Astronomy, Uppsala University, Box 516, SE-75120
Uppsala, Sweden  Spyridon D. Pappas  [  Department of Physics and Astronomy,
Uppsala University, Box 516, SE-75120 Uppsala, Sweden  Aaron Stein  Center of
Functional Nanomaterials, Brookhaven National Laboratory, P.O. Box 5000,
Upton, New York 11973, USA  Björgvin Hjörvarsson  Department of Physics and
Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden  Petra E.
Jönsson  Department of Physics and Astronomy, Uppsala University, Box 516,
SE-75120 Uppsala, Sweden  Vassilios Kapaklis<EMAIL_ADDRESS>Department of Physics and Astronomy, Uppsala University, Box 516, SE-75120
Uppsala, Sweden

###### Abstract

We report on the study of the thermal dynamics of square artificial spin ice,
probed by means of temperature and frequency dependent AC susceptibility.
Pronounced influence of the inter-island coupling strength was found on the
frequency response of the samples. Through the subsequent analysis of the
frequency- and coupling-dependent freezing temperatures, we discuss the
phenomenological parameters obtained in the framework of Vogel-Fulcher-Tammann
law in terms of the samples microscopic features. The high sensitivity and
robust signal to noise ratio of AC susceptibility validates the latter as a
promising and simple experimental technique for resolving the dynamics and
temperature driven dynamics crossovers for the case of artificial spin ice.

##  I Introduction

Artificial Spin Ice (ASI), i.e. arrays of magnetostatically coupled
ferromagnetic islands – mesospins Östman  _et al._ (2018a) – fabricated by
nanolithography Wang  _et al._ (2006); Nisoli  _et al._ (2013); Heyderman and
Stamps (2013); Rougemaille and Canals (2019), exhibit collective phenomena,
and importantly, their interaction strength and geometry can be tailored
almost at will Kapaklis  _et al._ (2014); Perrin  _et al._ (2016); Östman  _et
al._ (2018b); Nisoli  _et al._ (2017); Pohlit  _et al._ (2015, 2016). Properly
designed to support thermal fluctuations, ASI systems can serve as a platform
for the investigations of thermal magnetization dynamics and freezing
transitions in tailored nanostructures, which can also be used to mimic the
dynamical properties of frustrated, naturally-occurring magnetic spin systems
Kapaklis  _et al._ (2012); Morgan  _et al._ (2010); Farhan  _et al._ (2013);
Kapaklis  _et al._ (2014); Porro  _et al._ (2013); Östman  _et al._ (2018b).
Insights into the freezing transition and the nature of the frozen low
temperature states were obtained by investigations using magnetometry
Andersson  _et al._ (2016) and synchrotron-based scattering- and microscopy-
techniques Morley  _et al._ (2017); Sendetskyi  _et al._ (2019); Saccone  _et
al._ (2019). With the exception of early work based on temperature dependent
magneto-optical measurements Kapaklis  _et al._ (2012) and more recent works
using synchrotron-based magnetic microscopy Kapaklis  _et al._ (2014) and muon
relaxation Anghinolfi  _et al._ (2015); Leo  _et al._ (2018), experimental
studies of thermally induced transitions are scarce. Furthermore, experimental
techniques based on synchrotron radiation and muons impose limitations on
availability and accessible time-scales. To this end, AC susceptibility is a
well established and accessible technique for probing magnetization dynamics,
giving access to a wide frequency range Bedanta and Kleemann (2008); Topping
and Blundell (2018).

In this work, we report on AC susceptibility measurements of thermally active
extended square ASI arrays measured in a wide frequency and temperature range.
We study arrays that are composed of close to identical mesospins, with
different gaps between the elements, in order to explore the influence of
coupling strength on the collective dynamics. Exploring the frequency
dependence of the AC susceptibility signal we employ the VFT law, that can be
used for describing the low-field magnetic relaxation of weakly interacting
nanoparticle systems Landi (2013); Vernay  _et al._ (2014) but has recently
been applied also to ASI systemsAndersson  _et al._ (2016); Morley  _et al._
(2017), attempting to extract parameters that can be directly related to the
magnetostatic energies of the ASI arrays. We discuss the validity of this
simplified approach and address the limitations of such models, in the
framework of thermal ASI arrays.

##  II Experimental Details

The extended square ASI structures were produced by post-patterning of
$\delta$-doped Pd(Fe) thin films Pärnaste  _et al._ (2007), employing
electron-beam lithography (EBL) followed by argon-ion milling. The films,
consisting of $40\,\rm{nm}$ Palladium, 2.2 monolayers of Iron, and a
$2\,\rm{nm}$ Palladium capping layer, were all grown on a Vanadium seeding
layer on top of Magnesium oxide (MgO(001)) substrates by DC magnetron
sputtering. The total effective thickness of the magnetic layer (Fe and
magnetically polarized Pd) was previously estimated to be $2\,\rm{nm}$
Pärnaste  _et al._ (2007); Hase  _et al._ (2014). Vibrating Sample
Magnetometetry (VSM) revealed that the temperature dependence of the in-field
volume magnetization can be described by:
$M_{\rm{s}}(T)=M_{\rm{0}}(1-T/T_{\rm{0}})^{0.5}$, with a Curie temperature
$T_{\rm{0}}=410\,\rm{K}$ Andersson  _et al._ (2016). The size of the islands,
the lattice parameter and the distance between the islands were determined
after the EBL process using Scanning Electron Microscopy (SEM). A typical SEM
image is shown as an inset of Fig. 1. All the islands have a length of $310\pm
15\,\rm{nm}$ and a width of $140\pm 15\,\rm{nm}$. Arrays with different inter-
island lattice spacing, $d$, of 380, 420 and 460 nm were prepared which yields
different gaps between the islands ($g$ = 70, 110, and 150 nm). The difference
in distances between the islands results in different magnetostatic coupling
strengths of the mesopins.

The magnetic moment of a mesospin at $T=5\,\rm{K}$ was determined to be
$m_{\rm{0}}=M_{\rm{0}}V=6.5\times 10^{6}\mu_{\rm{B}}$, where V is the volume
of the magnetic material in the island.111The small difference between the
moment stated here compared to Andersson  _et al._ (2016) stems from a more
detailed analyses of the island dimensions from SEM imaging. Micromagnetic
calculations using $\rm{MuMax}^{3}$ Vansteenkiste  _et al._ (2014) revealed
non-collinearities of the moment within the islands. A reduced effective
moment of approximately $m_{\rm{0,eff}}=0.65\times m_{\rm{0}}=4.2\times
10^{6}\mu_{\rm{B}}$ at $T=0\,\rm{K}$, is used to compensate for the dynamic
non-collinear internal magnetic structure of the elements (see Andersson  _et
al._ (2016), Bessarab  _et al._ (2012), Bessarab  _et al._ (2013) and Gliga
_et al._ (2015)). Furthermore, the intrinsic effective moment $m_{\rm eff}(T)$
of the mesospins is temperature dependentKapaklis  _et al._ (2014); Andersson
_et al._ (2016); Östman  _et al._ (2018b):

$m_{\rm eff}(T)=m_{\rm 0,eff}(1-T/T_{0})^{0.5}$ (1)

Figure 1: Temperature dependence of the in-phase component of the AC
susceptibility $\chi^{\prime}$ (upper panel) and out-of-phase component
$\chi^{\prime\prime}$ (lower panel) for three square ASI arrays with different
inter-island coupling strengths. The inset depicts the geometry and a typical
SEM image of the nanostructures composed of $310\pm 15\,\rm{nm}\times 140\pm
15\,\rm{nm}$ magnetic islands. Only the gap $g$, differs between the three
investigated arrays, as the center-to-center distance $d$ is varied to be
$d\approx 380\,\rm{nm}$, $420\,\rm{nm}$, or $460\,\rm{nm}$ respectively. The
measurements were conducted with an AC magnetic field strength of
$\mu_{\rm{0}}H_{\rm{ac}}=0.1\,\rm{mT}$ and a frequency of $f=33\,\rm{Hz}$
aligned along the [10]-direction of the array. For direct comparison the
$\chi^{\prime}$ and $\chi^{\prime\prime}$ curves were normalized to the value
of the corresponding maximum of $\chi^{\prime}$ ($T_{\rm{m}}$, blue dots).

The dynamic response of the three extended square ASI arrays (2$\times$2 mm2
each) was investigated by AC susceptibility employing a magneto-optical Kerr
effect (MOKE) magnetometer in longitudinal modeAspelmeier  _et al._ (1995).
For this purpose the samples were mounted in a cryostat with optical access
($4\,\rm{K}\leq T\leq 300\,\rm{K}$). All experiments were performed with a 20
mW laser with a wavelength of 660 nm. A pair of Helmholtz coils was used to
generate a small sinusoidal magnetic field with a given frequency, aligned
along the [10]-direction of the nanostructured array (see inset of Fig. 1,
upper panel). For frequencies between $f=0.1-333\,\rm{Hz}$ an amplitude of
$\mu_{\rm{0}}H_{\rm{ac}}=0.1\,\rm{mT}$ was used, while for the higher
frequencies $f=1111\,\rm{Hz}$ and $f=3333\,\rm{Hz}$ the amplitude was reduced
to $\mu_{\rm{0}}H_{\rm{ac}}=0.01\,\rm{mT}$. A lock-in amplifier (Stanford
Research SR830) was used as a voltage source for generating the magnetic field
and measuring the AC susceptibility at the corresponding frequency. The sample
was shielded from the earth’s magnetic field by a double-wall mu-metal
cylinder, and demagnetization was performed prior to each cool-down.

Applying the magnetic field along the [10]-direction of the samples (see Fig.
1 ) in the longitudinal MOKE configuration probes primarily the mesospins with
long axes parallel to the direction of the magnetic field. Both in- and out-
of-phase components, $\chi^{\prime}$ and $\chi^{\prime\prime}$, were recorded
during warming up to the maximum temperature of $T=300\,\rm{K}$, starting
either at $T\approx 80\,\rm{K}$ ($f=33\,\rm{Hz}$) or $T=200\,\rm{K}$ (other
frequencies) and by stabilizing the cryostat in discrete $1\,\rm{K}$
temperature steps while keeping the amplitude and the frequency of the
magnetic field fixed. In every step, sufficient time (60 s) was taken to allow
the measurement system (cryostat) to stabilise before starting a measurement
222As temperature gradients between sample surface and temperature sensor of
the set-up cannot be fully excluded it is plausible to assume the presence of
temperature drifts of the sample of up to $1\,\rm{K}$.. The time for data
acquisition at each temperature step was increased for lower frequencies to
ensure an acceptable signal to noise ratio. In this respect it should be noted
that the magnetic signal of the samples is extremely small, as it originates
from a 2.2 ML magnetic Fe layer in Pd, with a coverage of about $\approx
37-54\%$ (depending on the geometry of the pattern). For the weakest coupled
sample ($d=460\,\rm{nm}$) this equals an effective magnetic coverage of a sub-
monolayer thick continuous Fe film in Pd.

##  III Results and discussion

Figure 2:  Temperature dependence of the in-phase component $\chi^{\prime}$ of
the AC susceptibility measured at different excitation field frequencies. The
magnetic susceptibility was recorded while heating the sample. Typical dataset
obtained on the array with the strongest inter island-interactions
($d=380\,\rm{nm}$) using excitation fields of
$\mu_{\rm{0}}H_{\rm{ac}}=0.1\,\rm{mT}$ ($f=0.1-333\,\rm{Hz}$) and
$\mu_{\rm{0}}H_{\rm{ac}}=0.01\,\rm{mT}$ ($f=1111-3333\,\rm{Hz}$) oriented
along [10]-direction. The curves are normalized to their respective magnetic
field amplitudes. The blue dots mark the position of the maximum $T_{\rm{m}}$.
The dotted blue line is a guide to the eye.

For each of the three studied arrays, a single peak in both $\chi^{\prime}$
and $\chi^{\prime\prime}$ is observed, as illustrated in Fig. 1. The shape of
the peaks is similar for the three arrays, while a shift in the peak position
$T_{\rm{m}}$ towards higher temperatures is observed with decreasing inter-
island distance, $i.e.$ increasing inter-island coupling strength. We
attribute the maximum in the AC susceptibility to the condition when the
average relaxation time is equal to the observation time window
$\tau_{\rm{m}}(T_{\rm{m}})=1/(2\pi f)$Souletie and Tholence (1985), where $f$
is the frequency of the applied magnetic field and $T_{\rm{m}}$ is the
corresponding temperature commonly referred to as the blocking temperature in
superparamagnetic samples. Consequently, by determining the peak positions
$T_{\rm{m}}$ for different observation time windows (frequencies), the average
relaxation time of the system as a function of temperature can be extracted
from the AC susceptibility data.

In order to investigate both the temperature dependence of the relaxation time
and the effect of inter-island interactions, $T_{\rm{m}}$ was determined for
different frequencies, effectively resulting in different observation time
windows. Typical results are illustrated for the array with the strongest
interactions ($d=380\,\rm{nm}$) in Fig. 2. $T_{\rm{m}}$ was determined by
fitting a parabola to a region of interest around the maximum of each curve.
For each of the three arrays, nine $T_{\rm{m}}$ values are extracted,
corresponding to the employed excitation frequencies $f$, as illustrated in
Fig. 2. The peaks are found to shift towards higher temperatures with
increasing frequency.

Figure 3: Compiled $T_{\rm{m}}$ datasets, as obtained from the in-phase
component $\chi^{\prime}$ at nine frequencies $f$ for the three arrays,
plotted against $ln(\tau_{\rm{m}})=ln(1/(2\pi f))$. The bold lines represent
fits based on the VFT law with a fixed $\tau_{\rm{0}}=1\cdot
10^{-11}\,\rm{s}$. The results of the fits are summarized in Table 1. Inset:
Systematic reduction of the extracted freezing Temperatures $T_{\rm{F}}$ with
increased island spacing i.e. reduced coupling strength.

###  III.1 Fitting the experimental data using the Vogel-Fulcher-Tammann law

We start the analysis by employing the empirical VFT law for the relaxation
time $\tau$ of weakly interacting magnetic particles Shtrikman and Wohlfarth
(1981), an approach which has previously been used to describe the relaxation
in artificial spin ice structures Andersson  _et al._ (2016); Morley  _et al._
(2017). Within this approach the relaxation time can be calculated by

$\tau=\tau_{\rm{0}}\cdot\exp{{\Big{(}\frac{E_{\rm{K}}}{k_{\rm{B}}(T-T_{\rm{F}})}\Big{)},}}$
(2)

where $\tau_{0}$, $E_{\rm{K}}$, $k_{\rm{B}}$, $T$ and $T_{\rm{F}}$ correspond
the inverse attempt frequency, the intrinsic energy barrier, the Boltzmann
constant, the temperature and the freezing or Fulcher temperature,
respectively. While the mesospin’s energy barrier $E_{\rm{K}}$ depends on the
shape and magnetisation of the mesospins, the Fulcher temperature $T_{\rm{F}}$
is indicative of the interaction strength of the elements. The intrinsic
energy barrier $E_{\rm{K}}$ in Eq. (1) is attributed to the shape anisotropy
and its temperature dependence can be captured by:

$E_{\rm{K}}=\mu_{\rm{0}}\Delta
N\frac{[m_{\rm{eff}}(T)]^{2}V}{2}=\mu_{\rm{0}}\Delta
N\frac{m_{\rm{0,eff}}^{2}V}{2}(1-\frac{T}{T_{0}}),$ (3)

where $\mu_{\rm{0}}$ is the vacuum magnetic permeability and $\Delta N$ is the
island differential demagnetizing factor according to Osborn (1945). Taking
into account the temperature dependence of the island saturation magnetization
$M_{\rm{S}}$ in the VFT law, the temperature dependence of the energy barrier
becomes $E_{\rm{K}}=E_{\rm{K}}^{\rm{0}}\cdot(1-T/T_{\rm{0}})$, where
$E_{\rm{K}}^{\rm{0}}=\mu_{\rm{0}}\Delta N\frac{m_{\rm{0,eff}}^{2}V}{2}$ is the
energy barrier at zero temperature.

Fig. 3 presents the extracted $T_{\rm{m}}$ values along with their
corresponding VFT fits. The fitting was performed by reversing Eq. (2), i.e.
solving for $T$, and using a Levenberg-Marquardt algorithm to find the best
fit for $E_{\rm{K}}^{\rm{0}}$ and $T_{\rm{F}}$. This procedure facilitates
weighting the data points with their respective uncertainties on the maxima
positions $T_{\rm{m}}$, while the uncertainty in the magnetic field frequency
given by the lock-in amplifier was considered negligible.333Note that, both
for the parabola as well as the nonlinear VFT-fits, the weighting is defined
by the inverse quadratic form of the individual errors for each data point,
i.e. $w_{\rm{i}}=1/\sigma_{\rm{i}}^{2}$. Furthermore, the standard errors of
the fitting parameters are scaled with the square root of the reduced chi
squared. The uncertainties of $T_{\rm{m}}$ were taken from the fits used for
determining the maximum positions, while the uncertainty in selecting the
individual temperature ranges used for peak finding were considered
negligible. This approach qualitatively captures the effect that the
uncertainty is larger for the highest and lowest frequencies measured. Since
all arrays are composed of the same size elements, the flipping time was
assumed to be constant, $\tau_{\rm{0}}=10^{-11}\,\rm{s}$, in accordance with a
previous relaxation study Andersson  _et al._ (2016). The same energy barrier,
$E_{\rm{K}}^{\rm{0}}$ was used for all three gap sizes studied. The summary of
the results from these analysis is found in Table 1.

For comparison, an independent magnetostatic estimation of the intrinsic
energy barrier was made, calculating $\Delta N$ using the Osborn methodology
Osborn (1945). The effective magnetic moment of an island $m_{\rm{0,eff}}$ was
used in order to take the non-collinearities of the moment within the islands
into account. By using this approach and the temperature dependence of
$m_{\rm{eff}}$, the mean value of the energy barrier was determined to be
$E_{\rm{K}}/k_{\rm{B}}\approx 3200\,\rm{K}$ at $250\,\rm{K}$. Using the
temperature-scaling for the fitted value of $E_{\rm{K}}^{\rm{0}}$, we obtain
$E_{\rm{K}}/k_{\rm{B}}(T=250\,\rm{K})=2510\,\rm{K}\pm 120\,\rm{K}$, a value
which, for the current case, is not so far from the magnetostatic estimation.

###  III.2 On the extraction of characteristic interaction energies

With a single energy barrier fitted for all three datasets, the impact of the
different gap sizes is inherently accounted for by the Fulcher temperature,
$T_{\rm{F}}$. This dependence is represented in the inset of Fig. 3,
highlighting a connection between $T_{\rm{F}}$ and the inter-island
magnetostatic interactions. This further raises the question of whether the
Fulcher temperature can be systematically and accurately used to extract the
characteristic interaction strength between mesospins.

In the framework of a weak-coupling regime in spin glasses, Shtrikman and
Wohlfarth (1981) offered a recipe for extracting the typical interaction
energy between the magnetic components by using the mean-field based formula:

$E_{\rm{i}}^{0}=\sqrt{k_{\rm{B}}\cdot T_{\rm{F}}\cdot E_{\rm{K}}^{0}}\\\ $ (4)

where $E_{\rm{K}}^{0}$ represents the average intrinsic energy barrier of a
single magnetic element, while $E_{\rm{i}}^{0}$ is the characteristic
interaction energy between the magnetic constituents, in turn defined as
$E_{\rm{i}}^{0}=\mu_{\rm{0}}m_{\rm 0,eff}\cdot H_{\rm{i}}^{0}$, with
$H_{\rm{i}}^{0}$ representing the characteristic mean interaction field. Note
that the latter quantity is determined by considering the ground state
manifold of the magnetic system, i.e. the configurations for which the mean
interaction field is the strongest.

We have extracted the characteristic interaction energies for the three
studied arrays, presented in Table 1, and compared them with the corresponding
ground-state energies of a square spin ice, i.e. type I tiling, given by
conventional interaction models based on the point-dipole approximation, the
dumbbell representation, as well as micromagnetic simulations 444The
micromagnetic simulations have been performed using the $\rm{MuMax}^{3}$ GPU-
accelerated micromagnetic simulation program Vansteenkiste  _et al._ (2014).
The simulations were performed using the following parameters:
$M=499289\,\rm{A/m}$, $A_{\rm{ex}}=3.36\cdot 10^{-12}\rm{J/m}$, length
$=310\,\rm{nm}$, width $=140\,\rm{nm}$, and a thickness of $2\,\rm{nm}$.. As
shown in Fig. 4, the VFT-law obtained values deviate considerably from the
corresponding estimations of all the models employed, with resulting energies
that can be more than 5 times smaller than expected from micromagnetic
estimates. Furthermore, although limited by the number of experimentally-
available points, the gap-dependent scaling of the extracted values appears to
be much less pronounced than the model predictions. This approach was
previously also employed by Morley  _et al._ (2017) \- attempting to extract
experimental values of the characteristic mesospin interaction energies - who
reported a significant underestimation of the extracted energies.

Table 1: Left: Energy barrier $E_{\rm{K}}^{\rm{0}}$ and Fulcher temperature $T_{\rm{F}}$ as extracted from the VFT fits of the experimental data plotted in Fig. 3, assuming a fixed $\tau_{\rm{0}}=1\cdot 10^{-11}\,\rm{s}$, as done in Andersson  _et al._ (2016) for the same samples. Right column: Characteristic magnetostatic interaction energies obtained from the VFT fits in the weak coupling formalism using Eq. (4). The given uncertainties correspond to the fitted parameter’s standard errors as provided from the weighted non-linear fitting procedure described in the text. $\rm{d}$ | $\rm{gap}$ | $E_{\rm{K}}^{\rm{0}}/k_{\rm{B}}$ | $T_{\rm{F}}$ |  | $E_{\rm{i}}^{\rm{0}}/k_{\rm{B}}$
---|---|---|---|---|---
(nm) | (nm) | (K) | (K) |  | (K)
380 | $70$ |  | $145\pm 6$ |  | $966\pm 30$
420 | $110$ | $6438\pm 308$ | $120\pm 6$ |  | $879\pm 32$
460 | $150$ |  | $105\pm 7$ |  | $822\pm 32$
Figure 4: (a) Comparison between the interaction energies per island
($E_{i}^{0}$) extracted for the three experimental gap sizes (black dots) and
the type I tilling ground-state interaction energies per island computed for
various inverse gap sizes ($1/g$) using micromagnetics (blue diamonds), the
dumbbell model (red line) and the point-dipole model (green line) for the
given sample geometries. The shaded areas account for the error originating
from a gap uncertainty of $\pm 15\,\rm{nm}$, extracted from the analysis of
the SEM images. (b) Normalized temperature dependence of the normalized mean
interaction field for various thermal dynamics ranges, calculated using the
Fulcher-based expression (dotted lines) and the statistical approach
(continuous lines). The black lines correspond to the case with temperature-
independent energy barriers and island moments, while the coloured lines
consider the temperature scaling for different normalized Curie temperatures
$T_{0}/T_{F}$. The vertical dashed lines mark the boundaries of each
considered range. The inset gives the normalized temperature dependence of the
ratio between the Fulcher-based and statistical-based interaction fields. The
experimental values provided in Table 1 for $E_{\rm{K}}^{\rm{0}}$ and
$T_{\rm{F}}$ are used in both figures.

While several factors can contribute to the mismatch between the experimental
values and the various models, we mainly attribute these discrepancies to the
incompatibility between the frameworks for Eq. (4) and artificial spin
systems. The basic assumption of the formalism is that the VFT law can be
treated as an interacting Arrhenius law, i.e. an Arrhenius-like expression in
which the intrinsic energy barrier is further biased by a local interaction
field:

$\ln{\Big{(}\frac{\tau}{\tau_{\rm{0}}}\Big{)}}=\frac{E_{\rm{K}}}{k_{\rm{B}}(T-T_{\rm{F}})}=\frac{E_{\rm{K}}+\mu_{\rm{0}}m_{\rm
0,eff}\cdot H^{\rm{F}}_{\rm{i}}(T)}{k_{\rm{B}}T}$ (5)

where $H^{\rm{F}}_{\rm{i}}(T)$ represents the temperature-dependent mean
interaction field and $E_{\rm{K}}$ is generally assumed to be a temperature-
independent quantity. This equality imposes a certain temperature dependence
on this variable, parameterized by $E_{\rm{K}}$ and $T_{\rm{F}}$:

$H^{\rm{F}}_{\rm{i}}(T)=\frac{1}{\mu_{\rm{0}}m_{\rm
0,eff}}\cdot\frac{E_{\rm{K}}\cdot T_{\rm{F}}}{T-T_{\rm{F}}}$ (6)

Notice that this interaction field presents a divergence around the Fulcher
temperature, a physically unrealistic feature. Taking now an at-equilibrium
statistical perspective, the mean interaction field should converge at low
temperatures towards a finite value, $H^{\rm{0}}_{\rm{i}}$, corresponding to
the ground-state manifold. We shall further consider the choice of Shtrikman
and Wohlfarth (1981) assuming a hyperbolic tangent behavior, akin to the case
of a paramagnetic system in an externally applied field:

$H^{\rm{S}}_{\rm{i}}(T)=H^{\rm{0}}_{\rm{i}}\cdot\tanh{\Big{(}\frac{\mu_{\rm{0}}m_{\rm
0,eff}\cdot H^{\rm{0}}_{\rm{i}}}{k_{\rm{B}}T}\Big{)}}$ (7)

This brings us to the remaining conditions for the validity of Eq. (4):

  * •

$T\gg T_{\rm{F}}$, i.e. the temperature should be sufficiently far away from
the Fulcher temperature.

  * •

$k_{\rm{B}}T\gg E^{\rm{0}}_{\rm{i}}$, i.e. the energy associated to the
thermal bath should be much larger than the characteristic interaction energy.

It should also be noted that the intrinsic energy barrier is also assumed to
be much higher than both the thermal bath and the interaction energy, i.e.
$E_{K}\gg k_{\rm{B}}T$ and $E_{K}\gg E^{0}_{i}$. With these in place, the
functions describing the two interaction fields, $H^{\rm{F}}_{\rm{i}}(T)$ and
$H^{\rm{S}}_{\rm{i}}(T)$, become compatible to a first order expansion. This
validity region is illustrated in Fig. 4(b). Notice how the two interaction
fields present almost identical temperature scaling once the temperature is
about two orders of magnitude higher than $T_{\rm{F}}$.

If we now consider the temperature dependence of the intrinsic energy barrier
and the magnetic moment, a key feature of mesoscopic spin systems, the
possibility of finding a compatibility region for the two interaction fields
is severely limited. Fig. 4(b) illustrates the impact of having a range in
thermal dynamics bound by the Fulcher and Curie temperatures. Here, the energy
barrier, $E_{\rm{K}}$, from Eq. (5) is replaced with its temperature dependent
form given by Eq. (3), while the characteristic interaction field from Eq.
(7), $H_{\rm{i}}^{0}$, is similarly replaced with a temperature dependent
expression, $H_{\rm{i}}^{0}\cdot(1-T/T_{0})^{0.5}$, thus accounting for the
scaling of the inter-island couplings. As it can be seen, only for a thermal
range spanning several orders of magnitude between $T_{0}$ and $T_{\rm{F}}$,
can one achieve a compatibility region that accommodates the aforementioned
formalism. This is particularly highlighted by the inset of Fig. 4(b), where
the ratio between the two interaction field expressions is plotted as a
function of temperature. The red lines correspond to the most weakly
interacting sample, with an average gap of 150 nm and characterized by a
$T_{0}/T_{\rm{F}}\cong 4$ ratio. Even for this scenario, there is no clear
overlap between the two mean interaction fields within the experimental
temperature window, which should therefore compromise the matching with the
interaction models considered in Fig. 4(a).

A microscopic description of the phenomenological parameter $T_{\rm{F}}$ was
provided in a study of the AC susceptibility of weakly interacting magnetic
nanoparticles by Vernay  _et al._ (2014). The VFT law, assuming the case of
both weak dipolar interactions and surface anisotropy for the magnetic
nanoparticles, can be transformed into a semi-analytical expression, linking
$T_{\rm{F}}$ to deviations from uniaxial anisotropy and dipole-dipole
interactionsVernay  _et al._ (2014). While this analysis assumes weak
interactions being valid also in our case, the modelling of the interactions
with discrete point dipoles is an imprecise description of our spatially
extensive thermal mesospins (see Fig. 4(b)). On the other hand this
approachVernay  _et al._ (2014) can serve as a stimulation for the development
of a new and revised formalism capable of capturing in detail the collective
temporalAndersson  _et al._ (2016); Östman  _et al._ (2018b) and thermal
dynamicsPohlit  _et al._ (2015, 2016); Kapaklis  _et al._ (2014); Östman  _et
al._ (2018b) of artificial spin ice.

##  IV Conclusions

We studied the AC susceptibility of thermally active square ASI arrays of
varying interaction strength. The freezing of the mesospin dynamics was
measured over a wide range of observation times and a systematic dependence of
the freezing temperature on the inter-island coupling strength was found.
Extracting magnetostatic energies from the frequency dependence using the VFT
law, which was recently applied to ASIAndersson  _et al._ (2016); Morley  _et
al._ (2017), revealed significant discrepancies of the obtained interaction
energies compared to theoretical estimates, similar to the case of Morley  _et
al._ (2017). Besides experimental uncertainties, we attribute this to the
violation of the requirements of weak coupling to extract energies in the VFT
formalism, along with the inability to obtain measurements at temperatures far
above the freezing transition, while staying well below the material’s Curie
temperature. We note that these requirements are generally difficult to meet
in thermally active mesospin systems. Therefore a more advanced model enabling
the extraction of microscopic variables and accounting for the details of our
mesospin systems, such as the internal magnetic structure, temperature
dependence of the energy barriers and interaction energies, is highly
desirable and will be developed in future works.

Nonetheless, AC susceptibility using a longitudinal MOKE setup was proven to
be a simple yet powerful technique for studying magnetization dynamics of
thermally active nanostructures. Its high sensitivity allows investigations of
minute changes in the mesospin dynamics in arrays with an intrinsically low
magnetic moment. The method lends itself to temperature and frequency
dependent studies on a laboratory scaleTopping and Blundell (2018), which
facilitates an evaluation of well-established models, or scaling laws for the
description of mesospin systems. The method can further serve as an excellent
tool for the characterization of thermal dynamics, collective behaviour Gliga
_et al._ (2013); Ciuciulkaite  _et al._ (2019) and thermodynamic phase
transitions in magnetic metamaterials Heyderman and Stamps (2013); Leo  _et
al._ (2018), such as ASIs Anghinolfi  _et al._ (2015); Östman  _et al._
(2018b).

The data that support this study are available via the Zenodo repository
Pohlit  _et al._ (2019).

###### Acknowledgements.

The authors thank Prof. Per Nordblad (Department of Engineering Sciences,
Uppsala University) for fruitful and stimulating discussions. The authors
acknowledge support from the Swedish Research Council, the Swedish Foundation
for International Cooperation in Research and Higher Education (STINT) project
KO2016-6889 and the Knut and Alice Wallenberg Foundation project “Harnessing
light and spins through plasmons at the nanoscale” (2015.0060). This research
used resources of the Center for Functional Nanomaterials, which is a U.S. DOE
Office of Science Facility, at Brookhaven National Laboratory under Contract
No. DE-SC0012704.

## References

  * Östman  _et al._ (2018a) E. Östman, U. B. Arnalds, V. Kapaklis, A. Taroni, and B. Hjörvarsson, Journal of Physics: Condensed Matter 30, 365301 (2018a).
  * 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 (London) 439, 303 (2006).
  * Nisoli  _et al._ (2013) C. Nisoli, R. Moessner, and P. Schiffer, Rev. Mod. Phys. 85, 1473 (2013).
  * Heyderman and Stamps (2013) L. J. Heyderman and R. L. Stamps, Journal of Physics: Condensed Matter 25, 363201 (2013).
  * Rougemaille and Canals (2019) N. Rougemaille and B. Canals, Eur. Phys. J. B 92, 62 (2019).
  * Kapaklis  _et al._ (2014) V. Kapaklis, U. B. Arnalds, A. Farhan, R. V. Chopdekar, A. Balan, A. Scholl, L. J. Heyderman, and B. Hjörvarsson, Nat. Nanotechnol. 9, 514–519 (2014).
  * Perrin  _et al._ (2016) Y. Perrin, B. Canals, and N. Rougemaille, Nature 540, 410 (2016).
  * Östman  _et al._ (2018b) E. Östman, H. Stopfel, I.-A. Chioar, U. B. Arnalds, A. Stein, V. Kapaklis, and B. Hjörvarsson, Nature Physics 14, 375 (2018b).
  * Nisoli  _et al._ (2017) C. Nisoli, V. Kapaklis, and P. Schiffer, Nature Physics 13, 200 (2017).
  * Pohlit  _et al._ (2015) M. Pohlit, F. Porrati, M. Huth, Y. Ohno, H. Ohno, and J. Müller, Journal of Applied Physics 117, 17C746 (2015).
  * Pohlit  _et al._ (2016) M. Pohlit, I. Stockem, F. Porrati, M. Huth, C. Schröder, and J. Müller, Journal of Applied Physics 120, 142103 (2016).
  * Kapaklis  _et al._ (2012) V. Kapaklis, U. B. Arnalds, A. Harman-Clarke, E. T. Papaioannou, M. Karimipour, P. Korelis, A. Taroni, P. C. W. Holdsworth, S. T. Bramwell, and B. Hjörvarsson, New Journal of Physics 14, 035009 (2012).
  * Morgan  _et al._ (2010) J. P. Morgan, A. Stein, S. Langridge, and C. H. Marrows, Nature Physics 7, 75 (2010).
  * Farhan  _et al._ (2013) A. Farhan, P. M. Derlet, A. Kleibert, A. Balan, R. V. Chopdekar, M. Wyss, L. Anghinolfi, F. Nolting, and L. J. Heyderman, Nature Physics 9, 375 (2013).
  * Porro  _et al._ (2013) J. M. Porro, A. Bedoya-Pinto, A. Berger, and P. Vavassori, New Journal of Physics 15, 055012 (2013).
  * Andersson  _et al._ (2016) M. S. Andersson, S. D. Pappas, H. Stopfel, E. Östman, A. Stein, P. Nordblad, R. Mathieu, B. Hjörvarsson, and V. Kapaklis, Scientific Reports 6, 37097 (2016).
  * Morley  _et al._ (2017) S. A. Morley, D. A. Venero, J. M. Porro, S. T. Riley, A. Stein, P. Steadman, R. L. Stamps, S. Langridge, and C. H. Marrows, Physical Review B 95, 104422 (2017).
  * Sendetskyi  _et al._ (2019) O. Sendetskyi, V. Scagnoli, N. Leo, L. Anghinolfi, A. Alberca, J. Lüning, U. Staub, P. M. Derlet, and L. J. Heyderman, Physical Review B 99, 214430 (2019).
  * Saccone  _et al._ (2019) M. Saccone, A. Scholl, S. Velten, S. Dhuey, K. Hofhuis, C. Wuth, Y.-L. Huang, Z. Chen, R. V. Chopdekar, and A. Farhan, Physical Review B 99, 224403 (2019).
  * Anghinolfi  _et al._ (2015) L. Anghinolfi, H. Luetkens, J. Perron, M. G. Flokstra, O. Sendetskyi, A. Suter, T. Prokscha, P. M. Derlet, S. L. Lee, and L. J. Heyderman, Nat. Commun. 6, 8278 (2015).
  * Leo  _et al._ (2018) N. Leo, S. Holenstein, D. Schildknecht, O. Sendetskyi, H. Luetkens, P. M. Derlet, V. Scagnoli, D. Lançon, J. R. L. Mardegan, T. Prokscha, A. Suter, Z. Salman, S. Lee, and L. J. Heyderman, Nature Communications 9, 2850 (2018).
  * Bedanta and Kleemann (2008) S. Bedanta and W. Kleemann, Journal of Physics D: Applied Physics 42, 013001 (2008).
  * Topping and Blundell (2018) C. V. Topping and S. J. Blundell, Journal of Physics: Condensed Matter 31, 013001 (2018).
  * Landi (2013) G. T. Landi, Journal of Applied Physics 113, 163908 (2013).
  * Vernay  _et al._ (2014) F. Vernay, Z. Sabsabi, and H. Kachkachi, Phys. Rev. B 90, 094416 (2014).
  * Pärnaste  _et al._ (2007) M. Pärnaste, M. Marcellini, E. Holmström, N. Bock, J. Fransson, O. Eriksson, and B. Hjörvarsson, Journal of Physics: Condensed Matter 19, 246213 (2007).
  * Hase  _et al._ (2014) T. P. A. Hase, M. S. Brewer, U. B. Arnalds, M. Ahlberg, V. Kapaklis, M. Björck, L. Bouchenoire, P. Thompson, D. Haskel, Y. Choi, J. Lang, C. Sánchez-Hanke, and B. Hjörvarsson, Physical Review B 90, 104403 (2014).
  * Note (1) The small difference between the moment stated here compared to Andersson  _et al._ (2016) stems from a more detailed analyses of the island dimensions from SEM imaging.
  * Vansteenkiste  _et al._ (2014) A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. V. Waeyenberge, AIP Advances 4, 107133 (2014).
  * Bessarab  _et al._ (2012) P. F. Bessarab, V. M. Uzdin, and H. Jónsson, Physical Review B 85, 184409 (2012).
  * Bessarab  _et al._ (2013) P. F. Bessarab, V. M. Uzdin, and H. Jónsson, Physical Review Letters 110, 020604 (2013).
  * Gliga  _et al._ (2015) S. Gliga, A. Kákay, L. J. Heyderman, R. Hertel, and O. G. Heinonen, Physical Review B 92, 060413 (2015).
  * Aspelmeier  _et al._ (1995) A. Aspelmeier, M. Tischer, M. Farle, M. Russo, K. Baberschke, and D. Arvanitis, Journal of Magnetism and Magnetic Materials 146, 256 (1995).
  * Note (2) As temperature gradients between sample surface and temperature sensor of the set-up cannot be fully excluded it is plausible to assume the presence of temperature drifts of the sample of up to $1\tmspace+{.1667em}\rm{K}$.
  * Souletie and Tholence (1985) J. Souletie and J. L. Tholence, Physical Review B 32, 516 (1985).
  * Shtrikman and Wohlfarth (1981) S. Shtrikman and E. Wohlfarth, Physics Letters A 85, 467 (1981).
  * Osborn (1945) J. A. Osborn, Phys. Rev. 67, 351 (1945).
  * Note (3) Note that, both for the parabola as well as the nonlinear VFT-fits, the weighting is defined by the inverse quadratic form of the individual errors for each data point, i.e. $w_{\rm{i}}=1/\sigma_{\rm{i}}^{2}$. Furthermore, the standard errors of the fitting parameters are scaled with the square root of the reduced chi squared.
  * Note (4) The micromagnetic simulations have been performed using the $\rm{MuMax}^{3}$ GPU-accelerated micromagnetic simulation program Vansteenkiste  _et al._ (2014). The simulations were performed using the following parameters: $M=499289\tmspace+{.1667em}\rm{A/m}$, $A_{\rm{ex}}=3.36\cdot 10^{-12}\rm{J/m}$, length $=310\tmspace+{.1667em}\rm{nm}$, width $=140\tmspace+{.1667em}\rm{nm}$, and a thickness of $2\tmspace+{.1667em}\rm{nm}$.
  * Gliga  _et al._ (2013) S. Gliga, A. Kákay, R. Hertel, and O. G. Heinonen, Physical Review Letters 110, 117205 (2013).
  * Ciuciulkaite  _et al._ (2019) A. Ciuciulkaite, E. Östman, R. Brucas, A. Kumar, M. A. Verschuuren, P. Svedlindh, B. Hjörvarsson, and V. Kapaklis, Phys. Rev. B 99, 184415 (2019).
  * Pohlit  _et al._ (2019) M. Pohlit, G. Muscas, I.-A. Chioar, H. Stopfel, A. Ciuciulkaite, E. Östman, S. Pappas, A. Stein, B. Hjörvarsson, P. E. Jönsson, and V. Kapaklis, “Collective magnetic dynamics in artificial spin-ice probed by AC susceptibility,”  (2019), Zenodo, https://doi.org/10.5281/zenodo.3607135.