Degeneracy and hysteresis in a bidisperse colloidal ice

We use numerical simulations to investigate the low-energy states of a bidisperse colloidal ice, realized by conﬁning two types of magnetic particles into double wells of different lengths. For this system, theoretical calculations predict a highly degenerate ground state where all the vertices with zero topological charge have equal energy. When raising the applied ﬁeld, we ﬁnd a re-entrant transition where the system passes from the initial disordered state to a low-energy one and then back to disorder for large interaction strengths. The transition is due to the particle localization on top of the central hill of the double wells, as revealed from the position distributions. When we decrease the applied ﬁeld, the system displays hysteresis in the fraction of low-energy vertices, and a small return point memory by cycling the applied ﬁeld.

In spin ice materials, the dipole moments carried out by rare-earth ions sit on a lattice of corner-sharing tetrahedra. At each vertex, these spins have equal distances, and at low temperature they tend to follow the "ice rules" [17], where two spins point towards the center of the tetrahedron and two away from it. The multiple configurations of these low-energy states give rise to the GS degeneracy. However, when the spin lattice is projected onto a two dimensional (2D) plane, as in lithographically designed artificial spin ice systems (ASIs) [18][19][20][21][22], some features may be lost due to the reduced dimensionality [23]. In a square ASI, the distance between the ferromagnetic islands at each vertex is not the same, since opposing spins are separated by greater distances than adjacent ones, and the corresponding interaction energies are not equivalent. As a result, the degeneracy is lost, and the GS becomes an antiferromagnetic order filled only by one type of vertices (type 3, see Fig. 1). A way to recover the GS degeneracy in such systems was suggested by Möller and Moessner [24], based on the idea of introducing a height displacement between the ferromagnetic islands. Changing the elevation of only two of the four nanoislands at each vertex effectively relieves degeneracy, but it requires a complex nanofabrication procedure, as demonstrated recently [25]. Another possibility relies on the use of bicomponent systems [26,27], characterized by ferromagnetic nanoislands with different size and distances, and thus interaction strengths.
The artificial colloidal ice is an alternative soft matter system to investigate the complex physics emerging from geometric frustration, by using confined colloids as analog of interacting spins [28][29][30][31][32][33][34][35][36]. For lattices characterized by a single coordination number c N , it was shown both by theory [37] and experiments [38] that the collective interactions between the particles lead to similar vertex energetics than ASIs. Thus the square colloidal ice, similar to the square ASI, also has its GS degeneracy lifted. However, the flexibility of the colloidal system in changing the structure and tuning the pair interactions enables us to search for alternative solutions to this problem.
Here we use numerical simulations to explore the lowenergy states and the degeneracy of a bidisperse colloidal ice, composed of two populations of particles and double wells. We choose to explore this particular geometric arrangement since single vertex calculations of the magnetic interactions between the particles predict a degenerate GS, where all vertices that satisfy the ice rule have equal energy. We find that, by raising the interaction strength, the fraction of GS vertices first increases but later decreases reaching again a disordered state. We thus uncover a novel re-entrant effect which distinguishes the colloidal system with mobile particles from the ASI [39], and where the applied field can be used to disorder the system by raising its amplitude. This effect can be used to induces hysteresis in the fraction of GS vertices. We note that memory effects in colloidal [40] and artificial [41,42] spin ice systems have been reported in a few works, however those systems where driven by an in plane field rather than having tunable repulsive interactions induced by a perpendicular field.

II. NUMERICAL SIMULATIONS
We perform Brownian dynamics simulations of a 2D lattice of double wells, with lattice constant a = 33 μm and each well featuring a central hill [43,44]. For the monodisperse case ( Fig. 1) the traps have length l 1 = 23 μm and they contain one paramagnetic colloid of diameter d = 2 μm and magnetic volume susceptibility χ 1 = 0.5. On the other hand, for the binary mixture we keep the values χ 1 and l 1 for half of the particles and traps while, for the rest we use χ 2 = 0.0675 and l 2 = 30.356 μm and keep constant the particle diameter d. In both cases, for each particle i at position r i we integrate the overdamped equation of motion: where γ = 0.032 pN s μm −1 is the viscous drag. The first term on the right-hand side of Eq. (1) is given by, μ 0 is the induced moment from the applied field B, V = π d 3 /6 the particle volume, μ 0 = 4π × 10 −7 H/m the magnetic permeability of the medium andr i j = (r i − r j )/|r i − r j | is the relative position unit vector between two particles (i, j). The second term in Eq. (1) is the force F T i that acts on a colloid i due to the potential well, and it is given by where r and r ⊥ are components of a vector r parallel and perpendicular resp. to the line of length λ that joins the two minima in the double well. The stiffness k trap = 6 × 10 −4 pN/nm keeps the particle confined to the elongated region around the center of the trap, and k hill = 5 × 10 −6 pN/nm creates a potential hill that pushes the particles away towards one of the two potential wells. Finally, η represents a random force due to thermal fluctuation, with zero mean, η = 0 and delta correlated, Here, T = 300 K is the ambient temperature. The simulations are performed using a constant time step of 10 ms. We use always a square sample with a number of vertices on each side equal to N = 20, except for the calculations in Fig. 2(b) where we vary N ∈ [1,40]. We increase B linearly up to B = 300mT, at a rate of α B = 0.035 mT/s for the initial simulations, and later change the rate when investigating hysteresis effects. Finally, we neglect hydrodynamic and electrostatic interactions between the particles since in the experiments such interactions are screened by the topographic double wells.

III. THE MONODISPERSE COLLOIDAL ICE
The essential features of the colloidal ice are shown in the schematic in Fig. 1(a). It is composed of a square lattice of double wells, each filled by one paramagnetic colloidal particle that is confined due to gravity. An external magnetic field B perpendicular to the particle plane induces in each particle a dipole moment m i ∼ χ i B, and pair of particles (i, j) located at a distance r = |r i − r j | experience a repulsive dipolar interaction: which depends on both, the field amplitude and particle susceptibility. Thus particles located in close traps on the same vertex will repel each other, but the interaction strength is such that while they can cross the central hill, they cannot leave the double well. Geometric frustration in the system arises from the arrangement of the double wells, that do not allow the simultaneous minimization of all pairwise interactions. The mapping between the colloidal and spin ice systems can be done by assigning a pseudospin to each double well, such that it points where the particle is located [28]. As shown in Fig. 1(b), this mapping allows to construct a set of vertex rules, and for the square lattice, six energetic configuration of the particles are possible. Also, to each vertex one can assign a topological charge Q = 2n − c N [45], where n is the number of pseudo spins that point toward the vertex center, and c N the coordination number of the lattice, which for the square is c N = 4. In this formalism, the vertices that satisfy the ice rules are those which cancel the topological charge, Q = 0, i.e., the type 3 and type 4, while low (type 1 with Q = −4, type 2 with Q = −2) and high (type 6 with Q = +4, type 5 with Q = +2) energetic vertices break the ice rule and represent topological defects. The topological nature of such defects arises from the fact that they can disappear only when annihilating with other defects of opposite charge, for example type 2 with −2 and so on. In Fig. 1(c), we show the simulation results for a square lattice composed of N = 20 vertices per side, and plot the fraction φ of different vertex types versus the amplitude of the applied field. We find that for the monodisperse ice, the system displays a clear tendency to follow the ice rules (Q = 0) at high interaction strength, and it reaches a GS filled only by type 3 vertices for B > 115 mT. This is also in agreement with the energy hierarchy showed in Fig. 1(b), and a similar situation is observed for artificial spin ice systems, where the GS is characterized by type 3 vertices that form local loop configurations with alternating chirality [46][47][48][49]. As a result, the degeneracy is lost, and the GS has a twofold antiferromagnetic order.

IV. THE BIDISPERSE COLLOIDAL ICE
To recover degeneracy, we designed a colloidal ice characterized by two types of traps and double wells such that the energetic weight of vertices of type 3 and 4 are the same. The idea is illustrated in Fig. 2(a), where we use two types of particles characterised by two different magnetic susceptibilities (χ 1 , χ 2 ) and, thus, two different induced moments (m 1 , m 2 ) under the same applied field. On each vertex, the particles with lower magnetic susceptibility (χ 2 ) are placed inside longer double wells (length l 2 ) that are located closer, i.e., at a distance d 2 < d 1 . Instead the double wells occupied by the particles with larger magnetic susceptibility (χ 1 > χ 2 ) are located far away when pointing toward the vertex (distance d 1 ).
For arbitrary values of susceptibilities and trap lengths, this conformation breaks the symmetry, and splits the vertex energy as shown in Fig. 2(b): type 2 gives rise to 2a and 2b, type 3 to 3a and 3b, and type 5 to 5a and 5b. To engineer a degenerate vertex hierarchy, we fix the values (l 1 , χ 1 ) and a, and we solve the system of equations, u 4 (χ 2 , l 2 ) =ū 3a (χ 2 , l 2 ) =ū 3b (χ 2 , l 2 ), beingū = u the sum of the dipolar interactions between the vertex elements [Eq. (3)], where (l 2 , χ 2 ) are used as free parameters. The solution ensures that the pair interactions between opposing and adjacent particles at a vertex are the same and thus also the energetic weight of the type 3a, 3b and type 4, Fig. 2(c). The previous calculations are performed on individual vertices, but do not reflect the collective nature of the colloidal ice and the fact that the paramagnetic colloids interact beyond the nearest neighbor level. To take into account these effects, we consider a lattice made of N vertices with fixed particles arranged to create a state with only type 3 or type 4. In these two configurations, we calculate the mean energy per vertex, u /N, to distinguish whether these two configurations are energetically equivalent. As shown in Fig. 2(d), the vertex energy grows with system size reaching a constant value of ∼69k B T for N > 1500 vertices. However, the energetic difference between the types 4 and 3 is non zero, but it raises to ∼0.4k B T, as seen in the Fig. 2(d). This effect results from the collective interactions between the particles, that lower the energy of the type 3 vertices with respect to the type 4. Thus, at high interaction strength, the bidisperse colloidal ice is expected to reach a GS filled only by type 3.
To confirm this hypothesis, we run numerical simulations of the bidisperse system, see Fig. 3. Instead, we find that even though the fraction of type 3 vertices first raises rapidly from φ = 0.1 to φ = 0.5 at B ∼ 100 mT, above B ∼ 110 mT, it starts decreasing until it reaches a relative low value of φ = 0.2 even for the highest applied field of B ∼ 300 mT. Further, also the type 4 vertices first decrease from φ = 0.25 to φ = 0.13, but later increase again to recover almost their initial disordered value of φ ∼ 0.25. Interestingly, intermediate energy vertices with charge Q = ±2, which are topologically bound, appear with nearly the same probability as the ice rule vertices. This is in contrast with high and low charge (Q = ±4) vertices, which completely disappear already above B = 30 mT, as in previous observations with the monodisperse systems ( Fig. 1 and Ref [43]). In addition, this transition occurs at a similar magnetic field value when the monodisperse system reaches the GS, Fig. 1(c).
To understand this re-entrant behavior, we have analysed the distribution of particles positions within the double wells, shown in Fig. 4(b), where representative snapshots of the simulation system are shown for different values of B in Fig. 4(a). For low field, both types of particles remain confined within one of the two wells, and this behavior lasts until B ∼ 80 mT. However, above this value we observe that the particles with high magnetic susceptibility and located in the short traps tend to localize close to the central hill rather than away from it, Fig. 4(b) second column. This effect is remarkable, and it raises much more the energy of the type 4, decreasing their statistical fraction within the ensemble. Above B > 110 mT the induced dipolar interactions are so strong that they force also the localization of the particles with smaller magnetic susceptibility close to the central hill, Fig. 4(b) third column. This localization becomes more pronounced at higher field strength as shown by the narrow particle distribution around the central hill for B = 200 mT. The system thus tends to disorder back to its initial configuration, with the only difference being the absence of type 1 and 6 vertices that causes an increase in the final fraction of type 3, while type 2 and 5 reaches practically the same initial value of φ.
The discovered re-entrant behavior is rather robust, as it was observed for a wide range of field rate α B . By lowering α B , the fraction of type 3 vertices display a maximum at lower field amplitude and later reaches higher fraction φ, whereas increasing the rate causes this re-entrant phenomenon to occur later in time. In the latter case, the system does not have time to reorganize into low-energy state, and the fraction of type 3 vertices decreases rapidly. However, the final vertex count remains the same regardless of the field change rate. Thus this phenomenon represents an alternative way to disorder the system by raising now the applied field, rather than resetting the colloidal ice by repeating the experiment and recollocating the particles randomly within the double wells.
We now explore the presence of hysteresis and memory in our system when cycling B. We measure this effect by considering the fraction GS vertices (type 3) and how it changes by first increasing and then decreasing the applied field. Some representative results are shown in Fig. 5(a) for different rates α B . In general, we find that when the field is cycled sufficiently fast (α B > 0.5 mTs −1 ) the fraction φ tends to retrace itself and there is a small detectable hysteresis. In contrast, for smaller rates (α B < 0.1 mT s −1 ), the difference between both curves is significant as shown in Fig. 5(a). To better understand the rate dependent behavior, we plot in Fig. 5(b), such difference A between the area of the two curves versus the magnetic field rate α B . Effectively, we find that the behavior is power law with an exponent ∼−0.1 [inset Fig. 5(b)], while A → 0 starting from α B ∼ 0.25 mT s −1 . This rate dependent behavior is the opposite than the one observed for other bistable systems, like ferromagnets, where the area of the hysteresis loop increases with the rate of the field switching [50][51][52].
Another phenomenon that could arise in our system is the presence of return point memory, namely the possibility that our bidisperse colloidal ice could reproduce a microstate, i.e., a given particle configuration, when the field is cycled after the first minor loop. To quantify this effect, we follow previous works on colloidal [40] and artificial spin ice [42], and calculate the spin overlap order parameter defined as [53] i . Form the colloidal position, we associate an effective spin s to each of the bistable traps composing the lattice. Thus s (n) i represents the magnetic moment associated to a bistable trap after a minor loop (n), while s (n−1) i is calculated one minor loop earlier. This order parameter quantifies the presence of reproducible particle microstates (q s = 1) after a certain number of minor loops starting from a virgin curve. In contrast, q s = 0 indicates the total absence of correlation between the configuration of two subsequent loops n and n − 1. As shown in Fig. 5(c), the spin overlap parameter has two dramatically different behaviours as the magnetic ramp rate α B is changed. For α B < 0.7 mTs −1 the overlap increases, without reaching the fully reversible state (q s = 1) after 7 minor loops. As the field rate increases, the overlap becomes smaller, which is expected since the system lives in more excited states, and has access to a larger phase space. However, surprisingly for α B = 0.1 mT s −1 , we find that the overlap suddenly increases in the first cycle, but then decreases until it almost reaches the same steady value as the case of α B = 0.069 mT s −1 . These results indicate that our colloidal ice does not present a strong memory of its previous state, and such effect could be attributed to the presence of noise as thermal fluctuations of the particles within the double well or disorder resulting from the different magnetic couplings between the two types of particles.

V. CONCLUSIONS
In this paper, we have investigated the low-energy states of an artificial colloidal ice composed of a binary mixture of particles and traps. Before this work, the colloidal ice had mainly been explored with monodisperse particles, and the emerging physics had been studied extensively for different types of lattices and interactions strength. However, novel effects may arise when changing the types of particles or double well traps. In particular, we uncover a re-entrant behavior which is induced at high interaction strength, and that can be used as an effective way to disorder the system. Such effect is absent in the monodisperse system, where the particles localize in one of the two double wells at all interaction strengths. Furthermore, it produces hysteresis and small memory effects when changing the magnetic field rate.
This work opens different future directions. For example, the idea of using a binary mixture of magnetic colloids to fill bistable traps suggests the exciting possibility of employing a polydisperse suspension which will modify the magnetic coupling between the particles at parity of the applied magnetic field. This is equivalent to varying the amount of disorder in the system, without altering its geometrical arrangement and thus, it could be used to uncover the emergence of glassiness on a periodic lattice [54][55][56].