Cite this Article
Matsunaga Daiki, Chodo Takumi, Kawai Takuma. Diffusion and collective motion of rotlets in 2D space. Chinese Physics B, 2020, 29(6): 064705  
Permissions
Diffusion and collective motion of rotlets in 2D space
Matsunaga Daiki †, Chodo Takumi, Kawai Takuma
Graduate School of Engineering Science, Osaka University, 5608531 Osaka, Japan

 

† Corresponding author. E-mail: daiki.matsunaga@me.es.osaka-u.ac.jp

Project supported by JST, ACT-T Grant No. JPMJAX190S Japan and Multidisciplinary Research Laboratory System for Future Developments (MIRAI LAB)

Abstract

We investigate the collective motion of rotlets that are placed in a single plane. Due to the hydrodynamic interactions, the particles move through the two-dimensional (2D) plane and we analyze these diffusive motions. By analyzing the scaling of the values, we predict that the diffusion coefficient scales with ϕ0.5, the average velocity with ϕ, and relaxation time of the velocity autocorrelation function with ϕ–1.5, where ϕ is the area fraction of the particles. In this paper, we find that the predicted scaling could be seen only when the initial particle position is homogeneous. The particle collective motions are different by starting the simulation from random initial positions, and the diffusion coefficient is the largest at a minimum volume fraction of our parameter range, ϕ = 0.05. The deviations based on two initial positions can be explained by the frequency of the collision events. The particles collide during their movements and the inter-particle distances gradually increase. When the area fraction is large, the particles will result in relatively homogeneous configurations regardless of the initial positions because of many collision events. When the area fraction is small (ϕ < 0.25), on the other hand, two initial positions would fall into different local solutions because the rare collision events would not modify the inter-particle distances drastically. By starting from the homogeneous initial positions, the particles show the maximum diffusion coefficient at ϕ ≈ 0.20. The diffusion coefficient starts to decrease from this area fraction because the particles start to collide and hinder each other from a critical fraction ∼ 23 %. We believe our current work contributes to a basic understanding of the collective motion of rotating units.

PACS: ;47.85.-g;;83.10.-y;
Keyword:;rotlet;;diffusion;;Stokes flow;
1. Introduction

The collective motion of the rotating unit can be seen in the artificial or natural systems, and a number of interesting phenomena have been reported. By applying a rotating external magnetic to magnetic colloids, the particles rotate and show complex collective motions due to the hydrodynamic and magnetic interactions between the colloids. Previous works applied in-plane[1–4] or out-of-plane[5–8] rotational magnetic field to magnetic colloids that locate in a single plane and analyse their collective motions. Recently, we reported that the magnetic rotors that are placed in a square grid show interesting rotational patterns.[9,10] We found a quarter pattern, in which the rotors undergo full rotations with different quadrants of the array turning in different directions when the external magnetic field is dominant. When the dipolar interaction between the rotors is dominant, a stripe pattern in which the rotors swing upwards or downwards in alternating stripes emerges. In recent years, we also proposed that this magnetic torque is useful to control a single magnetic object.[11–13]

In nature, a group of bacteria that locates close to a rigid wall shows a crystal pattern[14] or clustering phenomenon.[15,16] When the bacteria are swimming toward the bottom wall, they have attractive interactions between each other due to the hydrodynamic interactions of rotating doublets. Previous studies[16,17] reported that the clustering of magnetotactic bacteria can be controlled by changing the external magnetic field. There is also a report that rotating units, such as bacteria flagella, that are fixed to a wall would create a metachronal beating pattern.[18]

When there are a number of rotating units under small but finite Reynolds number, they form a cluster and rotate together as a group.[19,20] There are also interesting reports that a binary mixture, which is two groups of rotating units that have opposite rotational directions, shows a phase separation.[21–23] All these physics can be categorized to the active matter physics, which is a research area that is growing in recent years.

In these previous studies, they reported the collective motions of rotating units that have not only hydrodynamic interactions but also the magnetic interactions, inertial force, binary mixtures, or other factors. In this paper, we set a goal to understand the basic character of collective motions of rotating units, under a simple problem setup. We assume Stokes flow in the present work and investigate the collective motion of rotlets (torques), which is a flow field that is driven by a point torque,[24] that are placed in a single 2D plane with changing the area fraction of the particles ϕ. When a particle is under the two extremes conditions, a dilute limit ϕ ≪ 1 or maximum packing fraction, the particle has no movement because there is no flow that translates the particle. Therefore, the particle would obtain the maximum velocity and diffusion coefficient at an intermediate area fraction, and we discuss this point in this paper.

2. Method

Consider N particles with the radius a that are placed in a single monolayer xy-plane, with an area fraction ϕ as shown in Fig. 1. The particles are suspended in an infinite fluid domain with the viscosity η and density ρ. External torque T = (0, 0, T) is applied to all particles, and the particles move collectively due to the hydrodynamic interactions. If the Reynolds number satisfies Re = a2 ρ ω/η ≪ 1, where ω is the particle angular velocity, the system can be regarded as in a Stokes flow regime. There are a number of previous researches that investigate the collective motion of magnetic rotors placed in a monolayer.[1–8] The current problem setup can be regarded as a case when the magnetic dipolar interaction is negligibly small compared to the external magnetic torque.

Fig. 1. Schematic of the problem setup.

In an ideal hexagonal lattice, the centre-to-centre distance ℓ can be described as a function of ϕ as

In this paper, we consider this length as the average inter-particle distance under homogeneous distributions. Note that the maximum area fraction of spheres in 2D is .

2.1. Governing equations

In this paper, we use far-field hydrodynamic interactions of the Stokesian dynamics.[25,26] In a dilute regime ϕ ≪ 1, the i-th particle would rotate about z-axis with an angular velocity

where 8 π η a3 is the rotational friction coefficient for a sphere.[24] Because of the far-field hydrodynamic coupling, the angular velocity of the i-th particle, ωi, is altered because of the other particles’ existence as

where rij = xj – xi is the particle relative position vector, rij = |rij|, and xi is the position of the i-th particle. Since all particles are in the same plane, the equations reduce to

where ωi = (0, 0, ωi) and Ti = (0, 0, T). The rotlets T also contributes to the other particle translational velocity v[24,27] as

Once the particles are overlapped, we assume that there are repulsive steric interactions defined as

where k is the spring constant. Due to the repulsive force, the particle moves with a velocity

where

is the Rotne–Prager–Yamakawa tensor.[28,29] Now, the total velocity of the particle is .

Note that this model stands only for the dilute to semi-dilute regime, and it should be revised when the particles are in the touching distance. In this paper, we keep this far-field assumption in order to understand the basic character of the collective motions of rotlets.

2.2. Dimensionless form

By defining as the characteristic time, the equations are rewritten in the dimensionless form as follows:

Note that the superscript * denotes the dimensionless values.

2.3. Numerical method

There is only one dimensionless parameter ka2/T in this setup, which is the strength of the repulsive force compared to the rotlet interactions. In this paper, we kept this parameter as a constant and set it small ka2/T = 0.1 since our main focus is on the rotlet contribution. The number of particles N = 2000 is kept constant and the domain size L × L is varied to control the volume fraction ϕ, where . The periodic boundary condition is applied to x- and y-directions, and we only take into account the hydrodynamic interactions between particles that have distance rij smaller than L/2 to avoid the effect of the square periodic boundaries. The 4th-order Runge–Kutta method with a time step is used for the time integration, and graphics processing unit (GPU) is used to speed up the simulation. Note that we confirmed that the diffusion coefficient of the particles are nearly the same (errors: 3.42 % (ϕ = 0.05), 1.87 % (ϕ = 0.20)) even if we set the time step 5 times smaller, .

In this paper, we utilize two different initial positions of the particles: random or homogeneous configurations. The particle positions are initially random for the random configurations, while the particles are positioned in a square grid for the homogeneous configurations.

3. Results

Figure 2 shows snapshots from the final frame (dimensionless time ; random initial positions) of the simulation for different area fractions ϕ. There is no apparent aggregation or cluster in the system and the distributions are rather uniform. When the area fraction is sufficiently large ϕ ≥ 0.40 as in Fig. 2(d), the particles form a hexagonal crystal structure and would not change their configurations by time. The collective motions of the particles are also shown in movies 1, 2, 3, and 4.

Fig. 2. Snapshot of the simulation (random initial positions) for different area fractions ϕ at dimensionless time . Note that the image shows only the center 44a × 44a region of the simulation domain. The black lines on the particle indicate the instantaneous velocity direction.

Figure 3 shows the trajectories of randomly-selected eight particles under a dimensionless time . Note that the initial positions of all particles are set at (x0, y0) = (0, 0) for easier understanding. By comparing the trajectories, the persistent length of the trajectories is different by the area fraction ϕ. When the area fraction is small, the particles travel rather straight and longer distances, while the particles exhibit wiggling trajectories and travel shorter for larger area fractions. Among these four cases, the particles travel the longest distance at the smallest area fraction ϕ = 0.05. In order to characterise the trajectories, we show the curvature radius R/a of the particle trajectory in Fig. 4(a). The curvature radius R/a at time t is evaluated by computing a circle radius that goes through three points: x(t – Δ t), x(t), and x (t + Δ t), where for the current analysis. The probability distribution function (PDF) shows that the trajectories of larger area fractions give smaller R/a. The PDF curve shows a maximum at small R, and it quickly decays with the radius. The inset figure shows the ensemble average of the radius 〈 R 〉/a, and the value is larger for lower area fraction as expected.

Fig. 3. Particle trajectories of eight particles for dimensionless time (initial state: random). The starting point of the particle is set at the origin (x0, y0) = (0, 0) for easier understanding.}

Many previous papers reported that the rotating units, such as rotating disks[19] or bacteria swimming close to a surface,[14] form hexagonal lattice structures. In these previous studies, they form a crystal-like structure due to the existences of the attractive and repulsive forces. In the case of rotating disks under finite Reynolds number,[19] the Magnus effect contributes to repulsive interactions while hydrodynamic interactions contribute to attractive interactions. In the present study, the steric repulsive force F drives the crystallization for ϕ ≥ 0.40. Figure 4(b) shows a particle trajectory that is surrounded by six rotating particles that are fixed in space. Note that we set the inter-particle distance of the six outer particles at ℓ/a (ϕ = 0.50) = 2.69. Due to the collision events, the inter-particle distance gradually increases until there is no collision, rij/a = 2. In other words, the particle would stay in the same circular orbit without being pushed towards the center if there is no steric force. The repulsive force alters the particle trajectories, and as a result, the inter-particle distance gradually increases. Since there are many collision events for large area fractions, the particles would form the crystal regardless of the initial positions.

Fig. 4. (a) The probability distribution function of the curvature radius R/a of the particle trajectories. The inset shows the ensemble average of the radius as a function of the area fraction ϕ. (b) The mechanism of the increase in the inter-particle distances. The figure shows a particle trajectory that is surrounded by six rotating particles that are fixed in space. Due to the collision events, the inter-particle distance gradually increases and results in a distance rij/a = 2.
3.1. Prediction of scalings

Before further reporting the numerical results, we show here estimations of the scalings of three parameters in this subsection: average velocity 〈 v 〉 = 〈 |v| 〉, diffusion constant D, and relaxation time of the velocity auto-correlation function τ0. Note that we ignore the effect of the steric force F in this estimation.

Coming back to the equation of translational velocity that is generated by the rotlets vT, Eq. (5), the velocity is only a function of the distance rij in our system (note that the torque is a constant in this work). In other words, varying the area fraction ϕ only changes the velocity magnitude, and the particle trajectories would be the same for different ϕ if we ignore the steric effect vF.

We now consider a new dimensionless system: fixing the distance ℓ and varying the area fraction ϕ to change the particle size a(ϕ). If we rewrite the velocity (5) as

our dimensionless velocity would scale with ϕ using , because should be a constant for any particle radius a (ϕ). This scaling matches the prediction in a previous study.[30]

In the same matter, scaling of other values can be evaluated as

From the next sections, we find that these scalings work only when we start the simulation from the homogeneous initial positions. Note again that these estimations stand only for small area fractions because we ignore the steric effects.

3.2. Diffusion coefficient

Figure 5 shows the mean square displacement (MSD) and diffusion coefficient, which are respectively given by

The MSD curves show the ballistic movement (slope 2) for a short time and the diffusive motion (slope 1) for a long time scale. Note that we evaluate the variables in this paper after a long time () from the initial state.

Fig. 5. (a) Mean square displacement 〈 r2 (τ) 〉/a2 as a function of elapsed time τ. Note the simulations started from random initial positions. (b) Diffusion coefficient D as a function of the area fraction. The inset shows the same result in a log–log scale.

Figure 5(b) shows the diffusion coefficients for two different initial positions, random and homogeneous (aligned) configurations, and we find that the diffusion coefficients drastically differ by the initial position especially for ϕ < 0.20. By starting the simulations from random initial positions, the diffusion is larger for lower area fractions. Note we run the simulations with three different random positions for ϕ = 0.05–0.20 and show the error bars in the figure. Although we could not find the area fraction that shows the maximum diffusion because of the computational time limitation, the maximum diffusion coefficient might be located at ϕ < 0.05.

By starting from the homogeneous positions, on the other hand, the diffusion coefficient scales with ϕ0.5 for small ϕ as we predicted, and the coefficient reaches a maximum around ϕ ∼ 0.20. The deviation from the scaling ϕ0.5 at ϕ ∼ 0.20 can be explained by the steric effect. When a particle would go through space between two particles i and j, the center-to-center distance between two particles needs to be larger than rij/a = 4. Solving ℓ/a = 4 gives the critical area fraction ϕc ≈ 0.23, and the particles hinder each other if the area fraction is larger than this value. Note we evaluate the diffusion coefficient with other homogeneous initial position (hexagonal lattice with a small spacial inhomogeneity) just for a single case (ϕ = 0.15), and find that the diffusion is nearly the same ().

3.3. Radial distribution function

In order to investigate the difference between two initial configurations, we show the radial distribution functions in Fig. 6. Note the value n (r) is the particle number density that is in the distance r from a particle, and n0 = N/L2 is the average number density of the whole domain. By starting from homogeneous initial positions, the particles have spacial patterns for all area fractions. The distribution reaches a maximum at shorter distances for larger ϕ, and the first peaks are nearly the same as ℓ/a (Eq. (1)), which are shown with dash-dot lines.

Fig. 6. Radial distribution function of the particles. Note the value n (r) is the particle number density at a distance r, and n0 = N/L2 is the average number density of the whole domain. Dash-dot lines show the average distances between particles at the homogeneous state, Eq. (1).

These distributions are drastically different for small area fractions ϕ < 0.20 by starting the simulation from random initial positions. Compared to the homogeneous initial positions, firstly, the particles locate closer to each other. At ϕ = 0.05 (Fig. 6) for example, the minimum distance of the non-zero distribution is r/a ∼ 4 for the homogeneous position while it is r/a ∼ 2 for the random position. Secondly, there are less dominant spacial structures in the distribution when starting from the random positions. These differences explain the deviations in the diffusion coefficient. Note we confirmed that the distribution of random initial positions does not change even we simulate ten times longer, until dimensionless time .

Contrary to the case of small area fractions, the distributions are the same for large area fractions. This contrast can be explained by the frequency of collision events. As we discussed in the previous sections, the inter-particle distance gradually increases due to the collision events. The radial distribution function is the same regardless of the initial positions for large area fractions since many collision events would quickly modify the particle configurations. For small area fractions, on the other hand, the radial distribution function at small distances would remain from the initial random positions, because the rare collision events would not modify the inter-particle distances drastically. As a result, the radial distribution function of the random initial positions fall into a different local solution from the homogeneous one.

3.4. Velocity

The average velocity of the particle 〈v〉 also depends on the initial positions, as shown in Fig. 7(a). By starting from the random initial positions, the average velocity is nearly constant. By starting from the homogeneous positions, on the other hand, the velocity grows linearly with ϕ for small volume fractions, as we predicted in the previous section.

Fig. 7. (a) Average velocity 〈 v 〉 = 〈 |v| 〉 as a function of the area fraction. (b) The velocity auto-correlation function as a function of the elapsed time τ (initial state: random). The inset shows the relaxation time of the velocity auto-correlation function τ0 as a function of the area fraction.

If we assume that the translational velocity of a particle is mainly determined by a single neighboring particle at a distance ℓ/a, the translational velocity is evaluated by substituting Eq. (1) to Eq. (5) as

This prediction, , is also shown in Fig. 7(a), but the estimation does not match the simulation result well. We obtain by fitting the velocity at the range ϕ ≤ 0.20, and the prediction might be failed because we did not take in to account many particle interactions. The average velocity starts to deviate from the scaling ϕ at ϕ = 0.25 again because of the steric effect. Since the particles collide and hinder each other beyond the critical area fraction ϕc ≈ 0.23, the velocity saturates and would not grow with the area fraction.

Figure 7(b) shows the velocity auto-correlation function, which is given by

for random initial positions. The function monotonously decreases with time for small area fractions, while it fluctuates for larger area fractions. The inset figure shows the decay time scale of the function, which is obtained by fitting the function with an equation exp (–t/τ0). The figure shows that the relaxation time τ0 scales with ϕ–1.5 at small area fractions, as we predicted.

4. Conclusion

In this paper, we investigate the collective motion of far-field rotlets that are placed in a single plane. Due to the hydrodynamic interactions, the particles move through the 2D plane and we analyze these diffusive motions. By analyzing the scaling of the values, we predict that the diffusion coefficient scales with ϕ0.5, the average velocity with ϕ, and relaxation time of the velocity autocorrelation function with ϕ–1.5, where ϕ is the area fraction of the particles. In this paper, we find that the predicted scaling could be seen only when the initial particle position is homogeneous. The particle collective motions are different by starting the simulation from random initial positions, and the diffusion coefficient is the largest at a minimum volume fraction of our parameter range, ϕ = 0.05.

The deviations based on two initial positions can be explained by the frequency of the collision events. The particles collide during their movements and the particle–particle distance gradually increases. When the area fraction is large, the particles will result in relatively homogeneous configurations regardless of the initial positions because of many collision events. When the area fraction is small (ϕ < 0.25), on the other hand, two initial positions would fall into different local solutions because the rare collision events would not modify the inter-particle distances drastically.

By starting from the homogeneous initial positions, the particles show the maximum diffusion coefficient at ϕ ≈ 0.20. The diffusion coefficient starts to decrease from this area fraction because the particles start to collide and hinder each other from a critical fraction ∼ 23 %.

We believe our current work contributes to a basic understanding of the collective motion of rotating units.

Reference
[1] Tierno P Muruganathan R Fischer T M 2007 Phys. Rev. Lett. 98 028301
[2] Coughlan A C H Bevan M A 2016 Phys. Rev. 94 042613
[3] Pham A T Zhuang Y Detwiler P Socolar J E S Charbonneau P Yellen B B 2017 Phys. Rev. 95 052607
[4] Soni V Bililign E S Magkiriadou S Sacanna S Bartolo D Shelley M J Irvine W T 2019 Nat. Phys. 15 1188
[5] Driscoll M Delmotte B Youssef M Sacanna S Donev A Chaikin P 2017 Nat. Phys. 13 375
[6] Kaiser A Snezhko A Aranson I S 2017 Sci. Adv. 3 e1601469
[7] Kokot G Snezhko A 2018 Nat. Commun. 9 2344
[8] Massana-Cid H Meng F Matsunaga D Golestanian R Tierno P 2019 Nat. Comm. 10 2444
[9] Matsunaga D Hamilton J K Meng F Bukin N Martin E L Ogrin F Y Yeomans J M Golestanian R 2019 Nat. Comm. 10 4696
[10] Kawai T Matsunaga D Meng F Yeomans J M Golestanian R 2020 arXiv: 2003.05082
[11] Matsunaga D Meng F Zöttl A Golestanian R Yeomans J M 2017 Phys. Rev. Lett. 119 198002
[12] Matsunaga D Zöttl A Meng F Golestanian R Yeomans J M 2018 IMA J. Appl. Math. 83 767
[13] Meng F Matsunaga D Yeomans J M Golestanian R 2019 Soft Matter 15 3864
[14] Petroff A P Wu X L Libchaber A 2015 Phys. Rev. Lett. 114 158102
[15] Chen X Yang X Yang M Zhang H 2015 Europhys. Lett. 111 54002
[16] Pierce C Wijesinghe H Mumper E Lower B Lower S Sooryakumar R 2018 Phys. Rev Lett. 121 188001
[17] Meng F Matsunaga D Golestanian R 2018 Phys. Rev. Lett. 120 188101
[18] Uchida N Golestanian R 2010 Phys. Rev. Lett. 104 178103
[19] Goto Y Tanaka H 2015 Nat. Commun. 6 5994
[20] Shen Z Lintuvuori J S 2020 Phys. Rev. Res. 2 013358
[21] Nguyen N H Klotsa D Engel M Glotzer S C 2014 Phys. Rev. Lett. 112 075701
[22] Yeo K Lushi E Vlahovska P M 2015 Phys. Rev. Lett. 114 188301
[23] Ai B Q Shao Z Q Zhong W R 2018 Soft Matter 14 4388
[24] Kim S Karrila J S 1991 Microhydrodynamics - Principles and Selected Applications New York Dover Publications, Inc. 10.1016/C2013-0-04644-0
[25] Durlofsky L Brady J F Bossis G 1987 J. Fluid Mech. 180 21
[26] Brady J F Bossis G 1988 Annu. Rev. Fluid Mech. 20 111
[27] Lushi E Vlahovska P M 2015 J. Nonlinear Sci. 25 1111
[28] Rotne J Prager S 1969 J. Chem. Phys. 50 4831
[29] Yamakawa H 1970 J. Chem. Phys. 53 436
[30] Llopis I Pagonabarraga I 2008 Eur. Phys. J. 26 103