Degenerate states, emergent dynamics and fluid mixing by magnetic rotors

We investigate the collective motion of magnetic rotors suspended in a viscous fluid under an uniform rotating magnetic field. The rotors are positioned on a square lattice, and low Reynolds hydrodynamics is assumed. For a $3 \times 3$ array of magnets, we observe three characteristic dynamical patterns as the external field strength is varied: a synchronized pattern, an oscillating pattern, and a chessboard pattern. The relative stability of these depends on the competition between the energy due to the external magnetic field and the energy of the magnetic dipole-dipole interactions among the rotors. We argue that the chessboard pattern can be understood as an alternation in the stability of two degenerate states, characterized by striped and spin-ice configurations, as the applied magnetic field rotates. For larger arrays, we observe propagation of slip waves that are similar to metachronal waves. The rotor arrays have potential as microfluidic devices that can mix fluids and create vortices of different sizes.


I. INTRODUCTION
Magnetic forces provide a useful way of driving small units in viscous fluids. For example, in recent years magnetic driving has been extensively used in microfluidic devices, for applications such as magnetic swimmers [1][2][3][4], magnetic pumps [5,6] and cilia [7][8][9][10][11], and to facilitate particle sorting and segregation [12][13][14][15]. Many of the approaches are designed to control a single magnetic unit or to actuate many units in exactly the same manner. For further development of magnetically-actuated devices, it would be desirable to be able to control the dynamics of the magnetic units just by changing the external magnetic field. However, it is not immediately obvious how to achieve different behavior among different magnetic units when they are driven by the same external magnetic driving mechanism.
In a recent work [6], we analyzed the collective motion of a square array of magnetic rotors and found a surprisingly rich dynamics under an external magnetic field which oscillated along one axis of the array. By changing the relative strength of the external field and the dipolar interactions between the rotors, different collective rotational patterns emerged. When the dipole interaction was dominant the rotors swung back and forth, clockwise or counterclockwise in alternating stripes. When the external field dominated over the dipolar interactions, the rotors underwent full rotations with different quadrants of the array turning in different directions.
In the present paper, we extend these results to analyze the motion of the rotor array under a rotating magnetic field (see Fig. 1 for a schematic representation of our system). This enables us to identify three different collective rotational patterns as the balance of the torque due to the external field and the torque due to dipolar interactions is varied. In the first regime that occurs for strong fields, we observe synchronized patterns in which all the magnetic units rotate in the same direction as the external field. We observe oscillating patterns for weak field strengths, when the rotors cannot achieve net rotation. We also observe a most interesting chessboard pattern, which appears when the contributions of the two torques are comparable. In this pattern, the rotors rotate in a direction opposite to their closest neighbouring rotors. The rotors are positioned on a square grid in the xy-plane with spacing . Arrows in the rotors represent magnetic dipole moments. The magnets rotate about the z-axis driven by an external field B.
We examine system sizes up to 100 × 100 rotors, and report the observation of the propagation of novel slip waves that are similar to the metachronal waves that arise from coordinated motion of cilia [16][17][18][19][20][21][22]. We demonstrate that our system can be used as a microfluidic device for mixing and pumping, as it can drive flow fields with different vortex sizes and mixing length scales.

II. DESCRIPTION OF THE SYSTEM AND GOVERNING EQUATIONS
A. The system Our system consists of magnetic rotors positioned on a square lattice in xy-plane with a lattice spacing as shown in Fig. 1. The total number of rotors is N = N x N y , where N x and N y are the number of rotors in the x-and y-directions, respectively. The rotors are embedded in an unbounded fluid domain with viscosity η and density ρ. The magnetic dipoles have no translational degrees of freedom because they are fixed in space, but each has a single rotational degree of freedom about the z-axis. Every rotor is taken to be a sphere with radius a, which has a magnetic dipole moment where m is the magnitude and φ i is the polar angle characterizing the direction. The rotors are driven by a rotating magnetic field where B is the field strength, t is the time, and θ(t) = 2πf t is the polar angle of the external field that is rotating at a frequency f . We assume that hydrodynamic is in the viscous regime (low Reynolds number Re = a 2 ρf /η 1).

B. Governing equations
Each rotor experiences a magnetic torque that comprises contributions from the external field T ext , and the dipole-dipole interactions T dd . The torque contributions on the i-th rotor are where r i is the position of i-th rotor, r ij = r j − r i , r ij = |r ij |, n ij = (cos γ ij , sin γ ij , 0) = (r j − r i )/r ij , andê z is the unit vector along the z direction. Torque balance in the viscous regime gives the angular velocity ω i of a rotor as where T i = T ext i + T dd i and 8πηa 3 is the friction constant for the rotation of a sphere [23]. Note that we have ignored the hydrodynamic interaction between rotors by assuming that the rotor radius a is sufficiently small compared to the lattice distance a . The leading order effect of hydrodynamic coupling is discussed in Appendix A.
The potential energies due to the external magnetic field, U ext , and magnetic interaction, U dd , are The velocity field and the vorticity field at a position x are described by the rotlet [23] v (x) = 1 8πη When the voriticity observation point is in the z = 0 plane, it can be simply rewritten as

C. Dimensionless parameters
We introduce two dimensionless parameters: where τ r is the characteristic relaxation time of the magnetic dipoles. The Mason number Mn characterizes the relative strength of viscous stresses as compared to the magnetic ones [24], and τ r f characterizes the relaxation time of the system as compared to the period of the rotation of the external field. Equation (5) can then be rewritten in non-dimensional form: where r * ij = r ij / , T ext * i = T ext i /(mB), and T dd * . Similarly, we define the dimensionless potential energies as U ext * = U ext /(mB) and U dd * = U dd × 4π 3 /(µ 0 m 2 ).

D. Numerical method
We discretize Eq. (13) and follow the time evolution of the rotor angles using a first-order Euler method with a time step f ∆t = 10 −4 , and the initial orientation angles φ i are set randomly. Before the main simulation process starts, a strong external field corresponding to 1/(τ r f ) = 10 4 is applied in the +x-direction for a dimensionless time f t = 1 to align the rotors. When we analyze the average motion of the rotors or the flow field, we run the simulation for 50 cycles and take the average of the last 30 cycles.

A. Motion of a single rotor
For completeness, we first summarize the results for a single rotor.In the absence of dipole interaction, the rotational velocity is determined purely by the torque from the external field T ext . The rotational velocity is therefore simply and the time evolution of the phase difference ∆ = θ − φ is∆ For 1/(τ r f ) > 2π the stable solution has∆ = 0 and the rotor follows the external field with a constant angular velocity ω = 2πf . For 1/(τ r f ) < 2π, however, there is no solution with∆ = 0: the rotor does not reach a constant velocity but periodically slips to move backwards as it fails to keep up with the external field [2,11,25]. Figure 2(a) is a phase diagram showing the dynamical configurations of 3 × 3 arrays of magnetic rotors (see also Movie 1). The competition between the two torques T ext and T dd , which are controlled by the two dimensionless parameters 1/(τ r f ) and 1/Mn, respectively, determines the rotational patterns. We identify three different responses to the field: "synchronized" (Fig. 2(b)), "chessboard" (Fig. 2(c)), and "oscillating" (Fig. 2(d)) patterns.
Note that we introduce a parameter R to characterize the rotational pattern [6], defined as where ∆φ i and ∆θ are the total rotational angles of the rotors and the external field, respectively.
The synchronized pattern appears when the torque due to the external field is dominant over the dipolar interactions 1/(τ r f ) > 1/Mn, which corresponds to the top left corner of Fig. 2(a). In this regime, all rotors rotate in the direction of the external field as shown in Fig. 2 The oscillating pattern appears when the dipolar interaction is dominant 1/Mn 1/(τ r f ), or when the rotors cannot catch up with the external field 1/(τ r f ) < 2π as also seen in the motion of a single rotor. The rotors simply oscillate and cannot achieve a net rotation because the external field is not strong enough to destroy the configurations favoured by the magnetic interaction, such as the spin-ice pattern [6,26,27].
The chessboard pattern appears when the contributions of the two torques are comparable, namely, 1/Mn ∼ 1/(τ r f ). In this regime, all rotors rotate in a direction opposite to their closest neighbours. This is counterintuitive because half of the rotors are rotating in the opposite direction to the magnetic field. There have been reports of such chessboard-like behaviour in other active matter systems, such as bacterial suspensions [28][29][30], although they do not operate under the effect of an external driving mechanism unlike our system.

The mechanism behind the chessboard pattern
In order to further understand the rotational patterns, we next apply a static magnetic field, B(θ 0 ) = (B cos θ 0 , B sin θ 0 , 0) where 0 ≤ θ 0 < 2π, and analyze the equilibrium orientation configuration for each angle θ 0 . Since the rotors are in the equilibrium state, there is only one relevant dimensionless parameter in the system, namely which compares the strength of the external field to the typical dipole field [6]. When there is no applied external field (i.e. α = 0) there exist several equilibrium configurations for a square rotor array such as the stripe and spin-ice configurations shown in Fig. 3(a). In a 3 × 3 array, the magnetic dipoledipole interaction energy U dd * is in fact the same in these two configurations for regular alignment of the magnetic moments, i.e., φ i = nπ/4 with n as integers. Note, however, that the directions of the magnetic dipoles in finite systems can have minor deviations from the regular alignment defined by φ i = nπ/4, which corresponds to infinite systems. This leads to the stripe configuration having slightly lower dipolar potential energy (U dd * ≈ −16.51) than the spin-ice configuration (U dd * ≈ −16.45) in the 3 × 3 array. Starting from random rotor orientations without an external field, namely α = 0, the simulation results in one of these equilibrium configurations depend on the details of the initial angles. If we now impose a static external magnetic field that is comparable to the dipolar field, namely α ∼ 1, the system changes the pattern to that dictated by the external field direction θ. As shown in the chessboard A pattern in Fig. 3(b), the system prefers stripe-like configurations for external field directions θ 0 = nπ/2 where n is an integer, while it prefers spin-ice-like configurations for θ 0 = (2n + 1)π/4 (see Appendix B). Since the chessboard rotational pattern can be obtained by following the sequence of 8 configurations in Fig. 3(b), we can understand that it is a consequence of alternate switching between stripe and spin-ice ordering.
Note that we classified chessboard patterns into two different patterns, chessboard A (that corresponds to smaller α) and chessboard B (that corresponds to larger α), depending on the transient configurations in Fig. 3(b). In the chessboard B pattern, configurations at θ 0 = nπ/2 are not stripe-like patterns (such as chessboard A), as shown in Movie 2.
The chessboard pattern can be seen only for α ∼ 1. Figure 3(d) shows the potential of dipolar interactions U dd * . When the external field is weak (α 1), as for α = 0.05 in Fig. 3(d), the system cannot transit from the stripe pattern to the spin-ice pattern because the external field strength is not sufficiently strong for the rotors to jump across the energy barrier ∆U dd * between the two configurations. As the result, they stay on a single stripe pattern and just result in the oscillating pattern. When the external field is moderately strong, say for α = 0.2, the system can make the transition into all 8 patterns as shown in Fig. 3 (b) and (d). When the external field is strong (α 1) on the other hand, the chessboard pattern cannot be seen because the rotors all align to the external field direction and the stripe or the vortex configurations disappear.

Predicting the phase boundaries
Considering a simpler problem allows us to obtain an estimate for the boundaries between the different patterns. Figure 3(c) shows the rotational patterns for different values of α. Under a negligible viscous friction, namely 1/Mn 1 and 1/(τ r f ) 1, we expect each rotor angle φ i to follow the local magnetic field instantaneously. The simplified system shows all three patterns as α is varied. The numerical value for the threshold between the oscillating and the chessboard pattern is α 1 = 0.05 − 0.06 while the corresponding value for the transition between the chessboard and synchronized patterns is α 2 = 1.  the three rotational patterns can be estimated using 1 The second term of right hand side, 2π, appears because there is a threshold for slipping, namely 1/(τ r f ) = 2π, as shown in the single rotor analysis 1/Mn 1. The lines in Fig. 2(a) correspond to Eq. (18) with the parameters α 1 and α 2 . The predictions agree well with the simulation results.
Analytical estimates for the values of α 1 and α 2 can be obtained by considering the energies of the different patterns. For simplicity, we consider only rotor orientations that satisfy φ i = nπ/4 where n is an integer, thus approximating the stripe and spin-ice configurations as shown in Fig. 3(a). The corresponding dipole interaction energy U dd * is also shown in the figure. For a given orientation patterns φ i , the total dimensionless energy can be estimated as An estimate of α 1 is obtained by comparing the total energy of the stripe and spin-ice configurations for an external field direction θ 0 = π/4. This is based on the observation that the chessboard pattern appears when the total energy of spin-ice configurations is lower than that of stripe configurations. Using U ext * (φ stripe i , π/4) = −3 √ 2/2 and U ext * (φ ice i , π/4) = −3, and assuming that the energy varies monotonically as a function of α, we obtain α 1 = 0.068, which is in good agreement with the simulation result.
For α 2 , the estimate follows from comparing the total energy of the striped and aligned (Fig. 3(a)) configurations, for an external field at θ 0 = 0 noting that the synchronized pattern appears when the total energy of aligned configurations is lower than that of the striped ones. Using U ext * (φ stripe i , 0) = −3 and U ext * (φ aligned i , 0) = −9, we obtain α 2 = 1.36, which is again in good agreement with the numerical value obtained from our simulations.

Velocity field
The collective rotational patterns can be used to produce complex rotational flow fields [6]. Figure 4 shows time-averaged streamlines and vorticity patterns, evaluated by using equations (8) and (10), respectively. For the synchronized pattern, the flow field is a single large vortex as shown in Fig. 4(a). Both the rotors and the large vortex rotate in the counterclockwise direction. Note that the local vorticity direction is opposite to that of the torque, as shown in Eq. (10). For the chessboard pattern, each rotor creates a smaller, individual vortex as shown in Fig. 4(b). There is no flow between the rotors because the flow field generated by neighbours cancels out. In the oscillating pattern there is no time-averaged flow field because the net rotation angle of each rotor is zero.
It is well known that mixing at small Reynolds numbers is challenging [31][32][33], because of the absence of inertia. This device would allow the vortex size and the mixing length scale to be controlled by the external magnetic field. Other types of flow fields are also possible for oscillatory magnetic driving mechanisms [6].
C. The effect of lattice size Figure 5 shows the phase diagrams for different sizes of the rotor arrays. The value of R as defined in Eq. (16) is always in the range −1 ≤ R ≤ 1, and the color in Fig. 5 represents the rotational directions: red (R = 1) corresponds to rotation along with the external field, blue (R = −1) shows the opposite rotation and gray (R = 0) shows rotors with no net rotation.
Although the basic structure of the phase diagram does not depend on the array size, the rotational patterns are different around the regions where the two torque contributions are comparable, namely 1/(τ r f ) ∼ 1/Mn. The arrays that have odd values of N x and N y show the same chessboard pattern as the 3 × 3 array (see also Movie 3). By comparing arrays of 3 × 3, 5 × 5, and 7 × 7 (data not shown), we find that the parameter range over which the chessboard pattern is stable becomes smaller with increasing system size. For 4 × 4 arrays, there are equal numbers of clockwise and counterclockwise rotors for 1/(τ r f ) ∼ 1/Mn. However, the pattern differs from the chessboard pattern: each rotor takes three periods to complete a single full rotation. Moreover, a 2 × 2 array shows only synchronized or oscillating patterns.

Metachronal waves in a large rotor array
To examine the effect of the system size even further, we consider the case of N x = N y = 100. Figure 6 shows the rotational pattern under a slip condition 1/(τ r f ) = 5 < 2π and 1/Mn = 1. Although the magnets rotate in the same direction as the external field (ω > 0; red) most of the time, they sometimes slip (ω < 0; blue) when they cannot keep up with the field. Interestingly, the slipping motion exhibits a wave-like collective pattern in the large system. We observe that the direction of wave propagation depends on whether or not we incorporate hydrodynamic interactions.
When interactions mediated by the fluid are negligible a * = a/ 1, the slipping motion is first triggered at the outer edges of the system and slowly propagates inwards, as shown in Fig. 6(a) and Movie 4. This is because the alignment of the outer layers is strongly distorted due to the open boundary. If we now include the first-order hydrodynamic coupling between the rotors (see Appendix A) and increase the effective density by using the parameter a * = 0.45, the direction of the wave propagation reverses and the rotors first slip at the inner layers as shown in Fig. 6(b) and Movie 5. This occurs because the hydrodynamic effects are stronger at the centre of the array. As shown by Eq. (10), each rotor creates a flow field with vorticity that is opposite to the rotating direction of the rotor. Therefore the hydrodynamic interactions tend to slow down neighbouring rotors. Since the inner layers have a larger number of neighbours, the slip starts from this region.

IV. CONCLUSION
In this paper, we describe the collective motion of magnetic rotors under a rotational magnetic field. We mainly focus on a 3×3 array of magnets, and observe three characteristic dynamical rotor patterns as the external field strength is varied. In the synchronized regime, which ap- pears when the magnetic torque due to the external field is dominant, all magnets rotate in the same direction as the field. If, however, the dipolar interaction between the rotors is dominant, all rotors exhibit oscillations (rather than full rotations) during a cycle. When the contributions from the external field and the dipolar interaction are comparable, there is an unexpected chessboard pattern of rotations where the magnets rotate in opposite directions to their closest neighbours. We argue that the chessboard pattern appears as a consequence of the rotor array switching between stripe and spin-ice configurations: essentially the field lifts the degeneracy of these two states. For large system sizes, we observe the propagation of slip waves reminiscent of the metachronal waves observed in natural systems [16][17][18][19][20][21][22].
Such rotor arrays have potential as microfluidic devices that can create various flow fields and vortices of different sizes. Depending on the external field condition, the flow around the device considered here can vary from a single large vortex, to small alternating vortices around each rotor. We believe that our work sheds light on new collective dynamics of magnetic rotors and opens possibilities for future experiments and applications.
Varying the deviation ∆φ from 0 to 2π, we find that the magnetic dipole-dipole interaction energy U dd is constant for odd N = N x = N y , not only for the stripe and spin-ice configurations, but also for all other angles. On the other hand, the potential U dd for even arrays is not constant. The potentials due to the external magnetic field, U ext , are U ext * (φ p,q (∆φ), θ 0 ) = −N cos(θ 0 − ∆φ) (odd N ) 0 (even N ) (B2) and arrays with even N show a constant energy level. As the result, U ext is dominant when N is odd, and the orientational pattern prefers ∆φ = θ 0 . For arrays with even N , there is no preferred configuration because U ext is constant.