Virtual magnetic hills to unlock the inner phases of hexagonal colloidal ice

Abstract

We study the low energy states in a hexagonal colloidal ice realized by using repulsive paramagnetic colloids confined by gravity within a honeycomb lattice of traps. In contrast to similar systems featuring optical or topographic double wells, here we introduce field tunable “virtual” magnetic hills. These hills are created by placing pairs of fixed paramagnetic particles close to the semi-cylindrical traps that contain the interacting, mobile colloids. With this strategy, a single magnetic field can be used to simultaneously tune the particle pair-interactions and the hill elevation, without losing the trap bistability at any field strength. We use numerical simulations to explore the rich low energy states of the system. By varying both the relative distance and the magnetic content of the fixed particles, not only the effects of the first but also of the second nearest neighbors can be accessed, allowing the inner charge-ordered ice-II phase to be reached. Our strategy of controlling the vertex energetics via fixed, field tunable interstitial units may be extended to other geometrically frustrated systems on different length scales, including nanoscale spin ice and macroscopic magnetic metamaterials.

---

1 Introduction

In geometrically frustrated systems, interacting units are arranged on a lattice with a geometry that makes it impossible to simultaneously satisfy all their interaction energies.^1–3 The classical example is the Ising net, where one has to locate 3 spins with antiferromagnetic interactions at a vertex of a triangle. No matter how these spins are arranged, two will always be parallel, generating a high energy, frustrated bond.⁴ The impossibility of reaching low energy states in geometrically frustrated systems has led to many intriguing phenomena including the presence of residual entropy,^5,6 ground state degeneracy⁷ and disordered spin liquids^8,9 or glassy^10,11 phases.

In recent years, modern technological advances have enabled the realization of pre-designed artificial lattices where strongly interacting units can be mapped to Ising-like spins with associated topological charges and binary degrees of freedom. These include artificial spin ice,^12–14 array of macroscopic magnets,^15–19 confined microgel particles,^20–22 mechanical metamaterials,^23–26 or colloidal particles dispersed in topographic double wells.^27–29 In contrast to bulk magnetic materials,^30–32 such systems enable the direct visualization of the interacting magnetic moments while systematically tuning the geometry, interactions, or boundary conditions by external means. They thus emerged as powerful platforms for investigating the physics of geometric frustration and supporting potential applications in data storage or logic circuits based on the manipulation and motion of topological charges.

In all such systems, perhaps the simplest two-dimensional geometry which can produce a highly degenerate ground state is the honeycomb one, where at each vertex, only three spins meet. Despite its low coordination number, z = 3, the honeycomb lattice has been the subject of intense research activities especially with artificial spin ice (ASI), e.g. lattices of interacting magnetic nanobars.^33–43 Indeed, nanomagnets arranged along the honeycomb geometry, known in the literature also as kagome ASI, have shown a rich spectrum of low temperature phases. These start from a high-temperature, disordered paramagnetic phase followed by a first Ice I state, where the spins at the vertices obey the “pseudo ice rule” of 2-in 1-out or 2-out 1-in, that represent q = ±1 topological charges. Further lowering the temperature reveals an Ice II phase where these topological charges become ordered, however, the associated magnetic moments are not, therefore ending in a long-range ordered phase (LRO) characterized by ordered magnetic moments and vertex charges. However, reaching such a state is rather challenging due to the presence of dynamical freezing,^44,45 where large energy barriers have to be overcome to reach further low energy configurations.

A similar sequence of phases was theoretically predicted to occur in a particle ice,^46,47 a soft matter analogue of an ASI.⁴⁸ In such a system, repulsive colloidal particles are confined within a lattice of double wells at a one-to-one filling ratio, such that each particle settles in one of the two wells separated by a central hill. The pair interactions are controlled by an external field and, at high interaction strengths, the particles can switch positions surmounting the central hill. Thus, the double well become bistable, and one can associate to it a binary variable. In this system, the double wells can be arranged in a lattice that impedes the minimization of all interaction energies at a vertex, and geometric frustration sets in. In previous experimental realizations with paramagnetic colloids^28,49 the pair interactions were tuned via a magnetic field. But, for large field amplitudes, the colloids were able to localize on top of the central hill, therefore losing the trap bistability.^50,51

Here we investigate the low energy states of a hexagonal colloidal ice made of interacting paramagnetic colloids confined within a lattice of semi-cylindrical traps. Instead of relying on double wells characterized by a central optical²⁷ or topographic²⁸ hill, we introduce “virtual” magnetic hills that are created by placing a pair of fixed particles with equal symmetric distances from the trap center. Using this strategy, we ensure trap bistability and tunable confinement for all field strengths. In addition, by changing the relative magnetic susceptibility or distance between the fixed and mobile colloids, we show that it is possible to reach an Ice II phase at high interaction strength, where ±1 topological charges arrange in an antiferromagnetic-like order alternating their sign periodically along the lattice. Finally, we show the feasibility of the proposed approach by providing an experimental realization of the magnetic hill, by combining soft lithography and magnetic field manipulation techniques.

2 The particle ice and the virtual magnetic hill

Fig. 1(a) illustrates a particle ice where repulsive colloids are placed within a lattice of double wells at a one-to-one filling ratio. The analogy with other frustrated systems results by assigning a binary spin variable to each double well. As shown in Fig. 1(b), this is achieved by considering a vector, similar to a spin, that points towards the well occupied by the particle. In this way, one can classify the different vertex types in a manner similar to other frustrated systems such as ASIs. Thus, to maintain bistability, it is important that the central hill confines the particle in one of the two wells, while still permitting the particle to switch wells due to the repulsive interaction from neighboring particles. As shown in Fig. 1(b), in previous realizations this bistability was ensured by the presence of a central hill of elevation h. In this case, a particle of diameter σ has to overcome a gravitational potential U_g ∼ σ³h ≫ k_BT to switch state. An external field perpendicular to the particle plane, B = Bẑ, induces a magnetic moment m = VχB/μ₀ in each particle pointing along its direction, being χ the particle magnetic volume susceptibility and μ₀ the permeability of the medium. Since the dipoles are parallel to each other, the colloids repel due to a tunable pair potential, U_dd = μ₀m²/(4πr³), where r = |r_i − r_j| is the relative distance between a pair of particles (i,j). Thus, a particle can switch location within the double well when the dipolar force from all neighboring particles is such that, U_dd > U_g.

Fig. 1 (a) Schematic showing a hexagonal colloidal ice made of paramagnetic colloidal particles confined by gravity within topographic double wells. A magnetic field B applied perpendicular to the particle plane (x̂, ŷ) induces an isotropic repulsion between the particles. (b) and (c) Two vertices with 5 double wells illustrate the gravitational hill (b) and 5 semi-cylindrical traps with 10 fixed colloids forming the “virtual” magnetic hill (c). Schematics at the top show the traps in the (x̂, ŷ) plane, and schematics at the bottom show the profile of a single trap in the (x̂, ẑ) plane. (d) Calculated total magnetostatic potential U_m felt by the central particle (within the red circle) versus distance from the center of the trap r/σ, where σ is the particle diameter. The top inset shows the log–log plot of the magnetic hill elevation U₀ versus field amplitude B. (e) Three-dimensional graph showing the magnetostatic potential U_m calculated for a field amplitude B = 0.5 mT.

In previous works the topographic hill h was fixed and could not be tuned.²⁸ Thus, at high field strengths, the particles were capable of localizing on top of the hills in order to maximize their distance from each other.^50,51 This effect impeded the use of large field amplitudes and limited the exploration of low energy states. We overcome this problem by introducing “virtual” magnetic hills, where the energy barrier is tuned by the applied field, preserving the particle bistability for all field strengths. As shown in Fig. 1(c), the magnetic hill is realized by placing a pair of fixed paramagnetic colloids at the side of a semi-cylindrical trap that contains the spin-associated mobile particle. Once the field is applied, the dipolar repulsion from the fixed colloids generates a magnetic barrier similar to a double well, as shown by the calculation reported in Fig. 1(d). In particular, the virtual magnetic hill that acts on a particle located near two vertices, can be captured by the total magnetostatic energy: . Here N = 10 since, as shown in Fig. 1(c), U_m contains the contribution from all the dipolar interactions due to the 4 nearest mobile colloids and 10 fixed ones. Moreover, the paramagnetic nature of the colloidal particles ensures that the magnetic hill is field tunable, and its elevation grows as U₀ ∼ B², reaching U₀ = 100 k_BT for B = 1.8 mT (inset in Fig. 1(d)). The corresponding three-dimensional representation of such potential is also shown in Fig. 1(e), calculated for a field amplitude of B = 0.5 mT.

3 The hexagonal colloidal ice

Using numerical simulations (see Appendix for details), we explore the low energy states of a honeycomb colloidal ice. The effect of geometric frustration in such a system can be investigated by measuring the fraction ϕ of topological charges q associated with each vertex in the lattice. As shown in Fig. 2(a), the hexagonal ice, with coordination number z = 3, is characterized by four vertex types depending on the particle location. The corresponding topological charge is given by q = 2n − z, where n is the number of spins (or particles) pointing towards the vertex center. In contrast to the square geometry (z = 4) where low energy states are characterized by ice rule vertices with q = 0, in the honeycomb one, there are two low charged states with q = ±1 (K_III, K_II) that survive at high interaction strengths. Theoretical arguments⁴⁸ supported by subsequent experiments⁵² demonstrated that for lattices characterized by a single coordination number z, the particle ice displays an effective vertex energetics similar to that of an ASI. Thus, we expect to observe a sequence of states similar to ASIs when increasing the interaction strength, as shown in Fig. 2(b). Starting from a disordered state characterized by low (q = ±1) and high (q = ±3) charged vertices, the particle may organize in a first Ice I state filled only of disordered q = ±1 topological charges. Further lowering the energy will reveal an Ice II phase where the q = ±1 charges become ordered, finally ending in a long-range spin ordered phase (LRO).

Fig. 2 (a) The four vertex types in the hexagonal colloidal ice with their topological charge q. High charges (q = ±3) display a green dot at the vertex. (b) Schematics showing the particle arrangement in 7 plaquettes of the hexagonal colloidal ice respectively in the disordered phase, Ice-I, Ice-II and long range ordered or spin solid phase. (c)–(e) Simulation results of hexagonal colloidal ice made of topographic gravitational hills (left) and virtual magnetic hills (right). (c) Partial view of simulation snapshots showing the particle positions and the topological charges for both systems after annealing. Light brown disks on the right represent fixed interstitial particles. In both cases, we simulate N = 339 mobile colloids. (d) Spatial probability n of finding a particle (diameter σ) at distance Δr from the trap center, normalized by σ/2. Here n_c shows the central occupancy, while insets display representative examples of particle locations. (e) Vertex fraction ϕ including the fraction ϕ_h of q = ±3 topological defects in green (top), central occupancy n_c (middle) and the charge order ψ_CO (bottom), all versus field amplitude B.

We start by comparing the properties of colloidal ices made of gravitational topographic hills (h = σ/2) and virtual magnetic hills on the left and right sides of Fig. 2(c)–(e), respectively. Illustrative simulation snapshots of the two systems after annealing are shown in Fig. 2(c). Here, the system is annealed by increasing linearly the magnetic field from B = 0 to 15 mT in a time interval Δt = 1500 s. The most striking difference relies in the particle location within the traps. At large field amplitudes, most of the particles localize in the center of the double wells for the gravitational case, while they remain on one side of them for the magnetic case. This is quantified in Fig. 2(d) by measuring the distribution of particle location n within the double wells for both cases. Importantly, for the gravitational case, this loss of bistability prevents the system from reaching a low energy state, or “inner” phase, characterized by high interaction strengths. Thus, as an observable for this loss of bistability, we define the “central occupancy” n_c as the fraction of particles located within the ±0.1 relative distance from the trap center.

Other observables to characterize the low energy states are plotted in Fig. 2(e) as a function of the applied field. The vertex fraction ϕ indicates that both systems tend to eliminate the high charges K_I and K_IV, as shown by their fraction ϕ_h which decreases to 0. Moreover, we find that for B ∼ 2 mT, ϕ ∼ 0.5 for the K_II and K_III characterized by q = −1 and q = +1, respectively. However, the fraction of K_II is slightly larger than that of K_III, with ϕ_KII = 0.51 and ϕ_KIII = 0.48. In contrast, the system with topographic hills does not present the same degeneracy, showing an unequal fraction of such defects with even a small amount of q = ±3 at large field amplitudes. This fraction of high charges results from the loss of bistability in the gravitational case, becoming pronounced when the central occupancy n_c becomes high. Note that, in both cases, to avoid considering the accumulation of charges at the boundaries, we excluded from the statistics all vertices located at the edge of the lattice. Considering these vertices will mainly affect the fraction of K_II and K_III for the magnetic hill. In particular, these fractions will separate further, becoming ϕ_KII = 0.62 and ϕ_KIII = 0.37. The reason is the radial polarization effect⁵³ since vertices filled with repulsive colloids will tend to expel the particles towards the edge of the lattice favoring the formation of ϕ_KII over ϕ_KIII.

A further observable is the charge ordering, defined as:

where N_v is the number of vertices.⁴⁶ This parameter measures the charge–charge correlation between nearest neighbor vertices. It is ψ_CO = 1/3 for random spins, ψ_CO = 0 for an ice state with fully disordered charges, while it becomes ψ_CO = 1 when the q = ±1 topological charges perfectly crystallize in an antiferromagnetic way. This value thus describes well the charge crystallization in the low energy states. For the gravitational case, when ϕ_h vanishes (ice state) ψ_CO ∼ 0.1 and then saturates to ψ_CO ∼ 0.25. The saturation is likely due to the loss of bistability. For the magnetic case, when ϕ_h becomes zero, ψ_CO ∼ 0.1 and does not increase with the magnetic field since U₀ in the magnetic hills grows fast with the interactions between the mobile particles and this effect prevents the particles from rearranging, as shown by the constant values of all observables starting from B = 2 mT.

4 Tunable magnetic hills to reach the Ice II phase

In contrast to the double wells with fixed hills, the magnetic ones enable to access lower energy states by maintaining the trap bistability for all applied fields. However, as previously shown in Fig. 2(e), using magnetic hills containing the same type of fixed and mobile particles, did not allow for the independent tuning of the hill elevation U₀ from the pair interactions. This effect precludes to further induce particle rearrangement and to explore low energy states. One way to overcome this limitation is by changing the relative distance between the fixed and mobile particles, or the corresponding magnetic content. We introduce this concept in Fig. 3, where we explore the low energy states obtained by varying the ratios d₂/d₁ and χ₂/χ₁. As shown in the top scheme in Fig. 3(a), d₂ (χ₂) is the distance of the fixed particle (its magnetic volume susceptibility) from the center of its enclosing plaquette, and d₁ (χ₁) is the corresponding parameter for the mobile colloid. Note that in Fig. 2(c)–(e), we fixed d₂/d₁ = 0.5 and χ₂ = χ₁.

Fig. 3 (a) Heatmap showing the different stationary regimes observed for a hexagonal colloidal ice with magnetic hills, obtained by varying the ratio between the two distances d₂/d₁ and the particle magnetic susceptibilities χ₂/χ₁. As shown in the schematic at the top, d₂ (d₁) represents the radius of the circle where the fixed particles (center of the cylindrical trap) are located and χ₂ (χ₁) denotes the magnetic volume susceptibility of the fixed (mobile) particles. Charge order ψ_CO is color-coded where the system is found in an ice state (n_c < 1% and ϕ_h < 1%). (b) Sequence of simulation snapshots taken along the line of constant χ₂/χ₁ = 0.5 and increasing d₂/d₁, showing successive states characterized by (from left to right): central occupancy, charge-ordered ice-II, disordered ice-I, residual high charges. (c) The three observable ψ_CO, ϕ_h and n_c versus field amplitude B for (left panel) the ice-II charge-ordered state and (right panel) the ice-I charge-disordered state. In these simulations, we use N = 627 mobile colloids.

Fig. 3(a) displays the heat map of the different stationary regimes observed by varying the ratios d₂/d₁ and χ₂/χ₁ using a larger lattice than before (each side is composed of 20 hexagons rather than 10). We distinguish three different areas in the parameter space, with representative snapshots shown in Fig. 3(b). When both d₂/d₁ and χ₂/χ₁ are low (bottom-left corner), corresponding to high pair interactions compared to U₀, some particles localize in the trap center, similar to the gravitational hills. In contrast, when both parameters are high (top-right corners), corresponding to the opposite situation, residual high topological charges are found. In between lies the range of parameters that lead to the observation of an Ice state, where q = ±1 defects organize in an ordered way. In this case, the charge order ψ_CO varies continuously from ∼0 (Ice-I state), reaching more than 0.6 (Ice-II state) close to the boundary of the diagram. Fig. 3(c) shows the behaviour of corresponding observables for both ranges of parameters. In particular, we see that ψ_CO of Ice-I saturates around 3 mT, whereas, in the other case, it continues to grow for large fields, saturating around B = 10 mT at ψ_CO ∼ 0.7. Given the presence of boundary effects which prevent reaching a perfect charge crystal and the finite time of our annealing, we can actually reach the Ice-II state, in which the ±1 charges alternate in a checkerboard pattern, while the spins do not display a long-range order. Indeed our simulations were performed under open boundary conditions, a physical situation similar to previous experimental realizations.

5 Experimental realization

After exploring numerically the effect of the magnetic hills on the low energy states, we demonstrate experimentally their viability for the singular case. Specifically, we fabricate a honeycomb lattice of topographic, semi-cylindrical wells, and fill 5 of the wells with mobile particles surrounded by 10 fixed colloids to generate the magnetic hill. We use a water dispersion of paramagnetic colloidal particles with diameter σ = 2.8 µm and place them with optical tweezers within semi-cylindrical wells. The wells are imprinted on a thin layer of polydimethylsiloxane (PDMS) using soft lithography techniques (see Appendix for technical details).

These semi-cylindrical wells are arranged along a honeycomb lattice with a lattice constant a = 24 µm and each has a pair of circular wells on either side with 3 µm diameters in which fits exactly one particle. The experimental images at the right side of Fig. 4(a) and (b) show the two scenarios considered. In the “no hill” configuration, as shown in Fig. 4(a), only 5 particles were confined within the cylindrical wells, while in the “magnetic hill” case, shown in Fig. 4(b), we place additionally 10 paramagnetic colloids pinned within the circular wells located around them. We then apply an external magnetic field B ∈ [0,10] mT perpendicular to the sample plane to induce repulsive interactions between the particles and monitor the location Δr of the central particle.

Fig. 4 Experimental realization of the magnetic hill: (a) and (b) sequence of histograms showing the equilibrium position of the central particle for different field amplitudes and for 5 (a) and 15 (b) paramagnetic colloids located within a honeycomb structure. In (b), 5 particles are within semi-cylindrical traps, while 10 are held fixed within circular traps with similar size. The two snapshots at the right sides in (a) and (b) show optical microscope images for the two cases at their respective highest fields.

The sequences of histograms in Fig. 4(a) and (b) show how the particle positions change with the applied field in the two cases. For B = 0 mT, Δr displays a flat histogram since the particle can easily diffuse along the semi-cylindrical well independently from its neighbors. In the absence of the pinned particles (top row), the colloids tend to localize along the center of the trap, and the corresponding histogram becomes narrow at the maximum field strength of B = 12 mT. This central localization arises due to the repulsive interactions induced by the surrounding particles driving them to move apart and reach their maximum separation distance. In this situation, repulsive colloids replicate the honeycomb lattice, and the binary character of the colloidal ice is lost. In contrast, in the presence of the fixed particles, a tunable magnetic hill is created and the corresponding histogram of the particle position shows not a central peak but rather two well-defined peaks separated by a distance ∼σ. Now the mobile particles localize on one of the two sides of the semi-cylindrical trap, as shown in Fig. 4(b), which demonstrates the efficiency of the magnetic hill in confining the trapped particle. It should be noted that the data in the third panel of Fig. 4(b) are an average of two measurements where the mobile particle has randomly chosen distinct sides of the double well in each experiment.

Note that in contrast to the topographic double wells, the magnetic ones do not require a physical central hill, making the lithographic fabrication process relatively easier. Indeed in previous realizations the hill was obtained from the lower illumination of the central part of the trap during UV light exposure, which easily induces artifacts related to blurring and multiple interference that makes it difficult to control with precision its elevation. In contrast, an external field can now be used to tune the elevation of the magnetic hill at will, or to redistribute the particles homogeneously within the traps due to thermal fluctuations by switching off the magnetic field.

6 Conclusions

In this work, we have investigated the low energy states of a honeycomb colloidal ice where repulsive paramagnetic colloidal particles are trapped within lithographic semi-cylindrical wells. Geometric frustration arises due to the competition between the symmetry of the pair interactions, which favors the arrangement of the particles along a triangular lattice, and the imposed geometrical constraint. To reach high field values and to maintain the binary nature of the topographic hills, we have introduced virtual magnetic hills which are made of a fixed pair of paramagnetic colloids trapped on the two sides of the particle well. With this strategy, we generated field tunable magnetic hills that allow exploring the low energy state of the system at high field amplitudes. In addition, by tuning the relative distances and magnetic susceptibility of the fixed particles, we were able to reach an Ice II phase where low charged q = ±1 topological defects crystallize in an alternating way and the charge order reaches ψ_CO ∼ 0.7.

Our experimental results demonstrate the principle of the magnetic hill, however for now are limited to only two lattice vertices. In principle, they can be extended also to large systems as those explored via simulations. There are two main limitations of our experimental implementation in large lattice structures. First, the large number of particles needed to fill the fixed holes within the lithographic structure will require a relatively long time to manually locate these particles via optical tweezers. To overcome this problem, one could implement an automated program that could recognize and track a particle in real time, trap it via optical tweezers and deposit it in a given semi-cylindrical well in order to speed up the process. Another problem is related to the ratio of the magnetic volume susceptibility required to reach the Ice-II phase in numerical simulations. Our experiments make use of commercially available paramagnetic colloids with χ = 0.8 and the Ice-II state was reached for a narrow range of values of the ratio χ₂/χ₁. Finding paramagnetic colloids with a size similar to our particles but with smaller or larger magnetic susceptibility to match this ratio is rather difficult. One could avoid this problem by directly synthesizing these particles. However, this may require further time to achieve a monodisperse suspension with the desired magnetic properties. Solving both issues would allow us to explore experimentally the effect of magnetic hills on large lattice systems and to confirm the simulation results through the observed low energy states.

Our magnetic hills could be equally implemented in other geometrically frustrated settings including artificial spin ice or macroscopic interacting magnets. In particular, in relation to ASIs, a previous work demonstrated the possibility of introducing interstitial microdisks with free magnetization to modify the energy landscape in artificial spin ice states.⁵⁴ These interstitials were introduced at the center of the vertices, but they could equally be placed parallel to the magnetic nanobars to reduce or increase the corresponding magnetic interactions. In addition to frustrated magnetic systems, the magnetic hills can be used to investigate other interesting features, including the thermally induced transition of magnetic nanoparticles when tuning the hill elevation or introducing a bias. This could be done either via an external force, or making the magnetic double well asymmetric like using two near, fixed particles characterized by different magnetic contents.

Conflicts of interest

There are no conflicts to declare.

Data availability

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

Appendix

A.1 Numerical simulations

We perform Brownian dynamic simulations by considering a set of I = 1…N paramagnetic colloidal particles of diameter σ = 2.8 µm and magnetic volume susceptibility χ = 0.8, confined to move within a 2D array of semi-cylindrical double wells. Each particle at position r_i obeys the overdamped equation of motion:Here is the magnetic dipolar force acting on particle i due to the neighbouring ones and it is expressed as: , where m = χVB/(μ₀) is the dipole moment induced by the external field B, μ₀ = 4π × 10⁻⁷ H m⁻¹, and . Long-range effects due to magnetic dipolar interactions are considered using a large cutoff distance in the dipolar coupling equal to twice the lattice constant. The second term, F^T_i, represents the gravitational force of the lithographic double well exerted on the bead. This term is considered in the simulations in the left of Fig. 2(c)–(e) and it is given by:
where k_trap = 4 × 10⁻² pN nm⁻¹ is the spring constant modeling the wall of the trap. It brings back the particle into the trap both along the longitudinal and the transversal direction. The other spring constant k_hill models the gravitational hill in the middle of the double well whose value is given by mgh = k_hill σ²/8. It is k_hill = 9.85 × 10⁻⁵ pN nm⁻¹ for the example of topographic gravitational colloidal ice. Initially, the particles are located on one side of the double well and it is stable there in the absence of a magnetic field.

In contrast, when introducing the “virtual” magnetic hill, we set k_hill = 0 and include within F^dd_i the contribution of 2N additional paramagnetic colloids with susceptibility χ₂ placed in pairs and at fixed positions close to the main traps, as shown in Fig. 1(c). In this situation, the mobile colloids are confined within elongated semi-cylindrical wells (length 2σ, limiting the position of the bead between +σ/2 and −σ/2). In the absence of an external field, the particles are free to move along the trap by diffusion.

Finally, η in eqn (2) represents the random force due to thermal fluctuations, with zero mean, 〈η〉 = 0 and delta correlated, 〈η(t)η(t′)〉 = 2k_BTγδ(t – t′), with a temperature T = 300 K and the damping coefficient γ = 6πμσ/2 with μ the water viscosity leading to a diffusion coefficient D = k_BT/γ = 0.156 µm² s⁻¹.

The equation of motion for the N mobile particles and 2N fixed ones is integrated numerically using a small time step δt = 1 ms and by increasing linearly, i.e. at constant rate, the magnetic field from B = 0 mT up to B = 15 mT during the simulation duration of Δt = 1500 s.

A.2 Fabrication of the colloidal ice

The protocol used to fabricate our microstructure follows the method developed in a previous work.⁵¹ The honeycomb lattice was designed using AutoCAD for a commercial 5″ chrome photomask with class 4500k dpi resolution (JD Photodata). The following protocol was carried out in a 10 000 class clean room to reduce surface contamination.

First, a 10 cm silica wafer was dehydrated for 15 min at 200 °C followed by a plasma cleaning treatment (Harrick, PCD-002-CE). Then, a thin layer of SU-8 3005 epoxy-based photoresist (∼5 µm in height) was spin coated (Laurell Tech, WS-650) onto the Si wafer. A pre-exposure bake (JP Selecta Plactronic hotplate) was performed on the wafer for 2 min at 95° and then it was left to cool to room temperature. The mask aligner (SÜSS Microtec, Model MJB4) set to i-line configuration with wavelength λ = 365 nm, was fitted with both the Cr mask and Si wafer. The Si wafer was exposed to a UV light dose of 100 mJ cm⁻² during 5.4 s. A post-exposure bake was performed on the Si wafer for 1 min at 65 °C followed by 1 min at 95 °C. The Si wafer was then submerged and developed in propylene-glycol monomethyl ether acetate for 30 s and rinsed with isopropanol and dried with N₂ gas. Finally, the wafer was hard baked for 20 min at 95 °C followed by 10 min at 65 °C.

A.3 Realization of the topographic structure in PDMS

The wafer was primed for silanization treatment by first being plasma treated and then exposed to a drop of trichloro(1H,1H,2H,2H-perfluorooctyl)silane under vacuum for an hour to make the surface hydrophobic. A small quantity of polydimethylsiloxane (PDMS, SYLGARD 184 Silicone Elastomer Kit, Sigma-Aldrich) was prepared in a 10 : 1 ratio of the base and curing agent, mixed rigorously, and then degassed. A thin layer of the PDMS (∼25 µm in height) is then spin coated onto the Si wafer and left to settle on a levelled surface overnight. The PDMS was then cured the following day in a vacuum drying oven (Memmert, model V029) at 65 °C for 4 h. Glass coverslips (Menzel-Gläser) were washed with soap, water, ethanol and isopropanol and dried with N₂ gas and then plasma-treated along with the wafer. The plasma-treated surfaces of the glass coverslips were placed face down on top of the microstructures of the wafer, bonding the coverslip to the thin PDMS membrane. The coverslip, with the inverted colloidal ice microstructure embedded in the thin membrane of PDMS, was then very carefully peeled off from the Si wafer. The particles and the structure were enclosed within a 15 mm × 15 mm × 260 µm plastic square spacer (Gene Frame).

A.4 Experimental setup

Individual particles were then allocated to each cylindrical trap, or circular well, using holographic optical tweezers integrated with a Nikon Eclipse Ti inverted microscope with a 100× oil objective. The stage of the microscope was equipped with a magnetic coil situated beneath the sample that generates a field oriented in the ẑ-axis inducing a repulsive interaction between particles. The coil is powered by a bipolar operational power amplifier (KEPCO) and controlled by a LabVIEW program. Image capturing is done using a metal oxide semiconductor camera (MQ013CG-E2, Ximea).

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 811234). P. T. acknowledges support from the Generalitat de Catalunya (ICREA Académia).

References

1. J. Vannimenus and G. Toulouse, J. Phys. C: Solid State Phys., 1977, 10, 2680–2684.
2. J. Sadoc and R. Mosseri, Geometrical Frustration, Cambridge University Press, Cambridge, UK, 1999.
3. R. Moessner and A. P. Ramirez, Phys. Today, 2006, 59, 24–29.
4. G. H. Wannier, Phys. Rev., 1950, 79, 357–364.
5. L. Pauling, J. Am. Chem. Soc., 1935, 57, 2680–2684.
6. W. Giauque and J. Stout, J. Am. Chem. Soc., 1936, 58, 1144–1150.
7. Y. Perrin, B. Canals and N. Rougemaille, Nature, 2016, 540, 410–413.
8. T.-H. Han, J. S. Helton, S. Chu, D. G. Nocera, J. A. Rodriguez-Rivera, C. Broholm and Y. S. Lee, Nature, 2012, 492, 406–410.
9. D. Zhou, F. Wang, B. Li, X. Lou and Y. Han, Phys. Rev. X, 2017, 7, 021030.
10. S. Kirkpatrick, Phys. Rev. B, 1977, 16, 4630–4641.
11. D. Chowdhury, Spin Glasses and Other Frustrated Systems, Cambridge University Press, Princeton University Press, Princeton, USA, 2014.
12. R. F. Wang, C. Nisoli, R. S. Freitas, J. Li, W. McConville, B. J. Cooley, M. S. Lund, N. Samarth, C. Leighton, V. H. Crespi and P. Schiffer, Nature, 2006, 439, 303.
13. C. Nisoli, R. Moessner and P. Schiffer, Rev. Mod. Phys., 2013, 85, 1473–1490.
14. S. H. Skjærvø, C. H. Marrows, R. L. Stamps and L. J. Heyderman, Nat. Rev. Phys., 2020, 2, 13–28.
15. P. Mellado, A. Concha and L. Mahadevan, Phys. Rev. Lett., 2012, 109, 257203.
16. E. Y. Vedmedenko, Phys. Rev. Lett., 2016, 116, 077202.
17. K. Wang, X.-J. Liu, L.-M. Tu, J.-J. Zhang, V. N. Gladilin and J.-Y. Ge, Phys. Rev. B, 2025, 111, 224418.
18. J.-Y. Ge, J.-J. Zhang, L.-M. Tu and V. N. Gladilin, Newton, 2025, 1, 100143.
19. C. Reichhardt, C. Nisoli and C. J. O. Reichhardt, Newton, 2025, 1, 100174.
20. Y. Han, Y. Shokef, A. M. Alsayed, P. Yunker, T. C. Lubensky and A. G. Yodh, Nature, 2008, 456, 898.
21. Y. Shokef, A. Souslov and T. C. Lubensky, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 11804–11809.
22. Y. Zhou, K. Kanoda and T.-K. Ng, Rev. Mod. Phys., 2017, 89, 025003.
23. S. H. Kang, S. Shan, A. Košmrlj, W. L. Noorduin, S. Shian, J. C. Weaver, D. R. Clarke and K. Bertoldi, Phys. Rev. Lett., 2014, 112, 098701.
24. C. Coulais, E. Teomy, K. de Reus, Y. Shokef and M. van Hecke, Nature, 2017, 535, 529–532.
25. A. S. Meeussen, E. C. Oǧuz, Y. Shokef and M. van Hecke, Nat. Phys., 2020, 16, 307–311.
26. C. Sirote-Katz, D. Shohat, C. Merrigan, Y. Lahini, C. Nisoli and Y. Shokef, Nat. Commun., 2024, 15, 4008.
27. A. Libál, C. Reichhardt and C. J. O. Reichhardt, Phys. Rev. Lett., 2006, 97, 228302.
28. A. Ortiz-Ambriz and P. Tierno, Nat. Commun., 2016, 7, 10575.
29. A. Ortiz-Ambriz, C. Nisoli, C. Reichhardt, C. J. O. Reichhardt and P. Tierno, Rev. Mod. Phys., 2019, 91, 041003.
30. S. T. Bramwell and M. J. Gingras, Science, 2001, 294, 1495–1501.
31. T. Fennell, P. P. Deen, A. R. Wildes, K. Schmalzl, D. Prabhakaran, A. T. Boothroyd, R. J. Aldus, D. F. McMorrow and S. T. Bramwell, Science, 2009, 326, 415–417.
32. J. S. Gardner, M. J. P. Gingras and J. E. Greedan, Rev. Mod. Phys., 2010, 82, 53–107.
33. M. Tanaka, E. Saitoh, H. Miyajima, T. Yamaoka and Y. Iye, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 052411.
34. Y. Qi, T. Brintlinger and J. Cumings, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 094418.
35. G. Möller and R. Moessner, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 140409.
36. G.-W. Chern, P. Mellado and O. Tchernyshyov, Phys. Rev. Lett., 2011, 106, 207202.
37. N. Rougemaille, F. Montaigne, B. Canals, A. Duluard, D. Lacour, M. Hehn, R. Belkhou, O. Fruchart, S. El Moussaoui, A. Bendounan and F. Maccherozzi, Phys. Rev. Lett., 2011, 106, 057209.
38. W. R. Branford, S. Ladak, D. E. Read, K. Zeissler and L. F. Cohen, Science, 2012, 335, 1597.
39. I. A. Chioar, B. Canals, D. Lacour, M. Hehn, B. Santos Burgos, T. O. Menteş, A. Locatelli, F. Montaigne and N. Rougemaille, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 220407.
40. B. Canals, I.-A. Chioar, V.-D. Nguyen, M. Hehn, D. Lacour, F. Montaigne, A. Locatelli, T. O. Mentes, B. S. Burgos and N. Rougemaille, Nat. Commun., 2016, 7, 11446.
41. A. Farhan, P. M. Derlet, L. Anghinolfi, A. Kleibert and L. J. Heyderman, Phys. Rev. B, 2017, 96, 064409.
42. T. Dion, D. M. Arroo, K. Yamanoi, T. Kimura, J. C. Gartside, L. F. Cohen, H. Kurebayashi and W. R. Branford, Phys. Rev. B, 2019, 100, 054433.
43. W.-C. Yue, Z. Yuan, Y.-Y. Lyu, S. Dong, J. Zhou, Z.-L. Xiao, L. He, X. Tu, Y. Dong, H. Wang, W. Xu, L. Kang, P. Wu, C. Nisoli, W.-K. Kwok and Y.-L. Wang, Phys. Rev. Lett., 2022, 129, 057202.
44. R. Melko and M. Gingras, J. Phys.: Condens. Matter, 2004, 16, R1277.
45. V. Schánilec, B. Canals, V. Uhlíř, L. Flajšman, J. Sadílek, T. S. Šikola and N. Rougemaille, Phys. Rev. Lett., 2020, 125, 057203.
46. A. Libál, C. Nisoli, C. J. O. Reichhardt and C. Reichhardt, Phys. Rev. Lett., 2018, 120, 027204.
47. A. Le Cunuder, I. Frérot, A. Ortiz-Ambriz and P. Tierno, Phys. Rev. B, 2019, 99, 140405.
48. C. Nisoli, New J. Phys., 2014, 16, 113049.
49. J. Loehr, A. Ortiz-Ambriz and P. Tierno, Phys. Rev. Lett., 2016, 117, 168001.
50. C. Rodríguez-Gallo, A. Ortiz-Ambriz and P. Tierno, Phys. Rev. Res., 2021, 3, 043023.
51. C. Rodríguez-Gallo, A. Ortiz-Ambriz, C. Nisoli and P. Tierno, New J. Phys., 2023, 25, 103007.
52. A. Libál, D. Y. Lee, A. Ortiz-Ambriz, C. Reichhardt, C. J. O. Reichhardt, P. Tierno and C. Nisoli, Nat. Commun., 2018, 9, 4146.
53. C. Nisoli, Phys. Rev. Lett., 2018, 120, 167205.
54. E. Östman, H. Stopfel, I.-A. Chioar, U. B. Arnalds, A. Stein, V. Kapaklis and B. Hjörvarsson, Nat. Phys., 2017, 14, 375.

---