Article

Emergent order in hydrodynamic spinlattices
https://doi.org/10.1038/s41586-021-03682-1
Received: 6 July 2020
Accepted: 1 June 2021
Published online: 4 August 2021
Check for updates

Pedro J. Sáenz1,2 ✉, Giuseppe Pucci2,3, Sam E. Turton2, Alexis Goujon2,4, Rodolfo R. Rosales2,
Jörn Dunkel2 & John W. M. Bush2 ✉Macroscale analogues1–3 of microscopic spin systems offer direct insights into
fundamental physical principles, thereby advancing our understanding ofsynchronization phenomena4 and informing the design of novel classes of chiralmetamaterials5–7. Here we introduce hydrodynamic spin lattices (HSLs) of ‘walking’droplets as a class of active spin systems with particle–wave coupling. HSLs reveal
various non-equilibrium symmetry-breaking phenomena, including transitions from
antiferromagnetic to ferromagnetic order that can be controlled by varying the latticegeometry and system rotation8. Theoretical predictions based on a generalized
Kuramoto model4 derived from first principles rationalize our experimental
observations, establishing HSLs as a versatile platform for exploring active phaseoscillator dynamics. The tunability of HSLs suggests exciting directions for future
research, from active spin–wave dynamics to hydrodynamic analogue computationand droplet-based topological insulators.

Broken-mirror symmetry is an essential feature of many physical andbiological systems, from magnetic metals and composites9 to DNA10.
The most fundamental manifestation of a broken left–right symmetry
are the spins of elementary particles such as photons or electrons11.
Microscopic spin systems and their macroscopic analogues1,2 are
attracting increasing interest12,13 as building blocks of information
storage14 and metamaterials3,15 devices. Over the last decade, experimental realizations of ultracold fermionic lattice gases16,17 have spurredadvances in the understanding of charge and spin transport, promising novel technological platforms for quantum computation18,19 andspintronics20. In parallel, recent progress in the control of natural andartificial active matter has led to the discovery of emergent collective
spin states in chemical21,22, biological23–25 and synthetic26 suspensions,inspiring the invention of chiral active metamaterials7. Here we demonstrate a class of HSL systems27 that combines active self-propulsion
and particle-wave dynamics.
An HSL consists of an array of in-phase bouncing droplets28, eachdrop confined by a submerged circular well and orbiting on the surface
of a vertically vibrated liquid bath (Fig. 1a, b; Supplementary Video 1).
Despite being macroscopic classical objects, these ‘walking’ droplets display wave-mediated self and pair interactions reminiscent of
those arising in microscopic systems29. To demonstrate the tunability of
HSLs, we induce a range of collective ‘magnetic’ ordering phenomena
through variations in the vibrational forcing and lattice geometry, andby mimicking the effect of an applied magnetic field through imposition of system rotation8. We explain the experimentally observed
particle–wave dynamics by deriving a mapping onto a reduced oscillator network description4,30, generalizing models used to describe
neuronal networks31 and other synchronization phenomena32,33. This
mapping establishes HSLs as a versatile platform for investigatingchiral symmetry breaking and synchronization phenomena far fromequilibrium.
In our experiments, both the droplets and underlying liquid bathconsist of the same 20 centistokes (cSt) silicon oil (1 cSt = 1 mm2 s−1;
Methods section ‘Experiments’). An electromagnetic shaker vibrates
the bath vertically with forcing acceleration Γ(t) = γcos(2πft), where γ
is the maximum acceleration, f the oscillation frequency, and t time
(Fig. 1b). For sufficiently strong acceleration, the air layer between thedrop and bath prevents direct contact, thus enabling persistent bouncing34. When the drop bounces at half the driving frequency—and soachieves resonance with the bath’s most unstable wave mode—it excites
radially decaying circular waves of characteristic wavelength λF on thebath surface34–36. The amplitude and spatio-temporal decay of thesewaves can be controlled through the driving35. When the driving amplitude γ is increased beyond a critical walking threshold γw, vertical bouncing becomes unstable and transitions to horizontal walking29,34. Aswalking droplets land on the slopes of the underlying wave field η, theyexperience a wave-induced horizontal force Fi (xi, t) ∝ − ∇η|x = xi proportional to the wave slope ∇η at the point of impact xi(t). Free droplets
thus walk at constant speed along straight lines when γ > γw (ref. 34).
To transform linear motion into circular motion, we confined eachwalker to a submerged circular well37 (Fig. 1b, Methods section ‘Experiments’). The resulting clockwise or anticlockwise motions were characterized by measuring the specific spin S(t) = Lz(t)/m, where m is thedroplet’s mass and Lz denotes its angular momentum with respect to thewell centre. Spin is thus understood in classical terms as corresponding to orbital motion. We first quantified how the spin of an isolated
walker depends on the driving strength by progressively increasing theforcing acceleration γ towards the Faraday threshold γF, at which thebath surface becomes unstable in the absence of droplets35,36. These1
Department of Mathematics, University of North Carolina, Chapel Hill, NC, USA. 2Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA, USA. 3Institut dePhysique de Rennes, CNRS, Université de Rennes, UMR 6251, Rennes, France. 4School of Engineering, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland. ✉e-mail: saenz@unc.edu;

bush@math.mit.edu

58 | Nature | Vol 596 | 5 August 2021

R
3

c

Ii

Spin, S (mm2 s–1)

1
Submerged wells

20
15

a

ACW
1
Spin, SdWell number

1 cm

0

10

–1

5

CW

Spin correlation, F
Magnetization, M

e
M

A

0

h
–1

Phase difference,ΔI (rad)

f

1

F

π

ΔI – (counter-rotation)
ΔI+ (co-rotation)

π_
2
0
0

0.5

1.0
Time, t (h)

1.5

2.0

Average spin correlation, < F>

2

Fig. 1 | Spontaneous antiferromagnetic order in a 1D
HSL. a, b, Oblique view (a) and schematic cross section(b) of a ring (of radius R) of submerged circular wells(of diameter D, depth h, centre separation L), eachDhcontaining a single drop executing clockwise orH
L
anticlockwise motion on the surface of a vibrating fluidbath. A thin fluid layer (of depth H) between wells enables
30
5wave-mediated interactions between neighbouring
droplets (Supplementary Video 1). c, Dependence of the
4
spin S and orbital radius r of an isolated walker (D = 14 mm,
20
H = 1 mm, f = 80 Hz) on the dimensionless bathS3
acceleration γ/γF. The solid lines are fits resulting fromsmoothing the piecewise linear curve obtained by
102connecting the data points. Once γ exceeds a criticalJcrthreshold γc, the drop’s circular trajectory destabilizes
00
into a precessing trefoil, increasing spin variation and
78
82
86
90
susceptibility to perturbation-induced spin flips. For aBath acceleration, J/JF (%)
given driving acceleration, only one state is observed.
g
Insets, two typical drop trajectories, coloured by theirinstantaneous spin value, at low (top left) and high(bottom right) accelerations. Centrelines denote
2
interpolated mean values and envelopes indicate
1
FA = –0.71
standard deviations. d–f, Temporal evolution of the spin
20
values Si (d), spin pair correlation χ and magnetization M19(e), and phase difference Δϕ (f) for a 1D periodic HSL with
N = 20 wells driven at γ/γF = 82.0%. g, Snapshot showing
collective antiferromagnetic order at instant A in e
(Supplementary Video 2). χA is the spin−spin correlation χ
0.50
Experiments
at instant A in e. Dashed lines highlight domains of local
GK model0.25
antiferromagnetic order. White circles mark dropJc
positions, background colour shows the spin magnitude.
0h, Dependence of the average spin correlation ⟨χ⟩ on γ/γF.
Solid lines are fits resulting from smoothing the piece­wise–0.25
linear curve obtained by connecting the data points as inc. The average spin correlation ⟨χ⟩ indicates maximal–0.50
78 82 86 90 94 antiferromagnetic ordering for γ ≈ γ , high­lighting thatcBath acceleration, J/JF (%)
stable single-spin trajectories promote collective
ordering. Error bars denote the standard deviations(Methods section ‘Statistics’). Experimental parameters
in d–h: L = 17.7 mm, D = 14 mm, H = 1 mm, N = 20, f = 80 Hz,
γF = 4.780g (where g is the gravitational acceleration),
λF = 4.75 mm. ACW, anticlockwise; CW, clockwise.
Wave-mediatedinteraction

Γ (t) = J cos(2πft)bRadius, r (mm)asingle-spin measurements revealed an optimal signal-to-driving ratio
for S(t) at a vibrational acceleration γc ≈ 0.85γF (Fig. 1c). For suboptimal
driving γ < γc, the mean spin decreases with the walker’s characteristic speed as γ decreases36. For γ > γc, the circular orbits deform into
precessing trefoils (Fig. 1c) that exhibit a larger spin variance, making
them more susceptible to spin flips when perturbed. Armed with thisunderstanding of the single-spin states, we investigated the collective
spin dynamics in one-dimensional (1D) and two-dimensional (2D) HSLs.
The wave-mediated spin–spin coupling in HSLs is reminiscent of
spatially oscillating interactions9,38 in Ruderman–Kittel–Kasuya–
Yosida (RKKY)-type spin models. Our experimental setup allowed us
to tune the magnitude of the spin–spin coupling by varying the driving acceleration γ and the depth H of the fluid bath between adjacentwells (Fig. 1b; Methods section ‘Experiments’). When the pair couplingis sufficiently strong, nearest-neighbour interactions may cause spin
flips. To determine whether such flips can facilitate coherent collective
dynamics across the lattice, we measured the normalized effective‘magnetization’ M(t) = ∑iSi(t) and spin–spin correlation χ(t) = ∑i~jSi(t)Sj(t),where ∑i~j denotes a sum over adjacent pairs. To achieve statisticalsignificance, experiments were run for several hours: in 1 h a droplet
performs approximately 105 bounces and approximately 1,800 orbits(Methods section ‘Statistics’). Positive values of χ signal parallel alignment of neighbouring spins (‘ferromagnetic’ order) whereas negativevalues of χ indicate antiparallel alignment (‘antiferromagnetic’ order).

The oscillatory wave-mediated spin–spin coupling suggests that
HSLs can support different types of collective order depending on theratio between lattice spacing L and Faraday wavelength λF. To test thishypothesis, we studied a 1D periodic (circular) HSL with N = 20 equally
spaced wells and L/λF = 3.7 (Fig. 1a, d–h). Starting from a random initial
spin configuration (Methods section ‘Statistics’), we observed that pairinteractions can trigger multiple spin flips (Fig. 1d), leading to fluctuations in the instantaneous magnetization and spin–spin correlation(Fig. 1e). The magnetization vanishing on average, ⟨M⟩ ≈ 0, indicatedthat global mirror symmetry was preserved (Fig. 1e). However, thenegative pair correlation, ⟨χ⟩ < 0, revealed a bias towards local antiferromagnetic order (Fig. 1g; Supplementary Video 2). Holding L/λF
fixed, a similar antiferromagnetic bias occurred in experiments with
non-periodic boundary conditions, different N values, and different
lattice radii R (Methods section ‘Statistics’), confirming that antiferromagnetic ordering is selected by the lattice spacing. The strengthof the emergent antiferromagnetic order depends non-monotonicallyon the driving amplitude γ (Fig. 1h; Methods section ‘Statistics’). The
strongest collective antiferromagnetic response was observed forγ ≈ γc, indicating a correlation between the global collective orderingand the robustness of a single spin state (Fig. 1c).
The emergence of local magnetic ordering can be understood byanalysing the drops’ horizontal phase ϕi, which plays a role analogous
to the polar angle in XY-type spin models (Fig. 1a). To quantify pairwise
Nature | Vol 596 | 5 August 2021 | 59

Article
ACW

25
20
15
10
5

Spin correlation, F
Magnetization, M
Phase difference,ΔI (rad)e0Spin, S1
J0(kFx – xj(Ij))

–1
CW

b

c

d

1

Fi,j

j

FA = 0.81

r

χ

Ii

Ij

Zj

0

i
L

Zi

A
M

–1

f

π

ΔI – (counter-rotation)
ΔI+ (co-rotation)

π_
2
0
0

0.5

1.0

1.5
Time, t (h)

2.0

2.5

FM+

AFM –

FM –

AFM+

Wave height, K

Well number

a

13

3.0

14

15
16
17
Drop–drop mean distance (mm)18

Fig. 2 | Tuning collective order through the lattice parameters inexperiments and theory. When the lattice spacing is reduced to L = 13.2 mm,
the 1D HSL spontaneously settles into ferromagnetic order. Parameters:
D = 10 mm, H = 0.8 mm, N = 28, f = 78 Hz, γ/γF = 86.0%, γF = 5.280g, λF = 4.84 mm.
a–c, The preferred local ferromagnetic order is reflected in the instantaneous
spins S i (a), spin–spin correlation χ and magnetization M (b), and the mean
phase differences Δϕ± (c). d, Spin configuration at instant A in b showing fourdomains of ferromagnetic order (Supplementary Video 3). e, Geometry
considered in the theoretical wave model. Neighbouring drops i and j movealong circular trajectories of radius r at angular speed ±ω. The drop position isdetermined by the angular phase ϕ, defined with respect to the latticecentreline. Drop i is subject to a force Fij along its direction of motion, owing to
the wave field generated by drop j, as is approximated by a Bessel function J0
centred at its position. f, Depending on the drop separation, four dominant
drop-pair synchronization modes are favoured, leading to the emergence ofalternating ferromagnetic and antiferromagnetic local order as the latticespacing is varied (Methods section ‘Simulations’). ACW, anticlockwise; CW,
clockwise.

phase synchronization, we recorded the mean phase differences Δϕ± for
co-rotating (+) and counter-rotating (−) neighbouring pairs, respectively(see Methods section ‘Statistics’). These observables signal in-phase(Δϕ± → 0) or out-of-phase (Δϕ± → π) rotations of neighbouring walkers,
and so also indicate their relative distance. For antiferromagnetic HSLs,
we found a substantial bias towards in-phase rotations (Fig. 1f). In particular, both spin–spin correlation (Fig. 1h) and phase synchronization
are maximized at driving strength γ ≈ γc, implying a causal link between
coherent orbital synchronization and emergent spin order.
To demonstrate how the collective spin dynamics depends on latticegeometry, we performed 1D experiments for a reduced lattice spacing
L/λF = 2.8 (Fig. 2a–d). The resulting changes in the wave-mediated spin
interactions (Fig. 2e) support extended domains of ferromagnetic
ordering, as reflected by a positive pair correlation, ⟨χ⟩ > 0 (Fig. 2b, d;
Supplementary Video 3), with a bias towards in-phase rotation (Fig. 2c).
As in the antiferromagnetic HSL with L/λF = 3.7, we observed no global
magnetization, ⟨M⟩ ≈ 0, indicating that changes in lattice geometry
can be used to control local magnetic order but do not lead to global
mirror-symmetry breaking.
To rationalize the manner in which geometry dictates collective spinorder and phase synchronization in HSLs, we derived a generic phaseoscillator model from a detailed hydrodynamic description of ourexperimental system (Fig. 2e; Methods section ‘Theory’). In this model,
we focus on subcritical driving γ ≤ γc (Fig. 1c), so that each walker i follows a circular trajectory xi(t) = r0(cosϕi(t), sinϕi(t)), with angular frequency ωi = ϕi̇, where (r0, ϕi) are the cylindrical coordinates with respect
to the centre of each well. The walking speed v0 then defines a preferredorbital frequency ±ω0 = ±v0/r0 (ref. 36). Modelling this effect by a nonlinear Rayleigh-type friction, the phase dynamics of a walker is
described by

Equation (1) can be interpreted as an active XY-type model9, in which
perturbations to ω0 decay over the timescale τ, and Fij is the wavemediated force exerted on drop i by its nearest neighbour j (Fig. 2e).
For weak accelerations and small orbital radii r0 ≪ L, the interaction
force Fij(ϕi, ϕj) = −∂U/∂ϕi can be derived from the effective coupling11
potential U (ϕi , ϕj ) = 2 α(cos ϕi − cos ϕj ) 2 + 2 β(sin ϕi − sin ϕj ) 2 (Methods
section ‘Theory’). The coefficients are proportional to the curvature,α ∝ J ″0 (kF L), and gradient, β ∝ J ′0 (kF L), of the coupling wave field, which
is well approximated by a Bessel function of the first kind J0(kF|x – xj|)
centred at the neighbouring drop position. Here prime symbols denote
derivatives with respect to rj = |x – xj|, and kF = 2π/λF the Faraday wavenumber. HSLs characterized by purely circular orbits, as arise for γ ≤ γc,thus belong to a generalized Kuramoto universality class4.
Analysis of the phase-oscillator dynamics (equation (1)) confirmsthat the spin–spin correlation alternates between ferromagnetic andantiferromagnetic order as the lattice spacing L is varied (Extended
Data Fig. 4a). The associated mean phase differences (Extended Data
Fig. 4b) revealed four different modes of pairwise symmetry breaking.
Specifically, there are two ferromagnetic phases FM± and two antiferromagnetic phases AFM±, distinguished by preferential in-phase (+)
and out-of-phase (−) rotation, respectively (Fig. 2f). Depending on L,
the force parameters α and β can be either positive or negative, yielding four magnetic phases corresponding to minima of the coupling
potential. The model predictions agree well with our experimental
data: the antiferromagnetic (Fig. 1d–f; L = 17.7 mm) and ferromagnetic(Fig. 2a–c; L = 13.2 mm) HSL experiments fall within the predicted AFM+
and FM+ ranges, respectively. The theory also captures how the preferred order varies with vibrational acceleration (Fig. 1h, Extended Data
Fig. 5).
The spins in an antiferromagnetic material may be realigned into aferromagnetic state by imposing a constant magnetic field B (ref. 9).
We next sought to determine whether HSLs can undergo similar global
symmetry breaking. The Lorentz force FB = q(v × B) acting on a charge q

1  ω2 ωi̇ = 1 − i2  ωi + ∑ Fij(ϕi , ϕj ).
τ  ω0 
j
60 | Nature | Vol 596 | 5 August 2021

(1)

Spin correlation, F
Magnetization, M
Phase difference,ΔI (rad)010

FA = –0.79

–1

5

CW
1

M

1

Ω

0

0
A

B

C

F

–1

–1

Experiments
GK model

0.6

1.2
1.0

0.3
0

0.8

–0.3

Ωc

–0.6

0.6

0
0.25 0.50 0.75
Bath angular speed,Ω (rad s–1)

FB = 0.92f

Co-rotation
Counter-rotation

00.25 0.50 0.75 1.00
Bath angular speed,Ω (rad s–1)

g

X/X0

1.2

π

ΔI – (counter-rotation)
ΔI+ (co-rotation)

Fw

1.1

FΩ

FC = 0.99

π_
2
0
0

h

0.9

Radius, r/r0eSpin, S

Well number

15

b

c

d

1

Average spin correlation, < F>

ACW

20

Bath angular speed,Ω (rad s–1)

a

Zi
v

0.25

0.50

0.75
Time, t (h)

1.00

1.25

v

FΩ
Fw

Zi

Ω

1.0
0.9
0.8

1.50

Fig. 3 | Inducing global polarization through applied rotation. Rotating thecircular HSL from Fig. 1d–h induces a transition from effective antiferromagnetic
to ferromagnetic order. a–c, The experiment was run for 1 h without rotation
before the recording process started. During the first 0.5 h of data acquisition,
the experiment continued without rotation (Ω = 0) with the HSL exhibiting local
antiferromagnetic order with no net magnetization. At t = 0.5 h, the bath rotation
was switched on (Ω = 0.7 rad s−1, anticlockwise) inducing a transition to
ferromagnetic order and positive magnetization. At t = 1 h, the rotation direction
was reversed (Ω = −0.7 rad s−1, clockwise), leading to negative magnetization.
Lattice parameters and forcing frequency are the same as in Fig. 1 and γ ≈ γc.
d, Instantaneous spin configurations at time points A–C in b, showing
antiferromagnetic and ferromagnetic ordering (Supplementary Video 4).

a

Ω

c

χA, χB and χC are the spin−spin correlations at instants A, B and C. e, Average spincorrelation ⟨χ⟩ signals a transition to ferromagnetic order at Ωc ≈ 0.22 rad s−1(Methods section ‘Statistics’). f, g, Force balance and measured trajectory,colour-coded by speed, for an isolated spin in a frame co-rotating with the bath.
The drop’s trajectory is nearly circular when its angular velocity ω and the bathrotation Ω are aligned, but approximates a precessing trefoil when the signs areopposite. h, Mean orbital radius r versus bath angular speed Ω for a walker
co-rotating and counter-rotating with respect to the bath, as in f and g,
respectively; r0 is the trajectory radius for Ω = 0 shown in Fig. 1c. Solid lines are fits
resulting from smoothing the piecewise linear curve obtained by connecting thedata points as in Fig. 1. The centre line denotes the mean radius, and the envelope
indicates the standard deviation. ACW, anticlockwise; CW, clockwise.

d

Ω =0

Ω = 0.75 rad s–1

1 cm

FB = 0.91

Spin correlation, F

b

e
1

B

0

–1
0

Ω (rad s–1)
0
0.75
0.5

A
1.0
Time, t (h)Fig. 4 | Rotation-induced polarization transition from antiferromagnetic
to ferromagnetic order in 2D square lattices. a, Top view showing the wave
field of a square lattice with L = 17.3 mm, D = 13 mm, H = 0.96 mm, f = 80 Hz,
γ/γF = 86% and γF ≈ 4.940g. b, Evolution in time of the spin correlation χ reveals
preferred antiferromagnetic order for a stationary lattice (Ω = 0), and

1.5

2.0

Average spin correlation <F>

FA = –0.70
0.50
0.25
0
–0.25

0

0.25 0.50 0.75 1.00
Bath angular speed,Ω (rad s–1)

ferromagnetic order for a rotating lattice (Ω = 0.75 rad s−1). c, d, Instantaneous
spin arrangements at time points A and B in b, showing the emergence of
antiferromagnetic (c), and ferromagnetic (d) domains, respectively(Supplementary Video 5). e, A transition from antiferromagnetic toferromagnetic order is prompted by increasing the bath’s angular speed.

Nature | Vol 596 | 5 August 2021 | 61

Article
moving with velocity v has precisely the same form as the Coriolis forceFΩ = m(v × 2Ω) acting on a mass m moving with velocity v in a frame rotating at constant angular velocity Ω (ref. 8) . This suggests the possibilityof transforming antiferromagnetic HSLs into ferromagnetic states byrotating the vibrating bath about a vertical axis at a rate Ω (Methods
section ‘Experiments’). To test this hypothesis, we initialized the antiferromagnetic HSL described in Fig. 1 without bath rotation (Ω = 0) and let
it equilibrate for 1 h. Keeping Ω = 0 for another 0.5 h, the HSL remained in
the antiferromagnetic state as expected (Fig. 3a, b; t < 0.5 h). At t = 0.5 h,
we began to rotate the setup anticlockwise by gradually increasing therotation rate to Ω = 0.7 rad s−1. In a matter of seconds, most spins had
co-aligned with the bath rotation, leading to a strongly positive spincorrelation, ⟨χ⟩ ≈ 0.58. However, unlike the steady HSL in Fig. 2a–d, therotated HSL had a net magnetization, ⟨M⟩ ≈ 0.78. Moreover, when therotation direction was reversed to Ω = −0.7 rad s−1 (clockwise) at t = 1 h,
most of the spins flipped to realign with the inverted ‘magnetic’ field.
The effective magnetization thus became negative, ⟨M⟩ < 0, while thespin pairs remained ferromagnetically arranged, ⟨χ⟩ > 0 (Fig. 3b; Supplementary Video 4). The transition to ferromagnetic order requires
a supercritical rotation rate |Ω| > Ωc ≈ 0.22 rad s−1 (Fig. 3e). Rotation
also affects the pairwise phase synchronization: as the spin dynamics
become dominated by the Coriolis force, the phase differences become
uncorrelated, thus averaging to the mean ⟨Δϕ±⟩ ≈ π/2 (Fig. 3c).
Examining the single-spin dynamics in the rotating frame reveals themechanism responsible for the field-induced polarization (Fig. 3f, g).
For optimal driving, γ ≈ γc, a walker executes circular motion in theabsence of rotation (Fig. 1c). When rotation is switched on, a drop
co-orbiting with the bath experiences a Coriolis force FΩ that opposes
the confining force of the submerged well Fw (ref. 37) (Fig. 3f). As a result,the drop preserves its circular trajectory but its orbital radius increases
with Ω (Fig. 3h). Conversely, when the walker is counter-rotating withrespect to the bath, both Fw and FΩ point inwards (Fig. 3g), destabilizing
the circular orbits into trefoil-like trajectories similar to those observedfor Ω = 0 and γ > γc (Fig. 1c). Droplets on these trefoil-like trajectories
pass closer to the well centres (Fig. 3h), where their walking directions
can be more easily reversed by perturbations from their neighbours.
Hence, confinement favours co-rotating orbits in HSLs; conversely,
free walkers prefer counter-rotating inertial orbits8. We note that ageneralization of our phase-oscillator model captures the effects ofthe external field (Fig. 3e; Methods section ‘Theory’).
Two-dimensional classical and quantum spin lattices display featuresthat are absent in their 1D counterparts, including geometric frustrationand topological order9. HSLs provide a promising platform for exploring such effects at the macroscale. For example, square HSLs (Fig. 4a;
Supplementary Video 5) that promote antiferromagnetic order in theabsence of rotation (Ω = 0; Fig. 4b, c), undergo a polarization transition(Ω > 0; Fig. 4b, d) as the Coriolis force is increased (Fig. 4e). Additionally, simulations of larger lattices confirm the emergence of collective
magnetic order in square lattices (Extended Data Fig. 6).
To conclude, since lattice geometry (Extended Data Fig. 7) and theinteractions between lattice sites can be directly tailored in various ways
(Methods section ‘Tunability’), we expect HSLs to provide a versatile
testbed for designing wave-coupled autonomous metamaterials39,
exploring active spin-wave dynamics13, and studying Anderson-typewave localization phenomena40 far from equilibrium.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions

62 | Nature | Vol 596 | 5 August 2021

and competing interests; and statements of data and code availability
are available at https://doi.org/10.1038/s41586-021-03682-1.
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.
30.
31.
32.
33.
34.
35.
36.
37.
38.
39.
40.

Kane, C. L. & Lubensky, T. C. Topological boundary modes in isostatic lattices. Nat. Phys.
10, 39–45 (2014).
Chen, B. G.-g., Upadhyaya, N. & Vitelli, V. Nonlinear conduction via solitons in a
topological mechanical insulator. Proc. Natl Acad. Sci. USA 111, 13004–13009 (2014).
Nash, L. M. et al. Topological mechanics of gyroscopic metamaterials. Proc. Natl Acad.
Sci. USA 112, 14495–14500 (2015).
Strogatz, S. H. From Kuramoto to Crawford: exploring the onset of synchronization in
populations of coupled oscillators. Physica D 143, 1–20 (2000).
Rechtsman, M. C. et al. Optimized structures for photonic quasicrystals. Phys. Rev. Lett.
101, 073902 (2008).
van Zuiden, B. C. et al. Spatiotemporal order and emergent edge currents in active
spinner materials. Proc. Natl Acad. Sci. USA 113, 12919–12924 (2016).
Soni, V. et al. The odd free surface flows of a colloidal chiral fluid. Nat. Phys. 15, 1188–1194(2019).
Fort, E. et al. Path-memory induced quantization of classical orbits. Proc. Natl Acad. Sci.
USA 107, 17515–17520 (2010).
Coey, J. M. D. Magnetism and Magnetic Materials (Cambridge Univ. Press, 2010).
Watson, J. D. & Crick, F. H. C. The structure of DNA. Cold Spring Harb. Symp. Quant. Biol.
18, 123–131 (1953).
Weinberg, S. The Quantum Theory of Fields Vol. 2 (Cambridge Univ. Press, 1996).
Baltz, V. et al. Antiferromagnetic spintronics. Rev. Mod. Phys. 90, 015005 (2018).
Chumak, A. V. et al. Magnon spintronics. Nat. Phys. 11, 453–461 (2015).
Maurer, P. C. et al. Room-temperature quantum bit memory exceeding one second.
Science 336, 1283–1286 (2012).
Süsstrunk, R. & Huber, S. D. Observation of phononic helical edge states in a mechanicaltopological insulator. Science 349, 47–50 (2015).
Schneider, U. et al. Metallic and insulating phases of repulsively interacting fermions in a
3D optical lattice. Science 322, 1520–1525 (2008).
Cheuk, L. W. et al. Quantum-gas microscope for fermionic atoms. Phys. Rev. Lett. 114,
193001 (2015).
Buluta, I., Asshab, S. & Nori, F. Natural and artificial atoms for quantum computation.
Rep. Prog. Phys. 74, 104401 (2011).
Gross, C. & Bloch, I. Quantum simulations with ultracold atoms in optical lattices. Science
357, 995–1001 (2017).
Dieny, B. et al. Opportunities and challenges for spintronics in the microelectronics
industry. Nat. Electron. 3, 446–459 (2020).
Biktashev, V. N., Barkley, D. & Biktasheva, I. V. Orbital motion of spiral waves in excitable
media. Phys. Rev. Lett. 104, 058302 (2010).
Totz, J. F. et al. Spiral wave chimera states in large populations of coupled chemical
oscillators. Nat. Phys. 14, 282–285 (2018).
Wioland, H. et al. Ferromagnetic and antiferromagnetic order in bacterial vortex lattices.
Nat. Phys. 12, 341–345 (2016).
Nishiguchi, D. et al. Engineering bacterial vortex lattice via direct laser lithography.
Nat. Commun. 9, 4486 (2018).
Riedel, I. H., Kruse, K. & Howard, J. A self-organized vortex array of hydrodynamicallyentrained sperm cells. Science 309, 300–303 (2005).
Souslov, A. et al. Topological sound in active-liquid metamaterials. Nat. Phys. 13,
1091–1094 (2017).
Sáenz, P. J. et al. Spin lattices of walking droplets. Phys. Rev. Fluids 3, 100508 (2018).
Couder, Y. et al. Dynamical phenomena: walking and orbiting droplets. Nature 437, 208(2005).
Bush, J. W. M. Pilot-wave hydrodynamics. Ann. Rev. Fluid Mech. 47, 269–292 (2015).
Lindner, B. et al. Effects of noise in excitable systems. Phys. Rep. 392, 321–424 (2004).
Cumin, D. & Unsworth, C. P. Generalising the Kuramoto model for the study of neuronalsynchronization in the brain. Physica D 226, 181–196 (2007).
Matheny, M. H. et al. Exotic states in a simple network of nanoelectromechanicaloscillators. Science 363, eaav7932 (2019).
Dörfler, F. & Bullo, F. Synchronization in complex networks of phase oscillators: a survey.
Automatica 50, 1539–1564 (2014).
Protière, S., Boudaoud, A. & Couder, Y. Particle-wave association on a fluid interface.
J. Fluid Mech. 554, 85–108 (2006).
Eddi, A. et al. Information stored in Faraday waves: the origin of a path memory. J. Fluid
Mech. 674, 433–463 (2011).
Moláček, J. & Bush, J. W. M. Drops walking on a vibrating bath: towards a hydrodynamic
pilot-wave theory. J. Fluid Mech. 727, 612–647 (2013).
Sáenz, P. J., Cristea-Platon, T. & Bush, J. W. M. A hydrodynamic analog of Friedel
oscillations. Sci. Adv. 6, eaay9234 (2020).
Ruderman, M. A. & Kittel, C. Indirect exchange coupling of nuclear magnetic moments byconduction electrons. Phys. Rev. 96, 99–102 (1954).
Ronellenfitsch, H. et al. Inverse design of discrete mechanical metamaterials. Phys. Rev.
Mater. 3, 095201 (2019).
Lagendijk, A., Van Tiggelen, B. & Wiersma, D. S. Fifty years of Anderson localization.
Phys. Today 62, 24 (2009).

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.
© The Author(s), under exclusive licence to Springer Nature Limited 2021

Methods
Introduction
Walking droplets in free space self-propel in a rectilinear fashion28,34–36.
Hydrodynamic spin states41, in which a droplet executes circular orbits
confined by its own wave field are generally thought to be unstable in thelaboratory42,43, unless confined by an external force such as a linear springforce44 or the Coriolis force that acts in a rotating frame8,45–47. Recent
experiments have demonstrated that variable bottom topography may
be exploited to realize new hydrodynamic quantum analogues37,48; inparticular, the influence of a submerged well has led to analogues of theso-called quantum mirage48 and Friedel oscillations37. In this work, we
exploit variable bottom topography to stabilize a collection of interacting hydrodynamic spin states27, and characterize the emergent collective
order. We note that, although several works have considered the stabilityof stationary bouncing droplet lattices49–53, and other works54,55 examinedwave-induced correlations of walker pairs confined to adjacent 1D cavities, we present here a macroscopic analogue of electronic spin systems
with large collections of orbiting droplets.
Experiments
Stationary bath. A schematic of the experimental setup is presented
in Extended Data Fig. 1. A circular aluminium bath was filled with silicon
oil with density ρ = 950 kg m−3, viscosity v = 20.9 cSt, and surface tension
σ = 20.6 mN m−1. The lattices of submerged wells were constructed bylaser-cutting acrylic sheets of thickness h = 6.2 mm, which were thenbolted to the base of the bath. The well depth h was thus sufficiently
large with respect to the Faraday wavelength, λF < 5 mm, for the liquid
inside the wells to be in the deep-water wave regime56. Similarly, theminimum distance d between the bath’s outer wall and the outermost
well was sufficiently large, d > 6λF, to ensure that the meniscus waves57
that originated at the bath wall had no effect on the experiments. We
Hnote that the Faraday wavelengths corresponding to the shallow λFhand deep λF regions are only marginally different, specificallyh H
λF /λF ≈ 1.05. For the sake of simplicity of notation, we thus simply refer
to both wavelengths as λF.
The bath was mounted on an optical table (Newport SG-34-4 custom
breadboard, 3.0′ × 4.0′ × 4.32″) and vibrated vertically by an electromagnetic shaker (Data Physics, V55) with an external amplifier (DataPhysics, PA300E) at acceleration Γ(t) = γcos(ωt), where γ and f = ω/2π
are the maximum acceleration and frequency, respectively. The shaker
is connected to the bath by a thin steel rod coupled with a linear air
bearing (PI L.P., 4 × 4″ cross-section, 6.5″ long hollow bar) that ensuresa spatially uniform vibration to within 0.1%58. The forcing is monitored
through a data-acquisition system (NI, USB-6343) with two piezoelectric accelerometers (PCB, 352C65), attached to the base plate on
opposite sides of the drive shaft, and a closed-loop feedback ensuresa constant acceleration amplitude to within ±0.002g (ref. 58); where gis the acceleration due to gravity. We define the Faraday threshold γF
as the critical vibrational acceleration above which Faraday waves59
appear above the wells. The experimental setup was left vibrating atγ ≈ γF for at least 1 h before running any experiment to ensure a steadystate in the bath. The Faraday threshold was measured before and after
each experimental run in order to ensure the constancy of the prescribed driving, γ/γF, which may drift slowly owing to weak variations
in viscosity and surface tension resulting from ambient temperaturechanges. With this protocol, the standard day-to-day variation in γF
was limited to ±0.04g, and the variation in each experimental run was
limited to |ΔγF|/γF < 0.003, where ΔγF represents the variation in the
Faraday threshold over a 1-h segment.
Droplets of the same silicon oil with diameter Dd = 0.75 ± 0.01 mm
were generated with a piezoelectric droplet-on-demand generator60,
and placed manually above the submerged wells. The bath was illuminated with a LED light ring in order to increase the contrast between
the drops and the black background. The motion of the droplets wasrecorded at 19 frames per second with a CCD camera mounted directly
above the bath, and tracked with an in-house particle-tracking algorithm written in Matlab. The form of the wave field was captured byrecording the normal reflection of light at the free surface57. To that
end, a semi-reflective mirror at 45° was placed between the CCD camera
and the bath, and the light ring was replaced by a diffuse-light lamp
facing the mirror horizontally. One then observes images with brightregions corresponding to horizontal parts of the surface, extrema or
saddle points57. To ensure that ambient air currents did not affect theresults, the vibrating bath was enclosed within a transparent acrylicchamber. Experimental results obtained using the acrylic chamber werecompared to others in which the spin system was shielded by a transparent acrylic lid bolted directly to the top of the bath. We found that theeffective noise in the system was slightly reduced with the lid, but thatthe emergent collective behaviour was statistically indistinguishable.
The experiments were performed with walking drops bouncing
subharmonically36 in resonance with the underlying Faraday wave
field, each with the same vertical bouncing phase. To ensure in-phasebouncing, before each experiment the bath acceleration was reduced
to γ = 1.1g, so that the droplets transitioned to harmonic bouncers36,
synchronized with the bath vibration. By observing the light reflectedfrom the upper side of the droplets, which intensifies when the droplets
deform during impact, we were able to deduce the vertical phase of all
the drops instantaneously. We confirmed that the drops continued
to bounce in phase when the driving acceleration was progressively
increased beyond γ = 1.1g, and throughout the duration of the experiments. Switching in the vertical bouncing phase was observed only for
driving accelerations, γ/γF > 90%, beyond the acceleration range where
coherent collective behaviour emerged (Fig. 1h).
Rotating bath. To subject hydrodynamic spin lattices to bulk rotation,
we mounted the aluminium bath on a vertical shaft attached to theupper plate of the linear air bearing through a system of preloadedball bearings. The freely rotating shaft was connected with a vibration
damping coupling to a d.c. electric motor with a planetary gearbox(Phidgets, 12 V d.c. motor, 139:1 gearbox) housed inside the hollow air
bearing (see Extended Data Fig. 1). The rotation rate was monitoredwith an optical rotary encoder (Phidgets, HKT22) and the electric motor controlled with an in-house Python programme, yielding a rotation
rate Ω that was constant to within ±2%.
The centrifugal force on the rotating bath induces a parabolic deformation of the fluid interface of the form

H0(r ) = H −Ω 2 Rb2 Ω 2 2
+
r ,
4g
2g

(2)

where r is the radial coordinate with respect to the centre of the bath,
H is the depth of the fluid at rest, and Rb the bath radius, which in therotating experiments was Rb = 88 mm (Extended Data Fig. 1). To compensate for the parabolic shape of the fluid surface, the liquid volume
was readjusted in experiments with periodic 1D lattices (Fig. 3) to ensurethat H0(R) = 1 mm, where R is the radius of the annular ring (Fig. 1a). The
resulting variation in depth across the wells did not exceed 40 μm, andso had a negligible influence on the dynamics. While such an adjustment
was not possible in experiments with 2D lattices, we note that for the
relatively slow rotation rates required to induce the transition from
antiferromagnetic to ferromagnetic order (Ω < 0.7 rad s−1, Fig. 4e), theliquid height was reduced in the centre and increased in the peripheryby <10% with respect to the depth of fluid at rest.
Lattice parameters. The lattice parameters for the 1D antiferromagnetic lattice discussed in Fig. 1 are L = 17.7 mm, D = 14 mm, H = 1 mm,
N = 20, f = 80 Hz, γF = 4.780g, λF = 4.75 mm. The lattice parameters of
the 1D ferromagnetic lattice shown in Fig. 2 are L = 13.2 mm, D = 10 mm,
H = 0.8 mm, N = 28, f = 78 Hz, γF = 5.280g, λF = 4.84 mm. We note that

Article
for the 1D ferromagnetic lattice, we reduced the well diameter D andliquid height H to compensate for the heightened pair coupling resulting from locating the wells closer to each other. Such modification helped reduce the effective noise by both increasing the wave
damping in the shallow region, and stabilizing the droplets’ circular
trajectories. We note that the vibration frequency, which effectively sets
the characteristic wavenumber kF, was also fine-tuned to accentuate
the ferromagnetic response. The lattice parameters for the 2D squarelattice shown in Fig. 4 are L = 17.3 mm, D = 13 mm, H = 0.96 mm, N = 37,
f = 80 Hz, γF = 4.940g, λF = 4.75 mm. Experiments were also performed
with lattices with size N = 7 × 7.

Statistics
Order parameters. To characterize the collective behaviour in HSLs,
we define three instantaneous order parameters, specifically the magnetization M(t), spin–spin correlation χ(t), and average phase differenceΔϕ(t). The magnetization, which serves to quantify the global symmetry
breaking, is defined as
M (t ) =∑i Si (t),∑i |Si (t)|(3)

where ∑i denotes the sum over all N spins, and Si(t) = Lz,i(t)/m is the
specific spin of droplet i with mass m and angular momentum Lz,i withrespect to the well centre. The denominator is the normalizing factor
such that M ∈ [−1, 1]. Positive (negative) values of M(t) signal a larger
number of walkers rotating anticlockwise (clockwise). Pairwise symmetry breaking is quantified through the normalized spin–spin correlation,χ (t ) =∑i ~j Si (t)Sj (t)∑i ~j |Si (t)Sj (t)|,(4)

where ∑i~j denotes sums over pairs {i, j} of adjacent walkers. χ ∈ [−1, 1]
thus signals a preference for co-rotation (χ > 0, ferromagnetic order)
or counter-rotation (χ < 0, antiferromagnetic order) among neighbouring walkers.
To quantify the relationship between the emergent collective behaviour and the pairwise orbital phase synchronization of adjacent walkers, we compute separately the mean phase difference Δϕ± ∈ [0, π] forco-rotating (+) and counter-rotating (−) walkers as

Δϕ±(t) =1∑ |ϕ (t) ∓ ϕj (t)|,N p± i ~ j i(5)±where ϕi is the instantaneous horizontal phase with respect to the lattice centreline (Fig. 1a), and ∑i j± denotes the sum over N p± adjacent~co-rotating (Np+ and i~j+) or counter-rotating (Np– and i~j–) pairs of walkers. The phase difference thus allows us to identify if there is preference
for in-phase (Δϕ → 0) or out-of-phase (Δϕ → π) rotation among adjacent
walkers. We note that incoherent orbital synchronization corresponds
to Δϕ → π/2.
Data-acquisition protocol. To guarantee robust statistics, we develop
a systematic data-acquisition protocol. The goal is to repeat multiple
times a sufficiently long experiment to ensure that (i) the initial spin
distribution is random, and (ii) the spin system reaches a statisticallysteady state. To guarantee random initial conditions, we manuallystirred with a needle all the drops at the beginning of each run. By
checking that the spin correlation started from 0 at the beginning ofeach experimental run, we confirmed that our stirring procedure wasadequate. We also checked that the stirring did not disrupt the verticalsynchronization of the walkers.
To determine the minimum time required for each experimental
run, we performed preliminary tests to estimate the characteristic

autocorrelation time τi for the spin arrangement, which is computed
through the autocorrelation function,

R(t) = S(t = 0) ⋅ S(t)/N ,(6)

as the first time for which R(τi) = 0. Here, S(t) denotes the vector, withcomponents Si(t), that characterizes the lattice spin arrangement attime t. The autocorrelation time provides a measure of the characteristic time beyond which the collective order at t + τi is statistically independent of the collective order at time t. Since τi fluctuates dependingon the experimental realization, we define the length of the ‘statisticallyindependent interval’ asτI = ⟨τi ⟩ + σ (τi),(7)

where ⟨τi⟩ and σ(τi) denote, respectively, the average and standard
deviation resulting from 10 different realizations of a given lattice configuration. Our experiments show that τI is inversely related to the rate
at which spins switch direction (higher switching rates lead to shorter
τI), which in turn is a non-monotonic function of the forcing acceleration
γ/γF. We thus recomputed τI as in equation (7) for each driving forcing
in the range of interest (Fig. 1h). The largest τI,max ≈ 0.5 h (lowest switching rate) was observed at γ ≈ γc, the driving acceleration for which thesignal-to-driving ratio for S(t) is maximal (Fig. 1c), and the emergent
collective behaviour strongest (Fig. 1h). Based on this observation, we
selected 1 h (~1,800 circulations for each drop) as the recording lengthfor our experimental realizations, thus guaranteeing that the system
has more than enough time to reach a steady state. This typical run
duration also yielded an optimal frequency to measure γF before andafter each realization.
Following this protocol, we recorded as many 1-h videos as needed
to obtain a minimum of NI ≥ 8 statistically independent intervals foreach γ/γF level. In some cases the number of independent intervals was
as high as NI ≈ 50. For each independent interval we then computed its
associated spin–spin correlation χj(t). The data reported in Fig. 1h correspond to the resulting average, ⟨∑χj(t)/NI⟩, and standard deviation,σ(∑χj(t)/NI), of these time series. The same data-acquisition protocol
was followed for all lattices discussed in the main text, including those
involving system rotation. We note that exceptions were made to obtain
Figs. 1d, 2a, 3a, 4b, for which we conducted longer experiments to better illustrate the emergent collective phenomena.
Boundary conditions, ring curvature R, and lattice size N. Experiments were performed to determine the effects of the boundary conditions, lattice size N, and ring curvature R in 1D HSLs (Fig. 1a). To assess
the effect of periodic versus non-periodic boundary conditions, weconsidered the antiferromagnetic configuration with N = 20 shown inthe main text (Fig. 1), and repeated the experiment at γ/γF = 86% afterplugging one of the submerged wells. We note that in the resulting
non-periodic lattice with N = 19, the two end walkers were separated
by a relatively large distance (2L > 7λF) by the thin layer (H = 1 mm) overwhich the waves dissipated quickly. Hence, the walkers at the ends did
not interact with each other, as was checked by direct visualizationof the wave field. Comparison of our experimental results revealed⟨χ(N = 19)⟩/⟨χ(N = 20)⟩ ≈ 1.05, which is well within the expected standard
deviation, and so demonstrated a negligible effect of the boundary
conditions. We also note that for 1D periodic lattices tuned for preferred
antiferromagnetic order, it is important whether N is even or odd. An
odd number of spins may lead to geometrical frustration effects9 thathinder the preferred collective order. To avoid such effects—the investigation of which is left for future studies—we selected even N for all
experiments reported in the main text.
To assess the effect of lattice size, we compared experiments with
N = 10 and N = 20, while keeping both L and R constant. Both sets ofexperiments led to preferred antiferromagnetic order, thus supporting

our inference that the emergent order is not a finite-size effect. We also
compared these results with those of a 1D straight lattice (R = ∞) with
N = 8, where the same preferred order emerged. We may thus eliminate
the curvature of the circular lattice from consideration as an importantfactor for the emergent collective behaviour. Finally, we also checked
for finite-size effects in the 2D square lattices: experiments with N = 37
and N = 49 resulted in the same collective order.

TheoryRayleigh–Bessel oscillator model. To rationalize the emergent orderin HSLs, we extend the modelling framework for pilot-wave dynamics
introduced previously61. In the weak-acceleration limit, |˙|v TM /|v| ≪ 1,
these authors show that the trajectory of a walking droplet, x(t), may
be described by the equationd(m v) + Dw v = F ,
dt w

(8)

where v = ẋ is the droplet velocity, F an externally imposed force acting
on the drop, Dw(|v|) the nonlinear drag, and mw the effective walker
mass. In the low-speed limit, in which the wave decay time, TM, is muchsmaller than the time taken for the droplet to travel one Faraday wavelength, TM ≪ λF/u0 (where u0 is the free walking speed of a droplet), one
may further simplify equation (8) to a Rayleigh oscillator-type equation62, similar to that introduced in a previous work63, where the droplet trajectory is described bymw

 |v|2 dv= DR 1 − 2  v + F .
dt
u0 


(9)

Here, DR denotes the magnitude of the nonlinear Rayleigh-typedrag, which drives the droplet towards its free walking speed, u0. The
dependence of the model parameters mw, u0 and DR on the physical variables is summarized in Extended Data Table 1. Equation (9) adequately
describes the motion of weakly accelerating, relatively low-speed droplets. We note that both asymptotic limits are satisfied by orbiting droplets in the regime γ ≤ γc, where TM ≈ 0.1 s. Although these limits are notstrictly satisfied in the relatively short periods of time during whichspin flips occur, our analysis will demonstrate that details of the spin
flips are not essential to rationalize the emergent collective behaviour.
We consider N droplets self-propelling above an array of submerged circular wells with centre-to-centre separation L. Equation (9) determines the trajectory of each droplet, which is definedby a two-component displacement vector xi(t), where i = 1, …, N. The
force acting on each droplet may be decomposed into two components,
F = Fw + Fi,j, where Fw represents the effective confining force exerted
by the submerged well, and Fi,j the wave force exerted on droplet i by
its nearest neighbouring droplet j. Motivated by experiments64, we
assume that Fw is a central force that points towards the well’s centre.
As described in the main text, in the relatively low-memory regime(γ ≤ γc), the drops remain on circular trajectories, with a fixed orbital
radius r. Focusing on finite 1D lattices for the sake of simplicity, we may
thus express the droplet trajectory and velocity vectors as xi(t)=(r cos ϕi (t) + (i − 1)L, r sin ϕi (t)) = r ˆ
ri + (i − 1)Lx
ˆ , and x
˙i(t) = rϕ˙i (− sin ϕi (t),cos ϕi (t)) = rϕ˙i θ^i , respectively. Here, rˆi and θ̂i denote the radial and azimuthal unit vectors for the ith spin with respect to the centre of its well,and x̂ the unit vector along the lattice centre line (Fig. 2e). Writing ωi = ϕi̇and substituting xi(t) and ẋ i(t)into equation (9) yields the following equations for the normal and tangential force balances:− mw rωi2 = F ⋅ rˆ,i

mw r

 r 2ω 2 dωi
ˆi .
− DR 1 − 2 i  rωi = Fi , j ⋅ θ
dt
u0 


(10)

(11)

We note that Fw only appears in the radial force balance, where it
acts to select the orbital radius r. Since the experiments suggest thatemergent order in HSLs is driven primarily by synchronization of thedroplets’ relative orbital phase (Figs. 1f, 2c), we restrict our attention
to the azimuthal force balance (equation (11)) and simply infer thevalue of r from our experiments. Thus, to rationalize the preferredmean phase differences, Δϕ±, in terms of the droplet coupling force,
Fi,j, we may write equation (11) as a nonlinear differential equation for
the evolution of each droplet’s angular frequency, ωi(t),ˆiFi , j ⋅ θ
1  ω2 ωi̇ = 1 − i2  ωi +,mw r
τ  ω0 (12)where ω0 = u0/r is the preferred angular frequency for a given driving
acceleration γ, as would arise for a single drop orbiting in isolation,and τ = mw/DR is the relaxation timescale (see Extended Data Table 1).
To complete our spin model, we must define the wave-mediatedforce exerted by droplet j on droplet i, Fi,j. The stroboscopic model ofref. 42 predicts that the wave field generated by the jth walker may beapproximated by

ηj (x, t) =A t
∫ J (k |x − xj(s)|)e−(t − s)/TMds ,TF −∞ 0 F

(13)

where A, TF and TM are defined in Extended Data Table 1. Thus, the coupling force may be expressed as

Fi , j =F t J1 (kF|xi(t) − xj(s)|)(xi(t) − xj(s))e−(t − s)/TMds.
∫
TF −∞ |xi(t) − xj(s)|(14)

This expression was derived for droplets unperturbed by topography and so does not incorporate the additional wave damping
expected in the shallow inter-well regions65. Nevertheless, our visualization of the wave field produced by a single droplet in the latticeyielded a waveform notably similar in form to that expected on thebasis of equation (14) (Methods section ‘Wave coupling’). We also
note that there is a relatively fast-propagating wavefront generated
at impact35, which decays rapidly with distance and so has a negligible
influence on droplet pair interactions when the drops are separated
by more than approximately λF (refs. 66,67). Guided by previous studies
and experimental evidence, we thus neglect this capillary wavefrontand consider that the liquid layer between wells has a negligible influence of the form of the drops’ wave fields and so seek a simplificationof equation (14) that yields insight into the mechanism responsible
for the emergent order.
In the parameter regime of interest, γ ≤ γc, the wave memory time,
TM < 0.1 s, is much smaller than the orbital period, 2π/ω0 ≈ 2 s. We may
thus simplify equation (14) in the limit of TM ≪ 2π/ω0 to obtain

Fi , j ≈(xi(t) − xj(t))
FTM
J (k |x (t) − xj(t)|).
TF 1 F i|xi(t) − xj(t)|(15)

Substituting equation (15) into equation (12) yields a system of ODEs
for each droplet interacting with its nearest neighbours in a spin lattice(xi − xj)
1  ω2 ωi̇ = 1 − i2  ωi + ℱ ∑ J1 (kF|xi − xj|)
⋅ (− sin ϕi , cos ϕi ) . (16)τ  ω0 
|
−
|
x
x


i
j
j 

Here, ℱ is the interaction force coefficient, and ∑j denotes the sum
over nearest neighbours j = i ± 1. Owing to the form of the coupling
force, we refer to this as the ‘Bessel model’, which predicts the emergent
spin order and phase synchronization in terms of three parameters(τ, ω0, ℱ) that may be deduced from the system’s fluid properties, inaddition to the lattice spacing, L, and the Faraday wavenumber, kF.
Expressions for the parameters τ, ω0 and ℱ are provided in Extended

Article
Data Table 1, and their dependence on γ/γF is presented in Supplementary Fig. 1. Finally, we note that investigating precisely the collective
order emerging from excited spin states, such as trefoils arising forγ > γc (Fig. 1c), would require more involved models such as a full implementation of the stroboscopic model42, or the variable-topographymodel of a previous work68.
Generalized Kuramoto model. By making a series of further approximations to the form of the coupling in the Bessel model (equation (16)),
we now demonstrate that the emergent collective order in HSLs is governed by a generalization of the classical Kuramoto model for coupledphase oscillators4. From equation (16), the force on droplet i due to theadjacent droplet (i + 1) is

Fi , i +1 = − ℱJ ′0 (kF|xi(t) − xi +1(t)|)(xi(t) − xi +1(t)).
|xi(t) − xi +1(t)|(17)

We first make the approximation that the orbital radius is small relative to the centre-to-centre lattice separation, r/L ≪ 1. This approximation is well satisfied in the experiments, where r ≤ 3 mm and L > 14 mm,so r/L ≤ 0.2. We can thus approximate the distance between the two
droplets as |xi(ϕi) – xi+1(ϕi+1)| ≈ L + r(cosϕi+1 − cosϕi), and expand theazimuthal component of the wave-induced force:θ^i ⋅ Fi , i +1 ≈ − ℱ[sin ϕi + (r /L)(sin ϕi − sin ϕi +1)cos ϕi ]J 0′ (kFL − kFr (cos ϕi − cos ϕi +1)).

(18)

We then make the additional approximation that ϕ − ϕi +1   ϕi + ϕi +1 
kFr|(cos ϕi − cos ϕi +1)| = 2kFr sin  i ≪1 sin 
2
2

 


(19)

which may be written equivalently as

xi − xi+1 ≪

λF
,
2π

(20)

where (xi, yi) are the Cartesian coordinates for each drop relative tothe well centre. Although this limit will not be satisfied by all solutionsof equation (16), for the range of L considered in our experimental
study, the droplets exhibit in-phase synchronization (Δϕ± ≈ 0; Figs. 1f,2c); consequently, the observed dynamics are expected to satisfy theapproximation in equation (19).
Exploiting equation (19), we thus deduce thatJ ′0 (kFL − kFr (cos ϕi − cos ϕi +1))
≈ J ′0 (kFL) − kFr J ″0 (kFL)(cos ϕi − cos ϕi +1),

(21)

from which equation (18) indicates thatˆi ⋅ Fi , i +1 ≈ − ℱJ ′ (kFL)sin ϕ + ℱkFr J ″ (kFL)(cos ϕ − cos ϕ )sin ϕ
θ
0
0
i +1
i
i
i
−

ℱr
J ′ (k L)(sin ϕi − sin ϕi +1)cos ϕi .
L 0 F

(22)

Noting that the first term in equation (22) is cancelled when we sum
the contribution of the forces from the two neighbouring droplets,i − 1 and i + 1, we may now express the coupling force as the gradient ofa potential U, yielding a reduced model for the spin lattice system as∂Ui , j1  ω2 ωi̇ = 1 − i2  ωi − ∑,τ  ω0 
j ∂ϕi

(23)

where the sum is only over nearest neighbours within the lattice. The
generalized Kuramoto coupling potential, Ui,j, may be expressed asαβUi , j ≔ U (ϕi , ϕj ) = | cos ϕi − cos ϕj |2 + | sin ϕi − sin ϕj |222αβ= 2 |xi − xj |2 + 2 |yi − yj |2 ,
2r
2r

(24)

where

α = ℱkFr J ″0 (kFL),(25)

β = (ℱr /L)J ′0 (kFL).

(26)

We refer to this system as the ‘generalized Kuramoto’ (GK) model, andillustrate the dependence of the ratio β/α on lattice spacing, L, in Supplementary Fig. 2. We note that the form of the potential U is prescribed
only by the local wave curvature and gradient, through the signs andrelative magnitudes of α and β. Consequently, the amplitude of thelocal wave field does not substantially alter the form of the potential,
U, which further justifies the neglect of topography-induced wave
damping. We note that our GK model (equation (23)) includes inertia,
as do the Kuramoto-like models typically used to study synchronization phenomena in power grids69,70.
Notably, we found that the GK model succeeds in rationalizing the
emergent order even at intermediate values of lattice spacing, L, forwhich the approximation equation (19) is not necessarily satisfied. To
understand this unforeseen success of the GK model, we compare theform of the Bessel coupling force (equation (16)) with that expectedfrom the GK model (equation (23)) for all ϕi and ϕj and four representative lattice separations (Supplementary Fig. 3). In all four cases, wefind that the GK coupling corresponds to the two fundamental Fouriermodes of the Bessel coupling, which explains the agreement between
the predictions of the Bessel and GK models over the range of latticespacings of interest. For smaller values of L, we find that the agreement
holds until the spacing becomes comparable to the orbital radius,L ≈ 2r. In the limit of rkF ≪ 1, when the orbital radius is much smaller thanthe Faraday wavelength, the approximation (equation (19)) is strictly
satisfied. Finally, we note that the amplitude of the GK coupling force
computed via the predictions (equations (25), (26)) is a factor of twolarger than the amplitude of the dominant Fourier modes computedfrom the Bessel force (Supplementary Fig. 3). As discussed in Extended
Data Fig. 5, we compensate for this overestimation in order to obtain
the quantitative agreement between the experiments and the GK model(equation (23)).
Spin flips and wave coupling. Our theory indicates that, in the parameter regime considered in our experimental study, the fundamental
mechanism responsible for the emergent collective order does notdepend strongly on either the precise details of the spin flips or thesmall deformation of the walker wave field due to variations in bottom
topography. Nevertheless, we detail these effects here for the sake ofcompleteness.
Spin flips. Extended Data Fig. 2 (Supplementary Video 6) illustrates
a typical spin flip event as observed in our experiments with the 1D
antiferromagnetic HSL described in Fig. 1. Owing to transient rearrangements from the random initial conditions, at some time t during thisexperimental run, three neighbouring spins happen to be rotating inthe same direction (see Extended Data Fig. 2a). Since the lattice geometry promotes antiferromagnetic order, this local arrangement is not a
preferred state. The collective interactions are such that the trajectoryof the middle walker then starts to become elliptical (see Extended
Data Fig. 2b). The elongation of the ellipse continues until the walker
executes quasi-rectilinear motion across the well centre (see Extended
Data Fig. 2c). Subsequently, the process is reversed, the result being thewalker rotating in the opposite direction (see Extended Data Fig. 2d).
After the spin flip, which typically takes 5 ± 1 s, the three walkers rotate
in the preferred antiferromagnetic configuration (see Extended DataFig. 2e). When the collective order is strongest (γ/γF ≈ 86% in Fig. 1h),the duration of a flip, of the order of 5 s, is negligible relative to thecharacteristic flipping period, of the order of 5–10 min. At higher accelerations, when the walkers switch direction much more frequently andthe coherent order is lost (γ/γF > 90% in Fig. 1h), the duration of the flips
becomes comparable to the flipping period. We note that ferromagenetic lattices, such as those shown in Fig. 2, undergo similar spin flips.
Wave coupling. A rough estimate of the extent to which the variablebottom topography will alter the form of the wave field relative to thatpredicted by the stroboscopic model (equation (13)) may be obtainedby comparing the Faraday wavelengths given by the water-wave disperhHsion relation for the shallow λF and deep λF regions. As noted earlier
in the experimental description, these wavelengths are only marginallyh H
different for our choice of parameters, specifically λF /λF ≈ 1.05. One
thus expects that the walker wave field in our lattices will not be appreciably different from that previously predicted42. Indeed, as shown in
Extended Data Fig. 3a, b, direct visualization of waves generated by asingle walker reveals a wave form comparable to that expected from asuperposition of sub-harmonic zeroth-order Bessel functions, as inequation (13). We note that over the range of parameters consideredin our study, no exotic Faraday patterns were excited, such as thosearising in nonlinear regimes71,72 or excited by droplets walking in elliptical corrals48.
To provide a quantitative characterization of the wave coupling,including the spatial decay36, we perform additional simulations witha previously developed model68, wherein changes in topography areincorporated through their influence on the local wave speed. We notethat this model solves the wave problem explicitly, and has been benchmarked against experiments in a number of settings involving walkers
interacting with submerged topography37,73–75. In accordance with ourexperimental observations, we find that the simulated wave field has aform similar to that of a J0 Bessel function (see Extended Data Fig. 3c). By
taking into consideration spatial damping, it has been demonstrated36
that in the deep-fluid regime, the wave field generated by a bouncing
drop located at xj may be approximated byh(|x − xj|) = 2A b J0 (kF|x − xj|)|x − xj|K1(2 b |x − xj|)

(27)

where K1 is the modified Bessel function of the second kind. Our
simulations indicate that the fitting parameters are A ≈ 15.4 μm andb ≈ 0.002 mm−2, which is roughly consistent with previous free-surface
synthetic Schlieren experiments76. We note that the spatial damping
included in equation (27) does not play a crucial role in rationalizing
the mechanism responsible for the emergent collective order; thus, we
can safely neglect its effects. Finally, we note that experiments support
our assumption that the nearest-neighbour couplings are dominant:
the waves become negligible before reaching the second-nearest wells.
Comparison to experiment. We proceed by comparing the results ofour experimental study with the theoretical predictions of both thefull wave coupling Bessel model (equation (16)), and the generalized
Kuramoto model (equation (23)). Specifically, we characterize thedependence of the collective lattice order on both the well separation, L, and vibrational acceleration, γ, and demonstrate that the
emergent orbital phase synchronization is determined by the formof the coupling force.
We compare the predictions of our theoretical models to the experimental results through the average spin-correlation, ⟨χ⟩. In the simulations, for every data point, we perform Nt independent trials of a given
lattice. We initialize each simulation with random initial conditions(ϕ0, ω0) for each droplet, independently drawn from uniform distributions with ϕ0 ∈ [0, 2π] and ω0 ∈ [−3ω0/2, 3ω0/2] to ensure a large spread.
From each independent simulation, we compute a spin correlation timeseries χn(t), for n = 1, …, Nt, as defined in equation (4). We then computethe average of the instantaneous spin correlations across all trials asNχ (t ) =1 t∑ χ (t).
Nt n =1 n(28)

We compute the final average spin correlation, ⟨χ⟩, by time-averaging
over the final quarter of this time series, so as to eliminate any influenceof the initial transient regime. In our simulations of equations (16) and(23), this transient period is generally <100 s. To guarantee robust statistics, we perform all simulations up to a time tmax ≈ 300–600 s. We
note that, owing to the approximations made in our modelling of the
wave form, the preferred collective order is attained more rapidly insimulations than experiments.
We perform numerical tests to ensure that our results are not affectedby either an insufficient number of trials or finite-size effects. Our simulations indicate that the collective behaviour becomes largely independent of the number of trials for Nt = 200 (Supplementary Fig. 4a).
Similarly, our exploration of the effect of lattice size on the spin correlation (Supplementary Fig. 4b) reveals that the statistics become
independent of the number of spins for N ≈ 150. We thus select Nt = 200
and N = 150 for all subsequent simulations.
Varying lattice spacing. As demonstrated experimentally (Fig. 2),
there exists a transition from antiferromagnetic to ferromagnetic orderwhen the well separation, L, decreases from 17.7 to 13.2 mm. Continuously varying L is extremely prohibitive experimentally; thus, we proceed to do so using simulations.
We computed the average spin correlation, ⟨χ⟩, in terms of the latticespacing, L, via the Bessel model (equation (16)). For our simulations,
we select a relaxation time of τ = 0.4 s, an angular velocity of ω0 = 3.3 s−1,
and an interaction coefficient of ℱ = 70 s−2, values consistent with theexperimental range γ ≤ γc (Supplementary Fig. 1). We prescribe theorbital radius to be r = 1.2 mm, which is close to the average observed
experimentally with L = 13.2 mm. The Faraday wavelength in the simulation is fixed at λF = 2π/kF = 4.95 mm. We vary the well-to-well separation
L across one Faraday wavelength, between 13 and 18 mm.
We present the dependence of the time-averaged spin correlation,⟨χ⟩, on lattice spacing, L, in Extended Data Fig. 4a. Our results showthat the preferred collective behaviour oscillates between regions ofpositive and negative spin correlation as the lattice separation is varied.
To better characterize the emergent ferromagnetic and antiferromagnetic order, we also compute the mean droplet phase difference, Δϕ±(see Extended Data Fig. 4b), which reveals that the collective lattice
behaviour is dominated by not two but four synchronization modes(see Extended Data Fig. 4c, d).
In the first ferromagnetic region, denoted FM+, the model predicts
that ⟨Δϕ±⟩ ≈ 0, indicating that neighbouring pairs of spins are on average in phase in the azimuthal direction, regardless of whether they
co-rotate or counter-rotate. The bias towards ferromagnetic order(see Extended Data Fig. 4a) is thus associated with the dominance ofthe co-rotating synchronization mode over the counter-rotating mode(see Extended Data Fig. 4c, d). Similarly, in the antiferromagnetic region
denoted by AFM+, we again observe a vanishing mean phase difference, ⟨Δϕ±⟩ ≈ 0, indicating in-phase orbital motion, but now the relative
dominance of the two synchronization modes is reversed. By contrast,
the ferromagnetic and antiferromagnetic regions denoted by FM− andAFM−, respectively, result from two different synchronization modes
that are characterized by out-of-phase orbital motion, ⟨Δϕ±⟩ ≈ π (seeExtended Data Fig. 4c, d). The relative dominance of these modes thus
determines the bias towards positive or negative spin correlation (seeExtended Data Fig. 4a).
The origins of these four forms of collective order are rationalized
through the generalized Kuramoto model (equation (23)). By examining the dependence of the control parameters α and β, as defined inequations (25), (26), on lattice spacing L, we note that the four dynamical modes are those that precisely minimize the reduced coupling
potential U (equation (24)). Specifically, when α > 0, which happens
when J ″0 (kFL) > 0, corresponding to FM+ (see Extended Data Fig. 4a),

Article
the first term of the coupling potential U strongly favours motion forwhich |xi − xj| is minimized. In-phase rotation is thus preferred, but thiscondition is equally satisfied by both co- and counter-rotating spins(see Extended Data Fig. 4c, d, column A). The observed bias towards
ferromagnetic order thus comes from the minimization of the second
term in U, which is determined by the sign of β, and so by the local slope
J ′0 (kFL). In particular, β > 0 in the FM+ region (see Extended Data Fig. 4a),indicating that, between the two possible in-phase modes, the system
will favour the co-rotating mode because it also minimizes |yi − yj| (seeExtended Data Fig. 4c, d). The emergence of ferromagnetic order FM+
is thus rationalized. Similar arguments may provide rationale for theobserved order in the remaining regions. For instance, as arises for
FM+, α > 0 in AFM+ and so in-phase rotation is preferred. However, since
β < 0, the system promotes counter-rotating spins as those now minimize the potential U by maximizing |yi − yj|. We note the regions where
no preferred collective order emerges (see Extended Data Fig. 4a)
precisely correspond to the extrema and inflection points of J0(kFL), at
which either α = 0 or β = 0. Finally, we note that simulations of the Bessel model (equation (16)) with 50 × 50 lattices confirmed that a similar
dependence of collective order on lattice spacing emerges for substantially larger 2D square lattices (see Extended Data Fig. 6a).
Varying interaction strength. In Extended Data Fig. 5, we compare
the experimentally observed dependence of the average spin correlation, ⟨χ⟩, on the bath acceleration, γ/γF, with the predictions basedon both the Bessel model (equation (16)) and the generalized Kuramoto model (equation (23)). For values of γ/γF just beyond the walking
threshold, both models predict the relatively weak spin correlations
associated with relatively weak coupling. As the coupling strength isincreased progressively, both models predict the emergence of a more
pronounced antiferromagnetic response, up to a critical acceleration
γc, where min(χ) ≈ −0.4, in excellent agreement with the experiments.
Beyond this optimal driving, the magnitude of the spin correlation
decays rapidly, until the preferred collective order vanishes. Strictly
speaking, the validity of our theoretical models is restricted to γ ≤ γc,where the droplet motion is confined to circular orbits. Nevertheless,both models capture the decay of the antiferromagnetic order beyond
γc. Experimentally, such decay is related to the destabilization of circular orbits into trefoil-like trajectories (Fig. 1h) that in turn prompts
chaotic switching in orbital direction. In both theoretical models,similar chaotic switching emerges as the coupling parameters, ℱ andα, β, are increased to values corresponding to γ > γc. Similar simulations with both the Bessel (equation (16)) and generalized Kuramoto(equation (23)) models were conducted for a 50 × 50 lattice in order toconfirm that the same collective order emerges in 2D square lattices
for variable memory (see Extended Data Fig. 6d).
Bath rotation. We conclude by generalizing our theoretical framework
to rationalize the transition from antiferromagnetic to ferromagnetic
order induced by applying a constant rotation, Ω, to the vibrating bath(Fig. 3a–e).
We recall that the rotating bath surface is an equipotential; consequently the radially outward centrifugal force on the drop is preciselybalanced by the inward interfacial force resulting from the parabolic
bath surface and so need not be considered henceforth8. The effect ofbath rotation may thus be incorporated into our modelling throughthe addition of the Coriolis force, Fc = − 2mΩ × ẋ (ref. 47), which acts in
the inward radial direction for circular cyclonic/anticyclonic orbits.
To correctly model the HLS dynamics in a rotating frame, we should
thus include the radial force balance into our models. To proceed, we
model the force exerted by the well in the low-memory regime by aconfining radial potential of the form V(r) = kr4, so that Fw = −dV/dr^r,
where r = |x| (ref. 64). We find that the fourth-order degree of the potential captures adequately the abrupt influence of topography at the wellboundary. We select the strength of the potential, k, so that the predictions for the free walking speed u0 and experimentally observed orbital

radius r are consistent with our experimental observations in theabsence of rotation (Fig. 1c).
Focusing on the wave-mediated droplet coupling Fi,j as approximatedby the generalized Kuramoto force (equation (23)), we may thus write
the 2D version of the Rayleigh oscillator equation (9), with an applied
force F = Fc + Fw + Fi,j as

1  |ẋ |2 ẍ i = 1 − i2  ẋ i − 4k|xi|2 xi − 2Ω × ẋ i + ∑ Fi , j (xi, xj, L).
τ
u0 
j

(29)

Noting that the orbital radius is now a variable ri(t), we may express
the droplet position around the centre of the ith well in polar coordinates asxi(t) = ri(t)(cos ϕi (t), sin ϕi (t)).

(30)

Since full resolution of the droplet’s dynamics in the presence ofrotation is beyond the scope of the theory presented here, we simplify
the equations to rationalize the observed behaviour. We thus neglect
the contribution of ri ̇ to the Rayleigh drag term, which is valid provided
that r i2̇ ≪ r i2ωi2, and the contribution of the Rayleigh drag term in theradial equation, which is valid provided that |ri|̇ ≪ τri ωi2. Since the droplet–droplet interaction primarily drives the phase synchronization
between the droplets, we also neglect the wave-mediated interaction
force from the radial force balance. Finally, we assume that the interaction force does not vary with the orbital radius, so that α and β are held
constant, which is valid provided the droplet trajectories do not deviate substantially from circles. With these assumptions, we may thensubstitute equation (30) into equation (29) to yieldr¨i = ri ωi2 + 2Ωωi ri − 4kr i3,(31)∂Ui , jr
2˙
1  r 2ω 2 
˙i = 1 − i 2 i ωi − i (Ω + ωi ) − ∑
ω
.
τ
ri
u0 
j ∂ϕi

(32)

The numerical prediction of the model formed by equations (31) and(32) are compared with our experimental results in Fig. 3e, showing
excellent agreement given the level of theory reduction. Simulations
are performed with N = 150 spins and Nt = 200 independent trials. We
fix τ = 0.1 s, u0 = 6 mm s−1, r = 1.8 mm, k = 4.34 mm−2 s−2, α = 10 s−2 andβ = −1.5 s−2. We then vary the rotation rate, Ω, over the range 0 < Ω < 0.7 s−1
as in the experiments. Our results reveal a dependence of the average spin correlation, ⟨χ⟩, on the rotation rate, Ω, very similar to thatobserved in the laboratory.

Tunability
HSLs also offer the opportunity to investigate more exotic forms ofcollective behaviour and frustration effects induced by wave-mediatedinteractions. Here, we provide several proof-of-concept experiments
and simulations in order to illustrate their potential for future research.
As discussed in Methods section ‘Stationary bath’, it is possible
to deduce the vertical phase of the droplets by observing the lightreflected from the upper side of the droplets, which intensifies
when the droplets deform during impact. Since the walkers bounce
sub-harmonically with respect to the bath vibration, drops may bounce
vertically in-phase or out-of-phase with respect to their neighbours(see Extended Data Fig. 7a). Vertically out-of-phase walkers are coupled through a wave field equal to that of in-phase walkers but with the
opposite sign (see Extended Data Fig. 4). Consequently, for example,
the emergent order of an HSL, whose lattice spacing promotes orbital
in-phase antiferromagnetic order AFM+ for vertically in-phase walkers,
will shift to out-of-phase antiferromagnetic order AFM− for verticallyout-of-phase walkers (see Extended Data Fig. 7a). An HSL with walkers
with random vertical phases will thus produce an antiferromagnetic

state with competing orbital synchronization that effects the emergent
spin correlation (see Extended Data Fig. 7b). Other forms of frustration
can be investigated by changing the lattice geometry in 2D. For example,
a triangular lattice (see Extended Data Fig. 7c) in which the spacing istuned to promote antiferromagnetic order leads to degenerate groundstates, the investigation of which is expected to yield insights into how
wave-coupled spin lattices respond to geometric frustration77.
We also perform simulations using a previously published model68
that solves the wave field explicitly, in order to demonstrate how alterations of the basic experimental setup might allow one to tailor the
coupling between spins and so study more exotic collective phenomena. We note that in these simulations, all the walkers are bouncing
vertically in phase. We first confirm that orbital in-phase ferromagnetic FM+ and antiferromagnetic AFM+ orders are reproduced withthe model68 (see Extended Data Fig. 8). We proceed by simulating alattice where the horizontal spacing promotes FM+, while the verticalspacing promotes AFM+ (see Extended Data Fig. 9a, left). Drops along
the horizontal lines are thus rotating in phase, but their direction alternates vertically in an antiferromagnetic fashion (see Extended Data
Fig. 9a, centre and right).
Next, we take as reference the ferromagnetic HSL (see Extended
Data Fig. 8a) and repeat the simulation with a lattice where the centreof each well has been shifted by a random distance ±ε in the horizontal
and vertical directions (see Extended Data Fig. 9b, left). As the degreeof randomness increases, the walkers’ trajectories become less regular;thus, it becomes more difficult for the walkers to maintain precise
in-phase orbital motion (see Extended Data Fig. 9b, centre and right).
In our third demonstration, with the same FM+ as the reference (seeExtended Data Fig. 8a), we alternate the sizes of the bouncing droplets,
thereby mixing two ‘species of spins’ with different frequency andinertia (see Extended Data Fig. 9c, left). Although the drop size does
not change the preferred mode of synchronization, which is dictatedby the lattice spacing, the different rotation frequencies lead to frustration. Notably, the drops adjust their orbital radii (see Extended Data
Fig. 9c, centre) in order to match their rotation frequencies, as can be
deduced from the equal slope of the ϕ(t) signals (see Extended Data
Fig. 9c, right).
In our fourth demonstration, we reduce the local liquid depth along
alternating columns of a FM+ lattice geometry (see Extended Data
Fig. 9d, left). As a result, the waves above the new submerged barriers
dissipate waves more effectively (see Extended Data Fig. 9d, centre).
Walkers are thus subjected to lateral couplings of different strengths,
which leads to pairs whose phase is either strongly or weakly correlated(see Extended Data Fig. 9d, right), and so to the possibility of exploring
SSH-type couplings as studied in topological insulators78. We also notethat the driving frequency can be used to dynamically tune the wavelength, range, and amplitude of the oscillatory coupling. Moreover, one
might superimpose random fluctuations on the vibrational driving inorder to study the effects of imposed white or coloured noise.
HSLs offer the opportunity to explore other potentially fruitful directions. One of them is the analysis of the collective order emerging fromhigher excited spin states such as trefoils and lemniscates44, which
emerge with larger well diameters64, and at higher driving acceleration(Fig. 1c). Notably, the accompanying theoretical rationalization would
require more involved models, such as a full implementation of thestroboscopic model42 or the variable-topography model68. Similarly,
the liquid depth, lattice geometry, and vibrational acceleration could
be tuned to study the effects of longer-range oscillatory couplings
beyond the nearest-neighbour interactions considered in this study,
which might lead to more exotic dynamical states. HSLs thus offer theopportunity to investigate macroscopic dynamics analogous to those
involving the oscillatory Ruderman–Kittel–Kasuya–Yosida (RKKY)38
electronic spin coupling. Finally, investigating the similarities anddifferences between the collective synchronization arising in HSLs andother active systems might also prove fruitful. For example, it is known

that rotating vortices in reaction–diffusion systems may be coupled
by nonlinear spiral waves79–82.

Data availability
The data that support the findings of this study are available from thecorresponding authors upon reasonable request.
41.
42.
43.
44.
45.
46.
47.
48.
49.
50.
51.
52.
53.
54.
55.
56.
57.
58.
59.

60.
61.
62.
63.
64.
65.
66.
67.
68.
69.
70.
71.
72.
73.
74.
75.
76.
77.
78.
79.
80.

Oza, A. U., Rosales, R. R. & Bush, J. W. M. Hydrodynamic spin states. Chaos 28, 096106(2018).
Oza, A. U., Rosales, R. R. & Bush, J. W. M. A trajectory equation for walking droplets:
hydrodynamic pilotwave theory. J. Fluid Mech. 737, 552–570 (2013).
Labousse, M. et al. Self-attraction into spinning eigenstates of a mobile wave source by
its emission backreaction. Phys. Rev. E 94, 042224 (2016).
Perrard, S. et al. Self-organization into quantized eigenstates of a classical wave-driven
particle. Nat. Commun. 5, 3219 (2014).
Eddi, A. et al. Level splitting at macroscopic scale. Phys. Rev. Lett. 108, 264503 (2012).
Harris, D. M. & Bush, J. W. M. Droplets walking in a rotating frame: from quantized orbits tomultimodal statistics. J. Fluid Mech. 739, 444–464 (2014).
Oza, A. U. et al. Pilot-wave dynamics in a rotating frame: on the emergence of orbital
quantization. J. Fluid Mech. 744, 404–429 (2014).
Sáenz, P. J., Cristea-Platon, T. & Bush, J. W. M. Statistical projection effects in ahydrodynamic pilot-wave system. Nat. Phys. 14, 315–319 (2018).
Eddi, A. et al. Archimedean lattices in the bound states of wave interacting particles.
Europhys. Lett. 87, 56002 (2009).
Eddi, A., Boudaoud, A. & Couder, Y. Oscillating instability in bouncing droplet crystals.
Europhys. Lett. 94, 20004 (2011).
Thomson, S. J., Couchman, M. M. P. & Bush, J. W. M. Collective vibrations of confinedlevitating droplets. Phys. Rev. Fluids 5, 083601 (2020).
Thomson, S. J., Durey, M. & Rosales, R. R. Collective vibrations of a hydrodynamic active
lattice. Proc. R. Soc. A 476, 20200155 (2020).
Couchman, M. M. P. & Bush, J. W. M. Free rings of bouncing droplets: stability anddynamics. J. Fluid Mech. 903, A49 (2020).
Nachbin, A. Walking droplets correlated at a distance. Chaos 28, 096110 (2018).
Nachbin, A. Kuramoto-like synchronization mediated through Faraday surface waves.
Fluids 5, 226 (2020).
Whitham, G. B. Linear and Nonlinear Waves (Wiley, 1974).
Douady, S. Experimental study of the Faraday instability. J. Fluid Mech. 221, 383–409(1990).
Harris, D. M. & Bush, J. W. M. Generating uniaxial vibration with an electrodynamic shaker
and external air bearing. J. Sound Vibrat. 334, 255–269 (2015).
Faraday, M. On a peculiar class of acoustical figures and on certain forms assumed bygroups of particles upon vibrating elastic surfaces. Philos. Trans. R. Soc. Lond. 121,
299–340 (1831).
Harris, D. M., Liu, T. & Bush, J. W. M. A low-cost, precise piezoelectric droplet-on-demand
generator. Exp. Fluids 56, 83 (2015).
Bush, J. W. M., Oza, A. U. & Moláček, J. The wave-induced added mass of walking droplets.
J. Fluid Mech. 755, R7 (2014).
Turton, S. E. Theoretical Modeling of Pilot-wave Hydrodynamics. PhD thesis,
Massachusetts Institute of Technology (2020).
Labousse, M. & Perrard, S. Non-Hamiltonian features of a classical pilot-wave dynamics.
Phys. Rev. E 90, 022913 (2014).
Cristea-Platon, T., Sáenz, P. J. & Bush, J. W. M. Walking droplets in a circular corral:quantisation and chaos. Chaos 28, 096116 (2018).
Tadrist, L. et al. Faraday instability and subthreshold Faraday waves: surface waves
emitted by walkers. J. Fluid Mech. 848, 906–945 (2018).
Galeano-Rios, C. A. et al. Ratcheting droplet pairs. Chaos 28, 096112 (2018).
Couchman, M. M. P., Turton, S. E. & Bush, J. W. M. Bouncing phase variations in pilot-wavehydrodynamics and the stability of droplet pairs. J. Fluid Mech. 871, 212–243 (2019).
Faria, L. M. A model for Faraday pilot waves over variable topography. J. Fluid Mech. 811,
51–66 (2017).
Filatrella, G., Nielsen, A. H. & Pedersen, N. F. Analysis of a power grid using a
Kuramoto-like model. Eur. Phys. J. B 61, 485–491 (2008).
Rohden, M. et al. Self-organized synchronization in decentralized power grids. Phys. Rev.
Lett. 109, 064101 (2012).
Chen, P. & Wu, K.-A. Subcritical bifurcations and nonlinear balloons in Faraday waves.
Phys. Rev. Lett. 85, 3813–3816 (2000).
Périnet, N., Juric, D. & Tuckerman, L. S. Alternating hexagonal and striped patterns in
Faraday surface waves. Phys. Rev. Lett. 109, 164501 (2012).
Pucci, G. et al. Non-specular reflection of walking droplets. J. Fluid Mech. 804, R3(2016).
Pucci, G. et al. Walking droplets interacting with single and double slits. J. Fluid Mech.
835, 1136–1156 (2018).
Harris, D. M. et al. The interaction of a walking droplet and a submerged pillar: fromscattering to the logarithmic spiral. Chaos 28, 096105 (2018).
Damiano, A. P. et al. Surface topography measurements of the bouncing droplet
experiment. Exp. Fluids 57, 163 (2016).
Diep, H. T. Frustrated Spin Systems (World Scientific, 2004).
Qi, X.-L. & Zhang, S.-C. Topological insulators and superconductors. Rev. Mod. Phys. 83,
1057–1110 (2011).
Biktasheva, I. V. & Biktashev, V. N. Wave–particle dualism of spiral waves dynamics.
Phys. Rev. E 67, 026221 (2003).
Biktasheva, I. V. et al. Computation of the drift velocity of spiral waves using response
functions. Phys. Rev. E 81, 066202 (2010).

ArticleAzhand, A., Totz, J. F. & Engel, H. Three-dimensional autonomous pacemaker in thephotosensitive Belousov–Zhabotinsky medium. Europhys. Lett. 108, 10004 (2014).
82. Totz, J. F., Engel, H. & Steinbock, O. Spatial confinement causes lifetime enhancement
and expansion of vortex rings with positive filament tension. New J. Phys. 17, 093043(2015).

preliminary experiments. A.G. contributed to the preliminary experiments. J.D. contributed tothe theoretical modelling and the writing of the paper. J.W.M.B. contributed to the conceptionof the experiments and theory, and to writing the paper.

Acknowledgements We gratefully acknowledge financial support from the NSF throughCMMI-1727565 (J.W.M.B. and P.J.S.) and DMS-1719637 (R.R.R.), MIT Solomon Buchsbaum
Research Fund (J.D.), and CNRS Momentum programme (G.P.).

Additional information
Supplementary information The online version contains supplementary material available athttps://doi.org/10.1038/s41586-021-03682-1.
Correspondence and requests for materials should be addressed to P.J.S. or J.W.M.B.
Peer review information Nature thanks Michael Shats and the other, anonymous, reviewer(s)for their contribution to the peer review of this work.
Reprints and permissions information is available at http://www.nature.com/reprints.

81.

Author contributions P.J.S. conceived the study, led the experimental developments and thewriting of the paper, and contributed to the theoretical modelling. S.E.T. and R.R.R. contributed
to the theoretical modelling. G.P. contributed to the conception and execution of theCompeting interests The authors declare no competing interests.

Chamber

Rb

Bath

H

d
h
Acrylic plate

Accelerometer

Levelling plate

Optical
table

Motor

Air
Bearing

Coupling rod
Extended Data Fig. 1 | Schematic of the experimental set-up. The test
section was mounted on an optical table and vibrated vertically by anelectromagnetic shaker. The shaker was connected to the bath by a thin steel
rod coupled with a linear air bearing. The forcing acceleration was monitored

cos(2

)

by two piezoelectric accelerometers. The bath was enclosed within atransparent acrylic chamber to ensure that ambient air currents did not affectthe experiments. A d.c. motor housed inside the hollow air bearing enabled
system rotation.

Article
a

b

c

d

e

10
9.5
9
8.5
8
7.5

t

t + 9.03 s

t + 11.49 s

Extended Data Fig. 2 | Experimental spin flips. a−e, Snapshots at different
times illustrate a typical spin-flip event arising in a 1D spin lattice(Supplementary Video 6). At each time, the panels at left are colour-coded
according to the instantaneous spin S i (same colour map as in Fig. 1d), while
those on the right depict the recent walker trajectories, colour-coded
according to the local speed. Perturbed by the wave fields emitted by its

t + 12.72 s

t + 16.00 s

v (mm/s)neighbours, the middle walker initially follows an elliptical path. The length ofthe minor axis decreases until the walker trajectory essentially becomes astraight line across the well centre. Subsequently, the process is reversed,
resulting in the walker rotating in the opposite sense. The three walkers shown
are part of the 1D antiferromagnetic lattice described in Fig. 1h with
γ/γF = 86.0%.

acsubmerged wellsbExtended Data Fig. 3 | Wave coupling. a, Experimental visualization of the
wave field generated by a single walker in a 1D lattice with the same geometry as
in Fig. 1h and γ/γF = 92.0%. The submerged wells can be identified as the regionswith a different shade of grey. b, Superposition of the wave field shown in a andthe zeros of the drop-centred Bessel function J0 (kF|x ‐x i|). c, Wave field of abouncer computed with the theoretical model developed previously68 forwalkers over variable topography. The bouncer is located at (x, y) = (3D/8, 0) in a
2D square lattice with the same well diameter D and centre-to-centre separation
L as in a and γ/γF = 88.0%. Solid blue lines denote the submerged wells anddashed lines the zeros of a Bessel function J0 centred at the drop position.

Article
0.5

ave. phase diff. ave. spin corr.
(rad)
<

a

d

FM _

AFM+

A

B

C

D

_
2
0

_ (counter-rot)+ (co-rot)

13

14

B

C

17

18

D

dominant
mode

A

15
16
lattice spacing L (mm)sub-dominant
mode

c

AFM _

0

-0.5

b

FM+

Extended Data Fig. 4 | Emergent order for varying lattice spacing. a, b, The
dependence of the average spin correlation, ⟨χ⟩ (a) and the mean phasedifference, ⟨Δϕ⟩ (b) on the lattice spacing, L, as predicted by the Bessel model(equation (16)) with τ = 0.4 s, ω 0 = 3.3 s−1, and ℱ = 70 s−2. The average
droplet-induced wave field, which is approximated by a Bessel function J0(k FL),
centred in the neighbouring well at L = 0, is provided for comparison in a.
c, d, Dominant (c) and subdominant (d) synchronization modes for thecorresponding sections (A–D) in a, b.

ave. spin corr. <

0.5
0.25

experiments
Bessel model
GK model
c

0

-0.25
-0.5

78 82 86 90 94
bath acceleration F (%)

Extended Data Fig. 5 | Emergent order for varying bath acceleration.
Comparison of the experimentally observed average spin correlation, ⟨χ⟩, withthe predictions of the Bessel model (equation (16)) and the generalized
Kuramoto model (equation (23)). Solid lines are fits resulting from smoothing
the piecewise linear plot given by connecting the points. Bessel model
parameters: L = 17.7 mm, r = 1.8 mm, λF = 2π/k F = 5.1 mm. The interaction
parameter ℱ is varied across the range 70 < ℱ < 130 s−2 and transformed back toγ/γF using the relation from Extended Data Table 1. To simplify the simulations,
we fix the relaxation time to be τ = 0.1 s and the natural angular frequency to beω 0 = 3.3 s−1, values consistent with experimental observation (Fig. 1c). The
effective walker mass is set to mw = 1.65m, in line with prior work61. GK model
parameters: α is varied over the range 8.5 < α < 15 s−2, while maintaining β < 0 anda constant ratio |β/α| = 0.3. By dividing the expressions for α and β inequations (25), (26) by a factor of two, in accordance with the mismatch
between the GK and Bessel models discussed in Supplementary Fig. 3, theminimum ⟨χ⟩ predicted by the GK model emerged in the vicinity of γc.

ave. spin corr. <a 0.2

Bessel model

b

c

0

-0.2

17
18
19
lattice spacing, L (mm)d0ave. spin corr. <Article

-0.1
-0.2
-0.3
-0.4
76

Bessel model
GK model78
80
82
bath acceleration

84
(%)
F

b

L = 17.1 mmc

L = 18.5 mm
Extended Data Fig. 6 | Emergent order in large 2D square lattices.
Simulations of the Bessel model (equation (16)) and the generalized Kuramotomodel (equation (23)) for a 50 × 50 square lattice demonstrate the emergenceof antiferromagnetic and ferromagnetic order in 2D for various lattice spacingsand bath accelerations. a, The lattice spacing determines the emergent
antiferromagnetic (ADM+) or ferromagnetic (FM+) order in a manner predictedby our reduced theory (equation (24)). b, c, Specifically, preferred in-phaserotation between neighbouring pairs can be clearly observed in theantiferromagnetic AFM+ (b), and ferromagnetic FM+ (c) regimes. We note that b
and c correspond to simulations of the Bessel model with the spacings

indicated on a. d, The emergent in-phase antiferromagnetic order (AFM+) as afunction of bath acceleration. Bessel model parameters in a: τ = 0.1 s, ℱ = 72 s−2,
16.8 ≤ L ≤ 19 mm, ω 0 = 3.3 s−1, and λF = 4.95 mm. Bessel model parameters in d are
the same as in a, but with L = 17.1 mm and the interaction parameter varies
across the range 65 ≤ ℱ ≤ 85 s−2, which is transformed back to γ/γF using therelation from Extended Data Table 1. GK model parameters in d: α is varied over
the range 9.5 < α < 13 s−2 while maintaining β < 0 and a constant ratio |β/α| = 0.07,
τ = 0.2 s and ω 0 = 3.3 s−1. In all cases, each data point results from averaging 50
simulations of 600 s each, to ensure statistical significance.

avertically in-phase pairb

out-of-phase pair

0.5

d

0.25
ave. spin corr. <ecall out-of-phase

0

-0.25

randomall in-phase

-0.5

78

82
86
bath acceleration

F

(%)

90

94

Extended Data Fig. 7 | Emergent order for different vertical bouncing
synchronizations and lattice geometries. a, Oblique view of a 2D spin latticewhere the walker in the centre is bouncing vertically in-phase and out-of-phasewith its left and right neighbours, respectively. b, Average spin correlation for
lattices with the same geometry as those described in Fig. 1 when the walkers all
bounce in phase (blue, result presented in text), out of phase (green), or haverandomly distributed bouncing phases (red). Vertically in-phase pairs promote
in-phase orbital antiferromagnetic order (AFM+), whereas verticallyout-of-phase pairs promote out-of-phase orbital antiferromagnetic order(AFM−). A random distribution of vertical phases thus leads to competing
orbital synchronization modes, which has an effect on the emergent spincorrelation. Solid lines are fits resulting from smoothing the piecewise linear
plot given by connecting the points. c, Triangular HSLs with a lattice spacing
tuned to promote antiferromagnetic order can be used to investigate
frustration effects.

Article-10

wave height
-5
0

5

10

Periodic BCs
FM+

5

F

(rad)

FM+

FM+

0

phase

a

0

-5

AFM+

5

F

AFM+

(rad)

AFM+

0

phase

b

0

-5

-5

0

5
F

Extended Data Fig. 8 | Emergent collective order in simulations.
Simulations of square HSLs with the theoretical model developed previously68
is used to explore collective order in 2D. a, b, For appropriate lattice spacings,

300

350

400
450
t/TF

500

in-phase ferromagnetic order FM+ (a) and in-phase antiferromagnetic order
AFM+ (b) are simulated. Left, schematic; middle, wave field and walker’s
trajectories; right, time evolution of the orbital phases.

wave height
-5
0

-10

5

10

Periodic BCs
FM+

5

AFM+

F

FM+ horizontallybut
AFM+ vertically

0

phase

a

0

-5

b
5

0

phase

F

FM+
in well position

0

-5

c

14

FM+

5

FM+

12
10
0
8
6

-5

speed

drop sizes

0

phase

F

FM+

4
FM+

local couplings
strongcoupling

Cross section

coupling

0

phase

FM+

5

F

FM+

shallower liquid layer

d

0

-5

-5

Extended Data Fig. 9 | Model tunability. Proof-of-concept simulations
performed with the model developed previously68 illustrate the tunability andpotential for future research of HSLs. Left, schematic; middle, wave field anddrop trajectories; and right, orbital phase evolution. a, An HSL tuned topromote FM+ along the horizontal direction, but AFM+ across vertical pairs.

0

5
F

300

350

400
450
t/TF

500

b, FM+ lattice geometry with a random vertical and horizontal shift ±ε in theposition of each well. c, FM+ lattice geometry with two drop sizes (and so, twowalker speeds). d, FM+ lattice geometry with coupling strength controlled
locally through the thickness of the liquid layer.

Article
Extended Data Table 1 | Model parameters and physical variables for HSLs

Variableσ, µ, ρ, ν = µ/ρ, νe , µa

Description
Liquid surface tension, dynamic viscosity, density, kinematic viscosity,
effective kinematic viscosity, air viscosity

m, Dd , Rd

Droplet mass, droplet diameter, droplet radius

D, L, rWell diameter, center-center separation, orbital radiusγ, fBath forcing acceleration, frequency

kF , TF = 2/f , γF

Faraday wavenumber, period and threshold

TD = (νe kF2 )−1 ,TD
TM = (1−γ/γ
F)

C = 0.17, D = Cmg
F =2

νe TF (mgkF Rd )
2 +ρg
2π
3σkF

mw = m 1 +DR = DM

TM
F = TFF m
wrsin Φ, sin Φ2
F kF TM
2m

2
F kF TM
2TF D − 1

u20 = 3k28T 2
F

2

ρRd
σ + 6πRd µa

Viscous decay time, wave-decay memory-time [36]Drag coefficient, time-averaged horizontal drag coefficient [36]Wave-force coefficient, phase parameter [36]Effective walker mass

w
, τ=m
DR

2
F kF TM
2TF D − 1

, ω0 = ur0

Rayleigh-drag strength, Relaxation time
Free walking speed, angular velocity
Interaction force coefficient