Pedro J Sáenz
email: saenz@unc.edu

Giuseppe Pucci

Sam E Turton

Alexis Goujon

Rodolfo R Rosales

Jörn Dunkel

John W M Bush
email: bush@math.mit.edu

Emergent order in hydrodynamic spin lattices

 derived from first principles rationalize our experimental observations, establishing HSLs as a versatile platform for exploring active phase oscillator dynamics. The tunability of HSLs suggests exciting directions for future research, from active spin-wave dynamics to hydrodynamic analog computation and droplet-based topological insulators.

Broken mirror symmetry is an essential feature of many physical and biological systems, from magnetic metals and composites [START_REF] Coey | Magnetism and Magnetic Materials. 3rd printing[END_REF] to DNA [START_REF] Watson | The Structure of DNA[END_REF]. The most fundamental manifestation of a broken left-right symmetry are the spins of elementary particles such as photons or electrons [START_REF] Weinberg | The Quantum Theory of Fields[END_REF]. Microscopic spin systems and their macroscopic analogs [START_REF] Kane | Topological boundary modes in isostatic lattices[END_REF][START_REF] -G. Chen | Nonlinear conduction via solitons in a topological mechanical insulator[END_REF] are attracting increasing interest [START_REF] Baltz | Antiferromagnetic spintronics[END_REF][START_REF] Chumak | Magnon spintronics[END_REF] as building blocks of information storage [START_REF] Maurer | Room-Temperature Quantum Bit Memory Exceeding One Second[END_REF] and metamaterials [START_REF] Susstrunk | Observation of phononic helical edge states in a mechanical topological insulator[END_REF][START_REF] Nash | Topological mechanics of gyroscopic metamaterials[END_REF] devices. Over the last decade, experimental realizations of ultra-cold fermionic lattice gases [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF][START_REF] Cheuk | Quantum-Gas Microscope for Fermionic Atoms[END_REF] have spurred advances in the understanding of charge and spin transport, promising novel technological platforms for quantum computation [START_REF] Buluta | Natural and artificial atoms for quantum computation[END_REF][START_REF] Gross | Quantum simulations with ultracold atoms in optical lattices[END_REF] and spintronics [START_REF] Dieny | Opportunities and challenges for spintronics in the microelectronics industry[END_REF]. In parallel, recent progress in the control of natural and artificial active matter has led to the discovery of emergent collective spin states in chemical [START_REF] Biktashev | Orbital Motion of Spiral Waves in Excitable Media[END_REF][START_REF] Totz | Spiral wave chimera states in large populations of coupled chemical oscillators[END_REF], biological [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF][START_REF] Nishiguchi | Engineering bacterial vortex lattice via direct laser lithography[END_REF][START_REF] Riedel | A Self-Organized Vortex Array of Hydrodynamically Entrained Sperm Cells[END_REF] and synthetic [START_REF] Souslov | Topological sound in active-liquid metamaterials[END_REF] suspensions, inspiring the invention of chiral active metamaterials [START_REF] Soni | The odd free surface flows of a colloidal chiral fluid[END_REF]. Here, we demonstrate a new class of hydrodynamic spin lattice (HSL) systems [START_REF] Sáenz | Spin lattices of walking droplets[END_REF] that combines active self-propulsion and particle-wave dynamics.

A HSL consists of an array of in-phase bouncing droplets [START_REF] Couder | Dynamical phenomena: walking and orbiting droplets[END_REF], each drop confined by a submerged circular well and orbiting on the surface of a vertically vibrated liquid bath (Fig. 1a,b; Movie S1). Despite being macroscopic classical objects, these walking droplets display wave-mediated self and pair interactions reminiscent of those arising in microscopic systems [START_REF] Bush | Pilot-wave hydrodynamics[END_REF]. To demonstrate the tunability of HSLs, we induce a range of collective 'magnetic' ordering phenomena through variations in the vibrational forcing and lattice geometry, and by mimicking the effect of an applied magnetic field through imposition of system rotation [START_REF] Fort | Path-memory induced quantization of classical orbits[END_REF]. We explain the experimentally observed particlewave dynamics by deriving a mapping onto a reduced oscillator network description [START_REF] Strogatz | From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF][START_REF] Lindner | Effects of noise in excitable systems[END_REF], generalizing models used to describe neuronal networks [START_REF] Cumin | Generalising the Kuramoto model for the study of neuronal synchronisation in the brain[END_REF] and other synchronization phenomena [START_REF] Matheny | Exotic states in a simple network of nanoelectromechanical oscillators[END_REF][START_REF] Dörfler | Synchronization in complex networks of phase oscillators: A survey[END_REF]. This mapping establishes HSLs as a versatile platform for investigating chiral symmetry breaking and synchronization phenomena far from equilibrium.

In our experiments, both the droplets and underlying liquid bath consist of the same 20 cSt silicon oil (Methods, Experiments). An electromagnetic shaker vibrates the bath vertically with forcing acceleration Γ(t) = γ cos(2πf t), where γ is the maximum acceleration, f the oscillation frequency, and t time (Fig. 1b). For sufficiently strong acceleration, the air layer between the drop and bath prevents direct contact, thus enabling persistent bouncing [START_REF] Protière | Particle-wave association on a fluid interface[END_REF]. When the drop bounces at half the driving frequency, so achieves resonance with the bath's most unstable wave mode, it excites radially decaying circular waves of characteristic wavelength λ F on the bath surface [START_REF] Protière | Particle-wave association on a fluid interface[END_REF][START_REF] Eddi | Information stored in Faraday waves: the origin of a path memory[END_REF][START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF]. The amplitude and spatio-temporal decay of these waves can be controlled through the driving [START_REF] Eddi | Information stored in Faraday waves: the origin of a path memory[END_REF]. When the driving amplitude γ is increased beyond a critical walking threshold γ w , vertical bouncing becomes unstable to horizontal walking [START_REF] Protière | Particle-wave association on a fluid interface[END_REF][START_REF] Bush | Pilot-wave hydrodynamics[END_REF]. As walking droplets land on the slopes of the underlying wavefield η, they experience a wave-induced horizontal force F i (x i , t) ∼ -∇η| x=xi proportional to the wave slope ∇η at the point of impact x i (t). Free droplets thus walk at constant speed along straight lines when γ > γ w [START_REF] Protière | Particle-wave association on a fluid interface[END_REF].

To transform linear into circular motion, we confined each walker to a submerged circular well [START_REF] Sáenz | A hydrodynamic analog of Friedel oscillations[END_REF] (Fig. 1b, Methods, Experiments). The resulting clockwise (CW) or anti-clockwise (ACW) motions were characterized by measuring the specific spin S(t) = L z (t)/m, where m is the droplet's mass and L z denotes its angular momentum with respect to the well center. 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 the forcing acceleration γ towards the Faraday threshold γ F , at which the bath surface becomes unstable in the absence of droplets [START_REF] Eddi | Information stored in Faraday waves: the origin of a path memory[END_REF][START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF]. These single-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 γ decreases [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF]. 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 this understanding of the single-spin states, we investigated the collective spin dynamics in one-(1D) and two-dimensional (2D) HSLs.

The wave-mediated spin-spin coupling in HSLs is reminiscent of spatially oscillating interactions [START_REF] Ruderman | Indirect Exchange Coupling of Nuclear Magnetic Moments by Conduction Electrons[END_REF][START_REF] Coey | Magnetism and Magnetic Materials. 3rd printing[END_REF] in RKKYtype 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 adjacent wells (Fig. 1b; Methods, Experiments). When the pair coupling is sufficiently strong, nearest-neighbor 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) = i S i (t) and spin-spin correlation χ(t) = i∼j S i (t)S j (t), where i∼j denotes a sum over adjacent pairs. To achieve statistical significance, experiments were run for several hours: In 1 hour a droplet performs ∼ 10 5 bounces and ∼ 1, 800 orbits (Methods, Statistics). Positive values of χ signal parallel alignment of neighboring spins ('ferromagnetic' order) whereas negative values of χ indicate anti-parallel alignment ('anti-ferromagnetic' order).

The oscillatory wave-mediated spin-spin coupling suggests that HSLs can support different types of collective order depending on the ratio between lattice spacing L and Faraday wavelength λ F . To test this hypothesis, 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, Statistics), we observed that pair interactions 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, indicated that global mirror symmetry was preserved (Fig. 1e). However, the negative pair correlation, χ < 0, revealed a bias towards local anti-ferromagnetic order (Fig. 1g; Movie S2). Holding L/λ F fixed, a similar anti-ferromagnetic bias occurred in experiments with non-periodic boundary conditions, different N values, and different lattice radii R (Methods, Statistics), confirming that anti-ferromagnetic ordering is selected by the lattice spacing. The strength of the emergent anti-ferromagnetic order depends non-monotically on the driving amplitude γ (Fig. 1h; Methods, Statistics). The strongest collective anti-ferromagnetic response was observed for γ ≈ γ c , indicating a correlation between the global collective ordering and the robustness of a single spin state (Fig. 1c).

The emergence of local magnetic ordering can be understood by analyzing the drops' horizontal phase φ i , which plays a role analogous to the polar angle in XY -type spin models (Fig. 1a). To quantify pairwise phase synchronization, we recorded the mean phase differences ∆φ ± for co-rotating (+) and counter-rotating (-) neighboring pairs, respectively (see Methods, Statistics). These observables signal in-phase (∆φ ± → 0) or out-of-phase (∆φ ± → π) rotations of neighboring 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 lattice geometry, 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; Movie S3), with a bias towards in-phase rotation (Fig. 2c). As in the anti-ferromagnetic 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 spin order and phase synchronization in HSLs, we derived a generic phase oscillator model from a detailed hydrodynamic description of our experimental system (Fig. 2e; Methods, Theory). In this model, we focus on subcritical driving γ ≤ γ c (Fig. 1c), so that each walker i follows a circular trajectory x i (t) = r 0 (cos φ i (t), sin φ i (t)), with angular frequency ω i = φi , where (r 0 , φ) are the cylindrical coordinates with respect to the center of each well. The walking speed v 0 then defines a preferred orbital frequency ±ω 0 = ±v 0 /r 0 [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF]. Modeling this effect by a nonlinear Rayleigh-type friction, the phase dynamics of a walker is described by ωi

= 1 τ 1 - ω 2 i ω 2 0 ω i + j F ij (φ i , φ j ). (1)
Equation ( 1) can be interpreted as an active XY -type model [START_REF] Coey | Magnetism and Magnetic Materials. 3rd printing[END_REF], in which perturbations to ω 0 decay over the time-scale τ , and F ij is the wave-mediated force exerted on drop i by its nearest neighbor j (Fig. 2e). For weak accelerations and small orbital radii r 0 L, the interaction force F ij (φ i , φ j ) = -∂U/∂φ i can be derived from the effective coupling potential U (φ i , φ j ) = 1 2 α (cos φ i -cos φ j ) 2 + 1 2 β (sin φ i -sin φ j ) 2 (Methods, Theory). The coefficients are proportional to the curvature, α ∼ J 0 (k F L), and gradient, β ∼ J 0 (k F L), of the coupling wave field, which is well approximated by a Bessel function of the first kind J 0 (|x -x j |) centered at the neighboring drop position. HSLs characterized by purely circular orbits, as arise for γ ≤ γ c , thus belong to a generalized Kuramoto universality class [START_REF] Strogatz | From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF].

Analysis of the phase-oscillator dynamics (1) confirms that the spin-spin correlation alternates between ferroand anti-ferromagnetic 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 anti-ferromagnetic 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 anti-ferromagnetic (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, and Extended Data Fig. 5).

The spins in an anti-ferromagnetic material may be realigned into a ferromagnetic state by imposing a constant magnetic field B [START_REF] Coey | Magnetism and Magnetic Materials. 3rd printing[END_REF]. Can HSLs undergo similar global symmetry breaking? The Lorentz force F B = q(v × B) acting on a charge q moving with velocity v has precisely the same form as the Coriolis force F Ω = m(v × 2Ω) acting on a mass m moving with velocity v in a frame rotating at constant angular rate Ω [START_REF] Fort | Path-memory induced quantization of classical orbits[END_REF]. This suggests the possibility of transforming anti-ferromagnetic HSLs into ferromagnetic states by rotating the vibrating bath about a vertical axis (Methods, Experiments). To test this hypothesis, we initialized the anti-ferromagnetic 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 anti-ferromagnetic state as expected (Fig. 3a,b, t < 0.5 h). At t = 0.5 h, we began to rotate the setup anti-clockwise by gradually increasing the rotation rate to Ω = 0.7 rad/s. In a matter of seconds, most spins had co-aligned with the bath rotation leading to a strongly positive spin correlation, χ ≈ 0.58. However, unlike the steady HSL in Fig. 2ad, the rotated HSL had a net magnetization, M ≈ 0.78. Moreover, when the rotation direction was reversed to Ω = -0.7 rad/s (CW) 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 the spin pairs remained ferromagnetically arranged, χ > 0 (Fig. 3b; Movie S4). The transition to ferromagnetic order requires a super-critical rotation rate |Ω| > Ω c ≈ 0.22 rad/s (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 the mechanism responsible for the field-induced polarization (Fig. 3f-g). For optimal driving γ ≈ γ c , a walker executes circular motion in the absence of rotation (Fig. 1c). When rotation is switched on, a drop co-orbiting with the bath feels a Coriolis force F Ω that opposes the confining force of the submerged well F w [START_REF] Sáenz | A hydrodynamic analog of Friedel oscillations[END_REF] (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 with respect to the bath, both F w and F Ω point inwards (Fig. 3g), destabilizing the circular orbits into trefoil-like trajectories similar to those observed for Ω = 0 and γ > γ c (Fig. 1c). Droplets on these trefoil-like trajectories pass closer to the well centers (Fig. 3h), where their walking directions can be more easily reversed by perturbations from their neighbors. Hence, confinement favors co-rotating orbits in HSLs; conversely, free walkers prefer counter-rotating inertial orbits [START_REF] Fort | Path-memory induced quantization of classical orbits[END_REF]. We note that a generalization of our phase-oscillator model captures the effects of the external field (Fig. 3e; Methods, Theory).

Two-dimensional classical and quantum spin lattices display features that are absent in their 1D counterparts, including geometric frustration and topological order [START_REF] Coey | Magnetism and Magnetic Materials. 3rd printing[END_REF]. HSLs provide a promising platform for exploring such effects at the macroscale. For example, square HSLs (Fig. 4a; Movie S5) that promote anti-ferromagnetic order in the absence 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 the interactions between lattice sites can be directly tailored in various ways (Methods, Tunablity), we expect HSLs to provide a versatile testbed for designing S3). (e) Geometry considered in the theoretical wave model. Neighboring drops i and j move along circular trajectories of radius r at angular speed ±ω. The drop position is determined by the angular phase φ, defined with respect to the lattice centerline. Drop i is subject to a force F ij along its direction of motion due to the wave field generated by drop j, as is approximated by a Bessel function J 0 centered at its position.

(f) Depending on the drop separation, four dominant drop-pair synchronization modes are favored, leading to the emergence of alternating ferromagnetic and anti-ferromagnetic local order as the lattice spacing is varied (Methods, Simulations).

f g

Fig. 3: Inducing global polarization through applied rotation. Rotating the circular HSL from Fig. 1d-h induces a transition from effective anti-ferromagnetic 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 anti-ferromagnetic order with no net magnetization. At t = 0.5 h, the bath rotation was switched on (Ω = 0.7 rad/s, ACW) inducing a transition to ferromagnetic order and positive magnetization. At t = 1 h, the rotation direction was inverted (Ω = -0.7 rad/s, CW) leading to negative magnetization. Lattice parameters and forcing frequency are the same as in Fig. 1 (f,g) Force balance and measured trajectory, color-coded by speed, for an isolated spin in a frame corrotating with the bath. The drop's trajectory is nearly circular when its angular velocity ω and the bath rotation Ω are aligned, but approximates a precessing trefoil when the signs are opposite. (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; r 0 is the trajectory radius for Ω = 0 shown in Fig. 1c. The center line denotes the mean radius, and the envelope indicates the standard deviation.

Introduction

Walking droplets in free space self-propel in a rectilinear fashion [START_REF] Couder | Dynamical phenomena: walking and orbiting droplets[END_REF][START_REF] Protière | Particle-wave association on a fluid interface[END_REF][START_REF] Eddi | Information stored in Faraday waves: the origin of a path memory[END_REF][START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF]. Hydrodynamic spin states [START_REF] Anand | Hydrodynamic spin states[END_REF], in which a droplet executes circular orbits confined by its own wave field are unstable in the laboratory [START_REF] Oza | A trajectory equation for walking droplets: hydrodynamic pilotwave theory[END_REF][START_REF] Labousse | Self-attraction into spinning eigenstates of a mobile wave source by its emission backreaction[END_REF], unless confined by an external force such as a linear spring force [START_REF] Perrard | Self-organization into quantized eigenstates of a classical wave-driven particle[END_REF] or the Coriolis force that acts in a rotating frame [START_REF] Fort | Path-memory induced quantization of classical orbits[END_REF][START_REF] Eddi | Level Splitting at Macroscopic Scale[END_REF][START_REF] Harris | Droplets walking in a rotating frame: from quantized orbits to multimodal statistics[END_REF][START_REF] Oza | Pilot-wave dynamics in a rotating frame: on the emergence of orbital quantization[END_REF]. Recent experiments have demonstrated that variable bottom topography may be exploited to realize new hydrodynamic quantum analogs [START_REF] Sáenz | A hydrodynamic analog of Friedel oscillations[END_REF][START_REF] Saénz | Statistical projection effects in a hydrodynamic pilot-wave system[END_REF]; in particular, the influence of a submerged well has led to analogs of the so-called quantum mirage [START_REF] Saénz | Statistical projection effects in a hydrodynamic pilot-wave system[END_REF] and Friedel oscillations [START_REF] Sáenz | A hydrodynamic analog of Friedel oscillations[END_REF]. In this work, we exploit variable bottom topography to stabilize a collection of interacting hydrodynamic spin states [START_REF] Sáenz | Spin lattices of walking droplets[END_REF], and characterize the emergent collective order. We note that, while several workers have considered the stability of stationary bouncing droplet lattices [START_REF] Eddi | Archimedean lattices in the bound states of wave interacting particles[END_REF][START_REF] Eddi | Oscillating instability in bouncing droplet crystals[END_REF][START_REF] Thomson | Collective vibrations of confined levitating droplets[END_REF][START_REF] Thomson | Collective vibrations of a hydrodynamic active lattice[END_REF][START_REF] Couchman | Free rings of bouncing droplets: stability and dynamics[END_REF], and Nachbin [START_REF] Nachbin | Walking droplets correlated at a distance[END_REF][START_REF] Nachbin | Kuramoto-Like Synchronization Mediated through Faraday Surface Waves[END_REF] examined wave-induced correlations of walker pairs confined to adjacent one-dimensional cavities, we present here the first macroscopic analog of electronic spin systems with orbiting droplets.

Experiments

Stationary bath

A schematic of the experimental setup is presented in Extended Data Fig. 1. A circular aluminum bath was filled with silicon oil with density ρ = 950 kg m -3 , viscosity ν = 20.9 cSt, and surface tension σ = 20.6 mN m -1 . The lattices of submerged wells were constructed by laser cutting acrylic sheets of thickness h = 6.2 mm, which were then bolted 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 regime [START_REF] Whitham | Linear and Nonlinear Waves[END_REF]. Similarly, the minimum distance d between the bath's outer wall and the outermost well was sufficiently large, d > 6λ F , to ensure that the meniscus waves [START_REF] Douady | Experimental study of the Faraday instability[END_REF] originated at the bath wall had no effect on the experiments. We note that the Faraday wavelengths corresponding to the shallow λ H F and deep λ h F regions are only marginally different, specifically λ h F /λ H 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 (Data Physics, 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 ensures a spatially uniform vibration to within 0.1%. [START_REF] Harris | Generating uniaxial vibration with an electrodynamic shaker and external air bearing[END_REF] 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 ensures a constant acceleration amplitude to within ±0.002g. [START_REF] Harris | Generating uniaxial vibration with an electrodynamic shaker and external air bearing[END_REF] We define the Faraday threshold γ F as the critical vibrational acceleration above which Faraday waves [START_REF] Faraday | On a Peculiar Class of Acoustical Figures and on Certain Forms Assumed by Groups of Particles upon Vibrating Elastic Surfaces[END_REF] appear above the wells. The experimental setup was left vibrating at γ ≈ γ F for at least 1 h before running any experiment to ensure a steady state 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 temperature changes. With this protocol, the standard day-to-day variation in γ F was limited to ±0.04g, and the variation in each experimental run to |∆γ F |/γ F < 0.003 , where ∆γ F represents the variation in the Faraday threshold over a 1-hour segment.

Droplets of the same silicon oil and diameter D d = 0.75 ± 0.01 mm were generated with a piezoelectric dropleton-demand generator [START_REF] Harris | A low-cost, precise piezoelectric droplet-on-demand generator[END_REF], 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 was recorded 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 by recording the normal reflection of light at the free surface. [START_REF] Douady | Experimental study of the Faraday instability[END_REF] 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 bright regions corresponding to horizontal parts of the surface, extrema or saddle points [START_REF] Douady | Experimental study of the Faraday instability[END_REF]. To ensure that ambient air currents did not affect the results, the vibrating bath was enclosed within a transparent acrylic chamber. Experimental results obtained using the acrylic chamber were compared 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 the effective noise in the system was slightly reduced with the lid, but that the emergent collective behavior was statistically indistinguishable.

The experiments were performed with walking drops bouncing subharmonically [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF] in resonance with the underlying Faraday wave field, each with the same vertical bouncing phase. To ensure in-phase bouncing, before each experiment the bath acceleration was reduced to γ = 1.1g, so that the droplets transitioned to harmonic bouncers [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF], synchronized with the bath vibration. By observing the light reflected from 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 behavior emerged (Fig. 1h).

Rotating bath

To subject hydrodynamic spin lattices to bulk rotation, we mounted the aluminum bath on a vertical shaft attached to the upper plate of the linear air bearing through a system of preloaded ball bearings. The freely rotating shaft was connected with a vibration damping coupling to a DC electric motor with a planetary gearbox (Phidgets Inc., 12V DC motor, 139:1 gearbox) housed inside the hollow air bearing (See Extended Data Fig. 1). The rotation rate was monitored with an optical rotary encoder (Phidgets Inc., HKT22) and the electric motor controlled with an in-house Python program, 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

H 0 (r) = H - Ω 2 R 2 b 4g + Ω 2 2g r 2 , ( 2
)
where r is the radial coordinate with respect to the center of the bath, H is the depth of the fluid at rest, and R b the bath radius, which in the rotating experiments was R b = 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 ensure that H 0 (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 microns, and so 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 anti-ferromagnetic to ferromagnetic order (Ω < 0.7 rad/s, Fig. 4e), the liquid height was reduced in the center and increased in the periphery by < 10% with respect to the depth of fluid at rest.

Lattice parameters

The lattice parameters for the 1D anti-ferromagnetic lattice discussed in Fig. 1 We note that for the 1D ferromagnetic lattice, we reduced the well diameter D and liquid 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 k F , was also fine-tuned to accentuate the ferromagnetic response. The lattice parameters for the 2D square lattice 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 behavior 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 S i (t) i |S i (t)| , (3)
where i denotes the sum over all N spins, and S i (t) = L z,i (t)/m is the specific spin of droplet i with mass m and angular momentum L z,i with respect to the well center. The denominator is the normalizing factor such that M ∈ [-1, 1]. Positive (negative) values of M (t) signal a larger number of walkers rotating ACW (CW). Pairwisesymmetry breaking is quantified through the normalized spin-spin correlation,

χ(t) = i∼j S i (t)S j (t) i∼j |S i (t)S j (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, anti-ferromagnetic order) among neighboring walkers.

To quantify the relationship between the emergent collective behavior and the pairwise orbital phase synchronization of adjacent walkers, we compute separately the mean phase difference ∆φ ± ∈ [0, π] for co-rotating (+) and counter-rotating (-) walkers as

∆φ ± (t) = 1 N p± i∼j± |φ i (t) ∓ φ j (t)|, (5)
where φ i is the instantaneous horizontal phase with respect to the lattice centerline (Fig. 1a), and i∼j± denotes the sum over N p± adjacent co-rotating or counter-rotating 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 statistically steady-state. To guarantee random initial conditions, we manually stirred with a needle all the drops at the beginning of each run. By checking that the spin correlation started from 0 at the beginning of each experimental run, we confirmed that our stirring procedure was adequate. We also checked that the stirring did not disrupt the vertical synchronization 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, with components S i (t), that characterizes the lattice spin arrangement at time 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 depending on the experimental realization, we define the length of the 'statistically-independent 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 [START_REF] Soni | The odd free surface flows of a colloidal chiral fluid[END_REF] 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 the signal-to-driving ratio for S(t) is maximal (Fig. 1c), and the emergent collective behavior strongest (Fig. 1h). Based on this observation, we selected 1 h (∼ 1, 800 circulations for each drop) as the recording length for 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 and after each realization. Following this protocol, we recorded as many 1-h videos as needed to obtain a minimum of N I ≥ 8 statistically independent intervals for each γ/γ F level. In some cases the number of independent intervals was as high as N I ≈ 50.

For each independent interval we then computed its associated spin-spin correlation χ j (t). The data reported in Fig. 1h corresponds to the resulting average, χ j (t)/N I , and standard deviation, σ( χ j (t)/N I ), 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 and 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 BCs, we considered the anti-ferromagnetic configuration with N = 20 shown in the main text (Fig. 1), and repeated the experiment at γ/γ F = 86% after plugging 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) over which the waves dissipated quickly. Hence, the walkers at the ends did not interact with each other, as was checked by direct visualization of 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 BCs. We also note that for 1D periodic lattices tuned for preferred anti-ferromagnetic order, it is important whether N is even or odd. An odd number of spins may lead to geometrical frustration effects [START_REF] Coey | Magnetism and Magnetic Materials. 3rd printing[END_REF] that hinder the preferred collective order. To avoid such effects, whose investigation 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 of experiments led to preferred anti-ferromagnetic 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 important factor for the emergent collective behavior. 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.

Theory

Rayleigh-Bessel oscillator model

To rationalize the emergent order in HSLs, we extend the modeling framework for pilot-wave dynamics introduced by Bush et al. [START_REF] Bush | The wave-induced added mass of walking droplets[END_REF]. In the weak-acceleration limit, | v|T M /|v| 1, these authors show that the trajectory of a walking droplet, x(t), may be described by the equation

d dt (m w v) + D w v = F, (8)
where v = ẋ is the droplet velocity, F an externally imposed force acting on the drop, D w (|v|) the nonlinear drag, and m w the effective walker mass. In the low-velocity limit, in which the wave decay time, T M , is much smaller than the time taken for the droplet to travel one Faraday wavelength, T M λ F /u 0 , one may further simplify (8) to a Rayleigh-oscillator-type equation [START_REF] Turton | Theoretical modeling of pilot-wave hydrodynamics[END_REF], similar to that introduced by Labousse & Perrard [START_REF] Labousse | Non-Hamiltonian features of a classical pilot-wave dynamics[END_REF], where the droplet trajectory is described by

m w dv dt = D R 1 - |v| 2 u 2 0 v + F. (9)
Here, D R denotes the magnitude of the non-linear Rayleigh-type drag, which drives the droplet towards its free walking speed, u 0 . The dependence of the model parameters m w , u 0 and D R 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 T M ∼ 0.1 s. While these limits are not strictly satisfied in the relatively short periods of time during which spin flips occur, our analysis will demonstrate that details of the spin flips are not essential to rationalize the emergent collective behavior.

We consider N droplets self-propelling above an array of submerged circular wells with center-to-center separation L. Equation ( 9) determines the trajectory of each droplet, which is defined by a two-component displacement vector x i (t), where i = 1, ..., N . The force acting on each droplet may be decomposed into two components, F = F w + F i,j , where F w represents the effective confining force exerted by the submerged well, and F i,j the wave force exerted on droplet i by its nearest neighboring droplet j. Motivated by experiments [START_REF] Cristea-Platon | Walking Droplets in a Circular Corral: Quantisation and Chaos[END_REF], we assume that F w is a central force that points towards the well's center.

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 one-dimensional lattices for the sake of simplicity, we may thus express the droplet trajectory and velocity vectors as x i (t) = (r cos φ i (t) + (i -1)L, r sin φ i (t)) = rr i + (i -1)Lx, and ẋi (t) = r φi (-sin φ i (t), cos φ i (t)) = r φi θi , respectively. Here, ri and θi denote the radial and azimuthal unit vectors for the i-th spin with respect to the center of its well, and x the unit vector along the lattice center line (Fig. 2e). Writing ω i = φi and substituting x i (t) and ẋi (t) into (9) yields the following equations for the normal and tangential force balances:

-m w rω 2 i = F • ri , (10)
m w r dω i dt -D R 1 - r 2 ω 2 i u 2 0 rω i = F i,j • θi . (11)
We note that F w only appears in the radial force balance, where it acts to select the orbital radius r. Since the experiments suggest that emergent order in HSLs is driven primarily by synchronization of the droplets' relative orbital phase (Fig. 1f and2c), we restrict our attention to the azimuthal force balance [START_REF] Weinberg | The Quantum Theory of Fields[END_REF] and simply infer the value of r from our experiments. Thus, to rationalize preferred mean phase differences, ∆φ ± , in terms of the droplet coupling force, F i,j , we may write equation [START_REF] Weinberg | The Quantum Theory of Fields[END_REF] as a non-linear differential equation for the evolution of each droplet's angular frequency,

ω i (t), ωi = 1 τ 1 - ω 2 i ω 2 0 ω i + F i,j • θi m w r , (12)
where ω 0 = u 0 /r is the preferred angular frequency for a given driving acceleration γ, as would arise for a single drop orbiting in isolation, and τ = m w /D R is the relaxation timescale (Extended Data Table 1).

To complete our spin model, we must define the wave-mediated force exerted by droplet j on droplet i, F i,j . The stroboscopic model of Oza et al. [START_REF] Oza | A trajectory equation for walking droplets: hydrodynamic pilotwave theory[END_REF] predicts that the wave field generated by a walker may be approximated by

η(x j , t) = A T F t -∞ J 0 (k F |x -x j (s)|)e -(t-s)/T M ds, (13)
which indicates that the coupling force may be expressed as

F i,j = F T F t -∞ J 1 (k F |x i (t) -x j (s)|) |x i (t) -x j (s)| (x i (t) -x j (s)) e -(t-s)/T M ds. ( 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 regions [START_REF] Tadrist | Faraday instability and subthreshold Faraday waves: surface waves emitted by walkers[END_REF]. Nevertheless, our visualization of the wave field produced by a single droplet in the lattice yielded a waveform notably similar in form to that expected on the basis of ( 14) (Section 4.3.2). We also note that there is a relatively fast-propagating wave front generated at impact [START_REF] Eddi | Information stored in Faraday waves: the origin of a path memory[END_REF], 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 [START_REF] Galeano-Rios | Ratcheting droplet pairs[END_REF][START_REF] Couchman | Bouncing phase variations in pilot-wave hydrodynamics and the stability of droplet pairs[END_REF]. Guided by previous studies and experimental evidence, we thus neglect this capillary wave front and consider that the liquid layer between wells has a negligible influence of the form of the drops' wave fields and so seek a simplification of ( 14) that yields insight into the mechanism responsible for the emergent order.

In the parameter regime of interest, γ ≤ γ c , the wave memory time, T M < 0.1 s, is much smaller than the orbital period, 2π/ω 0 ∼ 2 s. We may thus simplify ( 14) in the limit of T M 2π/ω 0 to obtain

F i,j ∼ F T M T F J 1 (k F |x i (t) -x j (t)|) (x i (t) -x j (t)) |x i (t) -x j (t)| . (15)
Substituting [START_REF] Susstrunk | Observation of phononic helical edge states in a mechanical topological insulator[END_REF] into equation ( 12) yields a system of ODEs for each droplet interacting with its nearest neighbors in a spin lattice

ωi = 1 τ 1 - ω 2 i ω 2 0 ω i + F j J 1 (k F |x i -x j |) (x i -x j ) |x i -x j | • (-sin φ i , cos φ i ) . (16)
Here, F is the interaction force coefficient, and j denotes the sum over nearest neighbors 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 , F) that may be deduced from the system's fluid properties, in addition to the lattice spacing, L, and the Faraday wavenumber, k F . Expressions for the parameters τ , ω 0 and F are provided in Extended Data Table 1, and their dependence on γ/γ F is presented in Supplementary Information 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 model [START_REF] Oza | A trajectory equation for walking droplets: hydrodynamic pilotwave theory[END_REF], or the variable-topography model of Faria [START_REF] Faria | A model for Faraday pilot waves over variable topography[END_REF].

Generalized Kuramoto (GK) model

By making a series of further approximations to the form of the coupling in the Bessel model [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF], we now demonstrate that the emergent collective order in HSLs is governed by a generalization of the classical Kuramoto model for coupled phase oscillators [START_REF] Strogatz | From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators[END_REF]. From ( 16), the force on droplet i due to the adjacent droplet (i + 1) is

F i,i+1 = -FJ 0 (k F |x i (t) -x i+1 (t)|) (x i (t) -x i+1 (t)) |x i (t) -x i+1 (t)| . ( 17
)
We first make the approximation that the orbital radius is small relative to the center-to-center lattice separation, r/L 1. This approximation is well satisfied in the experiments, where r ≤ 3mm and L > 14mm, so r/L ≤ 0.2. We can thus approximate the distance between the two droplets as |x i (φ i ) -x i+1 (φ i+1 )| ≈ L + r (cos φ i+1 -cos φ i ), and expand the azimuthal component of the wave-induced force:

θi • F i,i+1 ≈ -F [sin φ i + r/L(sin φ i -sin φ i+1 ) cos φ i ] J 0 (k F L -k F r(cos φ i -cos φ i+1 )) . ( 18
)
We then make the additional approximation that

k F r |(cos φ i -cos φ i+1 )| = 2k F r sin φ i -φ i+1 2 sin φ i + φ i+1 2 1 ( 19
)
which may be written equivalently as

x i -x i+1 λ F 2π , (20)
where (x i , y i ) are the Cartesian coordinates for each drop relative to the well center. While this limit will not be satisfied by all solutions of equation ( 16), for the range of L considered in our experimental study, the droplets exhibit in-phase synchronization (∆φ ± ≈ 0; Fig. 1f and2c); consequently, the observed dynamics are expected to satisfy the approximation in equation [START_REF] Gross | Quantum simulations with ultracold atoms in optical lattices[END_REF]. Exploiting [START_REF] Gross | Quantum simulations with ultracold atoms in optical lattices[END_REF], we thus deduce that

J 0 (k F L -k F r(cos φ i -cos φ i+1 )) ≈ J 0 (k F L) -k F rJ 0 (k F L)(cos φ i -cos φ i+1 ), (21)
from which [START_REF] Buluta | Natural and artificial atoms for quantum computation[END_REF] indicates that θi

• F i,i+1 ≈ -FJ 0 (k F L) sin φ i +Fk F rJ 0 (k F L)(cos φ i -cos φ i+1 ) sin φ i - Fr L J 0 (k F L)(sin φ i -sin φ i+1 ) cos φ i . ( 22
)
Noting that the first term in ( 22) is cancelled when we sum the contribution of the forces from the two neighboring droplets, i -1 and i + 1, we may now express the coupling force as the gradient of a potential U , yielding a reduced model for the spin lattice system as

ωi = 1 τ 1 - ω 2 i ω 2 0 ω i - j ∂U i,j ∂φ i , ( 23
)
where the sum is only over nearest neighbors within the lattice. The generalized Kuramoto coupling potential, U i,j , may be expressed as

U i,j := U (φ i , φ j ) = α 2 | cos φ i -cos φ j | 2 + β 2 | sin φ i -sin φ j | 2 = α 2r 2 |x i -x j | 2 + β 2r 2 |y i -y j | 2 , ( 24
)
where

α = Fk F rJ 0 (k F L), ( 25
) β = (Fr/L)J 0 (k F L). (26)
We refer to this as the 'Generalized Kuramoto' (GK) model, and illustrate the dependence of the ratio β/α on lattice spacing, L, in Supplementary Information Fig. 2. We note that the form of the potential U is prescribed only by the local wave curvature and gradient, through the signs and relative magnitudes of α and β. Consequently, the amplitude of the local wavefield does not significantly alter the form of the potential, U , which further justifies the neglect of topography-induced wave damping. We note that our GK model [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF] includes inertia, as do the Kuramoto-like models typically used to study synchronization phenomena in power grids [START_REF] Filatrella | Analysis of a power grid using a Kuramoto-like model[END_REF][START_REF] Rohden | Self-Organized Synchronization in Decentralized Power Grids[END_REF]. Notably, we found that the GK model succeeds in rationalizing the emergent order even at intermediate values of lattice spacing, L, for which the approximation [START_REF] Gross | Quantum simulations with ultracold atoms in optical lattices[END_REF] is not necessarily satisfied. To understand this unforeseen success of the GK model, we compare the form of the Bessel coupling force [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF] with that expected from the GK model ( 23) for all φ i and φ j and four representative lattice separations (Supplementary Information Fig. 3). In all four cases, we find that the GK coupling corresponds to the two fundamental Fourier modes of the Bessel coupling, which explains the agreement between the predictions of the Bessel and GK models over the range of lattice spacings of interest. For even 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 rk F 1, when the orbital radius is much smaller than the Faraday wavelength, the approximation ( 19) is strictly satisfied. Finally, we note that the amplitude of the GK coupling force computed via the predictions ( 25)-( 26) is a factor of two larger than the amplitude of the dominant Fourier modes computed from the Bessel force (Supplementary Information 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 [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF].

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 not depend strongly on either the precise details of the spin flips or the small deformation of the walker wave field due to variations in bottom topography. Nevertheless, we detail these effects here for the sake of completeness.

Spin flips

Extended Data Fig. 2 (Movie S6) illustrates a typical spin flip event as observed in our experiments with the 1D anti-ferromagnetic HSL described in Fig. 1. Due to transient rearrangements from the random initial conditions, at some time t during this experimental run, three neighboring spins happen to be rotating in the same direction (Extended Data Fig. 2a). Since the lattice geometry promotes anti-ferromagnetic order, this local arrangement is not a preferred state. The collective interactions are such that the trajectory of the middle walker then starts to become elliptical (Extended Data Fig. 2b). The elongation of the ellipse continues until the walker executes quasi-rectilinear motion across the well center (Extended Data Fig. 2c). Subsequently, the process is reversed, the result being the walker rotating in the opposite direction (Extended Data Fig. 2d). After the spin flip, which typically takes 5 ± 1 s, the three walkers rotate in the preferred anti-ferromagnetic configuration (Extended Data Fig. 2e). When the collective order is strongest (γ/γ F ∼ 86% in Fig. 1h), the duration of a flip, O(5 s), is negligible relative to the characteristic flipping period, O(5 -10 min). At higher accelerations, when the walkers switch direction much more frequently and the 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 variable bottom topography will alter the form of the wave field relative to that predicted by the stroboscopic model ( 13) may be obtained by comparing the Faraday wavelengths given by the water-wave dispersion relation for the shallow λ H F and deep λ h F regions. As noted earlier in the experimental description, these wavelengths are only marginally different for our choice of parameters, specifically λ h F /λ H F ≈ 1.05. One thus expects that the walker wave field in our lattices will not be significantly different from that predicted by Oza et al. [START_REF] Oza | A trajectory equation for walking droplets: hydrodynamic pilotwave theory[END_REF]. Indeed, as shown in Extended Data Fig. 3a-b, direct visualization of waves generated by a single walker reveals a wave form comparable to that expected from a superposition of sub-harmonic zeroth-order Bessel functions, as in [START_REF] Chumak | Magnon spintronics[END_REF]. We note that over the range of parameters considered in our study, no exotic Faraday patterns were excited, such as those arising in nonlinear regimes [START_REF] Chen | Subcritical Bifurcations and Nonlinear Balloons in Faraday Waves[END_REF][START_REF] Périnet | Alternating Hexagonal and Striped Patterns in Faraday Surface Waves[END_REF] or excited by droplets walking in elliptical corrals [START_REF] Saénz | Statistical projection effects in a hydrodynamic pilot-wave system[END_REF].

To provide a quantitative characterization of the wave coupling, including the spatial decay [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF], we perform additional simulations with the model developed by Faria [START_REF] Faria | A model for Faraday pilot waves over variable topography[END_REF], wherein changes in topography are incorporated through their influence on the local wave speed. We note that this model solves the wave problem explicitly, and has been bench-marked against experiments in a number of settings involving walkers interacting with submerged topography [START_REF] Sáenz | A hydrodynamic analog of Friedel oscillations[END_REF][START_REF] Pucci | Non-specular reflection of walking droplets[END_REF][START_REF] Pucci | Walking droplets interacting with single and double slits[END_REF][START_REF] Harris | The interaction of a walking droplet and a submerged pillar: From scattering to the logarithmic spiral[END_REF]. In accordance with our experimental observations, we find that the simulated wave field has a form similar to that of a J 0 Bessel function (Extended Data Fig. 3c). By taking into consideration spatial damping, Molácek and Bush [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF] demonstrated that in the deep-fluid regime, the wave field generated by a bouncing drop located at x j may be approximated by

h(|x -x j |) = 2A √ bJ 0 (k F |x -x j |)|x -x j |K 1 (2 √ b|x -x j |) ( 27
)
where K 1 is the modified Bessel function of the second kind. Our simulations indicate that A ≈ 15.4 µm and b ≈ 0.002 mm -2 , which is roughly consistent with the free-surface synthetic Schlieren experiments of Damiano et al. [START_REF] Damiano | Surface topography measurements of the bouncing droplet experiment[END_REF]. We note that the spatial damping included in [START_REF] Sáenz | Spin lattices of walking droplets[END_REF] 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 of our experimental study with the theoretical predictions of both the full wave coupling Bessel model [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF], and the generalized Kuramoto model [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF]. Specifically, we characterize the dependence 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 form of the coupling force.

We compare the predictions of our theoretical models to the experimental results through the average spincorrelation, χ . In the simulations, for every data point, we perform N t 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 time series χ n (t), for n = 1, ..., N t , as defined in equation ( 4). We then compute the average of the instantaneous spin correlations across all trials as

χ(t) = 1 N t Nt n=1 χ n (t). (28)
We compute the final average spin correlation, χ , by time-averaging over the final quarter of this time series, so as to eliminate any influence of the initial transient regime. In our simulations of ( 16) and [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF], this transient period is generally < 100 s. To guarantee robust statistics, we perform all simulations up to a time t max ∼ 300 -600 s.

We note that, owing to the approximations made in our modeling of the wave form, the preferred collective order is attained more rapidly in simulations than experiments. We perform numerical tests to ensure that our results are not affected by either an insufficient number of trials or finite-size effects. Our simulations indicate that the collective behavior becomes largely independent of the number of trials for N t ≈ 200 (Supplementary Information Fig. 4a). Similarly, our exploration of the effect of lattice size on the spin correlation (Supplementary Information Fig. 4b) reveals that the statistics become independent of the number of spins for N ≈ 150. We thus select N t = 200 and N = 150 for all subsequent simulations.

Varying lattice spacing

As demonstrated experimentally (Fig. 2), there exists a transition from anti-ferromagnetic to ferro-magnetic order when the well separation, L, decreases from 17.7 to 13.2 mm. Continuously varying L is extremely prohibitive experimentally, we thus proceed to do so using simulations.

We computed the average spin correlation, χ , in terms of the lattice spacing, L, via the Bessel model [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF]. 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 F = 70 s -2 , values consistent with the experimental range γ ≤ γ c (Supplementary Information Fig. 1). We prescribe the orbital 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π/k F = 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 show that the preferred collective behavior oscillates between regions of positive 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, ∆φ ± (Extended Data Fig. 4b), which reveals that the collective lattice behavior is dominated by not two but four synchronization modes (Extended Data Fig. 4c-d).

In the first ferromagnetic region, denoted FM + , the model predicts that ∆φ ± ∼ 0, indicating that neighboring pairs of spins are on average in phase in the azimuthal direction, regardless of whether they co-rotate or counterrotate. The bias towards ferromagnetic order (Extended Data Fig. 4a) is thus associated with the dominance of the co-rotating synchronization mode over the counter-rotating mode (Extended Data Fig. 4c-d). Similarly, in the anti-ferromagnetic 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. In contrast, the ferromagnetic and anti-ferromagnetic regions denoted by FM -and AFM -, respectively, result from two different synchronization modes that are characterized by out-of-phase orbital motion, ∆φ ± ∼ π (Extended Data Fig. 4c-d). The relative dominance of these modes thus determines the bias towards positive or negative spin correlation (Extended Data Fig. 4a).

The origins of these four forms of collective order are rationalized through the generalized Kuramoto model [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF]. By examining the dependence of the control parameters α and β, as defined in ( 25)- [START_REF] Souslov | Topological sound in active-liquid metamaterials[END_REF], on lattice spacing L, we note that the four dynamical modes are those that precisely minimize the reduced coupling potential U (24). Specifically, when α > 0, which happens when J 0 (k F L) > 0, corresponding to FM + (Extended Data Fig. 4a), the first term of the coupling potential U strongly favors motion for which |x i -x j | is minimized. In-phase rotation is thus preferred, but this condition is equally satisfied by both co-and counter-rotating spins (Extended Data Fig. 4c-d,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 (k F L). In particular, β > 0 in the FM + region (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 |y i -y j | (Extended Data Fig. 4c-d). The emergence of ferromagnetic order FM + is thus rationalized. Similar arguments may provide rationale for the observed 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 |y i -y j |. We note the regions where no preferred collective order emerges (Extended Data Fig. 4a) precisely correspond to the extrema and inflection points of J 0 (k F L), at which either α = 0 or β = 0. Finally, we note that simulations of the Bessel model ( 16) with 50 × 50 lattices confirmed that a similar dependence of collective order on lattice spacing emerges for substantially larger 2D square lattices (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 based on both the Bessel model [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF] and the generalized Kuramoto model [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF]. 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 is increased progressively, both models predict the emergence of a more pronounced anti-ferromagnetic 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 anti-ferromagnetic 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, F and α, β, are increased to values corresponding to γ > γ c . Similar simulations with both the Bessel ( 16) and generalized Kuramoto (23) models were conducted for a 50 × 50 lattice in order to confirm that the same collective order emerges in 2D square lattices for variable memory (Extended Data Fig. 6d).

Bath rotation

We conclude by generalizing our theoretical framework to rationalize the transition from anti-ferromagnetic 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 precisely balanced by the inward interfacial force resulting from the parabolic bath surface and so need not be considered henceforth [START_REF] Fort | Path-memory induced quantization of classical orbits[END_REF]. The effect of bath rotation may thus be incorporated into our modeling through the addition of the Coriolis force, F c = -2Ω × ẋ [START_REF] Oza | Pilot-wave dynamics in a rotating frame: on the emergence of orbital quantization[END_REF], which acts in the inward radial direction for circular 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 a confining radial potential of the form V (x) = k|x| 4 , so that F w = -∂V ∂x [START_REF] Cristea-Platon | Walking Droplets in a Circular Corral: Quantisation and Chaos[END_REF]. We find that the fourth order degree of the potential captures adequately the abrupt influence of topography at the well boundary. We select the strength of the potential, k, so that the predictions for the free walking speed u 0 and experimentally observed orbital radius r are consistent with our experimental observations in the absence of rotation (Fig. 1c).

Focusing on the wave-mediated droplet coupling F i,j as approximated by the generalized Kuramoto force ( 23), we may thus write the 2D version of the Rayleigh oscillator-equation ( 9), with an applied force

F = F c + F w + F i,j as ẍi = 1 τ 1 - | ẋi | 2 u 2 0 ẋi -4k|x i | 2 x i -2Ω × ẋi + j F i,j (x i , x j , L). (29)
Noting that the orbital radius is now a variable r i (t), we may express the droplet position around the center of the i -th well in polar coordinates as

x i (t) = r i (t) (cos φ i (t), sin φ i (t)) . (30)
Since full resolution of the droplet's dynamics in the presence of rotation is beyond the scope of the theory presented here, we simplify the equations to rationalize the observed behavior. We thus neglect the contribution of ṙi to the Rayleigh drag term, which is valid provided that ṙ2 i r 2 i ω 2 i , and the contribution to the Rayleigh drag term in the radial equation, which is valid provided that | ṙi | τ r i ω 2 i . 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 significantly from circles. With these assumptions, we may then substitute [START_REF] Lindner | Effects of noise in excitable systems[END_REF] into [START_REF] Bush | Pilot-wave hydrodynamics[END_REF]

to yield ri = r i ω 2 i + 2Ωω i r i -4kr 3 i , (31) ωi
= 1 τ 1 - r 2 i ω 2 i u 2 0 ω i - 2 ṙi r i (Ω + ω i ) - j ∂U i,j ∂φ i . ( 32
)
The numerical prediction of the model formed by [START_REF] Cumin | Generalising the Kuramoto model for the study of neuronal synchronisation in the brain[END_REF] 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 N t = 200 independent trials. We fix τ = 0.1 s, u 0 = 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 that observed in the laboratory.

Tunability

HSLs also offer the opportunity to investigate more exotic forms of collective behavior and frustration effects induced by wave-mediated interactions. Here, we provide several proof-of-concept experiments and simulations in order to illustrate their potential for future research.

As discussed in 2.1, it is possible to deduce the vertical phase of the droplets by observing the light reflected 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 (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 (Extended Data Fig. 4). Consequently, for example, the emergent order of a HSL whose lattice spacing promotes orbital in-phase anti-ferromagnetic order AFM + for vertically in-phase walkers will shift to out-of-phase anti-ferromagnetic order AFM -for vertically out-ofphase walkers (Extended Data Fig. 7a). A HSL with walkers whose vertical phases are random will thus produce an anti-ferromagnetic state with competing orbital synchronization that effects the emergent spin correlation (Extended Data Fig. 7b). Other forms of frustration can be investigated by changing the lattice geometry in 2D. For example, a triangular lattice (Extended Data Fig. 7c) whose spacing is tuned to promote anti-ferromagnetic order leads to degenerate ground states, the investigation of which is expected to yield insights into how wave-coupled spin lattices respond to geometric frustration [START_REF] Diep | Frustrated Spin Systems[END_REF].

We also perform simulations with the model of Faria [START_REF] Faria | A model for Faraday pilot waves over variable topography[END_REF], which 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 anti-ferromagnetic AFM + orders are reproduced with Faria's model (Extended Data Fig. 8). We proceed by simulating a lattice where the horizontal spacing promotes FM + , while the vertical spacing promotes AFM + (Extended Data Fig. 9a, left). Drops along the horizontal lines are thus rotating in phase, but their direction alternates vertically in an anti-ferromagnetic fashion (Extended Data Fig. 9a, center & right).

Next, we take as reference the ferromagnetic HSL (Extended Data Fig. 8a) and repeat the simulation with a lattice where the center of each well has been shifted by a random distance ±ε in the horizontal and vertical directions (Extended Data Fig. 9b, left). As the degree of randomness increases, the walkers' trajectories become less regular; thus, it becomes more difficult for the walkers to maintain precise in-phase orbital motion (Extended Data Fig. 9b, center & right).

In our third demonstration, with the same FM + as the reference (Extended Data Fig. 8a), we alternate the sizes of the bouncing droplets, thereby mixing two species of 'spins' with different frequency and inertia (Extended Data Fig. 9c,left). While the drop size does not change the preferred mode of synchronization, which is dictated by the lattice spacing, the different rotation frequencies lead to frustration. Notably, the drops adjust their orbital radii (Extended Data Fig. 9c, center) in order to match their rotation frequencies, as can be deduced from the equal slope of the φ(t) signals (Extended Data Fig. 9c,right).

In our fourth demonstration, we reduce the local liquid depth along alternating columns of a FM + lattice geometry (Extended Data Fig. 9d,left). As a result, the waves above the new submerged barriers dissipate waves more effectively (Extended Data Fig. 9d, center). Walkers are thus subjected to lateral couplings of different strengths, which leads to pairs whose phase is either strongly or weakly correlated (Extended Data Fig. 9d, right), and so to possibility of exploring SSH-type couplings as studied in topological insulators [START_REF] Qi | Topological insulators and superconductors[END_REF]. We also note that 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 in order to study the effects of imposed white or colored noise.

HSLs offer the opportunity to explore other potentially fruitful directions. One of them is the analysis of the collective order emerging from higher excited spin states such as trefoils and lemniscates [START_REF] Perrard | Self-organization into quantized eigenstates of a classical wave-driven particle[END_REF], which emerge with larger well diameters [START_REF] Cristea-Platon | Walking Droplets in a Circular Corral: Quantisation and Chaos[END_REF], and at higher driving acceleration (Fig. 1c). Notably, the accompanying theoretical rationalization would require more involved models, such as a full implementation of the stroboscopic model [START_REF] Oza | A trajectory equation for walking droplets: hydrodynamic pilotwave theory[END_REF] or the variabletopography model of Faria [START_REF] Faria | A model for Faraday pilot waves over variable topography[END_REF]. 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 the opportunity to investigate macroscopic dynamics analogous to those involving the oscillatory Ruderman-Kittel-Kasuya-Yosida (RKKY) [START_REF] Ruderman | Indirect Exchange Coupling of Nuclear Magnetic Moments by Conduction Electrons[END_REF] electronic spin coupling. Finally, investigating the similarities and differences between the collective synchronization arising in HSLs and other active systems might also prove fruitful. For example, it is known that rotating vortices in reaction-diffusion systems may be coupled by nonlinear spiral waves [START_REF] Biktasheva | Wave-particle dualism of spiral waves dynamics[END_REF][START_REF] Biktasheva | Computation of the drift velocity of spiral waves using response functions[END_REF][START_REF] Azhand | Three-dimensional autonomous pacemaker in the photosensitive Belousov-Zhabotinsky medium[END_REF][START_REF] Totz | Spatial confinement causes lifetime enhancement and expansion of vortex rings with positive filament tension[END_REF]. Snapshots at different times illustrate a typical spin-flip event arising in a 1D spin lattice (Movie S6). At each time, the panels at left are color coded according to the instantaneous spin S i while those on the right depict the recent walker trajectories, color coded according to the local speed. Perturbed by the wave fields emitted by its neighbours, the middle walker initially follows an elliptical path. The length of the minor axis decreases until the walker trajectory essentially becomes a straight line across the well center. Subsequently, the process is reversed, resulting in the walker rotating in the opposite sense. The three walkers shown are part of the 1D anti-ferromagnetic lattice described in Fig. 1h with γ/γ F = 86.0%. Extended 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 regions with a different shade of gray. (b) Superposition of the wave field shown in a and the zeros of the drop-centered Bessel function J 0 (k F x). (c) Wave field of a bouncer computed with the theoretical model developed by Faria [START_REF] Faria | A model for Faraday pilot waves over variable topography[END_REF] for walkers 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 center-to-center separation L as in a and γ/γ F = 88.0%. Solid blue lines denote the submerged wells and dashed lines the zeros of a Bessel function J 0 centered at the drop position.

Motor

Extended Data Fig. 4: Emergent order for varying lattice spacing. The dependence of (a) the average spin correlation, χ , and (b) the mean phase difference, ∆φ , on the lattice spacing, L, as predicted by the Bessel model ( 16) with τ = 0.4 s, ω 0 = 3.3 s -1 , and F = 70 s -2 . The average droplet-induced wavefield, which is approximated by a Bessel function J 0 (k F L), centered in the neighboring well at L = 0, is provided for comparison in panel a. (c) Dominant and (d) subdominant synchronization modes for the corresponding sections in panels a-b.

Extended Data Fig. 5: Emergent order for varying bath acceleration. Comparison of the experimentally observed average spin correlation, χ , with the predictions of the Bessel model [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF] and the generalized Kuramoto model [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF]. Bessel model parameters: L = 17.7 mm, r = 1.8 mm, λ F = 2π/k F = 5.1 mm. The interaction parameter F is varied across the range 70 < F < 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 m w = 1.65m, in line with prior work [START_REF] Bush | The wave-induced added mass of walking droplets[END_REF]. GK model parameters: α is varied over the range 8.5 < α < 15 s -2 , while maintaining β < 0 and a constant ratio β/α = 0.3. By dividing the expressions for α and β in equation ( 25)-( 26) by a factor of two, in accordance with the mismatch between the GK and Bessel models discussed in Supplementary Information Fig. 3, the minimum χ predicted by the GK model emerged in the vicinity of γ c . Extended Data Fig. 6: Emergent order in large 2D square lattices. Simulations of the Bessel model ( 16) and the generalized Kuramoto model [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF] for a 50×50 square lattice demonstrate the emergence of anti-ferromagnetic and ferromagnetic order in 2D for various lattice spacings and bath accelerations. (a) The lattice spacing determines the emergent anti-ferromagnetic (AFM + ) or ferromagnetic (FM + ) order in a manner predicted by our reduced theory [START_REF] Nishiguchi | Engineering bacterial vortex lattice via direct laser lithography[END_REF]. Specifically, preferred in-phase rotation between neighboring pairs can be clearly observed in the (b) antiferromagnetic AFM + , and (c) ferromagnetic FM + regimes. We note that (b) and 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, in order to ensure statistical significance. Extended Data Fig. 7: Emergent order for different vertical bouncing synchronizations and lattice geometry. (a) Oblique view of a 2D spin lattices where the walker in the center is bouncing vertically in-phase and out-of-phase with its left and right neighbours, respectively. (b) Average spin correlation for lattices with the same geometry as those described in Fig. 1 in the main text when the walkers all bounce in phase (blue, result presented in main text), out of phase (green), or have randomly distributed bouncing phases (red). Vertically in-phase pairs promote in-phase orbital anti-ferromagnetic order (AFM + ) while vertically out-of-phase pairs promote out-of-phase orbital anti-ferromagnetic order (AFM -). A random distribution of vertical phases thus leads to competing orbital synchronization modes, which has an effect on the emergent spin correlation. (c) Triangular HSLs whose lattice spacing is tuned to promote anti-ferromagnetic order can be used to investigate frustration effects.

Extended Data Fig. 8: Emergent collective order in simulations. Simulations of square HSLs with the theoretical model developed by Faria [START_REF] Faria | A model for Faraday pilot waves over variable topography[END_REF] is used to explore collective order in 2D. For appropriate lattice spacings, (a) in-phase ferromagnetic order FM + and (b) in-phase anti-ferromagnetic order AFM + are simulated. (Left) Schematic, (middle) wave field and walker's trajectories, and (right) time evolution of the orbital phases. . Extended Data Fig. 9: Model tunability. Proof-of-concept simulations performed with the model developed by Faria [START_REF] Faria | A model for Faraday pilot waves over variable topography[END_REF] illustrate the tunability and potential for future research of HSLs. (Left) Schematic, (middle) wave field and drop trajectories, and (right) orbital phase evolution. (a) A HSL tuned to promote FM + along the horizontal direction, but AFM + across vertical pairs. (b) FM + lattice geometry with a random vertical and horizontal shift ±ε in the position of each well. (c) FM + lattice geometry with two drop sizes (and so, two walker speeds). (d) FM + lattice geometry with coupling strength controlled locally through the thickness of the liquid layer.

T D = (ν e k 2 F ) -1 , T M = T D (1-γ/γ F )
Viscous decay time, wave-decay memory-time [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF] C = 0.17, D = Cmg ρR d σ + 6πR d µ a Drag coefficient, time-averaged horizontal drag coefficient [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF]

F = νeT F 2π (mgk 2 F R d ) 2 3σk 2
F +ρg sin Φ, sin Φ Wave-force coefficient, phase parameter [START_REF] Moláček | Drops walking on a vibrating bath: towards a hydrodynamic pilot-wave theory[END_REF]

m w = m 1 + F k F T 2 M 2m
Effective walker mass

D R = D F k F T 2 M 2T F D -1 , τ = mw D R
Rayleigh-drag strength, Relaxation time

u 2 0 = 8 3k 2 F T 2 M F k F T 2 M 2T F D -1 , ω 0 = u0 r
Free walking speed, angular velocity

F = F T M T F mwr
Interaction force coefficient

Fig. 1 :Fig. 2 :

 12 Fig. 1: Spontaneous anti-ferromagnetic order in a 1D HSL. (a) Oblique view and (b) schematic cross section of a ring (of radius R) of submerged circular wells (of diameter D, depth h, center separation L), each containing a single drop executing clockwise (CW) or anti-clockwise (ACW) motion on the surface of a vibrating fluid bath. A thin fluid layer (of depth H) between wells enables wave-mediated interactions between neighboring droplets (Movie S1). (c) Dependence of the spin S and orbital radius r of an isolated walker (D = 14 mm, H = 1 mm, f = 80 Hz) on the dimensionless bath acceleration γ/γ F . Once γ exceeds a critical threshold γ c , the drop's circular trajectory destabilizes into a precessing trefoil, increasing spin variation and susceptibility to perturbation-induced spin flips. For a given driving acceleration, only one state is observed. Insets: Two typical drop trajectories, colored by instantaneous spin value, at low and high accelerations. Centerlines denote interpolated mean values and envelopes indicate standard deviations. (d-f) Temporal evolution of the (d) spin values S i , (e) spin pair correlation χ and magnetization M , and (f) phase difference ∆φ for a 1D periodic HSL with N = 20 wells driven at γ/γ F = 82.0 %. (g) Snapshot showing collective anti-ferromagnetic order at instant A in e (Movie S2). Dashed lines highlight domains of local anti-ferromagnetic order. White circles mark drop positions, background color shows the spin magnitude. (h) The average spin correlation χ indicates maximal anti-ferromagnetic ordering for γ ≈ γ c , highlighting that stable single-spin trajectories promote collective ordering. Error bars denote the standard deviations (Methods, Statistics). Experimental parameters in d-h: L = 17.7 mm, D = 14 mm, H = 1 mm, N = 20, f = 80 Hz, γ F = 4.780g, λ F = 4.75 mm.

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

  Fig.3: Inducing global polarization through applied rotation. Rotating the circular HSL from Fig.1d-hinduces a transition from effective anti-ferromagnetic 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 anti-ferromagnetic order with no net magnetization. At t = 0.5 h, the bath rotation was switched on (Ω = 0.7 rad/s, ACW) inducing a transition to ferromagnetic order and positive magnetization. At t = 1 h, the rotation direction was inverted (Ω = -0.7 rad/s, CW) leading to negative magnetization. Lattice parameters and forcing frequency are the same as in Fig.1and γ ≈ γ c . (d) Instantaneous spin configurations at time points A-C in b showing anti-ferromagnetic and ferromagnetic ordering (Movie S4). (e) Average spin correlation χ signals a transition to ferromagnetic order at Ω c ≈ 0.22 rad/s (Methods, Statistics).(f,g) Force balance and measured trajectory, color-coded by speed, for an isolated spin in a frame corrotating with the bath. The drop's trajectory is nearly circular when its angular velocity ω and the bath rotation Ω are aligned, but approximates a precessing trefoil when the signs are opposite. (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; r 0 is the trajectory radius for Ω = 0 shown in Fig.1c. The center line denotes the mean radius, and the envelope indicates the standard deviation.

Fig. 4 :

 4 Fig. 4: Rotation-induced polarization transition from anti-ferromagnetic 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 ferromagnetic order for a rotating lattice (Ω = 0.75 rad/s). Instantaneous spin arrangements at time points A and B in panel b showing the emergence of (c) anti-ferromagnetic, and (d) ferromagnetic domains, respectively (Movie S5). (e) A transition from anti-ferromagnetic to ferromagnetic order is prompted by increasing the bath's angular speed.

METHODS

  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.

. 1 : 2 :

 12 Schematic of the experimental set-up. Extended Data Fig. Experimental spin flips.

  (c) correspond to simulations of the Bessel model with the spacings indicated on (a). (d) The emergent in-phase anti-ferromagnetic order (AFM + ) as a function of bath acceleration. Bessel model parameters in (a): τ = 0.1 s, F = 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 ≤ F ≤ 85 s -2 , which is transformed back to γ/γ F using the relation from Extended Data Table

  Stationary bath . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 Rotating bath . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3 Lattice parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3 Statistics 3.1 Order parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2 Data acquisition protocol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3 Boundary conditions, ring curvature R, and lattice size N . . . . . . . . . . . . . . . . . . . . . . . . . 14 4 Theory 4.1 Rayleigh-Bessel oscillator model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.2 Generalized Kuramoto (GK) model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 4.3 Spin flips and wave coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Bath rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

	1 Introduction
	2 Experiments
	2.1 5 Tunability

4.3.1 Spin flips . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.3.2 Wave coupling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.4 Comparison to experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.4.1 Varying lattice spacing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.4.2 Varying interaction strength . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.5

Extended Data Table 1 :

 Data1 Model parameters and physical variables for HSLs. Description σ, µ, ρ, ν = µ/ρ, ν e , µ a Liquid surface tension, dynamic viscosity, density, kinematic viscosity, effective kinematic viscosity, air viscosity m, D d , R d Droplet mass, droplet diameter, droplet radius D, L, r Well diameter, center-center separation, orbital radius γ, f Bath forcing acceleration, frequency k F , T F = 2/f , γ F Faraday wavenumber, period and threshold

	Variable

The authors gratefully acknowledge financial support from the NSF through CMMI-1727565 (J.W.M.B. & P.J.S) and DMS-1719637 (R.R.R.), MIT Solomon Buchsbaum Research Fund (J.D.), and CNRS Momentum program (G.P.).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

SUPPLEMENTARY INFORMATION

Supplementary Information Fig. 1: Parameters of the Bessel model.

Dependence on the normalized bath acceleration γ/γ F , in the range of interest, of the parameters arising in the Bessel model ( 15): (a) amplitude of the coupling force, F, (b) stable frequency of circular orbits, ω 0 , and (c) relaxation timescale, τ . We recall that for γ > γ c ≈ 0.85γ F , droplets execute trefoil-like trajectories, demarcating the limits of our theoretical model. Supplementary Information Fig. 2: Parameters of the generalized Kuramoto model. Dependence of the ratio of β to α on the lattice separation, L, for the generalized Kuramoto model ( 23)- [START_REF] Nishiguchi | Engineering bacterial vortex lattice via direct laser lithography[END_REF]. A plot of J 0 (k F L) is included for the sake of comparison. Supplementary Information Fig. 3: Theoretical wave-coupling force. Computed wave-coupling force F acting on droplet i due to its two neighbors j = i ± 1, for center-center well separations of (a) L = 13.6 mm, (b) L = 14.8 mm, (c) L = 16 mm, (d) L = 17.2 mm. First column: the prediction of the Bessel coupling force in [START_REF] Schneider | Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice[END_REF]. Second column: the prediction of the generalized Kuramoto model in [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF]. Third column: the sum of the dominant and subdominant Fourier components of the first column. Fourth column: the normalized power spectrum of the first column. The orbital radius is fixed at r = 1.2 mm, and the Faraday wavelength is λ F = 2π/k F = 4.95 mm. Supplementary Information Fig. 4: Dependence of the emergent order in the simulations on the lattice size and number of trials. Dependence of spin correlation, χ , on (a) the number of independent trials, N t , and (b) the number of spins, N , using the generalized Kuramoto model [START_REF] Wioland | Ferromagnetic and antiferromagnetic order in bacterial vortex lattices[END_REF] with α = 10 s -2 , β = -3 s -2 , τ = 0.1 s, ω 0 = 3.3 s -1 . The number of spins is set to N = 150 in a, and the number of trials to N t = 200 in b.