Dynamic Clustering and Scaling Behavior of Active Particles under Confinement

A systematic investigation of the dynamic clustering behavior of active particles under confinement, including the effects of both particle density and active driving force, is presented based on a hybrid coarse-grained molecular dynamics simulation. First, a series of scaling laws are derived with power relationships for the dynamic clustering time as a function of both particle density and active driving force. Notably, the average number of clusters N¯ assembled from active particles in the simulation system exhibits a scaling relationship with clustering time t described by N¯∝t−m. Simultaneously, the scaling behavior of the average cluster size S¯ is characterized by S¯∝tm. Our findings reveal the presence of up to four distinct dynamic regions concerning clustering over time, with transitions contingent upon the particle density within the system. Furthermore, as the active driving force increases, the aggregation behavior also accelerates, while an increase in density of active particles induces alterations in the dynamic procession of the system.


Introduction
Active particles, distinguished by their autonomous operations at the nano-and microscales within fluids, can be differentiated from particles which passively react to their environment.Their unique ability to self-propel and impart energy to a system has recently garnered considerable attention.The applications of active particles span a diverse array, encompassing catalysis [1], toxin detection [2,3], cellular marking [4], wastewater treatment [5][6][7][8][9][10][11][12], CO 2 scrubbing [13], chemical and biological warfare agent neutralization [14][15][16][17], on-the-fly hydrogel polymerization [18], micropatterning [19], and even energy generation [20].Indeed, the lion's share of interest in active particles has been concentrated in the realm of medicine and chemistry, particularly in drug delivery and synthesis [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36].Magnetically controlled active particles, such as those made of iron oxide, have shown great promise in drug delivery as well as in the removal of blockages from blood vessels within the body [37][38][39].In applications where active particles are abundantly present, such as in microfluidics devices or the human body, systems can be considered confined fluid systems [40].However, depending on factors like driving force and particle density, active particles within the body may exhibit tendencies to aggregate and cluster together, which could raise potential risks yet to be investigated in biomedical applications [21,41].One notable feature of active particles is their enhanced diffusion within a fluid system compared to similar passive particles.Moreover, an increase in active particle density corelates with a greater effective diffusion coefficient [42].Active colloids operate by dissipating energy at the level of the active particles in order to cause motion with no equilibrium equivalent, resulting in clustering under low dissipation and collective motion under high dissipation [43].The coexistence of active and passive colloidal particles has been shown to prompt clustering and phase separation, forming chains or clusters [44].The collective behavior of active particles introduces numerous effects not observed in the behavior of individual particles, surpassing the simple summation of individual behaviors due to interactions arising between active particles [45][46][47][48][49][50].The tendency of colloids to aggregate can likewise demonstrate complexity when considering the autonomous motion of active particles.
Recent research has demonstrated that the directed assembly of colloidal dispersions can be used to synthesize a diverse range of functional materials in a controlled fashion [51][52][53].While self-assembly is inherent to passive colloids, systems composed of active particles exhibit an augmented capacity for self-assembly, due to the driven nature of these particles [54,55].This type of active aggregation has been shown to create dynamic clusters that evolve over time, propelled by both fluid velocity and the individual motion of active particles.This mechanism effectively overcomes the long-range repulsion commonly observed in colloidal systems [56,57].The aggregation of active colloidal particles has been shown to occur naturally in systems with colloidal volume fractions ranging from 3% to 50% [58].The complexity of active particle aggregation is noteworthy, with behavior resembling living cell aggregates emerging as active colloids self-organize into clusters [59].Motility-induced clustering and phase separation have been demonstrated in quasi-2D environments [60].Under weak confinement, autophoretic active particles exhibit both clustering and dispersion modes, dependent on the driving force and interparticle interactions [61].In many cases, the presence of geometrical confinement grants far different dynamics to active particles as compared to the more widely modeled unbounded cases, highlighting the importance of considering explicit confinement when modeling non-bulk systems [62,63].The presence of physical or other manner of confinement has a profound effect on the behavior of active particles, as the existence of a boundary introduces complexity to the motion of the driven motors, leading to phase separation, boundary clustering, and other emergent behaviors [64][65][66][67].Confinement also serve as a tool for predicting and controlling the behavior of active particles such as micromotors [68].Active colloidal particle aggregation has been shown in previous works in the literature to be heavily dependent on the size and surface properties of the particles [69,70].In addition, the application of external electric or magnetic fields can also drive the aggregation process [51,71].However, there remains a lack of research towards the confined collective behavior of active colloids, specifically covering a range of volume fractions and driving forces, especially as it relates to their aggregation [72][73][74].Consequently, the micrometer-scale behavior of active particles, especially within confined systems, warrants comprehensive investigation, in order to supplement the rapidly growing literature on computational methods for modeling active matter [75].
The focus of this work is to investigate the effects of the density and driving force of active colloidal particles on their clustering behavior within a confined fluid system.To achieve this, a hybrid coarse-grained molecular dynamics model is employed to simulate the clustering and aggregation dynamics of micrometer-sized colloidal particles.Given the crucial role of hydrodynamics in the clustering process, an explicit representation of the fluid is incorporated into the model.Stochastic rotation dynamics, also known as multiparticle collision dynamics, is used to well replicate both long-and short-range hydrodynamic behaviors at a low Reynolds number.The key variables of interest are the particle coverage ratio θ and the driving force F. Here, θ denotes the 2D projection area coverage ratio of the colloidal particles to the 2D size of the simulation box, while F indicates the driving force acting on the colloidal particles along their axial direction.Results reveal that both coverage and driving force have strong effects on the aggregation behavior of the colloidal particles.A set of empirical power laws were derived to describe the rate of cluster formation over time, in which the scaling exponents strongly depend on the particle coverage ratio and driving force applied on each particle.This comprehensive investigation sheds light on the nuanced interplay between particle density, driving force, and clustering behavior, providing valuable insights into the underlying dynamics of active colloidal systems within confined fluid environments.

Methods and Models
The system of interest involves passive or active particles at the micrometer scale, which can be termed as colloids [76][77][78].In these colloidal systems, the suspended particles may interact through hydrodynamic, interparticle, and Brownian (or thermal) forces [79].Our system directly modeled hydrodynamic and Brownian forces through stochastic rotation dynamics, and the interparticle interactions were modeled with the Derjaguin-Landau-Verwey-Overbeek (DLVO) interaction.A multiparticle collision-based method was used to model the solvent.An extended explanation on the choice of DLVO theory to model this system can be found in the Supplemental Materials.

SRD-Explicit Hydrodynamic Solvent
There is a tremendous difference in scale for both length and kinematic time between colloidal particles and solvent beads.For example, a typical colloidal particle (1 µm diameter) occupies a volume comparable to that of 10 3 water molecules.Thus, coarse graining of the solvent is necessary in order to reasonably reach significant scales of time and size for the simulated system.In this work, we used stochastic rotation dynamics (SRD), also called multiparticle collision dynamics, to establish a computationally efficient and hydrodynamically accurate fluid system capable of handling thermal noise [80].The coarsegraining length scale was chosen to be smaller than that of mesoscopic colloidal particles but much larger than the natural length scales of the solvent molecules (in this case, water).By adhering to local momentum conservation, the SRD method reproduces Navier-Stokes hydrodynamics at larger length scales.A detailed rationale for choosing SRD in this work is given in the Supplemental Materials (Section of stochastic rotation dynamics).
In the framework of SRD, the solvent is represented by a large number of point-like beads, each possessing a mass of 1 m s (an explanation of units is given below).These beads are termed as the fluid beads, and it is crucial to note that they are not merely composite particles or clusters composed of aggregates of water molecules.Instead, these beads can be considered as a convenient computational tool to facilitate the coarse graining of fluid properties.The SRD works as follows: first is the streaming step, during which the positions and velocities of the fluid beads are calculated directly by integrating Newton's equations of motion: The forces acting on the fluid beads are produced through collisions with the system walls or active colloidal particles.Notably, there are no direct forces considered between fluid beads, and the simplification of such pairwise calculations significantly contributes to the overall efficiency of the SRD method.Moving to the second step, the collision step simulates the collisions between fluid beads.The system is partitioned into cubic cells, each having a length of 1 σ.During this step, the velocities of all fluid beads are rotated relative to the center of mass velocity (V) of their respective cell, as given below: During this step, momentum is effectively transferred between the fluid beads, and the rotation procedure can, thus, be viewed as a coarse graining of particle collisions across both time and space.As a result of the local conservation of mass, momentum, and energy, accurate Navier-Stokes hydrodynamic effects can be captured, including those stemming from thermal noise.The SRD method's coarse graining of the fluid allows for easy control over viscosity and coupling properties, enabling the depiction of phase segregation and reactive hydrodynamics arising from complex solutes.As a validation study to prove the suitability of SRD for studying the behavior of the fluid beads under a confinement environment, Figure S1 illustrates the inverse relationship between velocity and size in Stokes flow for spherical beads of different radii traveling through the SRD beads, with the simulation results aligning well with theory.Figure S2 shows the mean squared displacement curves for beads subjected to different applied forces, along with their trajectories, showcasing the enhanced diffusion resulting from the driving force acting on the active particles.Notably, the absence of an applied driving force is equivalent to the Brownian motion observed in passive colloidal particles.All simulations were conducted using the LAMMPS software (2 August 2023 version) [81].

Coarse-Grained Molecular Dynamics (CGMD) Active Particles and DLVO Potential
The system of interest in this paper involves a group of attractive colloidal active particles immersed in a fluid under a quasi-2D confinement scheme, as shown in Figure 1.In this configuration, the active particles (schematically depicted as gray beads in Figure 1a) are subjected to various influences, including a driving force resulting from external fields (such as a magnetic field or electric field), DLVO interactions originating from neighboring active colloids, drag force due to hydrodynamic effects, and random forces arising from the thermal fluctuations within the system.In the simulation, colloidal spheres with a mass M are propagated through the velocity Verlet algorithm using a timestep of 0.01 τ.These colloids are immersed in the fluid and interact with the fluid beads through a repulsiveonly Weeks-Chandler-Andersen potential (van der Waals force).While the fluid-fluid interactions are coarse-grained using SRD, the colloid-colloid and colloid-fluid interactions are integrated using a normal molecular dynamics procedure.
study to prove the suitability of SRD for studying the behavior of the fluid beads under a confinement environment, Figure S1 illustrates the inverse relationship between velocity and size in Stokes flow for spherical beads of different radii traveling through the SRD beads, with the simulation results aligning well with theory.Figure S2 shows the mean squared displacement curves for beads subjected to different applied forces, along with their trajectories, showcasing the enhanced diffusion resulting from the driving force acting on the active particles.Notably, the absence of an applied driving force is equivalent to the Brownian motion observed in passive colloidal particles.All simulations were conducted using the LAMMPS software (2 August 2023 version) [81].

Coarse-Grained Molecular Dynamics (CGMD) Active Particles and DLVO Potential
The system of interest in this paper involves a group of attractive colloidal active particles immersed in a fluid under a quasi-2D confinement scheme, as shown in Figure 1.In this configuration, the active particles (schematically depicted as gray beads in Figure 1a) are subjected to various influences, including a driving force resulting from external fields (such as a magnetic field or electric field), DLVO interactions originating from neighboring active colloids, drag force due to hydrodynamic effects, and random forces arising from the thermal fluctuations within the system.In the simulation, colloidal spheres with a mass M are propagated through the velocity Verlet algorithm using a timestep of 0.01 .These colloids are immersed in the fluid and interact with the fluid beads through a repulsive-only Weeks-Chandler-Andersen potential (van der Waals force).While the fluid-fluid interactions are coarse-grained using SRD, the colloid-colloid and colloid-fluid interactions are integrated using a normal molecular dynamics procedure.The CGMD colloid-colloid interaction potential is described by the DLVO potential specifically developed for colloidal interactions [82]: Here, the total potential energy between two colloidal particles is the summation of the van der Waals (attractive) force, Coulombic (screened electrostatic) force, and Hertz (hard repulsive contact) force.The interparticle colloidal forces are given as: where, for the van der Waals interaction, A H is the Hamaker constant, d is the colloidal particle diameter (1 µm), and r is the distance between the colloidal particle centers.For the Coulombic interaction, r is the electric permittivity of the fluid, 0 is the vacuum permittivity (8.854 × 10 −12 F/m), q is the colloidal particle charge, Ψ 0 denotes the effective surface potential of the colloid, and κ is the inverse Debye screening length.To model steric repulsion between colloidal particles and prevent overlap or penetration of the particles, a Hertz contact force is included: The SRD fluid beads interact with the CGMD particles via the Weeks-Chandler-Anderson potential: The hydrodynamic volume fractions of interest range from 5% to 30%; at higher volume fractions, steric interactions between the colloids start to dominate and frustrated systems can occur [83].All quantities are described in dimensionless units with mass in units of m s (the mass of the solvent particles), energy in units of , length in units of σ, force in units of σ , and time in units of τ = m s σ 2 / .
Figures 1b and 1c shows a representative snapshot of the CGMD-SRD hybrid model used for the simulations with CGMD active colloids and without active colloids in SRD solvent, respectively.In these snapshots, solvent beads are colored based on their relative velocities to illustrate their distribution.The SRD beads effectively capture both longand short-range hydrodynamic behaviors, while the CGMD particles serve as models for colloidal micron-sized active particles (micromotors).The CGMD active particles are finite-sized spheres with a radius of 1σ and a density of 1m s /σ 3 .The simulation box spans 256 σ × 256 σ in the x-y plane and 10σ in the z dimension.The confinement size was specifically chosen, designed to have an energy minimum width equivalent to that of a single CGMD bead.This selection allows beads to pass by one another, preventing geometric frustration.Simultaneously, the confinement force ensures that there is no appreciable cluster size in the z direction; thus, the 2D x-y spectral analysis is directly indicative of the size of the cluster.
A harmonic restraining force, applied perpendicular to the z direction, acts solely on the CGMD colloidal particles.The simulation box is periodic in all dimensions.To maintain a quasi-2D setup with a geometric confinement environment, a collision force between CGMD active particles and the system wall is implemented.This ensures that the CGMD active particles move within a plane of approximate thickness 3σ at the center, discouraging particle stacking in the z dimension while permitting particles to cross over and pass each other.
Initially, all CGMD active particles are regularly placed on a square lattice, as illustrated in Figure S3, with an areal density or particle coverage ratio (the 2D projection area coverage ratio of the colloidal particles to the 2D size of the simulation box) denoted as θ = 30%.In the simulation, θ = 5%, 10%, 15%, 20%, 25%, and 30%, each value corresponding to the number of active particles N = 1089, 2116, 3136, 4225, 5184, and 6241, respectively.Upon particle creation, each particle has its orientation randomly set, initiating the driving force to act in a random direction at the beginning of the simulation.The driving force applied to the active particles is consistently directed along the axis of each particle.Due to the reorientation of spherical particles in a fluid, there is no net directional movement resulting from the driving force.Each active particle imparts a constant driving force F through its geometric center, with F = 0, 0.5, 1, 5, and 10.As a specific case study, Figure S4 demonstrates the clustering of particles over time and depicts their nearest neighbors.

Dynamic Clustering Behavior
The fundamental quantities of interest for this system are the number and distribution of clusters, which change over time due to the diffusion-limited aggregation of individual particles and clusters.In the context of this confined setup with a predetermined number of particles, the system tends to converge towards fewer and larger clusters.Individual particles aggregate to form clusters, and as time progresses, these clusters undergo further aggregation.

Time Evolution of Cluster Number and Cluster Size
Figure 2 gives a visual representation of how both time and driving force affect the clustering behavior for a given coverage (in this case, θ = 15%).An increase in the active driving force causes more rapid clustering behavior, especially in the initial stages of the simulation.This leads to quicker formation and aggregation of clusters compared to the passive colloidal particles (those with no driving force; F = 0).The dynamics of clustering are intrinsically changed by the addition of a driving force, exhibiting more complex effects than a simple acceleration of the observed clustering.As evidenced in Figure 3, both coverage and driving force significantly contribute to the time evolution of clustering behavior, with snapshot times at t = 1000τ.Higher coverages promote faster clustering, particularly at smaller times, due to the decreased average inter-bead distances resulting from the geometrical realities of increased particle coverage.As particle aggregation into clusters proceeds, the cluster size becomes a crucial metric used in conjunction with total cluster number to determine the uniformity of clustering behavior.For this section, cluster size is calculated as the total number of particles within each cluster.
Figure 4 illustrates the cluster size distribution P S and the mean cluster size S, where S represents the number of particles within an individual cluster.Fitting lines are provided to determine the time evolution of S and the peak of the probability distribution, demonstrating how the exponents can be used to scale the cluster size distributions to a single curve.The figure shows that, for a given section of the curve, clustering behavior can be effectively represented by a simple power law.This form of scaling law for diffusionlimited cluster aggregation and growth has been shown to fit quite well with findings in previous works in the literature [84,85].This power law allows for the documentation of not only instantaneous but also dynamic behavior in such a clustering process, providing a comprehensive characterization for future reference.This scaling is based upon the stochastic nature of clustering, which tends to smooth out for large sample sizes.Figure S5 further demonstrates power-law behavior and collapse to a single curve for θ = 20%, F = 0.5, 465τ < t < 6000τ.The lower number of clusters in this case may result in more obvious scattering from the fit observed at lower time intervals.The probability distribution of the cluster size can be fitted well by the gamma distribution: The probability distribution of the cluster size can be fitted well by the gamma distribution: The gamma distribution has a shape parameter ( > 0), a scale parameter ( > 0), a cluster size ( > 0), and Γ() is the gamma function.As shape parameter () increases, the distribution tends to approach a normal distribution.The gamma distribution is useful in applications where there is a physical lower bound but no non-statistical upper bound.Figure S6 demonstrates a fitting of the gamma distribution to the cluster size over time for various  and  .The fitting parameters,  and  , can be plotted to visualize the skewness and mean of the distribution.Specifically, the mean of the distribution  ̅ is calculated as  * , and the skewness is inversely proportional to the shape, 2 √ ⁄ .The plot of the scale parameter, closely linked with cluster size, reveals a flattening or reduced clustering in the middle region, as shown in Figure S6.Additionally, the mean cluster size The gamma distribution has a shape parameter (a > 0), a scale parameter (b > 0), a cluster size (S > 0), and Γ(a) is the gamma function.As shape parameter (a) increases, the distribution tends to approach a normal distribution.The gamma distribution is useful in applications where there is a physical lower bound but no non-statistical upper bound.Figure S6 demonstrates a fitting of the gamma distribution to the cluster size over time for various θ and F. The fitting parameters, a and b, can be plotted to visualize the skewness and mean of the distribution.Specifically, the mean of the distribution S is calculated as a * b, and the skewness is inversely proportional to the shape, 2/ √ a.The plot of the scale parameter, closely linked with cluster size, reveals a flattening or reduced clustering in the middle region, as shown in Figure S6.Additionally, the mean cluster size is observed to increase as a power function of elapsed time, following the form S ∝ t m , where m is dependent on coverage and force.

Cluster Dynamic Spatial Distribution
Another aspect of clustering behavior, along with cluster size and number, is the physical distribution of the clusters in space [86].The periodic system can be analyzed through fast Fourier transformation (FFT) of snapshots captured in the x-y plane, as the quasi-2D nature of this system facilitates spatial analysis.Figure S7 illustrates a representative FFT transformation of a clustering image, including the circularly averaged power profile for θ = 15%, F = 0, t = 5000τ.For each snapshot used for analysis, an image is generated as seen in Figure S7a.Images are composed of 1024 × 1024 pixels for the 256 × 256 µm system (16 pixels/µm 2 resolution) and generated as completely black and white.Figure S7b shows a representative 2D fast Fourier transform (FFT) of a clustering image where white pixels (background) are treated as zeros, while black pixels (colloidal particles) are treated as ones.The circular average of the power of the FFT transform is calculated, and the k-space power spectra are plotted as seen in Figure S7c,d.To find the representative length k 0 , a Gaussian curve was fitted to the first peak for every snapshot.The values of k 0 and PSD(k 0 ) were determined and plotted, as shown in Figure 4c,d.Figure 5 demonstrates that, for each given dynamic region, a fitting power law allows for scaling of the FFT spectra by the power-law exponent.This scaling allows all snapshot spectra over the given region to collapse onto a single curve after scaling.More information on the FFT image analysis is provided in the Supplemental Material.

Phase Diagram of Dynamic Characteristics
This section presents definitions and descriptions of the different dynamic regions observed during clustering, considering the evolution of the number of clusters, cluster size distribution, as well as the spatial analyses which are discussed in detail in the preceding sections.Figure 6 shows the power-law curves fitting the different dynamic regions, and their intersection is termed as the critical time for the distinction of regions.The pre-clustering, individually separated particles prior to clustering were not suitable for FFT calculation, as demonstrated by Figure S8.In this figure, the initial square-lattice geometrical configuration of particles at setup causes crystalline FFT spectra, and the crystalline peaks gradually disappear as the system evolves to amorphous clustering.Notably, there is little to no clustering happening during this period, and consequently, no scaling laws for cluster distribution or spacing are presented or discussed for the preclustering behavior in this work.Figure 3 shows representative snapshots of low, medium, and high coverage and force to demonstrate how the snapshots used for FFT analysis change as a function of these parameters.
Figure S9 shows the same cases as Figure 3, but at a timestep of 5000τ.Over this duration, clustering proceeds towards fewer, larger clusters across each of the cases.However, the spatial distribution of clusters remains relatively uniform, as observed consistently across the five independent runs for every variation of F and θ.Due to this spatial uniformity, the FFT analysis used can reasonably be relied upon to determine the average intercluster spacing distance, as portrayed by Figure S7.This intercluster spacing distance is represented by its inverse, k 0 .It can be seen that higher coverage and driving force both accelerate the rate of clustering.At higher θ, there are more particles in solution and, thus, the rate of collision and clustering is increased.For higher F, the effective diffusion constant for the colloidal particles is higher and, thus, clustering proceeds following a scaling law with a higher exponent.Both FFT spatial distribution and clustering number can be well represented by their power-law scaling exponents.

Phase Diagram of Dynamic Characteristics
This section presents definitions and descriptions of the different dynamic regions observed during clustering, considering the evolution of the number of clusters, cluster size distribution, as well as the spatial analyses which are discussed in detail in the preceding sections.Figure 6 shows the power-law curves fitting the different dynamic regions, and their intersection is termed as the critical time for the distinction of regions.As seen in Figure 6, there are only two distinct dynamic regions for θ = 5% and 10%, but for θ = 15%-30%, four distinct dynamic regions occur.Table 1 briefly describes the different dynamic regions, with a detailed description as follows: The first dynamic region, Region I, is pre-clustering, where the average cluster size is approximately one (monomer region).As the coverage density grows higher, the time for this region becomes shorter due to the decrease in inter-bead spacing in the initial configuration.The length of Region I also decreases over time with an increase in driving force F, due to the increased average velocity of the active particles.In Region I, virtually no clustering can be seen; this pre-clustering region can be considered as a preliminary stage where the average cluster size is effectively zero.This region is visible in Figure 4a at t < 100 τ or in the initial stages of Figure 7c.Region I becomes vanishingly small as the coverage or the driving force increases, due to the accelerated initiation of clustering at higher particle densities and velocities (Figure 7).Region II can be considered as the initial clustering behavior, and this region dominates for θ = 5-10%.As seen in Figure 7d, the driving force F has minimal effect on the clustering behavior of this region at higher coverage densities.Region III demonstrates a marked slowing of clustering, especially prominent at higher θ and lower F.This region is characterized by larger clusters and a near-complete disappearance of the monomer phase, as seen in the 500 τ snapshots given in Figure 2. At exceptionally large driving forces, this region becomes far less prominent.The final region is Region IV, representing the long-time clustering behavior.The clustering here is comparable to that observed in Region II.Clustering in this region is characterized by large clusters flocculating towards the limit of one single cluster within the system.Figure 8 displays a 3D representation of the different dynamic regions, highlighting the critical time separating them as a function of both driving force F and coverage θ.It can be seen that the transition between Region I and Region II is smooth and monotonic across all F and θ, while the transition between Region III and Region IV exhibits greater variation.Table 2 offers a detailed breakdown of the critical times for transitions between dynamic regions.Figure S10 demonstrates that, for each given dynamic region (represented by Region IV, θ = 15% and F = 0), a fitting power law allows for scaling of the FFT spectra by the power-law exponent, which allows for all snapshot FFT spectra over the given region to collapse to a single curve.Figures S11-S15 show the representative clustering behavior evolution over time for each coverage θ and driving force F, providing a comprehensive overview of the procession of the clustering process.
Nanomaterials 2024, 14, x FOR PEER REVIEW 11 of 21 force increases, due to the accelerated initiation of clustering at higher particle densities and velocities (Figure 7).Region II can be considered as the initial clustering behavior, and this region dominates for  = 5% − 10% .As seen in Figure 7d, the driving force F has minimal effect on the clustering behavior of this region at higher coverage densities.Region III demonstrates a marked slowing of clustering, especially prominent at higher θ and lower F.This region is characterized by larger clusters and a near-complete disappearance of the monomer phase, as seen in the 500 τ snapshots given in Figure 2. At exceptionally large driving forces, this region becomes far less prominent.The final region is Region IV, representing the long-time clustering behavior.The clustering here is comparable to that observed in Region II.Clustering in this region is characterized by large clusters flocculating towards the limit of one single cluster within the system.Figure 8 displays a 3D representation of the different dynamic regions, highlighting the critical time separating them as a function of both driving force  and coverage .It can be seen that the transition between Region I and Region II is smooth and monotonic across all  and , while the transition between Region III and Region IV exhibits greater variation.Table 2 offers a detailed breakdown of the critical times for transitions between dynamic regions.Figure S10 demonstrates that, for each given dynamic region (represented by Region IV,  = 15% and  = 0), a fitting power law allows for scaling of the FFT spectra by the power-law exponent, which allows for all snapshot FFT spectra over the given region to collapse to a single curve.Figures S11-S15 show the representative clustering     One of the most important parameters influencing the clustering of colloidal particles is their packing density, represented in this work by the coverage, θ.The distinct behaviors observed in Figure 7 at low and high coverages are primarily influenced by the increased particle proximity at the beginning of the simulation and the overall larger and more numerous clusters as time proceeds.Herein, we leverage the qualitative behaviors discussed in the previous section and use the fitting exponents to quantitatively describe and scale the behavior.For each given dynamic region, the power-law fitting exponent is given in Table 3. Figure 7 demonstrates different clustering evolution behaviors for various F and θ.As can be seen, driving force F tends to flatten and accelerate the clustering behavior, in line with previous reports of driving force causing an increased effective diffusion coefficient for active particles.The clustering coverage θ, however, has a pronounced effect on the shape of the cluster number curve, with simple clustering behavior observed for low coverage (θ = 5%, 10%) and a more complex time evolution of the cluster number for higher coverage (θ > 10%).Figure S16 shows all of the five independent runs for θ = 5% at different driving forces.The similarity between cases remains consistent across different θ and F, while at F = 10, occasional temporary detachment of a bead from a cluster is observed due to the similar magnitude of the driving force and inter-bead attractive potential.From Figure S16, we can see that the average cluster number and cluster size can be taken as a good approximation for all runs across each variation of driving force F for a given coverage density θ. Figure S17 shows that only θ matters in terms of cluster number over time, and absolute bead number is not significant.For systems with 16x more or fewer beads, there was no big difference in clustering behavior over time, although there were minor size effects towards the end of the run when the number of beads was small (N 0 < 500).This can be attributed to the discrete behavior of very few (N < 10) clusters when compared with the more continuous behavior of many (N ≥ 100) clusters.Figure S18 shows the absolute cluster number over time for various θ and F (rather than the normalized number shown in Figure 7).This visualization helps to highlight the rapid clustering behavior seen at higher θ and also the similarities in long-time behavior (t > 3000 τ).
Figure 9 compiles different scaling exponents for the cluster size distribution and 2D cluster spatial distribution power laws across all coverages and driving forces for Region II.Within a given region, clear variations in clustering behavior as a function of θ are observed, with more pronounced differences in Region III and Region IV. Figure 10 demonstrates the differences between cluster size scaling exponents for Region II.The scaling exponent for Region II is not heavily affected by the coverage percentage, suggesting that early clustering behavior can be attributed to diffusion-driven or enhanced diffusion for F > 0, with less dependence on the density of colloidal particles or micromotors.Figure 11 illustrates the differences between cluster scaling exponents (m k ) and cluster spatial distribution scaling exponents (h k ) for Regions II, III, and IV for θ = 15-30%.As θ becomes larger, the scaling exponent for Region III decreases, indicating a marked slowing of the clustering behavior.However, as F increases, the scaling exponent for Region III increases rapidly, reaching a near-constant value for F = 10 (near 0.8 as seen in Table 3).The different behaviors observed at low and high driving forces highlight the distinction between passive colloidal clustering and the clustering of active colloidal micromotors.Active motion not only increases the effective diffusion coefficient of Brownian particles but influences the scaling of clustering times as well [87][88][89].
Nanomaterials 2024, 14, x FOR PEER REVIEW 15 of 21 or micromotors.Figure 11 illustrates the differences between cluster scaling exponents ( ) and cluster spatial distribution scaling exponents (ℎ ) for Regions II, III, and IV for  = 15% − 30% .As  becomes larger, the scaling exponent for Region III decreases, indicating a marked slowing of the clustering behavior.However, as  increases, the scaling exponent for Region III increases rapidly, reaching a near-constant value for  = 10 (near 0.8 as seen in Table 3).The different behaviors observed at low and high driving forces highlight the distinction between passive colloidal clustering and the clustering of active colloidal micromotors.Active motion not only increases the effective diffusion coefficient of Brownian particles but influences the scaling of clustering times as well [87][88][89].

Conclusions
In summary, we have established a hybrid coarse-grained molecular dynamics simulation model to investigate the dynamic clustering behavior of active particles under confinement.Several scaling laws were identified, capturing both cluster behavior and overall morphological changes.For the cluster behavior, the average cluster size S was found to scale with time (t) as S ∝ t m (0 < m < 1), and the amplitude of the cluster distribution probability P S followed P S ∝ t h (−1 < h < 0).For the morphological behavior, the representative wavevector k 0 , characterizing the quasi-periodic spatial distribution of clusters, scaled as k 0 ∝ t m k (−1 < m k < 0), and the power spectral density value at k 0 , PSD(k 0 ), obeyed PSD(k 0 ) ∝ t h k (0 < h k < 1).Detailed study shows that, depending on the colloidal particle coverage θ, the dynamic clustering process can be divided into four different regions: Region I (pre-clustering), Region II (initial clustering), Region III (initial cluster-cluster merging), and Region IV (final cluster-cluster merging).The behaviors in these regions were shown to strongly depend on θ and F: at low θ and F, diffusion-induced aggregation dominates the clustering dynamics; and cluster-cluster merging plays the most important role at high θ.The role of F was primarily observed in shortening the periods of Regions I and II and accelerating the cluster-cluster merging process.Due to limitations of computational resources, it is determined that a 2D system would provide valuable insights into clustering, as may occur on films or liquid surfaces or two liquid interfaces, at a larger scale than would be possible with a 3D system.However, future experiments may pursue whether the scaling laws observed in this work also present in cluster formation in 3D systems.
The study of active particles is irrevocably tied to the future of various scientific fields, especially those pertaining to medicine or health.Active particles, such as micromotors, hold promise in performing controllable tasks at the microscale, akin to natural bacteria or cells.As the field of active matter advances, active particles may become integral in scientific research work and potentially find applications in our daily lives.This work contributes to the understanding of the aggregation and possible removal of micromotors in medical and waste removal fields, particularly within laminar blood flow at biologically relevant Reynolds numbers.

Figure 1 .
Figure 1.(a) Schematic of the spherical active particles, with on-axis driving force, under confinement.The particles can attract each other (cluster) at close distances due to the DLVO interactions.The explicit SRD fluid grants hydrodynamic effects such as drag forces in the opposite

Figure 1 .
Figure 1.(a) Schematic of the spherical active particles, with on-axis driving force, under confinement.The particles can attract each other (cluster) at close distances due to the DLVO interactions.The explicit SRD fluid grants hydrodynamic effects such as drag forces in the opposite direction of motion and random thermal forces.Model setup: (b) The CGMD beads (red) in SRD solvent, and (c) the solvent with transparent beads to show solvent density.Solvent beads are colored as a function of their relative velocity, to demonstrate the distribution.

Nanomaterials 2024 ,
14,  x FOR PEER REVIEW 7 of 21 large sample sizes.FigureS5further demonstrates power-law behavior and collapse to a single curve for  = 20%,  = 0.5, 465 <  < 6000.The lower number of clusters in this case may result in more obvious scattering from the fit observed at lower time intervals.

Figure 3 .
Figure 3. Clustering behaviors at different densities and driving forces.Time is t = 1000 τ.

Figure 3 .
Figure 3. Clustering behaviors at different densities and driving forces.Time is t = 1000 τ.Figure 3. Clustering behaviors at different densities and driving forces.Time is t = 1000 τ.

Figure 3 .
Figure 3. Clustering behaviors at different densities and driving forces.Time is t = 1000 τ.Figure 3. Clustering behaviors at different densities and driving forces.Time is t = 1000 τ.

Figure 4 .
Figure 4. Cluster size distributions over time for  = 15% ,  = 5 .(a) The mean of the cluster number distribution, with fitting line for the region  = 150  to  = 3000 ; (b) the peak height of the cluster number distribution, with fitting line.(c) The circularly averaged FFT peak center for k0, with fitting line for the region of interest; (d) the peak intensity of k0, with fitting line for the region of interest.Dotted lines show the beginning and end of the region of interest.

Figure 4 .
Figure 4. Cluster size distributions over time for θ = 15%, F = 5.(a) The mean of the cluster number distribution, with fitting line for the region t = 150 τ to t = 3000 τ; (b) the peak height of the cluster number distribution, with fitting line.(c) The circularly averaged FFT peak center for k 0 , with fitting line for the region of interest; (d) the peak intensity of k 0 , with fitting line for the region of interest.Dotted lines show the beginning and end of the region of interest.

Nanomaterials 2024 , 21 Figure 5 .
Figure 5. Cluster dynamics over time for  = 15%,  = 5.(a) The cluster size distributions for  = 150  to  = 3000 .(b) The same cluster size distribution after being scaled by the scaling function for the pertinent region.(c) The circularly averaged FFT spectra for  = 150  to  = 3000 .(d) The same spectra after being scaled by the scaling function for the pertinent region.

Figure 5 .
Figure 5. Cluster dynamics over time for θ = 15%, F = 5.(a) The cluster size distributions for t = 150 τ to t = 3000 τ.(b) The same cluster size distribution after being scaled by the scaling function for the pertinent region.(c) The circularly averaged FFT spectra for t = 150 τ to t = 3000 τ.(d) The same spectra after being scaled by the scaling function for the pertinent region.

Figure 6 .
Figure 6.Demonstration of fitting curves for the average cluster size  ̅ and the associated different dynamic regions for (a)  = 5%, (b)  = 10%, (c)  = 15%, and (d)  = 30%. = 0 for all cases shown.The intersection points of fitting curves are denoted by asterisks.

Figure 8 .
Figure 8. Regions I, II, III, and IV (plotted as filled colored volumes), and the transition times between them (plotted as gridded surfaces), as a function of both coverage () and driving force ().The surfaces have been smoothed between discrete points using a modified cubic spline interpolation.

Figure 8 .
Figure 8. Regions I, II, III, and IV (plotted as filled colored volumes), and the transition times between them (plotted as gridded surfaces), as a function of both coverage (θ) and driving force (F).The surfaces have been smoothed between discrete points using a modified cubic spline interpolation.

Figure 9 .
Figure 9. Scaling exponents (a) m, (b) h, (c) mk, and (d) hk for Region II across all θ and F tested.

Table 3 .
Dynamic region scaling exponent m for Regions II-IV.