Dynamical self-assembly of dipolar active Brownian particles in two dimensions

Abstract

Based on Brownian Dynamics (BD) simulations, we study the dynamical self-assembly of active Brownian particles with dipole–dipole interactions, stemming from a permanent point dipole at the particle center. The propulsion direction of each particle is chosen to be parallel to its dipole moment. We explore a wide range of motilities and dipolar coupling strengths and characterize the corresponding behavior based on several order parameters. At low densities and low motilities, the most important structural phenomenon is the aggregation of the dipolar particles into chains. Upon increasing the particle motility, these chain-like structures break, and the system transforms into a weakly correlated isotropic fluid. At high densities, we observe that the motility-induced phase separation is strongly suppressed by the dipolar coupling. Once the dipolar coupling dominates the thermal energy, the phase separation disappears, and the system rather displays a flocking state, where particles form giant clusters and move collective along one direction. We provide arguments for the emergence of the flocking behavior, which is absent in the passive dipolar system.

---

1 Introduction

Self-propelled particles are capable of converting energy from an internal source or the surroundings into their own active motion.^1,2 It is now well established that active motion of individual particles, combined with different types of interactions among these particles, can lead to remarkable collective behavior. For example, repulsive interactions alone can generate the so-called motility-induced phase separation, where large, densely packed clusters coexist with freely moving particles.³ A standard model describing this non-equilibrium phase transition is that of active Brownian particles, where each particle moves with a constant propulsion speed along a direction subject to white noise, and the interactions are purely steric and isotropic. Clearly, one could expect more complex behavior in active systems with anisotropic interactions. This situation has been investigated by a large number of theoretical and simulation studies.^4–7 A famous example in this direction is the Vicsek model, where active, point-like particles interact such that they tend to align their velocities with those of their neighbors. The resulting collective behavior includes traveling bands⁶ and flocking.⁸ A further representative system is a suspension of active Brownian particles with polar interactions favoring parallel orientations of the velocities irrespective of the spatial configuration. This system displays a state with moving lanes and bands even at low densities, where conventional active Brownian particles do not show significant behavior.⁹ Further, at high densities, the polar coupling favors motility-induced phase separation compared to the non-polar case.¹⁰

In most studies of active systems with anisotropic interactions, they are assumed to be of short range.¹¹ Less effort has been spent on active systems involving long-range anisotropic interactions, such as dipole–dipole interactions stemming from intrinsic (permanent) dipole moments or dipoles induced by an electric or magnetic field.¹² A notable feature of dipolar interactions is that they depend not only on the orientations of the two involved particles, but also on their configuration in space. Specifically, two particles with point dipoles at their center tend to align head-to-tail, where the arrowhead of one dipole moment is directed toward the arrowtail of the other dipole. In contrast, particles configured side by side favor antiparallel alignment. Due to this complexity, one would expect a system of dipolar active particles to show macroscopic structures significantly different from active models with simple polar interactions, such as the Vicsek model⁸ and its variants.

There are some recent studies, which have focused particularly on the dynamics of dipolar active colloids. For example, ref. 13 and 14 have studied the structural transformation of small dipolar clusters under the impact of activity starting from different (meta)stable configurations, either neglecting¹³ or including¹⁴ hydrodynamic interactions. Another example is a system of spherical, ferromagnetic rollers confined to a fluctuating surface. Here it has been shown, both, experimentally and in simulations, that these particles can exhibit swarming or vortex patterns when energized by a vertical alternating field.¹⁵ Owing to the interplay between self-propulsion and dipolar interactions, these vortices persist even on a flat surface.¹⁶ In nature, magnetotactic bacteria are known to sense the earth's magnetic field and move along or against the field direction.^17–19 Experiments have observed that these bacteria, when confined in a microfluidic channel and placed under an external magnetic field, display clustering behavior along the channel.²⁰ The underlying mechanism of this clustering instability was also investigated analytically and in simulations.²¹ We also mention dipole-like, active Janus particles with two screened electric charges in each hemisphere. It is found that these particles can self-organize into fingerprint-like patterns at high densities.²² Moreover, an external electric or magnetic field can induce two point dipoles in the respective hemispheres of a Janus particle and thereby further complicate the interactions between particles. By fine-tuning the interactions, it has been reported that the external field can be used to control the collective behavior of Janus particles.^23,24 However, these latter models involve even more complex anisotropic interactions as compared to the purely dipolar case. Indeed, the collective behavior of large ensembles of the simplest dipolar active model, that is, dipolar active Brownian particles, is not yet explored. This is the motivation for the present work.

Specifically, we present Brownian Dynamics simulation results for the dynamical self-assembly of a two-dimensional system of active particles, where each particle has a permanent point dipole moment oriented in the plane. In addition, each particle is subject to a self-propulsion force, which is directed along the dipole, as well as to thermal noise. We investigate the collective behavior of our model system for three density regimes, considering a wide range of motilities and dipolar coupling strengths. Our model can be realized, e.g., by Janus particles with a magnetic material coated on one of the hemispheres.^25,26 In that case, both the dipole moment and the propulsion force of the Janus particle are directed along its symmetric axis.

For all three densities we present state diagrams illustrating the complex interplay between essentially three phenomena: chain formation (which already occurs in the passive case), motility-induced phase separation, and polar ordering.

The remainder of this paper is organized as follows. In Section 2 we present our model of dipolar active particles and the simulation details, as well as the target quantities and order parameters investigated. Based on analysis of these quantities for a range of parameters, we discuss the collective behavior of the model system in Section 3. Finally, in Section 4 we summarize our findings.

2 Model and methods of investigation

2.1 Model system

Our dipolar active system consists of N disk-shaped Brownian particles with diameter σ dispersed in a monolayer in the xy-plane. Each particle carries a fluctuating point dipole moment μ_i (i = 1,…,N) located at its center. For passive monolayers of dipolar disks, it is well established that the particles tend to orient along in-plane directions to form chains and rings,^27–31 or dense ordered states,^27,32–34 if the dipolar interactions are sufficiently strong. Having this in mind, we assume beforehand that μ_i lies in the xy-plane, i.e., fluctuations in z-direction are neglected. To model the self-propulsion, we assume that each particle is subject to a force F_i, which has a constant magnitude and is directed along μ_i at each instant of time.

The pair potential between particles located at positions r_i and r_j (i ≠ j) is of the form

u_pair(r_ij,μ_i,μ_j) = u_sr(r_ij) + u_dd(r_ij,μ_i,μ_j), (1)

where the first term on the right-hand side stands for the short-range steric repulsion (sr), which only depends on the distance r_ij = |r_ij| = |r_j − r_i|. We assume that the steric repulsion can be described by the Weeks–Chandler–Anderson potential

In this study, we set the length unit to be σ and fix the repulsive strength ε* = βε = 10, where the thermal energy, β⁻¹ = k_BT, is set to be the energy unit (with k_B being Boltzmann's constant and T being the temperature). This potential is truncated at the cut-off radius r_c = 2^1/6σ, such that eqn (2) and its derivative vanish to zero continuously at the truncation point. The last term in eqn (1) represents the (long-range) dipole–dipole interaction. Its functional form is given by

We define the dipolar coupling strength as λ = βμ²σ⁻³ with μ being the magnitude of a dipole moment, i.e., μ_i = μ[small mu, Greek, circumflex]_i.

2.2 Brownian dynamics simulation

To investigate the system's dynamical behavior, we perform conventional Brownian Dynamics (BD) simulations without hydrodynamic interactions. The motion of the ith particle is then determined by the coupled Langevin equations³⁵ for its position r_i and orientation ê_i = (cos ψ_i,sin ψ_i)^T,

ṙ_i = β[Doublestruck D][F₀ê_i − ∇_riU_i + ξ_i(t)], (4)

[small psi, Greek, dot above]_i = βD_r[−∂_ψiU_i + Γ_i(t)], (5)

where the dots denote time derivatives and ψ_i is the polar angle. In eqn (4) and (5), the potential energy for the ith particle is given bywith u_pair being defined in eqn (1). Since the shape of each particle is modeled as a disk, we set the translational diffusion tensor [Doublestruck D] = D_t[Doublestruck I], where D_t is the (isotropic) translational diffusion constant and [Doublestruck I] is the 2 × 2 identity matrix. Correspondingly, D_r denotes the rotational diffusion constant. As described in ref. 36, we find that the relationship between these two diffusion constants of a hard sphere in the low Reynolds number regime is given via D_r = 3D_t/σ_h². We also define the diameter of an effective (eff) hard sphere via. With the repulsive strength ε* = 10 considered in the present system, σ_eff ≈ 1.07851σ. By choosing σ_h = σ_eff, we obtain D_r = 2.57914D_t/σ².

In eqn (4), the effective propulsion force, which drives the active motion, is given by F_i = F₀ê_i with ê_i = [small mu, Greek, circumflex]_i. For simplicity, in the remainder of this work we present the impact of the effective propulsion force via the motility v₀ = βD_tF₀. Finally, the random force ξ_i(t) and torque Γ_i(t) for the ith particle are zero-mean Gaussian white noise, which satisfy and . The angle brackets 〈⋯〉 denote ensemble average, and the symbol ⊗ represents dyadic product.

Eqn (4) and (5) are solved via the Euler–Maruyama method³⁷ with the discrete time step Δt = 2 × 10⁻⁵τ, where τ is the Brownian diffusion time, given by τ = σ²/D_t. We choose a quadratic box with side length L_x = L_y = L, and we use periodic boundary conditions in both directions. In order to treat the long-range dipole–dipole interactions, we employ a two-dimensional (2D) Ewald summation, as outlined in Appendix A.1. The particle number is set to N = 1156, unless otherwise stated. All simulations are started with randomly oriented particles located on a square lattice. A typical run then consists of at least 5 × 10⁵ steps for reaching a steady state, followed by a production period of 5 × 10⁵ steps. The statistical properties of the dipolar particles (see Section 2.3) are measured every 500 steps. The simulations are carried out at three values of the mean area fraction Φ = Nπσ_eff²/(4L²), where the area of a single particle is defined as πσ_eff²/4. Specifically, we consider the values Φ = 0.12, Φ = 0.23, and Φ = 0.58 (for the exact values, see note ref. 38). To investigate the impact of activity, we perform simulations at different dimensionless propulsion speeds v₀* = v₀σ/D_t and various dipolar coupling strengths λ. We note that v₀* is indeed the same as the commonly used Péclet number Pe.^3,39,40

2.3 Target quantities

In this section we introduce the target quantities which will later be used to characterize the system's behavior. The choice of these quantities is inspired by earlier research on non-dipolar active systems on the one hand, and passive dipolar systems on the other hand.

2.3.1 Clustering behavior. It is well established that active particles with purely repulsive interactions have a tendency to form large clusters and even phase-separate if the motility is sufficiently high (motility-induced phase separation).^39,41–43 We would therefore expect that the present system displays similar behavior at least at low dipolar coupling strengths, i.e., small values of λ. To characterize the clustering behavior of the dipolar active particles, we perform a cluster size analysis based on a simple distance criterion: Two particles are regarded as being in contact if their center-to-center distance is smaller than a distance r_L. A cluster is then a set of particles that are in contact with each other. We quantify the cluster formation by

ϕ_c = 〈n_lcl〉/N, (7)

where n_lcl denotes the number of particles in the largest cluster. The quantity ϕ_c approaches 0 in a state with particles being “essentially uncorrelated”, where not only the one-particle density is a constant, but also the particle correlations (measured by g(r)) are weak. Further, ϕ_c is close to zero in a state with finite-sized chains or clusters with a size much smaller than N. In contrast, ϕ_c ≈1 when the active particles form “giant” clusters with a size comparable to N.³ We choose r_L as the distance corresponding to the first peak in the radial distribution function of a corresponding “reference system” (v₀* = 0, λ = 0) at the given density. In this way we can systematically investigate the impact of the motility and dipolar coupling without changing the cluster criteria. As a result, we obtain r_L ≈ 1.16σ for Φ = 0.12, r_L ≈ 1.13σ for Φ = 0.23, and r_L ≈ 1.12σ for Φ = 0.58 (for the exact values, see note ref. 44).

2.3.2 Global orientational ordering. The appearance of global orientational order in passive systems of dipolar particles has a long history. In three dimensions, the occurrence of ferromagnetic liquid and solid states has been confirmed both by computer simulations⁴⁵ and by theory.⁴⁶ In spatially confined systems, the situation is less clear. A ferromagnetic state has indeed been observed in slab-like systems composed of three and more layers of dipolar particles.^32,47 Systems with less than three layers have been considered in ref. 47. There, one did not find clear hints for the ferromagnetic order in very thin films, consistent with other studies.^48,49

The behavior of the global orientational ordering has also been investigated in models of self-propelled particles with “velocity-alignment” interactions, such as the Vicsek model⁸ and its variants (see, e.g., ref. 9). In these systems, one observes so-called flocking states, where particles gather together and move collectively toward a certain direction. In the present system, the orientational (dipolar) interaction is more complicated (than the Heisenberg-like interactions in the Vicsek model); still, one could imagine that activity-induced flocking states also emerge in the dipolar active system. To this end we consider (as is common in the literature on flocking states⁹) the parameter

which corresponds to the magnitude of the average orientation. This order parameter is unity if all particles self-propel toward the same direction, and zero if the particle orientations are uncorrelated.

2.3.3 Chain formation. At low densities, passive particles with strong dipole–dipole interactions (λ > 1) are known to self-assemble into chains and rings.^{28,29,45,50–57} The reason is easily seen from eqn (3): Considering two dipolar spheres separated by a distance σ, the configuration with the lowest energy is the head-to-tail configuration (with u_dd = −2μ²σ⁻³). In contrast, if the particles are arranged side by side, the energy reaches a maximum for parallel orientation (with u_dd = μ²σ⁻³), suggesting that such a configuration is energetically unfavorable.

To characterize chain formation in the dipolar active system, we use different strategies depending on the density regime considered. The first strategy is adequate for low densities (Φ ≲ 0.23), where the chains can essentially be considered as isolated objects. In this situation, we consider a chain as a set of, at least, three particles which are mutually “bonded.” Specifically, two particles i, j are regarded as “bonded” in the chain, if the following criteria are fulfilled: |r_ij| ≤ r_p, [small mu, Greek, circumflex]_i·[small mu, Greek, circumflex]_j > 0, and ([small mu, Greek, circumflex]_i·r_ij)([small mu, Greek, circumflex]_j·r_ij) > 0. Here, r_p = 1.25σ is set to a distance between the location of the first peak and the first valley of the pair correlation function.^58–60 Based on these rules, we quantify the low-density chain formation via the parameter

ϕ_p = 〈N_p〉/N, (9)

where N_p is the total number of particles which reside in chains. In the context of aggregating molecular systems, ϕ_p is often called the degree of polymerization.

At high densities (e.g., Φ = 0.58), the particles are obviously closer to each other, and the string-like structures are no longer isolated. This causes the degree of polymerization ϕ_p to no longer be appropriate to characterize the self-assembly. To overcome this difficulty, we propose a different order parameter to describe the chain formation quantitatively. The starting point is the dipole–dipole correlation function

where the longitudinal and transverse displacement of particle j relative to i is r^‖_ij = r_ij·[small mu, Greek, circumflex]_i and , respectively. To extract the angular dependence at short distances, we transform the Cartesian coordinates employed in eqn (10) to polar coordinates by using and θ = atan2(−r_⊥,r_‖), and compute the function

Here, the distance r_s is set to r_s = L/8 ≈ 5σ such that r_s is much smaller than the half of the box size (which corresponds to the length scale of the periodic boundary conditions). In the presence of chain-like structures with head-to-tail alignment of neighboring particles, we expect the angular correlation function, g̃_μ(θ), to display positive maxima both, in front of and behind the reference particle (i.e., at θ = 0 or π). In contrast, the dipole moments of the particles on the right- and left-hand side of the reference particle should remain rather uncorrelated (i.e., g̃_μ(θ) ≈ 0 for θ = π/2 or −π/2). With this picture in mind, we measure the quantity

Z = g̃_μ(θ_max) − g̃_μ(θ_min), (12)

where g̃_μ(θ) reaches its maximum and minimum, respectively, at θ = θ_max and θ = θ_min. By properly choosing a threshold value Z_thres, we define the dipolar particles as exhibiting chain-like structure when Z ≥ Z_thres. As will be discussed in more detail in Section 2.3.3, a reasonable value for Z_thres at Φ = 0.58 is Z_thres = 0.17.

3 Simulation results

Based on the target quantities described in Section 2.3, we can classify the states observed in our BD simulations performed at three densities and at various values of the motility v₀* and the dipolar coupling strength λ. In the following Section 3.1–3.3 we discuss, for each density, the state diagram in the (v₀*, λ) plane. Specifically, we identify five states whose characteristics are summarized in Table 1.

Table 1 Characterization of the states of dipolar active particles according to the order parameters defined in Section 2.3

State	Clustering	Orientational ordering	Chain formation
			Φ ≲ 0.23	Φ = 0.58
Homogeneous, isotropic fluids	ϕ c ≤ 0.5	ϕ e ≤ 0.5	ϕ p ≤ 0.5	Z ≤ 0.17
Chain-like structures	ϕ c ≤ 0.5	ϕ e ≤ 0.5	ϕ p > 0.5	Z > 0.17
Micro-flocking	ϕ c ≤ 0.5	ϕ e > 0.5	—
Macro-flocking	ϕ c > 0.5	ϕ e > 0.5	—
Motility-induced clustering	ϕ c > 0.5	ϕ e ≤ 0.5	—

---

3.1 The low density regime (Φ = 0.12)

We start by considering the state diagram in the low-density regime, taking Φ = 0.12 as a representative example, see Fig. 1. In the absence of dipolar coupling (λ = 0), non-dipolar active particles at such a low density typically display a homogeneous isotropic fluid state, without significant translational correlations or orientational ordering.³ As long as the dipolar coupling strength remains relatively small (λ ≲ 4), the system still shows the same homogeneous isotropic fluid state for all motilities considered. In contrast, when λ exceeds a value of about 6, the dipolar active particles form chain-like structures provided that the motility is below a critical motility v_0,c*(λ), whose value depends on λ. Above this critical motility v_0,c*(λ), the chain-like structures found at smaller motilities break, and the system becomes homogeneous and isotropic. As Fig. 1 reveals, v_0,c*(λ) increases with increasing λ. As a visualization of the impact of motility on dipolar active particles at a high coupling strength, such as λ = 10, we plot in Fig. 2 representative snapshots. At v₀* = 0, nearly all particles are bound into chains and rings, as it is typical for passive dipolar particles in 2D.²⁷ With increasing v₀* the ring structures are seemingly more abundant, as can be seen in Fig. 2(a–c). Further increase in v₀* eventually causes the chain-like structures to break into fragments of short chains and single particles, as shown in Fig. 2(d).

Fig. 1 State diagram of dipolar active particles in the (v₀*, λ) plane at Φ = 0.12. The points on the diagram indicate the parameter combinations used in the simulations. At Φ = 0.12, we have observed homogeneous, isotropic fluid states (black dots) and chain-like structures (blue crosses).

Fig. 2 Representative simulation snapshots at Φ = 0.12. Particles are colored according to their dipole orientations as indicated by the color ring in the inset.

To quantitatively characterize the motility-induced destruction of chain-like structures, we plot in Fig. 3 the degree of polymerization ϕ_p as a function of the motility v₀* for various coupling strengths λ. In the absence of dipole–dipole interactions (λ = 0), particles do not self-assembly into chain-like structures, still we observe a slight increase in ϕ_p [see eqn (9)] upon increasing motilities v₀*. This can be attributed to the dynamical clustering of active Brownian particles into small (finite-sized) aggregates.^3,36,61 Indeed, if the particles inside these finite-sized clusters fulfill our criteria of being bonded, the order parameter ϕ_p will be non-zero (yet small), despite the fact that dipole–dipole interactions are absent. Once the dipole–dipole interactions are introduced, pronounced chain-like structures appear for coupling strengths λ ≳ 6 and low motilities v₀*, giving rise to large values of ϕ_p (ϕ_p ≈ 1). Upon an increase in the motility, the degree of polymerization ϕ_p gradually decreases, reflecting that the self-propulsion opposes the formation of chain-like structures. This behavior resembles that seen in passive, dilute systems of self-assembled dipolar particles, when the temperature is increased.⁶² In this sense, the particle motility in the active dipolar system may be viewed as an analog to the temperature in the passive system.

Fig. 3 Degree of polymerization ϕ_p as a function of the motility v₀* at Φ = 0.12 for the coupling strength λ = 0 (black circles), 4 (red squares), 6 (green diamonds), and 10 (blue triangles). The dashed horizontal line at ϕ_p = 0.5 marks our criterion for a state with chain-like structures (see Table 1). The solid lines are guides to the eye.

The degree of polymerization ϕ_p only describes the chaining behavior globally. To supplement our analysis of chain formation, we measure the chain size distribution function P(n), which represents the probability that a randomly selected chain consists of n particles. In the presence of chain formation, the simulated systems usually contain several long chains and many short chains, such that P(n) becomes very small at large n. Therefore, we plot in Fig. 4 the chain size distribution weighted by n, where nP(n) is proportional to the probability that a randomly selected particle belongs to a chain with a size n. As seen in Fig. 4(a), the weighted distribution curves of non-dipolar particles vanish at around a chain size n ≲ 10, and the values slightly increase with increasing motilities. This is consistent with the aforementioned dynamical clustering. The corresponding functions at large dipolar coupling strengths (λ = 10) look very different, as shown in Fig. 4(b). Here, we observe a broad peak between 10 ≲ n ≲ 100 at v₀* = 0, reflecting formation of long chains. Upon an increase in the motility, this peak is gradually shifted to a smaller chain size. Finally, this peak disappears once v₀* ≳ 60. This suggests a vanishing of chains with the most probable size at large values of v₀*.

Fig. 4 Weighted distribution of the chain size at Φ = 0.12 for motilities v₀* = 0 (black circles), 20 (red squares), 40 (green diamonds), 60 (blue triangles up), 80 (orange triangles down), and 100 (brown crosses) with the coupling strength λ = 0 (a) and λ = 10 (b).

Inspecting again Fig. 2(d) we see that even at the largest coupling strength considered, the dipole moments do not align on a length scale comparable to that of the simulation box, regardless of the values of v₀*. In other words, there is no pronounced global orientational order. Indeed, through measuring the orientational order parameter ϕ_e, we confirm that ϕ_e remains small (ϕ_e < 0.5) for the parameter combinations explored in Fig. 1, indicating that the systems are globally isotropic at low densities.

3.2 The intermediate density regime (Φ = 0.23)

At the density Φ = 0.23, the state diagram (see Fig. 5) at small and intermediate motilities (v₀* ≲ 20) is very similar to that at Φ = 0.12 (see Fig. 1): For weakly coupled systems (λ ≲ 3), we observe a homogeneous and isotropic fluid state, whereas strong dipolar coupling leads to chain-like structures. This is visualized by the simulation snapshots presented in Fig. 6(a) and (c). For these “nearly passive” systems, the main difference compared to the dilute systems in Fig. 1 is that the change from the homogeneous isotropic fluid state into the chain-like state occurs at a somewhat smaller value of λ. This is plausible, because the higher density leads to a smaller mean separation between the particles and thus, to a higher probability for chain formation.

Fig. 5 State diagram of dipolar active particles in the (v₀*, λ) plane at Φ = 0.23. The points on the diagram indicate the parameter combinations used in the simulations. At Φ = 0.23, we have observed homogeneous, isotropic fluid states (black dots), micro-flocking (green triangles), and chain-like structures (blue crosses).

Fig. 6 Representative simulation snapshots at Φ = 0.23. Particles are colored according to their orientations as indicated by the color ring in the inset.

We now turn to the impact of the motility. At small values of λ, an increase in v₀* toward high values leads to the formation of finite-sized clusters [see Fig. 6(b)] which are, however, not stable. Thus, the system remains homogeneous and isotropic on average. This finding is consistent with earlier research on active Brownian particles. In particular, the phenomenon of motility-induced phase separation, accompanied by formation of “giant” clusters (ϕ_c ≈ 1), is known to occur only at higher densities.³ As Fig. 5 reveals, the picture changes if the dipolar coupling becomes larger. Consider, e.g., the case λ = 6. Increasing the motility from zero, the chain-like structures observed at small v₀* break [see Fig. 6(c)], and the particles instead form small clusters characterized by the same orientation [see Fig. 6(d)]. More quantitatively, the clustering parameter ϕ_c remains small (ϕ_c < 0.5) but the orientational order parameter ϕ_e reaches values above 0.5.

The full behavior of the orientational order parameter ϕ_e as function of v₀* is shown in Fig. 7. For strongly coupled systems (λ ≥ 3), the order parameter abruptly changes from (essentially) zero to large values at a “critical” motility v_0,c*(λ) ≈ 20. We take this as an indication for a motility-induced formation of ordered, yet small, clusters. According to Table 1, this is characteristic of a “micro-flocking” state. We recall that, for active particles, the development of average alignment between neighboring particles implies that they move along the same direction. To explore whether these small clusters are indeed “microscopic” structures, we measured the cluster size distribution. It turns out that the characteristic cluster size n₀ does not scale up with the system size N (see the detailed analysis in Appendix A.2).

Fig. 7 Global polarization ϕ_e as a function of the motility v₀* at Φ = 0.23 for the coupling strength λ = 0 (black circles), 2 (red squares), 3 (green diamonds), 4 (blue triangles up), and 6 (orange triangles down). The dashed line indicates the value corresponding to the emergence of flocking states (see Table 1). The solid lines are guides to the eye.

Inspecting Fig. 7 again, it is found that for systems at an intermediate coupling strength λ = 2, the order parameter ϕ_e gradually increases from zero as v₀* ≳ 100, suggesting that the “micro-flocking” state might still appear. However, further increase in v₀* requires to employ even smaller simulation time steps Δt < 2 × 10⁻⁵τ, so as to prevent the numerical instability. Combining this aspect with the time-consuming Ewald summation, it turns out that the simulations become computationally unfeasible for investigating the flocking state transition at small dipolar coupling strengths (λ < 3).

A behavior similar to the “micro-flocking” state has been previously observed in systems of active Brownian disks with polar alignment.⁹ Indeed, the density Φ = 0.23 considered here is close to the effective packing fraction ϕ* = 0.256 defined in ref. 9, where a comparable state with ‘microscopic’ polar clusters has been discovered at high motilities and intermediate polar coupling. Upon further increase in polar coupling, the system in ref. 9 displays macroscopic structures, such as moving patterns of bands or lanes. Here, we did not find such patterns in the parameter range explored. In this context, we also mention a state observed in the Vicsek model at low densities and low noise (yielding strong alignment).⁸

In this state, the point-like particles form small groups within which the particles move together along one (random) direction. However, the system does not show global ordering. In both of these models, i.e., active polar disks and the Vicsek model, the particles favor alignment of the propulsion direction with their neighbors, independent of their relative positions. This is different from the particles considered here whose interaction depends not only on the relative alignment but also on their configuration.

In order to gain a deeper insight into the flocking behavior of the present system, we draw in Fig. 8(a–d) sketches of the motion of two highly motile, strongly coupled particles starting from four representative configurations. The initial configurations (related to time t₀) in Fig. 8(a) and (b) are those with the lowest dipolar interaction energy. These configurations are therefore stable (up to thermal fluctuations) in the passive case. Upon switching on the motility (v₀* > 0), the particles in the head-to-tail configuration [see Fig. 8(a)] will move along a straight line in the same direction (as already argued in ref. 13), until thermal fluctuations set in. We thus consider this configuration as relatively stable not only in the passive, but also in the active case. This is different in the antiparallel side-by-side configuration depicted in Fig. 8(b) at t₀. Here, activity leads to motion in opposite direction, yielding this configuration unstable. In Fig. 8(c), the initial configuration at t₀ is energetically unstable. Still, at least for short time, one expects the particles to move together along the same direction until the dipolar repulsion drives them away from one another. Finally, in Fig. 8(d) the initial configuration is a slightly distorted (induced by rotational noise) head-to-head configuration. This leads to a strong dipolar torque which, combined with strong activity, pushes the particles apart. Combining these arguments, we conclude that the most stable configuration for two dipolar active particles is the head-to-tail alignment. While this is similar to the passive case, a major difference is that the antiparallel side-by-side configuration [Fig. 8(b)] breaks immediately apart when the particles become active, yielding this configuration unstable. At low densities, we therefore expect to find short straight chains moving along their individual long axis, consistent with the observation in the Fig. 2(d). Upon an increase in the density, short chains collide with each other more frequently than at low densities. As explained above, antiparallel alignment of the short chains is not stable. Therefore, once two short chains collide, they tend to align along the same direction and, thus, form a dynamical polar cluster, as illustrated in Fig. 8(e). At a sufficiently high density (e.g., Φ ≳ 0.23), the frequent collisions between short chains finally lead to a “micro-flocking” state, in which the dipolar active particles form finite-sized clusters with a polar order [see Fig. 6(d)]. This is fundamentally different from the alignment mechanism in other models such as the ones with ferromagnetic-like interactions, where the configuration at t₀ in Fig. 8(c) is kinetically stable.^9,10

Fig. 8 (a–d) Sketches of the motion of two dipolar active particles in three time steps (t₀ < t₁ < t₂) starting from different representative configurations. The dashed lines indicate the trajectories of the particles. (e) A scheme of a collision event of two short chains, consisting of aligned dipolar active particles.

3.3 The high density regime (Φ = 0.58)

At high densities such as Φ = 0.58, there are three major phenomena interfering with one another: motility-induced phase separation, polar ordering, and chain formation. In the following Sections 3.3.1–3.3.3, we discuss these issues in detail.

3.3.1 Motility-induced clustering. Fig. 9 shows the state diagram at Φ = 0.58 in the plane spanned by motilities v₀* and dipolar coupling strengths λ. We start by investigating the regime of small λ. Specifically, we are interested in the impact of dipolar interactions on the motility-induced phase separation known from conventional (non-dipolar) active particles. Indeed, as the motility v₀* increases from zero, non-dipolar active Brownian particles (λ = 0) undergo a transition from a homogeneous, isotropic fluid state into a state with large, dense clusters coexisting with freely moving colloids in the dilute region, as shown in Fig. 9, 10(a) and (b). This behavior, generally known as motility-induced phase separation, occurs even if attractive interactions are absent in the model.^39,42

Fig. 9 State diagram of dipolar active particles in the (v₀*, λ) plane at Φ = 0.58. The points on the diagram indicate the parameter combinations used in the simulations. At Φ = 0.58, we have observed homogeneous, isotropic fluid states (black dots), motility-induced clustering (red squares), macro-flocking (orange diamonds), micro-flocking (green triangles), and chain-like structures (blue crosses). The region surrounded by the dashed lines indicates a parameter regime where the simulations did not reach a steady state.

Fig. 10 Representative simulation snapshots at Φ = 0.58. Particles are colored according to their orientations as indicated by the color ring in the inset.

The motility-induced phase separation can be explained by a self-trapping mechanism: When a highly motile particle enters into a dense region occupied with other particles, this particle is temporally slowed down due to frequent collisions. Such a slowing-down effect makes this region even denser, thus creating a positive feedback loop which leads to the formation of giant clusters.³⁹ We note that the particle orientations in the two phases are essentially uncorrelated for non-dipolar particles (λ = 0), as seen in Fig. 10(b).

To get a first impression of the impact of dipole–dipole interactions on the clustering behavior, we plot in Fig. 11 the fraction of the largest cluster, ϕ_c [see eqn (7)], as a function of the motility v₀* for various coupling strengths λ. Consistent with the motility-induced phase separation, the curves for non-dipolar and weakly coupled active particles (λ = 0–1) display a sharp increase from zero at v₀* ≈ 20 and reach large values (ϕ_c ≈ 0.8) at v₀* = 100. The curves reach ϕ_c = 0.5 at larger v₀* with increasing λ, which is also reflected by the boundary between black dots and red squares in Fig. 9. For strong dipolar coupling (λ ≳ 2), the values of ϕ_c at high motilities remain in the same range (ϕ_c ≈ 0.8, indicating again formation of large clusters), but the increase in ϕ_c at small v₀* is much less pronounced. This already indicates a strong impact of dipolar interactions. For more detailed analysis of the clustering behavior, we discuss the cluster size distribution in Appendix A.2.

Fig. 11 Fraction of the largest cluster ϕ_c as a function of the motility v₀* for the coupling strength λ = 0 (black circles), 0.5 (red squares), 1 (green diamonds), 2 (blue triangles up), and 5 (orange triangles down). The dashed line ϕ_c = 0.5 marks the value corresponding to cluster formation.

To gain further information, we calculate position-resolved local area fractions, ϕ, based on a Voronoi tessellation.^36,40 (We note that, contrary to ref. 36, we did not perform a short time average of ϕ, since the polar clusters characterizing the dense state migrate over time and, therefore, are not stationary within the short-time interval. The flocking behavior will be discussed later in detail in Section 3.3.2.) Fig. 12 shows the probability distribution of the local area fractions, P(ϕ), for different coupling strengths at v₀* = 100. For non-dipolar active particles (λ = 0), P(ϕ) reveals a clear double-peak structure that reflects the separation between the dilute and the dense phase. The coexisting densities correspond to the location of the two local maxima of P(ϕ). Upon an increase in the coupling strength up to λ = 2, the double-peak structure gradually disappears. Instead, we observe the emergence of a single peak located at a density slightly larger than the mean density (Φ = 0.58), as well as a broad shoulder on the left. This suggests the disappearance of the phase separation observed at zero and small λ.

Fig. 12 Probability distribution function of local area fractions, P(ϕ), for the coupling strength λ = 0 (black circles), 0.5 (red squares), 0.75 (green diamonds), 1 (blue triangles up), and 2 (orange triangles down) with the motility v₀* = 100.

The same conclusion can be drawn from Fig. 13, where we plot the coexisting densities in the (v₀*, ϕ) plane for various λ. The location of the high-density branch at λ = 0 can be explained as follows: At high motilities v₀*, the non-dipolar system (λ = 0) exhibits “giant” clusters, composed of randomly oriented particles [see Fig. 10(b)]. These particles are separated by a distance close to the effective hard sphere diameter σ_eff defined in Section 2.2. As a result, the local area fraction in the densely packed region should be close to the close-packing fraction, . Upon a slight increase in dipolar coupling, the coexistence branches for the dilute region, ϕ_gas, are shifted toward higher densities, while the branches for the dense region, ϕ_den, move toward lower area fractions. As a result, the area surrounded by the curves of ϕ_den and ϕ_gas in Fig. 13 significantly shrinks with increasing λ, indicating that motility-induced phase separation is generally suppressed by the dipolar interactions. Within the present simulations, the phase separation and the corresponding coexistence curves disappear once λ > 1. We note that this observation is in contrast with the findings in a recent study of active Brownian particles with additional interactions of velocity alignment at high densities.¹⁰ This shows that different types of orientational interactions may have entirely different impacts on motility-induced phase separation.

Fig. 13 Coexisting densities in the (v₀*,ϕ) plane for the coupling strength λ = 0 (black circles), 0.5 (red squares), 0.75 (green diamonds), and 1 (blue triangle up). The black dashed line displays the effective close-packing fraction, .

Yet another perspective on the disappearance of the phase separation due to dipolar interactions emerges when we consider the self-trapping mechanism (which plays a key role for clustering at λ = 0). To this end, we compute the normalized speed of the ith particle, v_i*/v₀*, versus the local density ϕ_i = ϕ(r_i), where the speed is given by v_i* = |(Δr_i/σ)/(Δt_s/τ)| with Δt_s = 10⁻²τ. After averaging over N particles and over at least 1000 snapshots, we plot in Fig. 14 the normalized particle speed versus the local density. For non-coupled (λ = 0) and weakly coupled (0 < λ ≲ 1) systems, we observe a linear decay,

v*(ϕ)/v₀* = 1 − a_dϕ, (13)

consistent with the prediction from a phenomenological approach:^63,64 Particles move slower when traveling through a crowded area. The fitting parameter, a_d, represents the decay amplitude, and, therefore, describes the significance of the self-trapping mechanism. Upon increasing λ in the range λ ≲ 1, the decay becomes less pronounced. Consequently, the fitting parameter, a_d, monotonically decreases, as seen in the inset of Fig. 14. In contrast, for λ ≳ 2, the normalized speed as a function of ϕ is essentially constant, that is, a_d tends to zero. In other words, there is no self-trapping anymore. This may be attributed to the fact that the strongly-coupled dipolar particles tend to form polar clusters with local head-to-tail alignment, causing the particles to move along the same direction [see Fig. 10(d)]. As a result, the motility-induced phase separation is replaced by a macro-flocking state, as shown in Fig. 9. We will come back to the flocking behavior in Section 3.3.2.

Fig. 14 Normalized particle speed v*(ϕ)/v₀* as a function of the local area fraction ϕ for the coupling strength λ = 0 (black circles), 0.5 (red squares), 1 (green diamonds), 2 (blue triangles up), and 5 (orange triangles down) with the motility v₀* = 100. The solid lines represent fits to eqn (13). Inset: Fitted decay amplitude a_d as a function of λ. The line is drawn as a guide to the eye.

To obtain a more complete (yet qualitative) picture on the role of the function a_d(λ) for the appearance of the motility-induced phase separation, we consider the effective free energy proposed in ref. 41,

f(ϕ) = f₀(ϕ) + f_rep(ϕ), (14)

where the bulk contribution is given bywith density-dependent swim speed v(ρ) = 1 − a_dρ.^63,64 Further, the contribution from the excluded volume interactions between particles is written as

f_rep(ϕ) = k_repΘ(ϕ − ϕ_t)(ϕ − ϕ_t)⁴ (16)

with Θ(x) being the Heaviside step function, k_rep the repulsive strength, and ϕ_t the threshold area fraction.

For conventional active Brownian particles (λ = 0), the decay amplitude, as a crucial parameter describing the degree of self-trapping, is typically chosen to be a_d = 1 for 2D and a_d = 1.3 for 3D.⁶⁴ With this in mind, in Fig. 15 we plot the effective free energy f(ϕ) for the decay amplitude a_d = 1.08, 0.95, and 0.8. At low densities, f(ϕ) is nearly independent of a_d, while at high densities, f(ϕ) increases with decreasing a_d. From f(ϕ), we determine the coexisting densities through the common tangent construction. To compare the influence of the decay amplitude, a_d, on the phase separation from the theoretical perspective and the simulation results, we plot in the inset of Fig. 15 the coexisting densities in the (a_d,ϕ) plane for the effective free energy (blue circles) and simulations at v₀* = 100 and 0 ≤ λ ≤ 1 (orange squares). More specifically, for the above simulations we plot the coexisting densities determined from Fig. 12versus the associated a_d plotted in the inset of Fig. 14. The coexisting densities obtained from both methods describe the same trend, that is, a decrease in the decay amplitude a_d causes a reduction in ϕ_den and a growth in ϕ_gas. Further, the phase separation disappears once v(ϕ) decays sufficiently slow (a_d ≲ 0.6 for the effective free energy, and a_d ≲ 0.5 for simulations). This clearly indicates that the motility-induced phase separation is suppressed when the function v(ϕ) decreases too slowly, i.e., when the dipolar coupling becomes too strong.

Fig. 15 Effective free energy f(ϕ) with k_rep = 5000 and ϕ_t = 0.84 for various a_d = 1.08 (solid black line), 0.95 (solid red line), and 0.8 (solid green line). The dashed lines display the common tangent construction. Terms linear in ϕ are irrelevant for the common tangent construction and have been subtracted for clarity. Inset: Dependence of coexisting densities on a_d obtained from the effective free energy (blue circles) and simulations at v₀* = 100 and 0 ≤ λ ≤ 1 (orange squares).

3.3.2 Flocking. We now come back to the emergence of flocking and polar ordering (see green triangles and orange diamonds in Fig. 9). As argued in Section 3.2, once the density is sufficiently high (Φ ≳ 0.23), the interplay of dipolar interactions and activity allows for the formation of polar clusters. To quantify this behavior at high densities, Fig. 16 shows the magnitude of the average orientation, ϕ_e, as a function of λ for various values of v₀*. In the passive case (v₀* = 0), the order parameter remains small for all λ considered. Thus, there is no clear hint regarding whether global order appears, which is consistent with earlier simulations of monolayers of dipolar particles.^48,49,65 From the simulations, it is known that dense two-dimensional passive systems of dipolar particles tend to develop domain-like structures, which are highly frustrated and characterized by the polar order only locally. Coming back to the active case, we see from Fig. 16 that already at an intermediate value of the motility (v₀* = 20), the order parameter reaches significant values (ϕ_e ≳ 0.5) when λ exceeds a value of about two. Combining this finding with the cluster analysis shown in Fig. 11 and using the criteria in Table 1, we classify this behavior as a micro-flocking state (3 ≲ λ ≲ 6, 10 ≲ v₀* ≲ 40), see green triangles in Fig. 9. Finally, in the regime of high motilities (v₀* = 60–100), the data curve ϕ_e(λ) is reminiscent of a (polar) phase transition: ϕ_e rises suddenly from zero to values greater than 0.5 at a “critical” coupling strength λ_c ≈ 0.75. This critical value slightly decreases upon an increase in v₀*. In addition to the large values of ϕ_e, the fraction of the largest cluster, ϕ_c, is greater than 0.5 for v₀* ≳ 60, as shown in Fig. 11. Thus the systems at v₀* ≳ 60 and λ ≳ 1.25 are in a macro-flocking state. We note that the existence of the polar order is also reflected by the positive time correlations of individual dipole moments at a time interval much longer than the Brownian diffusion time τ (not shown here).

Fig. 16 Magnitude of the average orientation ϕ_e as a function of the coupling strength λ at Φ = 0.58 for v₀* = 0 (black circles), 20 (red squares), 60 (green diamonds), 80 (blue triangles up), and 100 (orange triangles down). The dashed line ϕ_e = 0.5 indicates the value above which the particles form a flocking state (according to Table 1). The solid lines are drawn as a guide to the eye.

3.3.3 Chain-like structures. Having discussed the emergence and interplay of phase separation and global polar order, we finally consider the fate of the chain-like structures characterizing passive, dense dipolar system upon increasing v₀* from zero. The question is how the system transforms from a state with chain-like structures into a macro-flocking state [see Fig. 9, 10(c), and (d)]. As discussed in Section 2.3.3, the order parameter ϕ_p (degree of polymerization) is no longer adequate to quantify such a transition at high densities. Therefore, an alternative method is proposed: Fig. 17(a–d) shows the spatial correlation function of dipole moments, g_μ(r_⊥,r_‖) [defined in eqn (10)], for v₀* = 0–100 and λ = 5. For the passive and nearly passive case [see Fig. 17(a and b)], the orange (curved) strips clearly indicate positive correlations in front of and behind the reference particle, while the purple regions suggest that the dipole moments in the equatorial zone are rather uncorrelated. This observation reflects chain formation along the direction of the reference dipole moment. As the motility v₀* increases [see Fig. 17(c and d)], the purple regions disappear. Instead, the whole figure turns orange, indicating positive correlations in all directions. This shows a rather uniform alignment of dipole moments regardless of the relative positions between particles. In other words, the chains disappear.

To better evaluate the dependence of orientational correlations on the direction of the connecting vector r_ij, Fig. 17(e) shows the angular correlation function, g̃_μ(θ) [defined in eqn (11)], for the strong dipolar coupling strength λ = 5. As mentioned in Section 2.3.3, the difference between the maximum g̃_μ(θ_max) and the minimum g̃_μ(θ_min) of the curve, Z [defined in eqn (12)], can be used as a measure of chain formation. Upon increasing v₀*, Z decreases and, at the same time, the curves in Fig. 17(e) are shifted toward lager positive correlations. An overview of the behavior of Z as a function of v₀* (for various λ) is given in Fig. 18. The data curves bear close resemblance to those for ϕ_p in Fig. 3. For a non-dipolar system (λ = 0), where chains are absent, the order parameter Z does not show any significant behavior within the range of motilities explored. In contrast, at large coupling strengths (λ = 4 and λ = 5), the decrease in Z with increasing v₀* reflects the disappearance of a state with chain-like structures. Specifically, the curves for λ ≥ 4 show a sharp decrease, accompanied by a point of inflection Z ≈ 0.17. We choose this as a threshold value. We note that the simulations of nearly passive, strongly-coupled systems (e.g., λ = 5 and 6) are severely plagued by large fluctuations of the order parameter Z over time. In these cases, we were unable to determine whether or not the simulations had reached a steady state even after a very long simulation time (t > 300τ). The corresponding region in the state diagram in Fig. 9 is marked by dashed lines.

Fig. 17 Spatial correlation function of dipole moments, g_μ(r_⊥, r_‖) defined in eqn (10), at Φ = 0.58 and λ = 5 for the motility v₀* = 0 (a), 8 (b), 16 (c), and 100 (d). The black circle represents the reference particle with the red arrow indicating the particle's orientation as well as the dipole moment. The area is drawn in white color to reflect that the center-to-center distance between particles can not be smaller than its diameter due to steric repulsion. (e) The angular correlation function of dipole moments, g̃_μ(θ) defined in eqn (11), as a function of θ = atan2(−r_⊥, r_‖) with the coupling strength λ = 5. The graph at the bottom-right corner of (e) illustrates the definition of θ.

Fig. 18 Order parameter Z as a function of the motility v₀* at Φ = 0.58 for the coupling strength λ = 0 (black circles), 4 (red squares), 5 (green diamonds), and 6 (blue triangles). The dashed line indicates the threshold value for the state transition Z_thres = 0.17. The solid lines are drawn as a guide to the eye.

4 Conclusions

Using Brownian dynamics simulations, we have studied how dipolar interactions and self-propulsion combine to influence the dynamical self-assembly of a monolayer of dipolar active Brownian particles. To this end, we have presented state diagrams in the plane spanned by the dipolar coupling and the motility for three representative densities.

When the motility is small, the state diagrams are similar for all densities considered. Specifically, homogeneous and isotropic fluids are observed for nearly passive particles with weak dipolar coupling, whereas strong dipolar coupling leads to chain-like structures. At high motilities, the state diagrams strongly depend on the mean area fraction Φ. At low densities (Φ = 0.12) and strong dipolar coupling, an increase in the motility v₀* from zero causes chain-like structures to break into fragments of short chains and individual beads. At intermediate densities, passive, strongly-coupled dipolar particles self-assemble into a state with chain-like structures. With increasing v₀* from zero, chains start to break and the system displays a micro-flocking state, where particles form finite-sized clusters with polar order. Finally, at a high densities and strong dipolar interactions, we observe a motility-induced transition from a state with chain-like structures into a micro-flocking state and finally into a macro-flocking state, where particles show global orientational ordering and form “giant” clusters.

To provide a simple argument for the emergence of polar order (which is absent in the passive 2D case), we have considered the time evolution of two dipolar particles starting from four representative initial configurations. As a result of the interplay between dipolar coupling and self-propulsion, the head-to-tail configuration remains the most stable one (same as the passive case), while the antiparallel side-by-side configuration is destabilized. With increasing particle number, the head-to-tail alignment mechanism can finally lead to a flocking state if the density is sufficiently high (Φ ≳ 0.23).

It has been reported earlier that finite size effects may hide crucial features of the flocking behavior in active systems, such as the first order nature of the flocking transition followed by the formation of traveling bands.⁶ We note that the sizes of the present simulations are limited to the order of 10³ particles, due to the long-range character of the dipolar interactions and, subsequently, the expensive computational cost. Besides the qualitative behavior, it would be very interesting to study the scaling behavior and explore whether the dipolar active particles belong to any of the existing universality classes, such as the one of the Vicsek model.^4,66,67 To this end, however, it would be necessary to perform extensive simulations with system sizes much larger than N ≈ 10³, which is, again, limited by the computational resources. Nevertheless, it is worth mentioning a fundamental difference regarding the flocking transition: in the Vicsek model the onset of the flocking state is characterized by simultaneous spontaneous appearance of density and orientational inhomogeneities.^6–8 In contrast, in our model the emergence of global orientational order is decoupled from the formation of large-scale structures, such as “giant” clusters.

At Φ = 0.58 and large values of v₀*, the system displays motility-induced phase separation, where “giant” clusters composed of densely-packed particles with random orientations coexist with freely moving particles in the dilute region. The phase separation persists as long as the dipolar coupling λ is negligible against thermal fluctuations (0 < λ ≲ 1). Once the dipolar interactions dominate (λ > 1), the orientations of the active particles are no longer uncorrelated and the particles tend to align with their neighbors, leading to a break down of the self-trapping mechanism and a subsequent suppression of the phase separation.

We note that our model does not account for hydrodynamic interactions between the particles. In the absence of dipolar coupling, simulation studies have reported that hydrodynamic interactions tend to suppress motility-induced phase separation due to, either, the near-field interactions,^68,69 or the rapid decorrelation of the particle orientations.^70–72 Through investigating the pair distribution function, it has been shown that hydrodynamic interactions generally damp out the translational structure of active particles at high densities.⁷³ Moreover, hydrodynamically interacting particles with a certain range of force dipole strengths can spontaneously form a state with global polar order,^68,69,74,75 which can be attributed to either the actively induced rotation-translation coupling,⁷⁶ or the near-field lubrication forces.^68,69 Therefore, one may expect that dipolar coupling and hydrodynamic interactions combine to further suppress the phase separation and promote the polar ordering. The details of the resulting collective behavior are a topic of future studies.

In the real world, active particles are often asymmetric.^77,78 As a result, the effective propulsion force does not coincide with the particle's center of mass, thus generating a propulsion torque which induces chiral active motion.^35,79 Such a mechanism is also expected for the magnetic or dielectric Janus particles, whose dipole moments (either permanent or induced) are mostly shifted from the center of the individual particles.^80–83 Therefore, it would be very interesting to investigate models accounting for the chiral motion and the shifted dipole moment, so as to obtain more comprehensive understanding of dipolar active systems.

Conflicts of interest

There are no conflicts to declare.

A Appendix

A.1 Ewald summation

To deal with the long-range character of the dipole–dipole interactions, we employ 2D Ewald summation with the “tinfoil” boundary condition.⁸⁴ Within this method, the total dipole–dipole energy is separated into different contributions,The first term on the right-hand side of eqn (17) corresponds to the real-space contribution, where the functions B(r,α) and C(r,α) are defined by³²andwith erfc(x) being the complementary error function. In eqn (17), the real-space contribution is formulated under the assumption that the convergence parameter α is large enough, such that we can consider only the interactions within the central simulation box. This is achieved by choosing α = 7/L (with L being the box size) and evaluating r_ij with the minimum-image convention.³² The second term in eqn (17) represents the reciprocal-space contribution withThe wave vectors in the reciprocal (square) lattice are given by k = (2π/L)m, where m = (m_x,m_y)^T with m_x and m_y being integers. The magnitude of the wave vector is denoted as . We evaluate the reciprocal-space summation for the k-vectors within the range m² = m_x² + m_y² ≤ 15².⁸⁵ The third term represents the correction due to unphysical self-interaction of dipole moments.⁸⁵

The Langevin eqn (4) and (5) involve the forces and torques due to the dipolar interactions. From eqn (17), the force acting on the ith particle due to dipole–dipole interactions is given by

F_i,dd = −∇_riU_dd = F^R_i,dd + F^k≠0_i,dd, (21)

where the real-space contribution isThe function C(r,α) is defined in eqn (19), and D(r,α) is given byFurther, the contribution from the Fourier-form is written aswhere and [see eqn (20)].

In the present work, every dipole moment is restricted to point in the xy-plane with ψ_i characterizing its direction relative to the x-axis. Therefore, the torque acting on the ith particle is always directed along the z-axis. It is given by

T_i,dd = −∂_ψiU_dd ẑ = T^R_i + T^k≠0_i, (25)

where U_dd is given in eqn (17), and ẑ denotes the unit vector along the positive z-axis. The real-space component of the torque is given by.Finally, the reciprocal-space contribution is given by

A.2 The cluster size distribution

The purpose of this paragraph is to show that the clusters observed in the “micro-flocking” state at intermediate densities are indeed microscopic structures. To this end, we performed a finite size analysis of the cluster size distribution. If these clusters are microscopic patterns, they should not grow as the total number of particles N (i.e., the system size) becomes larger. Fig. 19 shows the weighted cluster size distribution nP(n) at Φ = 0.23, v₀* = 100 and λ = 6 for various system sizes N, where the cluster size n represents the number of particles within a cluster. For N < 1000, the distribution decays faster as N decreases, indicating that the systems are influenced by the finite sizes. However, as we increase the particle number up to N ≳ 1000, the data points collapse onto a single curve. To obtain a more quantitative description, we fit the data points for each system size via the function,

nP(n) = P(1)n⁻¹e^{−(n−1)/n₀}, (28)

where n₀ denotes the characteristic cluster size. Similar fitting functions of the cluster size distribution have also been considered in systems of conventional active Brownian particles⁸⁶ and polar active disks.⁹ The inset of Fig. 19 shows that n₀ reaches a plateau once N ≳ 1000. This suggests that the cluster size indeed remains constant as the system size N increases. By this, we can confirm that the clusters observed in the “micro-flocking” state at Φ = 0.23 are microscopic structures.

Fig. 19 Weighted distribution of cluster size at Φ = 0.23 for the particle number N = 100 (black circles), 225 (red squares), 400 (green diamonds), 1156 (blue triangles up), and 2500 (orange triangles down). The orange solid line indicates the curve fitted to N = 2500 according to eqn (28). Inset: Characteristic size of clusters, n₀, as a function of particle number N.

To provide additional information on the clustering behavior at high densities (Φ = 0.58), we plot in Fig. 20 the weighted cluster size distribution for non-dipolar and strong coupled dipolar active particles at Φ = 0.58. At the coupling strength λ = 0, our model reduces to the limiting case of active Brownian particles.^3,63 Here, the correspondent weighted cluster size distribution function undergoes a transition from an exponential decay (v₀* ≲ 20) into a curve with a power law decay at small cluster sizes n and a peak at large n (v₀* ≲ 20), as shown in Fig. 20(a). Interestingly, the distribution functions of strongly-coupled dipolar particles (λ = 5) bear resemblance to the case of active Brownian particles, as can be seen in Fig. 20(b). For passive dipolar particles (v₀* = 0), the weighted distribution function vanishes at n ≈ 30, which is three times as large as in the non-dipolar case. This is due to the fact that the strongly-coupled particles tend to form head-to-tail configurations and therefore display chain-like structures. With increasing motilities, such structures disappear. Further, the dipolar particles start to align with their neighbors and form into clusters, as shown in Fig. 10(d). As indicated by the peaks at n ≈ 10³ in Fig. 20(b), these particles form “giant” clusters when v₀* ≳ 20.

Fig. 20 Weighted cluster size distribution at Φ = 0.58 for the motility v₀* = 0 (black circles), 20 (red squares), 60 (green diamonds), and 100 (blue triangles) with the coupling strength λ = 0 (a) and λ = 5 (b). We denote n as the cluster size.

Acknowledgements

The authors would like to thank Deutsche Forschungsgemeinschaft for the financial support from GRK 1524 (DFG No. 599982). This work was also supported by NSF's Research Triangle MRSEC under grant number DMR 1121107.

Notes and references

1. P. Romanczuk, M. Bär, W. Ebeling, B. Lindner and L. Schimansky-Geier, Eur. Phys. J.-Spec. Top., 2012, 202, 1.
2. A. Zöttl and H. Stark, J. Phys.: Condens. Matter, 2016, 28, 253001.
3. I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.
4. J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326.
5. E. Bertin, M. Droz and G. Grégoire, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 022101.
6. H. Chaté, F. Ginelli, G. Grégoire and F. Raynaud, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 046113.
7. A. P. Solon, H. Chaté and J. Tailleur, Phys. Rev. Lett., 2015, 114, 068101.
8. T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett., 1995, 75, 1226.
9. A. Martín-Gómez, D. Levis, A. Díaz-Guilera and I. Pagonabarraga, Soft Matter, 2018, 14, 2610.
10. E. Sesé-Sansa, I. Pagonabarraga and D. Levis, EPL, 2018, 124, 30004.
11. C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2016, 88, 045006.
12. F. Kogler and S. H. L. Klapp, EPL, 2015, 110, 10004.
13. A. Kaiser, K. Popowa and H. Löwen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 012301.
14. F. Guzmán-Lastra, A. Kaiser and H. Löwen, Nat. Commun., 2016, 7, 13519.
15. A. Kaiser, A. Snezhko and I. S. Aranson, Sci. Adv., 2017, 3, e1601469.
16. G. Kokot and A. Snezhko, Nat. Commun., 2018, 9, 2344.
17. R. P. Blakemore, Annu. Rev. Microbiol., 1982, 36, 217.
18. R. B. Frankel, Annu. Rev. Biophys. Bioeng., 1984, 13, 85.
19. S. Klumpp, C. T. Lefèvre, M. Bennet and D. Faivre, Phys. Rep., 2018, 789, 1.
20. N. Waisbord, C. T. Lefèvre, L. Bocquet, C. Ybert and C. Cottin-Bizonne, Phys. Rev. Fluids, 2016, 1, 053203.
21. F. Meng, D. Matsunaga and R. Golestanian, Phys. Rev. Lett., 2018, 120, 188101.
22. J. Harder and A. Cacciuto, Phys. Rev. E, 2018, 97, 022603.
23. J. Yan, M. Han, J. Zhang, C. Xu, E. Luijten and S. Granick, Nat. Mater., 2016, 15, 1095.
24. M. Han, J. Yan, S. Granick and E. Luijten, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 7513.
25. L. Baraban, D. Makarov, O. G. Schmidt, G. Cuniberti, P. Leiderer and A. Erbe, Nanoscale, 2013, 5, 1332.
26. L. Baraban, R. Streubel, D. Makarov, L. Han, D. Karnaushenko, O. G. Schmidt and G. Cuniberti, ACS Nano, 2013, 7, 1360.
27. J. M. Tavares, J. J. Weis and M. M. Telo da Gama, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 061201.
28. P. D. Duncan and P. J. Camp, J. Chem. Phys., 2004, 121, 11322.
29. P. D. Duncan and P. J. Camp, Phys. Rev. Lett., 2006, 97, 107202.
30. S. Kantorovich, J. J. Cerdà and C. Holm, Phys. Chem. Chem. Phys., 2008, 10, 1883.
31. J. J. Cerdà, S. Kantorovich and C. Holm, J. Phys.: Condens. Matter, 2008, 20, 204125.
32. S. H. L. Klapp and M. Schoen, J. Chem. Phys., 2002, 117, 8050.
33. W.-Z. Ouyang, S.-H. Xu and Z.-W. Sun, J. Chem. Phys., 2011, 134, 014901.
34. R. Geiger and S. H. L. Klapp, J. Mod. Phys., 2013, 04, 401.
35. S. van Teeffelen and H. Löwen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 020101.
36. G.-J. Liao and S. H. L. Klapp, Soft Matter, 2018, 14, 7873.
37. P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, Heidelberg, 1992.
38. The alternative definition of the area fraction is given by [capital Phi, Greek, tilde] = Nπσ²/(4L²). Therefore, Φ = [capital Phi, Greek, tilde](σ_eff/σ)². In this work, we consider three different density values [capital Phi, Greek, tilde] = 0.1, 0.2, and 0.5. Consequently, the exact area fractions to the fifth decimal place are Φ = 0.11632, 0.23264, and 0.58160.
39. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219.
40. J. Blaschke, M. Maurer, K. Menon, A. Zöttl and H. Stark, Soft Matter, 2016, 12, 9821.
41. J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2008, 100, 218103.
42. J. Bialké, H. Löwen and T. Speck, EPL, 2013, 103, 30008.
43. P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella and I. Pagonabarraga, Phys. Rev. Lett., 2018, 121, 098003.
44. The pair correlation function has been obtained with a spatial resolution Δr = 0.0125σ. The resulting exact values for r_L are r_L = 1.15625σ for Φ = 0.12, r_L = 1.13125σ for Φ = 0.23, and r_L = 1.11875σ for Φ = 0.58.
45. J.-J. Weis and D. Levesque, Advanced Computer Simulation Approaches for Soft Matter Sciences II, Springer, Berlin, Heidelberg, 2005, pp. 163–225.
46. J.-J. Weis, J. Chem. Phys., 2005, 123, 044503.
47. R. A. Trasca and S. H. L. Klapp, J. Chem. Phys., 2008, 129, 084702.
48. J. J. Weis, Mol. Phys., 2002, 100, 579.
49. J.-J. Weis, J. Phys.: Condens. Matter, 2003, 15, S1471.
50. J. J. Weis and D. Levesque, Phys. Rev. Lett., 1993, 71, 2729.
51. L. Rovigatti, J. Russo and F. Sciortino, Phys. Rev. Lett., 2011, 107, 237801.
52. L. Rovigatti, J. Russo and F. Sciortino, Soft Matter, 2012, 8, 6310.
53. L. Rovigatti, S. Kantorovich, A. O. Ivanov, J. M. Tavares and F. Sciortino, J. Chem. Phys., 2013, 139, 134901.
54. S. Kantorovich, A. O. Ivanov, L. Rovigatti, J. M. Tavares and F. Sciortino, Phys. Rev. Lett., 2013, 110, 148306.
55. S. S. Kantorovich and A. O. Ivanov, in Soft Matter Self-Assembly, ed. C. N. Likos, F. Sciortino, P. Ziherl and E. Zaccarelli, IOS Press, 2016, pp. 137–163.
56. M. Ronti, L. Rovigatti, J. M. Tavares, A. O. Ivanov, S. S. Kantorovich and F. Sciortino, Soft Matter, 2017, 13, 7870.
57. P. J. Camp, in Modern Problems of Molecular Physics, ed. L. A. Bulavin and A. V. Chalyi, Springer, Cham, 2018, pp. 185–204.
58. S. D. Peroukidis and S. H. L. Klapp, Soft Matter, 2016, 12, 6841.
59. H. Schmidle, C. K. Hall, O. D. Velev and S. H. L. Klapp, Soft Matter, 2012, 8, 1521.
60. H. Schmidle, S. Jäger, C. K. Hall, O. D. Velev and S. H. L. Klapp, Soft Matter, 2013, 9, 2518.
61. I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 268303.
62. S. S. Kantorovich, A. O. Ivanov, L. Rovigatti, J. M. Tavares and F. Sciortino, Phys. Chem. Chem. Phys., 2015, 17, 16601.
63. J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2013, 111, 145702.
64. J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489.
65. L. Luo and S. H. L. Klapp, J. Chem. Phys., 2009, 131, 034709.
66. G. Baglietto, E. V. Albano and J. Candia, Interface Focus, 2012, 2, 708.
67. B. Mahault, F. Ginelli and H. Chaté, 2019, arXiv: 1908.03794v1[cond-mat.stat-mech].
68. N. Yoshinaga and T. B. Liverpool, Phys. Rev. E, 2017, 96, 020603.
69. N. Yoshinaga and T. B. Liverpool, Eur. Phys. J. E, 2018, 41, 76.
70. R. Matas Navarro, R. Golestanian, T. B. Liverpool and S. M. Fielding, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 032304.
71. R. Matas Navarro and S. M. Fielding, Soft Matter, 2015, 11, 7525.
72. M. Theers, E. Westphal, K. Qi, R. G. Winkler and G. Gompper, Soft Matter, 2018, 14, 8590–8603.
73. F. J. Schwarzendahl and M. G. Mazza, J. Chem. Phys., 2019, 150, 184902.
74. F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56.
75. B. Delmotte, E. E. Keaveny, F. Plouraboué and E. Climent, J. Comput. Phys., 2015, 302, 524.
76. C. Hoell, H. Löwen and A. M. Menzel, J. Chem. Phys., 2018, 149, 144902.
77. E. Lauga, W. R. DiLuzio, G. M. Whitesides and H. A. Stone, Biophys. J., 2006, 90, 400.
78. M. Loose and T. J. Mitchison, Nat. Cell Biol., 2014, 16, 38.
79. F. Kümmel, B. ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen and C. Bechinger, Phys. Rev. Lett., 2013, 110, 198302.
80. L. Baraban, D. Makarov, M. Albrecht, N. Rivier, P. Leiderer and A. Erbe, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 031407.
81. G. Steinbach, S. Gemming and A. Erbe, Eur. Phys. J. E, 2016, 39, 69.
82. A. B. Yener and S. H. L. Klapp, Soft Matter, 2016, 12, 2066.
83. S. H. Klapp, Curr. Opin. Colloid Interface Sci., 2016, 21, 76.
84. M. Mazars, Phys. Rep., 2011, 500, 43.
85. M. Schoen and S. H. L. Klapp, Reviews in Computational Chemistry, John Wiley & Sons, Inc., 2007, pp. 301–340.
86. Y. Fily, S. Henkes and M. Cristina Marchetti, Soft Matter, 2014, 10, 2132.

---