## Abstract

Assemblages of microscopic colloidal particles exhibit fascinating collective motion when energized by electric or magnetic fields. The behaviors range from coherent vortical motion to phase separation and dynamic self-assembly. Although colloidal systems are relatively simple, understanding their collective response, especially under out-of-equilibrium conditions, remains elusive. We report on the emergence of flocking and global rotation in the system of rolling ferromagnetic microparticles energized by a vertical alternating magnetic field. By combing experiments and discrete particle simulations, we have identified primary physical mechanisms, leading to the emergence of large-scale collective motion: spontaneous symmetry breaking of the clockwise/counterclockwise particle rotation, collisional alignment of particle velocities, and random particle reorientations due to shape imperfections. We have also shown that hydrodynamic interactions between the particles do not have a qualitative effect on the collective dynamics. Our findings shed light on the onset of spatial and temporal coherence in a large class of active systems, both synthetic (colloids, swarms of robots, and biopolymers) and living (suspensions of bacteria, cell colonies, and bird flocks).

- self-assembly
- colloids
- active matter
- flocking
- Collective Behavior
- Synchronization

## INTRODUCTION

Colloids, that is, suspensions of microscopic particles, play an important role in everyday life and are crucial for many industries, from pharmaceutical to medicine and nanotechnology. A multitude of forces govern their behavior, from steric repulsion, electrostatic, magnetic, and van der Waals forces to gravity, hydrodynamics, etc. (*1*). Despite their seeming simplicity, colloidal systems may exhibit fascinating collective behavior when they are driven out of equilibrium by external magnetic or electric fields (*2*). Recent studies have led to the discovery of self-assembled swimmers, microrobots, and spinners in the system of driven ferromagnetic colloids (*3*–*5*), self-assembly of dynamic rings and vortices in electrostatically driven metallic particles (*6*), reconfigurable active patterns in metal-dielectric Janus colloids (*7*), formation of complex vorticity patterns of paramagnetic colloids in triaxial alternating fields (*8*, *9*), emergence of large vortices in the system of colloidal dielectric rollers energized by an electric field (*10*, *11*), ionic crystals and lane formation in bidisperse colloidal suspensions (*12*), crystal formation of light-induced colloidal surfers (*13*), and many others.

Studies of driven colloids have highlighted a close relation to a broad class of systems termed active matter: assemblies of self-propelled particles capable of transducing energy stored in the environment into mechanical motion. Active matter has a strong propensity toward the onset of large-scale collective behavior stemming from a simple alignment interaction between the self-propelled agents. In living systems, this collective behavior is exemplified by bacterial swarms, bird flocks, fish schools, etc. (*14*). Many aspects of collective behavior, at least on a qualitative level, are captured by a paradigmatic Vicsek model for the interacting self-propelled point particles (*15*). Although significant progress has been made in the theoretical understanding of model active systems [for example, see recent research (*16*–*18*)] in many experimental realizations of active matter (*19*–*24*), a multitude of long-range interactions (for example, hydrodynamic, magnetic, electrostatic, elastic, etc.) often obscure and complicate the simple phenomenology exhibited by the Vicsek model.

The emergence of large vortices and propagating densification fronts (so-called Vicsek bands) has been reported in the system of rolling colloids energized by the dc electric field (*10*, *11*). Self-propulsion is an outcome of the spontaneous rotation of a dielectric sphere in a conducting fluid when it is exposed to a dc electric field [Quincke effect (*25*)]. Although many experimental findings from the study of Bricard *et al*. (*10*, *11*) were successfully reproduced by coarse-grained continuum theory and simplified particle dynamics model, individual interactions between the particles are hard to quantify, especially the effects of the dielectric polarization and hydrodynamic forces.

Here, we report on the discovery of ferromagnetic flocking colloids. The experimental system is realized by ferromagnetic colloidal spheres energized by a vertical alternating (ac) magnetic field. Particles are settled on a slightly concave bottom surface of a liquid-filled container (see Materials and Methods). The rolling motion appears due to a spontaneous symmetry breaking of the clockwise/counterclockwise rotation experienced by a magnetic sphere in a uniaxial field (*5*). Depending on the frequency of the external magnetic field, a sequence of transitions can be observed: from gas-like motion of individual particles at low frequencies to the onset of flocking and global rotation, followed by a reentrant flocking and gas-like state for increasing frequency. Like in living systems, flocking behavior is characterized by a spontaneous onset of coherently moving groups of many particles. However, the flocks do not keep their identity; they often change their direction, break up, and reassemble. Our discrete particle simulations have provided insight into the behavior of our experimental system. We have established that the flocking is an outcome of mostly collisional and magnetic interactions between the particles, whereas the long-range hydrodynamic forces do not have a qualitative effect on the collective behavior. We have also observed a strong correlation between the onset of flocking and the synchronization of particle motion by the applied ac magnetic field. Our work provides fundamental insights into a broad class of active systems, both synthetic and living, where collective motion is caused by a subtle interplay of long- and short-range interactions. Control and prediction of collective behavior in out-of-equilibrium colloidal systems also lead to a better understanding of fundamental aspects of dynamic self-assembly in general and novel properties of functional colloidal materials.

## RESULTS

In our experiment, a dispersion of ferromagnetic nickel spheres was deposited on the bottom of a liquid-filled container. The bottom surface of the container was formed by a concave glass lens. The lens was used to drive the particles toward the center and to prevent them from accumulating at the container walls. A uniform vertical alternating magnetic field was used to energize the particles (see Materials and Methods for details). Typically, up to 400 particles were used in each experiment. The schematics of our experiment are shown in Fig. 1A.

The main experimental results are summarized in Figs. 1 to 4. Figure 1B illustrates the trajectories of individual particles. Whereas some particles roll persistently, other particles may stop and exhibit almost no rolling for extended periods of time. This variability in particle behaviors is likely due to the dispersion of the particle sizes and, correspondingly, magnetic moments as well as shapes and roughness. A sequence of experimentally observed transitions is illustrated in Fig. 2. The transitions between the states are quantified by the corresponding order parameter for in-plane rotational collective motion, φ_{R}, defined as (*26*)(1)where are the in-plane unit vectors in angular direction, *v*_{i} are the in-plane unit velocities of the individual particles, and *N* is the total number of particles. The order parameter attains its maximal values ±1 for pure rotation. Time evolution of the order parameters for the different regimes after reaching a steady state is shown in Fig. 3B.

The experimental transitions can be summarized as follows: For low frequencies *f* < 30 Hz, a gas-like state of randomly moving particles was observed (see Fig. 2, A and E, and movie S1). While individual particles roll, the order parameter for rotational collective motion ϕ_{R} was close to 0 (see Fig. 3B), and their individual velocities seem uncorrelated (see Fig. 3A). Correspondingly, for this state, the time-averaged rotational velocity is close to 0 (see Fig. 2E).

Increase in the frequency of the magnetic field resulted in a markedly new phenomenon: the onset of a behavior, which is reminiscent of bird flocks and fish schools, where large groups of particles start to move coherently and form well-defined flocks. These flocks are not static: they break up and reassemble in different locations (see Fig. 2, B and F, and movies S2 and S3). Correspondingly, the order parameter ϕ_{R} fluctuates over time at around a value of 0.2.

Further increase in the field frequency results in the increase in size and persistence time of the flocks. It eventually leads to the emergence of a vortex spanning the entire system and continuously rotating in a certain direction (see Fig. 2, C and G, and movies S4 and S5). The order parameter attains a maximal value at about 0.5 and shows fluctuations around the mean value (Fig. 3B). The sense of rotation is determined by an initial particle configuration and changes from experiment to experiment. The core of the vortex is characterized by an almost linear azimuthal velocity profile (Fig. 3C). Further increase in the frequency results in the breakdown of the coherent vortex motion and formation of a reentrant flocking state (see Fig. 2, D and H, and movies S6 and S7) and, finally, a gas-like state (see movie S8).

Figure 3A displays the spatial correlation functions for different values of the driving frequency *f*. Correlation length provides valuable information on the changes in the system order; however, the absolute value of it depends on the particle number density. As anticipated, the correlation length is minimal in the gas phase (*f* = 20 Hz). Then, the correlation length increases with the frequency and attains a value of 15 to 20 particle sizes in the flocking state (*f* = 35 Hz). It becomes comparable with the system size in the vortex state (*f* = 40 Hz). Moreover, the negative values in the correlation function indicate a large-scale rotational motion leading to anticorrelations. Then, with further increase in *f*, the correlation length decreases again. The inset in Fig. 3A shows that the correlation length peaks at the frequency *f* ≈ 39 Hz. In this case, the correlation length is possibly limited by the system size.

## DISCUSSION

### Individual particle dynamics

To obtain insights into the onset of collective behavior, we first investigated the dynamics of individual particles for very low concentrations (25 to 40 particles). The results are summarized in Fig. 4. Although individual particles roll and exhibit self-propelled motion, the trajectories resemble a random walk due to shape imperfections. The averaged mean square displacements of individual particles at a 40-Hz excitation frequency are shown in Fig. 4A. The behavior is consistent with the linear law characterized by translation diffusion coefficient *D*_{T}. For *f* = 40 Hz, we obtained *D*_{T} ≈ 1 mm^{2} s^{–1}. From the measurements of the translation diffusion, we infer the effective rotational diffusion coefficient *D*_{R} using the well-known relation for the self-propelled particles(2)where *V* is the typical particle velocity (*27*). The estimate yields the following value for the rotational diffusion *D*_{R} ≈ 3 to 4 s^{–1}, with a corresponding persistent length *l* = *V*/*D*_{R} ≈ 1 mm. We also estimated the rotational diffusion directly from the variance 〈δϕ^{2}〉 of in-plane particle velocity orientation angle φ(*t*), with δϕ = ϕ(t) − ϕ(0). This variance behaves as 〈δϕ^{2}〉 ≈ 2*D*_{R}*t* for time *t* ≫ 1 s (see inset in Fig. 4A). Using this procedure, we obtained *D*_{R} ≃ 3.5 s^{–1} for the excitation frequency *f* = 40 Hz. As the frequency of the excitation field increases, the rotational diffusion coefficient systematically decreases; however, variation of the magnetic field amplitude has not produced significant changes in the rotational diffusion coefficient.

Figure 4B shows the particle velocity distribution function *P*(*V*). The distribution is peaked at around a value of *V*_{p} ≈ 4 mm s^{–1}, which depends on the frequency itself. There is a significant population of particles that almost does not roll (likely due to lower values of the magnetic moment). This population is manifested by a plateau near *V* = 0.

Because of the contact with the surface, a rotating particle will roll. In the limit of no-slip contact with the surface, the rolling speed *V*_{0} would be *V*_{0} = ω*a*, where ω = 2π*f*. However, we have found from our experiments that this theoretical speed value is not attainable (see Fig. 4C). The particle rolling speeds are distributed around mean value , which is of the order in the range of our frequencies. Two main factors may contribute to the speed reduction. First, the fluid may work as a lubricant, leading to partial slip between the particle and the bottom. Second, a moving rotating particle near a solid surface will experience a lift force *F*_{L}, which is of the order of the gravity for our conditions. In the Stokes limit, the lift force on a rotating sphere is *F*_{L} = πρ_{f}*a*^{3}ω*V*, where ρ_{f} is the fluid density, and *V* is the particle translation velocity (*28*). A particle loses contact with the substrate when the lift force *F*_{L} becomes comparable with the gravity force *F*_{g} = *mg* = 4/3πρ_{p}*a*^{3}, where *m* is the particle mass, ρ_{p} is the particle density, and *g* is the gravity. By using *V* = *V*_{0} = ω*a*, we obtain the lift-off condition . For nickel particles (ρ_{p}/ρ_{f} ≈ 9) with the radius *a* ≈ 60 to 100 μm, we find that the lift force will exceed gravity for frequencies in the range of *f* ≈ 60 to 100 Hz. This implies that the particle cannot be in a protracted contact with the surface above these frequencies. Correspondingly, it should move parallel to the surface with a speed smaller than *V*_{0}.

### The onset of rolling

The onset of rolling is associated with the spontaneous symmetry breaking of the clockwise/counterclockwise rotation of an individual particle in an ac vertical magnetic field. Non-negligible particle inertia is a main cause of the directed rotation (*5*). By neglecting fluctuations, the evolution of the particle angle θ with respect to vertical direction (see Fig. 1A) as a function of time in a vertical uniaxial ac magnetic field *H*_{z} = *H*_{0} sin(ω*t*) is described by the following equation of motion(3)where *H*_{0} is the magnitude of the applied field, ω = 2π*f*, μ is the magnetic moment of the particle, and *I* = 2/5*ma*^{2}, and α_{r} = 8πη*a*^{3} are the moment of inertia and the rotational drag coefficient, respectively (*a* is the particle radius, *m* is the particle mass, and η is the fluid dynamic viscosity). If the particle inertia is negligible compared to the viscous forces, the first term in Eq. 3 can be dropped, and it is reduced to(4)

This equation can be solved exactly, leading to the family of periodic solutions(5)with ζ = μ*H*_{0}/α_{r}ω and an integration constant *C*_{0} determined by the initial conditions. Thus, these solutions would correspond to oscillations around some mean value of the angle θ. Any nonzero particle mass would lead to the breakdown of the family and locking of the angle to the values of either θ = 0 or θ = π. The effect of a small inertia can be demonstrated analytically by perturbing the solution given by Eq. 5. It is practical to consider the ε = *I*/α_{r} ≪ 1 limit. For simplicity, we set in Eq. 5 *C*_{0} ≪ 1, leading to θ = *C*_{0} exp(ζcos(ω*t*)). For small ε, we use the solution in the form θ = *C*_{0}(ε*t*)exp(ζcos(ω*t*)). Plugging this into Eq. 3, in the first order in ε, we find, after a simple algebra, that *C*_{0} ~ exp(−λ*t*), where . Thus, the particle will align in a vertical direction with the rate ~*I* for any small particle inertia.

However, the steady states θ = 0, π may become unstable with the increase of the particle inertia and/or field magnitude. This can be seen by linearizing Eq. 3 around θ = 0, π. Linearizing Eq. 3 after a proper scaling assumes the form(6)where are the normalized friction and field magnitude, respectively, and . Equation 6 can be reduced to the Mathieu equation whose solutions are well studied and tabulated (*29*). Correspondingly, Eq. 6 can be solved in terms of the Mathieu special functions (*29*). The time behavior of the Mathieu functions is determined by the Mathieu characteristic exponent *v*, which is the function of the parameters *p* and *q*. Thus, by transforming Eq. 5 to a canonical Mathieu form, one derives a condition for the instability of the steady states θ = 0, π(7)

The condition provides a critical value of the friction coefficient when locking with the field direction becomes unstable. In the limit of small inertia *I*, (that is, *p* = α_{r}/ω*I* ≫ 1), one finds from Eq. 7 that Λ < 0, consistent with our asymptotic analysis above. Thus, the locked states are stable. The threshold value of the normalized field *q* versus friction *p* obtained from the condition Λ = 0 is shown in Fig. 5A. One sees that the critical field value increases with the increase of the friction coefficient *p*. However, even for *p* → 0, there is a critical value of the field *q* ≈ 0.454 needed for instability of the locked state.

Furthermore, by decomposing the driving term into two counterrotating waves and neglecting one of the terms in this expression (for example, the second), we obtain a damped pendulum equation(8)

After making a transformation to a corotating frame, ψ = θ − ω*t*, one obtains an equation without explicit time dependence(9)

The rotating state corresponds to a fixed point of Eq. 9. The corresponding fixed points exist for μ*H*_{0} > 2α_{r}ω (or *q* = 2*p*), that is, for relatively small friction. Thus, Eq. 3 describes the rotation of the particle in clockwise direction with the frequency ω as long as *q* > 2*p*. A similar consideration is valid for the opposite direction of rotation. The direction of rotation is determined by the initial conditions and may change due to interactions with other particles. The rotating state existence boundary obtained from the condition *q* = 2*p* is shown in Fig. 5A by a thin dashed line.

Moreover, even a single particle may exhibit stochastic switching of its rotation direction if the field magnitude *H*_{0} becomes large enough. The approximate boundary for the chaotic regimes for very small friction α_{r} → 0 is given by the resonance overlap criterion (*30*), which is *q** = μ*H*_{0}/ω^{2}*I* ⪆ 1/2 in our notation. However, this condition for *q** underestimates the numerically obtained transition values somewhat by a factor of 3 (see Fig. 5A). For increased values of the friction coefficient α_{r}, stronger magnetic fields are required for the onset of these chaotic regimes. The existence boundary of chaotic regimes, obtained by numerical solution of Eq. 3, is shown in Fig. 5A by black diamonds.

Although our system of ferromagnetic rollers appears seemingly related to the Quincke colloidal rollers energized by the dc electric field (*10*, *11*), we want to emphasize a fundamental difference. The propulsion speed of the Quincke colloidal rollers depends on the magnitude of the electric field and vanishes at a critical field value. In contrast, for the ferromagnetic rollers energized by the ac magnetic field, the rotation speed does not depend on the field magnitude *H*_{0} and depends only on the field frequency *f*. Moreover, even individual ferromagnetic rollers exhibit complex dynamics characterized by spontaneous direction reversals. This particle behavior is not present in the electrostatic system of Quincke rollers. These intrinsic chaotic dynamics lead to the complexity of the phase diagram for ferromagnetic rollers.

Even in the absence of noise, the analysis of individual particle dynamics for ferromagnetic rollers provides nontrivial insight into the collective behavior of many interacting particles. We show below (see Fig. 5) that the transitions between different collective states are closely related to the boundaries in the phase diagram for individual particle (Fig. 5A).

### Computational modeling of collective dynamics

To obtain insights into the onset of collective motion and to identify the most important particle interactions, we performed computational modeling of a colloidal system driven by a vertical magnetic field (see Materials and Methods for details). A phase diagram spanned by the amplitude and the frequency of the magnetic field is shown in Fig. 5B. In good agreement with the experiment, we reproduced the three main states (gas, flocking, and vortex) and their correct sequence of transitions for increasing frequency of the applied field (see Fig. 6 and movies S9 to S12). Moreover, the temporal evolution of “computational flocks” closely resembles the experimental ones (as can be seen in Fig. 7A and movie S10), showing a polar order parameter fluctuating around ϕ_{R} ≈ 0.3. In agreement with the experimental results, we observed that in the vortex state, all particles move in a single flock, with a rotational order parameter ϕ_{R} close to 1 and a “void” with a radius *r*_{v} ≈10a, (see movie S8). Similar “hollow” vortices were previously observed in the model of interacting self-propelled particles (*31*). The individual particle behavior (see Fig. 5A) helps to understand the observed collective phenomena. The vortex with the maximal correlation length exists very close to the stability limit of the locked state for an individual particle. In this regime, the individual particles exhibit the most robust rotation, which leads to the emergence of a single vortex. We also observed vortices slightly below the stability limit. Whereas in the absence of fluctuations, individual particles do not rotate, interactions with other particles and the noise can cause the particle to spin. However, their propagation velocity decays with the increase of the frequency. This nonsteady rotation of an individual particle for high frequency leads to transient flocks in the case of collective dynamics. Accordingly, we find a flocking state for low frequencies near the transition from chaotic rotation toward persistent rolling.

We also emphasize a fundamental difference between the two flocking states observed for low and high frequencies. According to Fig. 5A, low-frequency flocking is observed in the regime when an individual particle exhibits a chaotic behavior manifested by spontaneous direction reversals. This implies that external noise does not play a significant role and the emergence of the vortex state is due to the suppression of the intrinsic stochasticity of individual particle motion. We checked through simulations that in this frequency range, the flocking-like motion exists even in the absence of external noise. In contrast, for high-frequency flocking, individual particles would be in a steady state, with their orientation along the field direction. Figure 5C displays the system phase diagram without external noise (*D*_{R} = 0). Without noise, a stationary phase with hexagonal order emerges (see inset in Fig. 5C). The transition occurs immediately at the boundary between the stationary and rotating states for individual particles (see again Fig. 5A). This observation implies that the high-frequency flocking is a noise-activated process and highlights again that the individual particle dynamics allows one to draw conclusions on the collective phenomena. In the absence of noise, the vortex state emerges only if an individual particle exhibits a robust rotation.

For the results shown in Figs. 5 and 6, we neglected long-range hydrodynamic forces exerted on the fluid by moving rotating particles. The fluid effect was only included in the translational and rotational drag cofficients α_{t} and α_{r}. A rolling particle will impose a net force on the fluid (force monopole or stokeslet). In bulk, the stokeslet velocity would decay as 1/*r*, where *r* is the distance from the origin (*32*). Near a nonslip boundary, the stokeslet will be screened, and in the far field, it will decay as 1/*r*^{2} due to an image stokeslet opposite in sign (*33*). We implemented the hydrodynamic interactions between the particles in the far-field approximation (for the distances large compared to the particle radius *a*). Our simulations produced no qualitative effect of the hydrodynamic interactions on the overall dynamics of the system (see Fig. 7 and movie S13). This implies that for typical experimental conditions, and except for the immediate vicinity of the particle, the hydrodynamic velocity is small compared to the rolling velocity ω*a* and can be neglected. All numerically obtained characteristics of the dynamic states (rotational order parameter, tangential velocity profiles, correlation functions, and correlation length dependencies; see Fig. 7) are in qualitative agreement with the experiments (see Fig. 3). However, there are some quantitative differences. First, in the simulations, the obtained correlation length is larger than that in the experiments because all particles are engaged in a vortical motion. In the experiment, even in the case of a vortex, some particles at the periphery move randomly. Second, there is some difference between the vortex tangential velocity profiles. In the experiment, the velocity profile is almost linear within the vortex core, whereas in the simulations, we have a void (see Fig. 7B). We found that the void size shrinks to zero as the number of particles in the system increases.

### The onset of flocking and synchronization phenomenon

A fundamental question in the context of the onset of large-scale collective behavior is how the increase of spatial coherence relates to the synchronization of particle orientations through the applied ac magnetic field. To address this question, we extracted the particle orientation angles θ_{i} at each field period (so-called stroboscopic plots) from the simulation data (see Fig. 8, A to D). The corresponding distribution functions *P*(θ), plotted after a transitional phase (about 200 periods), for different frequencies are shown in Fig. 8E. Figure 8 reveals a close relation of the phenomenon of phase synchronization and the onset of large-scale collective behavior. As anticipated, in the gas phase (Fig. 8A), the particles execute random-like motion, and their orientations with respect to the field are essentially uncorrelated. Correspondingly, the angle distribution is almost uniform (see Fig. 8). With the onset of flocking, a marked transition occurs (see Fig. 8B). After a short transitional phase, all particle angles converge to two well-separated narrow bands, and the distribution function exhibits two narrow peaks of almost equal heights. The width of the peaks is determined by the noise magnitude *D*_{R}. We verified that without noise (*D*_{R} = 0), the peaks become δ-functions. These two peaks in the distribution function *P*(θ) imply that the flocking phase is composed of two almost equal fractions of clockwise/counterclockwise rolling colloids. Because we do not see transitions between the bands over the entire duration of simulations, it means that the reversal of individual particle directions in low-frequency flocking is very rare. With the transition to a vortex, one of the bands disappears (see Fig. 8C). Correspondingly, the distribution function acquires only one narrow peak. Thus, vortex motion is characterized by the complete synchronization of particle angles with respect to the applied ac field. A further increase of the frequency leads to the breakdown of the synchronized vortex motion and the onset of noise flocking (Fig. 8D). In contrast to the low-frequency flocking shown in Fig. 8B, the high-frequency flocking is noise-activated. It is manifested by a significant amount of unsynchronized particles and relatively small broader peaks in the distribution function. Overall, the observed behavior is reminiscent of the synchronization and clustering phenomena in coupled noisy oscillator systems described by various versions of the Kuramoto model (*34*–*36*).

## CONCLUSION

We have demonstrated that a seemingly simple assemblage of rolling ferromagnetic colloids exhibits marked flocking behavior reminiscent of that in living systems. We have identified and characterized three main dynamic phases: gas, flocking, and vortex, as well as the transitions between them. The main advantage of our system compared to other experimental realizations of active matter is that the particle interactions are well characterized and can be tuned by the frequency and magnitude of the magnetic field, solvent viscosity, surface curvature, etc. Moreover, each and every particle can be tracked, and its trajectory can be analyzed. Thus, our work provides insights into a variety of synthetic and living active matter realizations. It also sheds light upon the relation between seemingly unrelated phenomena: the onset of collective behavior and phase synchronization in periodically driven systems. We have demonstrated that the emergence of flocking and vortex motion is manifested by the appearance of sharp peaks in the particle orientations and their synchronization with the ac applied field. Our results also emphasize a subtle role of the rotational noise: whereas one of the flocking phases appears to be noise-insensitive, the reentrant flocking happens to be noise-activated. Thus, our work uncovers new relations between collective motion, synchronization, and self-assembly (*37*). From the materials perspective, quantification and characterization of the main physical mechanisms governing the dynamics of out-of-equilibrium colloidal systems are crucial for reaching a better understanding of the dynamic self-assembly and self-organization at the microscale.

Our computational model for ferromagnetic rollers accurately reproduced the main experimental observations. We have found that the long-range hydrodynamic forces do not have a qualitative effect on the overall system behavior. However, this does not exclude the importance of near-field hydrodynamic interactions on the dynamics of particles that are almost in contact. In this respect, more accurate characterization of the hydrodynamic forces, for example, by the lattice Boltzmann algorithm (*38*) or by multiparticle collision dynamics (*39*), would be highly desirable but is computationally prohibitive.

## MATERIALS AND METHODS

### Experimental setup and image analysis

Nickel spherical particles with an average radius of 60 μm and an average magnetic moment μ = 2 × 10^{−5} emu (Alfa Aesar Company) were used for the experiments. Particles (up to 400) were dispersed at the bottom of a cylindrical container and filled with isopropanol. The bottom of the container was formed by a concave glass lens with a radius of curvature of 52 mm and a diameter of 50 mm. The particles were energized by a uniform uniaxial alternating magnetic field, *H* = *H*_{0} sin(2π*ft*), created by a custom-built Helmholtz air-core solenoid and applied parallel to the axis of the cylindrical container. The amplitude of the ac magnetic field *H*_{0}, was in the range of 0 to 75 Oe, and the frequency *f*, was varied from 20 to 100 Hz.

The container was mounted on a microscope stage (Leica MZ9.5). The trajectories of the particles were monitored by a fast charge-coupled device camera (Redlake MotionPro). We recorded image sequences (resolution, 1280 × 1024 pixels) at a frame rate of 50 frames/s. Image and data analysis of the time sequences were performed by ImageJ, MatPIV, and custom scripts. Particle tracking was carried out by a MATLAB script based on the Crocker and Grier algorithm (*40*).

The direction of a particle motion was color-coded in visualizations using the following prescription: [red, green, blue] = [(1 + sin(ϕ))/2, (1 + cos(ϕ))/2,0], where φ is the angle of the in-plane velocity vector with respect to the *x*- axis. The correlation length was determined as an integral of the absolute value of the velocity correlation function *C*(*r*).

### Discrete particle simulations

Positions and orientations of the particles are described by the following Newton’s equation of motion [compare to models of Kokot *et al*. (*5*) and Belkin *et al*. (*41*)](10)(11)where *m* is the mass of the particle with radius *a*, is the magnetic moment directed along the unit vector, *I* is the moment of inertia, α_{t} = 6πη*a* = 0.18 is the translation and α_{r} = 0.25 is the rotation drag coefficient, α_{s} = 0.25 is the slipping parameter, is the angular velocity of the particle, is the vector normal to the surface at the particle position, and is the applied external magnetic field with amplitude *H*_{0} normalized by *a*^{3}/μ. To account for the fluctuations due to the shape imperfections of the particles, we included a Gaussian distributed noise term with zero mean and variances , where *I* is an identity tensor. and are the forces and torques due to the interactions, respectively. Here, we include three contributions, the magnetic dipole-dipole interaction(12)with the particle distance *r*_{ij} and steric interactions modeled by a short-range hard-core repulsion(13)

Furthermore, the particles are confined in a harmonic potential , with . The particles are assumed to roll down the potential well, which causes an additional torque on the particles. For the unit of length, we chose the radius of the particle *a*, energy is measured in μ^{2}/*a*^{3}, and time is measured in α_{t}*a*^{5}/μ^{2}. Hydrodynamic interactions between the particles are included in the far-field approximation (see the Supplementary Materials for details) (*42*).

The particles are initially placed on a square lattice with random orientations. We equilibrated the system during a time interval of 15,000 field periods and gathered statistics over an interval of at least 30,000 periods of an applied magnetic field (corresponding to 5 to 10 min of the experiment). These simulations were implemented on graphical processing units, and the typical simulation time was about 30 min.

### Implementation of hydrodynamic interactions

Hydrodynamic interactions between the particles were included in the far-field approximation by considering the stokeslet and rotlet of spherical particles. Therefore, we added the following hydrodynamic velocity, whose *k – th* component can be written as(14)where is the Oseen tensor in Eq. 10 and the hydrodynamic torque in Eq. 11(15)

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/2/e1601469/DC1

movie S1. Experimental movie for the gas-like state at frequency *f* = 20 Hz.

movie S2. Experimental movie for the flocking state at frequency *f* = 30 Hz.

movie S3. Animation of the experimentally obtained flocking state at frequency *f* = 30 Hz, where the direction of motion is indicated by the color code.

movie S4. Experimental movie for the vortex state at frequency *f* = 40 Hz.

movie S5. Animation of the experimentally obtained vortex state at frequency *f* = 40 Hz, where the direction of motion is indicated by the color code.

movie S6. Experimental movie for the reentrant flocking state at frequency *f* = 50 Hz.

movie S7. Animation of the experimentally obtained reentrant flocking state at frequency *f* = 50 Hz, where the direction of motion is indicated by the color code.

movie S8. Experimental movie for the gas-like state at frequency *f* = 60 Hz.

movie S9. Numerically obtained gas-like state at frequency *f* = 6 Hz, indicating the direction of motion by the color code.

movie S10. Numerically obtained flocking state at frequency *f* = 20 Hz, indicating the direction of motion by the color code.

movie S11. Numerically obtained vortex state at frequency *f* = 47 Hz, indicating the direction of motion by the color code.

movie S12. Numerically obtained reentrant flocking state at frequency *f* = 88 Hz, indicating the direction of motion by the color code.

movie S13. Numerically obtained vortex state, considering hydrodynamic interactions, at frequency *f* = 47 Hz, indicating the direction of motion by the color code.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, Division of Materials Science and Engineering. A.K. was supported through a Postdoctoral Research Fellowship (KA 4255/1-1) from the Deutsche Forschungsgemeinschaft (DFG).

**Author contributions:**A.S. and I.S.A. conceived the research, A.S. performed the experiment and analysis of the experimental data, and A.K. and I.S.A. conducted the numerical simulations. All authors wrote the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2017, The Authors