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

TL;DR

This study uses Brownian Dynamics simulations to investigate how active dipolar particles self-assemble in two dimensions, revealing phenomena like chain formation, phase separation suppression, and flocking behavior depending on motility and interaction strength.

Contribution

It introduces a comprehensive analysis of active dipolar particles, highlighting the emergence of flocking and the suppression of phase separation due to dipolar interactions.

Findings

Dipolar particles form chains at low motility and density.

High motility breaks chains, leading to isotropic fluids.

Strong dipolar coupling suppresses phase separation and induces flocking.

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.

Tables

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

State	Clustering	Orientational	Chain formation
		ordering	$\Phi \lesssim 0.23$	$\Phi =0.58$
Homogeneous, isotropic fluids	${\varphi }_c\le 0.5$	${\varphi }_{𝒆}\le 0.5$	${\varphi }_p\le 0.5$	$𝒵\le 0.17$
Chain-like structures	${\varphi }_c\le 0.5$	${\varphi }_{𝒆}\le 0.5$	${\varphi }_p>0.5$	$𝒵>0.17$
Micro-flocking	${\varphi }_c\le 0.5$	${\varphi }_{𝒆}>0.5$	-
Macro-flocking	${\varphi }_c>0.5$	${\varphi }_{𝒆}>0.5$	-
Motility-induced clustering	${\varphi }_c>0.5$	${\varphi }_{𝒆}\le 0.5$	-

Equations

$$u_{pair}\big{(}\bm{r}_{ij},\bm{\mu}_{i},\bm{\mu}_{j}\big{)}=u_{sr}\left(r_{ij}\right)+u_{dd}\big{(}\bm{r}_{ij},\bm{\mu}_{i},\bm{\mu}_{j}\big{)}\text{,}$$

$$u_{sr}(r_{ij})=\{\begin{array}{ll}4ϵ\left[{\left(\frac{\sigma }{r_{ij}}\right)}^{12}-{\left(\frac{\sigma }{r_{ij}}\right)}^6+\frac{1}{4}\right]\text{,}& \text{if rij<rc,}\\ 0\text{,}& \text{else.}\end{array}$$

$$u_{dd}\big{(}\bm{r}_{ij},\bm{\mu}_{i},\bm{\mu}_{j}\big{)}=\dfrac{\bm{\mu}_{i}\cdot\bm{\mu}_{j}}{r_{ij}^{3}}-3\dfrac{\big{(}\bm{\mu}_{i}\cdot\bm{r}_{ij}\big{)}\big{(}\bm{\mu}_{j}\cdot\bm{r}_{ij}\big{)}}{r_{ij}^{5}}\text{.}$$

$${\dot{\mathbf{r}}}_i$$

$${\dot{\psi }}_i$$

$$U_{i}=\sum_{j=1,j\neq i}^{N}u_{pair}\big{(}\bm{r}_{ij},\bm{\mu}_{i},\bm{\mu}_{j}\big{)}\text{,}$$

$${\varphi }_c=⟨n_{lcl}⟩/N\text{, }$$

$${\varphi }_{\mathbf{e}}=\frac{1}{N}⟨\left|\sum _{i=1}^N{\hat{\mathbf{e}}}_i\right|⟩\text{,}$$

$${\varphi }_p=⟨N_p⟩/N\text{, }$$

$$g_{{}_{\bm{\mu}}}(r_{\bot},r_{\|})=\dfrac{\left\langle\sum_{i\neq j}\hat{\bm{\mu}}_{i}\cdot\hat{\bm{\mu}}_{j}\delta\big{(}r_{\bot}-r_{ij}^{\bot}\big{)}\delta\big{(}r_{\|}-r_{ij}^{\|}\big{)}\right\rangle}{\left\langle\sum_{i\neq j}\delta\big{(}r_{\bot}-r_{ij}^{\bot}\big{)}\delta\big{(}r_{\|}-r_{ij}^{\|}\big{)}\right\rangle}\text{, }$$

$${\tilde{g}}_{{}_{\mathbf{\mu }}}\left(\theta \right)=\frac{{\int }_0^{r_s}{g}_{{}_{\mathbf{\mu }}}(r,\theta )r\text{d}r}{{\int }_0^{r_s}r\text{d}r}\text{. }$$

$$\text{mathzapf}Z={\tilde{g}}_{{}_{\mathbf{\mu }}}\left({\theta }_{max}\right)-{\tilde{g}}_{{}_{\mathbf{\mu }}}\left({\theta }_{min}\right)\text{, }$$

$$v^{\ast }\left(\varphi \right)/v_0^{\ast }=1-a_d\varphi \text{,}$$

$$f\left(\varphi \right)=f_0\left(\varphi \right)+f_{rep}\left(\varphi \right)\text{,}$$

$$f_0\left(\varphi \right)=\varphi \left(\text{ln}\varphi -1\right)+{\int }_0^{\varphi }\text{ln}\left[v\left(\rho \right)\right]d\rho$$

$$f_{rep}\left(\varphi \right)=k_{rep}\Theta \left(\varphi -{\varphi }_t\right){\left(\varphi -{\varphi }_t\right)}^4$$

$$U_{dd}=$$

$$+\frac{1}{L^2}\sum _{\mathbf{k}\ne 0}\frac{\pi }{k}\text{erfc}\left(\frac{k}{2\alpha }\right)|\text{mathzapf}M\left(\mathbf{k}\right){|}^2$$

$$-\frac{2{\alpha }^3}{3\sqrt{\pi }}N{\mu }^2\text{.}$$

$$B\left(r,\alpha \right)\equiv$$

$$C\left(r,\alpha \right)\equiv$$

$$\text{mathzapf}M\left(\mathbf{k}\right)$$

$$\equiv \text{Re}\left\{\text{mathzapf}M\right\}+\text{mathzapf}i\text{Im}\left\{\text{mathzapf}M\right\}\text{.}$$

$${\mathbf{F}}_{i,dd}$$

$$={\mathbf{F}}_{i,dd}^R+{\mathbf{F}}_{i,dd}^{\mathbf{k}\ne 0}\text{,}$$

$$\displaystyle\bm{F}_{i,dd}^{R}=\sum_{i\neq j}\bigg{[}$$

$$\displaystyle-\big{(}\bm{\mu}_{i}\cdot\bm{r}_{ij}\big{)}\big{(}\bm{\mu}_{j}\cdot\bm{r}_{ij}\big{)}\bm{r}_{ij}D\left(r_{ij},\alpha\right)\bigg{]}\text{.}$$

$$\displaystyle D\left(r,\alpha\right)\equiv\dfrac{1}{r^{7}}\bigg{[}$$

$$\displaystyle+15\;\text{erfc}\left(\alpha r\right)\bigg{]}\text{.}$$

$${\mathbf{F}}_{i,dd}^{\mathbf{k}\ne 0}=\frac{2\pi }{L^2}\sum _{\mathbf{k}\ne 0}$$

$$\displaystyle\Big{(}\text{sin}\left(\bm{k}\cdot\bm{r}_{i}\right)\text{Re}\big{\{}\mathzapf{M}\left(\bm{k}\right)\big{\}}+\text{cos}\left(\bm{k}\cdot\bm{r}_{i}\right)\text{Im}\big{\{}\mathzapf{M}\left(\bm{k}\right)\big{\}}\Big{)}\text{,}$$

$${\mathbf{T}}_{i,dd}$$

$$={\mathbf{T}}_i^R+{\mathbf{T}}_i^{\mathbf{k}\ne 0}\text{,}$$

$$\bm{T}_{i}^{R}=-\sum_{j\neq i}\bigg{[}\big{(}\bm{\mu}_{i}\times\bm{\mu}_{j}\big{)}B\left(r_{ij},\alpha\right)-\big{(}\bm{\mu}_{i}\times\bm{r}_{ij}\big{)}\big{(}\bm{\mu}_{j}\cdot\bm{r}_{ij}\big{)}C\left(r_{ij},\alpha\right)\bigg{]}\text{.}$$

$${\mathbf{T}}_i^{\mathbf{k}\ne 0}=\frac{2\pi }{L^2}\sum _{\mathbf{k}\ne 0}$$

Full text

\makeFNbottom

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

Guo-Jun Liao,*∗a Carol K. Hall,b and Sabine H. L. Klapp∗*a

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.

††footnotetext: a*  Institut für Theoretische Physik, Technische Universität Berlin, Hardenbergstr. 36, D-10623 Berlin, Germany. E-mail: guo-jun.liao@campus.tu-berlin.de, klapp@physik.tu-berlin.de*††footnotetext: b* Department of Chemical & Biomolecular Engineering, North Carolina State University, Raleigh, NC 27695, USA. *

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. 3 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, 5, 6, 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 6 and flocking. 8 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. 9 Further, at high densities, the polar coupling favors motility-induced phase separation compared to the non-polar case. 10

In most studies of active systems with anisotropic interactions, they are assumed to be of short range. 11 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. 12 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 8 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 13 or including 14 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. 15 Owing to the interplay between self-propulsion and dipolar interactions, these vortices persist even on a flat surface. 16 In nature, magnetotactic bacteria are known to sense the earth’s magnetic field and move along or against the field direction. 17, 18, 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. 20 The underlying mechanism of this clustering instability was also investigated analytically and in simulations. 21 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. 22 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 Sec. 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 Sec. 3. Finally, in Sec. 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 $\sigma$ dispersed in a monolayer in the $xy$-plane. Each particle carries a fluctuating point dipole moment $\bm{\mu}_{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, 28, 29, 30, 31 or dense ordered states, 27, 32, 33, 34 if the dipolar interactions are sufficiently strong. Having this in mind, we assume beforehand that $\bm{\mu}_{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 $\bm{F}_{i}$, which has a constant magnitude and is directed along $\bm{\mu}_{i}$ at each instant of time.

The pair potential between particles located at positions $\bm{r}_{i}$ and $\bm{r}_{j}$ ($i\neq j$) is of the form

[TABLE]

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}=|\bm{r}_{ij}|=|\bm{r}_{j}-\bm{r}_{i}|$. We assume that the steric repulsion can be described by the Weeks-Chandler-Anderson potential

[TABLE]

In this study, we set the length unit to be $\sigma$ and fix the repulsive strength $\epsilon^{*}=\beta\epsilon=10$, where the thermal energy, $\beta^{-1}=k_{B}T$, 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}\sigma$, 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

[TABLE]

We define the dipolar coupling strength as $\lambda=\beta\mu^{2}\sigma^{-3}$ with $\mu$ being the magnitude of a dipole moment, i.e., $\bm{\mu}_{i}=\mu\hat{\bm{\mu}}_{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 $i$th particle is then determined by the coupled Langevin equations 35 for its position $\bm{r}_{i}$ and orientation $\hat{\bm{e}}_{i}=\left(\text{cos}\psi_{i},\text{sin}\psi_{i}\right)^{T}$,

[TABLE]

where the dots denote time derivatives and $\psi_{i}$ is the polar angle. In Eqs. (4)-(5), the potential energy for the $i$th particle is given by

[TABLE]

with $u_{pair}$ being defined in eqn (1). Since the shape of each particle is modeled as a disk, we set the translational diffusion tensor $\mathbb{D}=D_{t}\mathbb{I}$, where $D_{t}$ is the (isotropic) translational diffusion constant and $\mathbb{I}$ is the $2\times 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}/\sigma_{h}^{2}$. We also define the diameter of an effective (eff) hard sphere via $\sigma_{eff}=\int_{0}^{\infty}\left(1-\text{exp}\left[-\beta u_{sr}(r)\right]\right)\text{d}r$. With the repulsive strength $\epsilon^{*}=10$ considered in the present system, $\sigma_{eff}\approx 1.07851\sigma$. By choosing $\sigma_{h}=\sigma_{eff}$, we obtain $D_{r}=2.57914D_{t}/\sigma^{2}$.

In eqn (4), the effective propulsion force, which drives the active motion, is given by $\bm{F}_{i}=F_{0}\hat{\bm{e}}_{i}$ with $\hat{\bm{e}}_{i}=\hat{\bm{\mu}}_{i}$. For simplicity, in the remainder of this work we present the impact of the effective propulsion force via the motility $v_{0}=\beta D_{t}F_{0}$. Finally, the random force $\bm{\xi}_{i}\left(t\right)$ and torque $\Gamma_{i}\left(t\right)$ for the $i$th particle are zero-mean Gaussian white noise, which satisfy $\langle\bm{\xi}_{i}(t)\otimes\bm{\xi}_{j}(t^{\prime})\rangle=2\delta_{ij}\delta(t-t^{\prime})\mathbb{I}/(D_{t}\beta^{2})$ and $\langle\Gamma_{i}(t)\Gamma_{j}(t^{\prime})\rangle=2\delta_{ij}\delta(t-t^{\prime})/(D_{r}\beta^{2})$. The angle brackets $\langle\cdots\rangle$ denote ensemble average, and the symbol $\otimes$ represents dyadic product.

Equations (4) and (5) are solved via the Euler-Maruyama method 37 with the discrete time step $\Delta t=2\times 10^{-5}\tau$, where $\tau$ is the Brownian diffusion time, given by $\tau=\sigma^{2}/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\times 10^{5}$ steps for reaching a steady state, followed by a production period of $5\times 10^{5}$ steps. The statistical properties of the dipolar particles (see Sec. 2.3) are measured every $500$ steps. The simulations are carried out at three values of the mean area fraction $\Phi=N\pi\sigma_{eff}^{2}/\left(4L^{2}\right)$, where the area of a single particle is defined as $\pi\sigma_{eff}^{2}/4$. Specifically, we consider the values $\Phi=0.12$, $\Phi=0.23$, and $\Phi=0.58$ (for the exact values, see note 38). To investigate the impact of activity, we perform simulations at different dimensionless propulsion speeds $v_{0}^{*}=v_{0}\sigma/D_{t}$ and various dipolar coupling strengths $\lambda$. We note that $v_{0}^{*}$ 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). 41, 42, 39, 43 We would therefore expect that the present system displays similar behavior at least at low dipolar coupling strengths, i.e., small values of $\lambda$. 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

[TABLE]

where $n_{lcl}$ denotes the number of particles in the largest cluster. The quantity $\phi_{c}$ approaches [math] 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, $\phi_{c}$ is close to zero in a state with finite-sized chains or clusters with a size much smaller than $N$. In contrast, $\phi_{c}\approx 1$ when the active particles form “giant” clusters with a size comparable to $N$. 3 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$, $\lambda=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}\approx 1.16\sigma$ for $\Phi=0.12$, $r_{L}\approx 1.13\sigma$ for $\Phi=0.23$, and $r_{L}\approx 1.12\sigma$ for $\Phi=0.58$ (for the exact values, see note 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 45 and by theory. 46 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 8 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 9) the parameter

[TABLE]

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 ($\lambda>1$) are known to self-assemble into chains and rings. 50, 28, 45, 29, 51, 52, 53, 54, 55, 56, 57 The reason is easily seen from eqn (3): Considering two dipolar spheres separated by a distance $\sigma$, the configuration with the lowest energy is the head-to-tail configuration (with $u_{dd}=-2\mu^{2}\sigma^{-3}$). In contrast, if the particles are arranged side by side, the energy reaches a maximum for parallel orientation (with $u_{dd}=\mu^{2}\sigma^{-3}$), 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 ($\Phi\lesssim 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: $|\bm{r}_{ij}|\leq r_{p}$, $\hat{\bm{\mu}}_{i}\cdot\hat{\bm{\mu}}_{j}>0$, and \big{(}\hat{\bm{\mu}}_{i}\cdot\bm{r}_{ij}\big{)}\big{(}\hat{\bm{\mu}}_{j}\cdot\bm{r}_{ij}\big{)}>0. Here, $r_{p}=1.25\sigma$ is set to a distance between the location of the first peak and the first valley of the pair correlation function. 58, 59, 60 Based on these rules, we quantify the low-density chain formation via the parameter

[TABLE]

where $N_{p}$ is the total number of particles which reside in chains. In the context of aggregating molecular systems, $\phi_{p}$ is often called the degree of polymerization.

At high densities (e.g., $\Phi=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 $\phi_{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

[TABLE]

where the transverse and longitudinal displacement of particle $j$ relative to $i$ is $r_{ij}^{\|}=\bm{r}_{ij}\cdot\hat{\bm{\mu}}_{i}$ and $r_{ij}^{\bot}=\sqrt{|\bm{r}_{ij}|^{2}-\left(\bm{r}_{ij}\cdot\hat{\bm{\mu}}_{i}\right)^{2}}$, respectively. To extract the angular dependence at short distances, we transform the Cartesian coordinates employed in eqn (10) to polar coordinates by using $r=\sqrt{r_{\bot}^{2}+r_{\|}^{2}}$ and \theta=\text{atan2}\big{(}-r_{\bot},\>r_{\|}\big{)}, and compute the function

[TABLE]

Here, the distance $r_{s}$ is set to $r_{s}=L/8\approx 5\sigma$ 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, $\tilde{g}_{\bm{\mu}}\left(\theta\right)$, to display positive maxima both, in front of and behind the reference particle (i.e., at $\theta=0$ or $\pi$). 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., $\tilde{g}_{{}_{\bm{\mu}}}\left(\theta\right)\approx 0$ for $\theta=\pi/2$ or $-\pi/2$). With this picture in mind, we measure the quantity

[TABLE]

where $\tilde{g}_{{}_{\bm{\mu}}}\left(\theta\right)$ reaches its maximum and minimum, respectively, at $\theta=\theta_{max}$ and $\theta=\theta_{min}$. By properly choosing a threshold value $\mathzapf{Z}_{thres}$, we define the dipolar particles as exhibiting chain-like structure when $\mathzapf{Z}\geq\mathzapf{Z}_{thres}$. As will be discussed in more detail in Sec. 2.3.3, a reasonable value for $\mathzapf{Z}_{thres}$ at $\Phi=0.58$ is $\mathzapf{Z}_{thres}=0.17$.

3 Simulation results

Based on the target quantities described in Sec. 2.3, we can classify the states observed in our BD simulations performed at three densities and at various values of the motility $v_{0}^{*}$ and the dipolar coupling strength $\lambda$. In the following Sec. 3.1 $-$ 3.3 we discuss, for each density, the state diagram in the ($v_{0}^{*}$, $\lambda$) plane. Specifically, we identify five states whose characteristics are summarized in Table 1.

3.1 The low density regime ($\Phi=0.12$)

We start by considering the state diagram in the low-density regime, taking $\Phi=0.12$ as a representative example, see Fig. 1. In the absence of dipolar coupling ($\lambda=0$), non-dipolar active particles at such a low density typically display a homogeneous isotropic fluid state, without significant translational correlations or orientational ordering. 3 As long as the dipolar coupling strength remains relatively small ($\lambda\lesssim 4$), the system still shows the same homogeneous isotropic fluid state for all motilities considered. In contrast, when $\lambda$ 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}^{*}\left(\lambda\right)$, whose value depends on $\lambda$. Above this critical motility $v_{0,c}^{*}\left(\lambda\right)$, the chain-like structures found at smaller motilities break, and the system becomes homogeneous and isotropic. As Fig. 1 reveals, $v_{0,c}^{*}\left(\lambda\right)$ increases with increasing $\lambda$. As a visualization of the impact of motility on dipolar active particles at a high coupling strength, such as $\lambda=10$, we plot in Fig. 2 representative snapshots. At $v_{0}^{*}=0$, nearly all particles are bound into chains and rings, as it is typical for passive dipolar particles in 2D. 27 With increasing $v_{0}^{*}$ the ring structures are seemingly more abundant, as can be seen in Fig. 2(a-c). Further increase in $v_{0}^{*}$ eventually causes the chain-like structures to break into fragments of short chains and single particles, as shown in Fig. 2(d).

To quantitatively characterize the motility-induced destruction of chain-like structures, we plot in Fig. 3 the degree of polymerization $\phi_{p}$ as a function of the motility $v_{0}^{*}$ for various coupling strengths $\lambda$. In the absence of dipole-dipole interactions ($\lambda=0$), particles do not self-assembly into chain-like structures, still we observe a slight increase in $\phi_{p}$ [see eqn (9)] upon increasing motilities $v_{0}^{*}$. This can be attributed to the dynamical clustering of active Brownian particles into small (finite-sized) aggregates. 36, 61, 3

Indeed, if the particles inside these finite-sized clusters fulfill our criteria of being bonded, the order parameter $\phi_{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 $\lambda\gtrsim 6$ and low motilities $v_{0}^{*}$, giving rise to large values of $\phi_{p}$ ($\phi_{p}\approx 1$). Upon an increase in the motility, the degree of polymerization $\phi_{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. 62 In this sense, the particle motility in the active dipolar system may be viewed as an analog to the temperature in the passive system.

The degree of polymerization $\phi_{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\lesssim 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 ($\lambda=10$) look very different, as shown in Fig. 4(b). Here, we observe a broad peak between $10\lesssim n\lesssim 100$ at $v_{0}^{*}=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_{0}^{*}\gtrsim 60$. This suggests a vanishing of chains with the most probable size at large values of $v_{0}^{*}$.

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_{0}^{*}$. In other words, there is no pronounced global orientational order. Indeed, through measuring the orientational order parameter $\phi_{\bm{e}}$, we confirm that $\phi_{\bm{e}}$ remains small ($\phi_{\bm{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 ($\Phi=0.23$)

At the density $\Phi=0.23$, the state diagram (see Fig. 5) at small and intermediate motilities ($v_{0}^{*}\lesssim 20$) is very similar to that at $\Phi=0.12$ (see Fig. 1): For weakly coupled systems ($\lambda\lesssim 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 $\lambda$. 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.

We now turn to the impact of the motility. At small values of $\lambda$, an increase in $v_{0}^{*}$ 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 ($\phi_{c}\approx 1$), is known to occur only at higher densities. 3 As Fig. 5 reveals, the picture changes if the dipolar coupling becomes larger. Consider, e.g., the case $\lambda=6$. Increasing the motility from zero, the chain-like structures observed at small $v_{0}^{*}$ 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 $\phi_{c}$ remains small ($\phi_{c}<0.5$) but the orientational order parameter $\phi_{\bm{e}}$ reaches values above $0.5$.

The full behavior of the orientational order parameter $\phi_{\bm{e}}$ as function of $v_{0}^{*}$ is shown in Fig. 7. For strongly coupled systems ($\lambda\geq 3$), the order parameter abruptly changes from (essentially) zero to large values at a “critical” motility $v_{0,c}^{*}\left(\lambda\right)\approx 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_{0}$ does not scale up with the system size $N$ (see the detailed analysis in Appendix A.2).

Inspecting Fig. 7 again, it is found that for systems at an intermediate coupling strength $\lambda=2$, the order parameter $\phi_{\bm{e}}$ gradually increases from zero as $v_{0}^{*}\gtrsim 100$, suggesting that the “micro-flocking” state might still appear. However, further increase in $v_{0}^{*}$ requires to employ even smaller simulation time steps $\Delta t<2\times 10^{-5}\tau$, 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 ($\lambda<3$).

A behavior similar to the “micro-flocking” state has been previously observed in systems of active Brownian disks with polar alignment. 9 Indeed, the density $\Phi=0.23$ considered here is close to the effective packing fraction $\phi^{*}=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). 8

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_{0}$) 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}^{*}>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_{0}$. Here, activity leads to motion in opposite direction, yielding this configuration unstable. In Fig. 8(c), the initial configuration at $t_{0}$ 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., $\Phi\gtrsim 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_{0}$ in Fig. 8(c) is kinetically stable. 9, 10

3.3 The high density regime ($\Phi=0.58$)

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

3.3.1 Motility-induced clustering

Figure 9 shows the state diagram at $\Phi=0.58$ in the plane spanned by motilities $v_{0}^{*}$ and dipolar coupling strengths $\lambda$. We start by investigating the regime of small $\lambda$. 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_{0}^{*}$ increases from zero, non-dipolar active Brownian particles ($\lambda=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. 42, 39 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. 39 We note that the particle orientations in the two phases are essentially uncorrelated for non-dipolar particles ($\lambda=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, $\phi_{c}$ [see eqn (7)], as a function of the motility $v_{0}^{*}$ for various coupling strengths $\lambda$. Consistent with the motility-induced phase separation, the curves for non-dipolar and weakly coupled active particles ($\lambda=0-1$) display a sharp increase from zero at $v_{0}^{*}\approx 20$ and reach large values ($\phi_{c}\approx 0.8$) at $v_{0}^{*}=100$. The curves reach $\phi_{c}=0.5$ at larger $v_{0}^{*}$ with increasing $\lambda$, which is also reflected by the boundary between black circles and red squares in Fig 9. For strong dipolar coupling ($\lambda\gtrsim 2$), the values of $\phi_{c}$ at high motilities remain in the same range ($\phi_{c}\approx 0.8$, indicating again formation of large clusters), but the increase in $\phi_{c}$ at small $v_{0}^{*}$ 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.

To gain further information, we calculate position-resolved local area fractions, $\phi$, based on a Voronoi tessellation. 40, 36 (We note that, contrary to ref. 36, we did not perform a short time average of $\phi$, 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 Sec. 3.3.2.) Figure 12 shows the probability distribution of the local area fractions, $P\left(\phi\right)$, for different coupling strengths at $v_{0}^{*}=100$. For non-dipolar active particles ($\lambda=0$), $P\left(\phi\right)$ 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\left(\phi\right)$. Upon an increase in the coupling strength up to $\lambda=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 ($\Phi=0.58$), as well as a broad shoulder on the left. This suggests the disappearance of the phase separation observed at zero and small $\lambda$.

The same conclusion can be drawn from Fig. 13, where we plot the coexisting densities in the ($v_{0}^{*}$, $\phi$) plane for various $\lambda$. The location of the high-density branch at $\lambda=0$ can be explained as follows: At high motilities $v_{0}^{*}$, the non-dipolar system ($\lambda=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 $\sigma_{eff}$ defined in Sec. 2.2. As a result, the local area fraction in the densely packed region should be close to the close-packing fraction, $\phi_{cp}=\pi/\left(2\sqrt{3}\right)\approx 0.91$. Upon a slight increase in dipolar coupling, the coexistence branches for the dilute region, $\phi_{gas}$, are shifted toward higher densities, while the branches for the dense region, $\phi_{den}$, move toward lower area fractions. As a result, the area surrounded by the curves of $\phi_{den}$ and $\phi_{gas}$ in Fig. 13 significantly shrinks with increasing $\lambda$, 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 $\lambda>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. 10 This shows that different types of orientational interactions may have entirely different impacts on motility-induced phase separation.

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 $\lambda=0$). To this end, we compute the normalized speed of the $i$th particle, $v_{i}^{*}/v_{0}^{*}$, versus the local density $\phi_{i}=\phi\left(\bm{r}_{i}\right)$, where the speed is given by $v_{i}^{*}=|\left(\Delta\bm{r}_{i}/\sigma\right)/\left(\Delta t_{s}/\tau\right)|$ with $\Delta t_{s}=10^{-2}\tau$. 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 ($\lambda=0$) and weakly coupled ($0<\lambda\lesssim 1$) systems, we observe a linear decay,

[TABLE]

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 $\lambda$ in the range $\lambda\lesssim 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 $\lambda\gtrsim 2$, the normalized speed as a function of $\phi$ 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 Sec. 3.3.2.

To obtain a more complete (yet qualitative) picture on the role of the function $a_{d}\left(\lambda\right)$ for the appearance of the motility-induced phase separation, we consider the effective free energy proposed in ref. 41,

[TABLE]

where the bulk contribution is given by

[TABLE]

with density-dependent swim speed $v\left(\rho\right)=1-a_{d}\rho$. 63, 64 Further, the contribution from the excluded volume interactions between particles is written as

[TABLE]

with $\Theta\left(x\right)$ being the Heaviside step function, $k_{rep}$ the repulsive strength, and $\phi_{t}$ the threshold area fraction.

For conventional active Brownian particles ($\lambda=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. 64 With this in mind, in Fig. 15 we plot the effective free energy $f\left(\phi\right)$ for the decay amplitude $a_{d}=1.08$, $0.95$, and $0.8$. At low densities, $f\left(\phi\right)$ is nearly independent of $a_{d}$, while at high densities, $f\left(\phi\right)$ increases with decreasing $a_{d}$. From $f\left(\phi\right)$, 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 $\left(a_{d},\phi\right)$ plane for the effective free energy (blue circles) and simulations at $v_{0}^{*}=100$ and $0\leq\lambda\leq 1$ (orange squares). More specifically, for the above simulations we plot the coexisting densities determined from Fig. 12 versus 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 $\phi_{den}$ and a growth in $\phi_{gas}$. Further, the phase separation disappears once $v\left(\phi\right)$ decays sufficiently slow ($a_{d}\lesssim 0.6$ for the effective free energy, and $a_{d}\lesssim 0.5$ for simulations). This clearly indicates that the motility-induced phase separation is suppressed when the function $v\left(\phi\right)$ decreases too slowly, i.e., when the dipolar coupling becomes too strong.

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 Sec. 3.2, once the density is sufficiently high ($\Phi\gtrsim 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, $\phi_{\bm{e}}$, as a function of $\lambda$ for various values of $v_{0}^{*}$. In the passive case ($v_{0}^{*}=0$), the order parameter remains small for all $\lambda$ 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_{0}^{*}=20$), the order parameter reaches significant values ($\phi_{e}\gtrsim 0.5$) when $\lambda$ 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\lesssim\lambda\lesssim 6$, $10\lesssim v_{0}^{*}\lesssim 40$), see green triangles in Fig. 9. Finally, in the regime of high motilities ($v_{0}^{*}=60-100$), the data curve $\phi_{\bm{e}}\left(\lambda\right)$ is reminiscent of a (polar) phase transition: $\phi_{\bm{e}}$ rises suddenly from zero to values greater than $0.5$ at a “critical” coupling strength $\lambda_{c}\approx 0.75$. This critical value slightly decreases upon an increase in $v_{0}^{*}$.

In addition to the large values of $\phi_{\bm{e}}$, the fraction of the largest cluster, $\phi_{c}$, is greater than $0.5$ for $v_{0}^{*}\gtrsim 60$, as shown in Fig. 11. Thus the systems at $v_{0}^{*}\gtrsim 60$ and $\lambda\gtrsim 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 $\tau$ (not shown here).

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_{0}^{*}$ 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 10(d)]. As discussed in Sec. 2.3.3, the order parameter $\phi_{p}$ (degree of polymerization) is no longer adequate to quantify such a transition at high densities.

Therefore, an alternative method is proposed: Figure 17(a-d) shows the spatial correlation function of dipole moments, g_{{}_{\bm{\mu}}}\big{(}r_{\bot},r_{\|}\big{)} [defined in eqn (10)], for $v_{0}^{*}=0-100$ and $\lambda=5$. For the passive and nearly passive case [see Fig. 17(a-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_{0}^{*}$ increases [see Fig. 17(c-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 $\bm{r}_{ij}$, Fig. 17(e) shows the angular correlation function, $\tilde{g}_{{}_{\bm{\mu}}}\left(\theta\right)$ [defined in eqn (11)], for the strong dipolar coupling strength $\lambda=5$. As mentioned in Sec. 2.3.3, the difference between the maximum $\tilde{g}_{{}_{\bm{\mu}}}\left(\theta_{max}\right)$ and the minimum $\tilde{g}_{{}_{\bm{\mu}}}\left(\theta_{min}\right)$ of the curve, $\mathzapf{Z}$ [defined in eqn (12)], can be used as a measure of chain formation. Upon increasing $v_{0}^{*}$, $\mathzapf{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 $\mathzapf{Z}$ as a function of $v_{0}^{*}$ (for various $\lambda$) is given in Fig. 18. The data curves bear close resemblance to those for $\phi_{p}$ in Fig. 3. For a non-dipolar system ($\lambda=0$), where chains are absent, the order parameter $\mathzapf{Z}$ does not show any significant behavior within the range of motilities explored. In contrast, at large coupling strengths ($\lambda=4$ and $\lambda=5$), the decrease in $\mathzapf{Z}$ with increasing $v_{0}^{*}$ reflects the disappearance of a state with chain-like structures. Specifically, the curves for $\lambda\geq 4$ show a sharp decrease, accompanied by a point of inflection $\mathzapf{Z}\approx 0.17$. We choose this as a threshold value. We note that the simulations of nearly passive, strongly-coupled systems (e.g., $\lambda=5$ and $6$) are severely plagued by large fluctuations of the order parameter $\mathzapf{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\tau$). The corresponding region in the state diagram in Fig. 9 is marked by dashed lines.

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 $\Phi$. At low densities ($\Phi=0.12$) and strong dipolar coupling, an increase in the motility $v_{0}^{*}$ 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_{0}^{*}$ 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 ($\Phi\gtrsim 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. 6 We note that the sizes of the present simulations are limited to the order of $10^{3}$ 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\approx 10^{3}$, 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. 8, 6, 7 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 $\Phi=0.58$ and large values of $v_{0}^{*}$, 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 $\lambda$ is negligible against thermal fluctuations ($0<\lambda\lesssim 1$). Once the dipolar interactions dominate ($\lambda>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, 71, 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. 73 Moreover, hydrodynamically interacting particles with a certain range of force dipole strengths can spontaneously form a state with global polar order, 74, 75, 68, 69 which can be attributed to either the actively induced rotation-translation coupling, 76 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, 81, 82, 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.

Appendix 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. 84 Within this method, the total dipole-dipole energy is separated into different contributions,

[TABLE]

The first term on the right-hand side of eqn (A.1) corresponds to the real-space contribution, where the functions $B\left(r,\alpha\right)$ and $C\left(r,\alpha\right)$ are defined by 32

[TABLE]

with erfc($x$) being the complementary error function. In eqn (A.1), the real-space contribution is formulated under the assumption that the convergence parameter $\alpha$ is large enough, such that we can consider only the interactions within the central simulation box. This is achieved by choosing $\alpha=7/L$ (with $L$ being the box size) and evaluating $\bm{r}_{ij}$ with the minimum-image convention. 32 The second term in eqn (A.1) represents the reciprocal-space contribution with

[TABLE]

The wave vectors in the reciprocal (square) lattice are given by $\bm{k}=\left(2\pi/L\right)\bm{m}$, where $\bm{m}=\left(m_{x},m_{y}\right)^{T}$ with $m_{x}$ and $m_{y}$ being integers. The magnitude of the wave vector is denoted as $k=2\pi\sqrt{m_{x}^{2}+m_{y}^{2}}/L$. We evaluate the reciprocal-space summation for the $\bm{k}$-vectors within the range $\bm{m}^{2}=m_{x}^{2}+m_{y}^{2}\leq 15^{2}$. 85 The third term represents the correction due to unphysical self-interaction of dipole moments. 85

The Langevin equations (4) and (5) involve the forces and torques due to the dipolar interactions. From eqn (A.1), the force acting on the $i$th particle due to dipole-dipole interactions is given by

[TABLE]

where the real-space contribution is

[TABLE]

The function $C\left(r,\alpha\right)$ is defined in eqn (19), and $D\left(r,\alpha\right)$ is given by

[TABLE]

Further, the contribution from the Fourier-form is written as

[TABLE]

where \text{Re}\left\{\mathzapf{M}\left(\bm{k}\right)\right\}=\sum_{j=1}^{N}\big{(}\bm{k}\cdot\bm{\mu}_{j}\big{)}\text{cos}\left(\bm{k}\cdot\bm{r}_{j}\right) and \text{Im}\left\{\mathzapf{M}\left(\bm{k}\right)\right\}=-\sum_{j=1}^{N}\big{(}\bm{k}\cdot\bm{\mu}_{j}\big{)}\text{sin}\left(\bm{k}\cdot\bm{r}_{j}\right) [see eqn (20)].

In the present work, every dipole moment is restricted to point in the $xy$-plane with $\psi_{i}$ characterizing its direction relative to the $x$-axis. Therefore, the torque acting on the $i$th particle is always directed along the $z$-axis. It is given by

[TABLE]

where $U_{dd}$ is given in eqn (A.1), and $\hat{\bm{z}}$ denotes the unit vector along the positive $z$-axis. The real-space component of the torque is given by

[TABLE]

Finally, the reciprocal-space contribution is given by

[TABLE]

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. Figure 19 shows the weighted cluster size distribution $nP\left(n\right)$ at $\Phi=0.23$, $v_{0}^{*}=100$ and $\lambda=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\gtrsim 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,

[TABLE]

where $n_{0}$ denotes the characteristic cluster size. Similar fitting functions of the cluster size distribution have also been considered in systems of conventional active Brownian particles 86 and polar active disks. 9 The inset of Fig. 19 shows that $n_{0}$ reaches a plateau once $N\gtrsim 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 $\Phi=0.23$ are microscopic structures.

To provide additional information on the clustering behavior at high densities ($\Phi=0.58$), we plot in Fig. 20 the weighted cluster size distribution for non-dipolar and strong coupled dipolar active particles at $\Phi=0.58$. At the coupling strength $\lambda=0$, our model reduces to the limiting case of active Brownian particles. 63, 3 Here, the correspondent weighted cluster size distribution function undergoes a transition from an exponential decay ($v_{0}^{*}\lesssim 20$) into a curve with a power law decay at small cluster sizes $n$ and a peak at large $n$ ($v_{0}^{*}\gtrsim 20$), as shown in Fig. 20(a). Interestingly, the distribution functions of strongly-coupled dipolar particles ($\lambda=5$) bear resemblance to the case of active Brownian particles, as can be seen in Fig. 20(b). For passive dipolar particles ($v_{0}^{*}=0$), the weighted distribution function vanishes at $n\approx 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\approx 10^{3}$ in Fig. 20(b), these particles form “giant” clusters when $v_{0}^{*}\gtrsim 20$.

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

• Romanczuk et al. 2012

P. Romanczuk, M. Bär, W. Ebeling, B. Lindner and L. Schimansky-Geier, Eur. Phys. J. Spec. Top., 2012, 202, 1.

• Zöttl and Stark 2016

A. Zöttl and H. Stark, J. Phys.: Condens. Matter, 2016, 28, 253001.

• Buttinoni et al. 2013

I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger and T. Speck, Phys. Rev. Lett., 2013, 110, 238301.

• Toner and Tu 1995

J. Toner and Y. Tu, Phys. Rev. Lett., 1995, 75, 4326.

• Bertin et al. 2006

E. Bertin, M. Droz and G. Grégoire, Phys. Rev. E, 2006, 74, 022101.

• Chaté et al. 2008

H. Chaté, F. Ginelli, G. Grégoire and F. Raynaud, Phys. Rev. E, 2008, 77, 046113.

• Solon et al. 2015

A. P. Solon, H. Chaté and J. Tailleur, Phys. Rev. Lett., 2015, 114, 068101.

• Vicsek et al. 1995

T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen and O. Shochet, Phys. Rev. Lett., 1995, 75, 1226.

• Martín-Gómez et al. 2018

A. Martín-Gómez, D. Levis, A. Díaz-Guilera and I. Pagonabarraga, Soft Matter, 2018, 14, 2610.

• Sesé-Sansa et al. 2018

E. Sesé-Sansa, I. Pagonabarraga and D. Levis, EPL, 2018, 124, 30004.

• Bechinger et al. 2016

C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys, 2016, 88, 045006.

• Kogler and Klapp 2015

F. Kogler and S. H. L. Klapp, EPL, 2015, 110, 10004.

• Kaiser et al. 2015

A. Kaiser, K. Popowa and H. Löwen, Phys. Rev. E, 2015, 92, 012301.

• Guzmán-Lastra et al. 2016

F. Guzmán-Lastra, A. Kaiser and H. Löwen, Nat. Commun., 2016, 7, 13519.

• Kaiser et al. 2017

A. Kaiser, A. Snezhko and I. S. Aranson, Sci. Adv., 2017, 3, e1601469.

• Kokot and Snezhko 2018

G. Kokot and A. Snezhko, Nat. Commun., 2018, 9, 2344.

• Blakemore 1982

R. P. Blakemore, Ann. Rev. Microbiol., 1982, 36, 217.

• Frankel 1984

R. B. Frankel, Annu. Rev. Biophys. Bioeng., 1984, 13, 85.

• Klumpp et al. 2018

S. Klumpp, C. T. Lefèvre, M. Bennet and D. Faivre, Phys. Rep., 2018, 789, 1.

• Waisbord et al. 2016

N. Waisbord, C. T. Lefèvre, L. Bocquet, C. Ybert and C. Cottin-Bizonne, Phys. Rev. Fluids, 2016, 1, 053203.

• Meng et al. 2018

F. Meng, D. Matsunaga and R. Golestanian, Phys. Rev. Lett., 2018, 120, 188101.

• Harder and Cacciuto 2018

J. Harder and A. Cacciuto, Phys. Rev. E, 2018, 97, 022603.

• Yan et al. 2016

J. Yan, M. Han, J. Zhang, C. Xu, E. Luijten and S. Granick, Nat. Mater., 2016, 15, 1095.

• Han et al. 2017

M. Han, J. Yan, S. Granick and E. Luijten, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 7513.

• Baraban et al. 2013

L. Baraban, D. Makarov, O. G. Schmidt, G. Cuniberti, P. Leiderer and A. Erbe, Nanoscale, 2013, 5, 1332.

• Baraban et al. 2013

L. Baraban, R. Streubel, D. Makarov, L. Han, D. Karnaushenko, O. G. Schmidt and G. Cuniberti, ACS Nano, 2013, 7, 1360.

• Tavares et al. 2002

J. M. Tavares, J. J. Weis and M. M. Telo da Gama, Phys. Rev. E, 2002, 65, 061201.

• Duncan and Camp 2004

P. D. Duncan and P. J. Camp, J. Chem. Phys., 2004, 121, 11322.

• Duncan and Camp 2006

P. D. Duncan and P. J. Camp, Phys. Rev. Lett., 2006, 97, 107202.

• Kantorovich et al. 2008

S. Kantorovich, J. J. Cerdà and C. Holm, Phys. Chem. Chem. Phys., 2008, 10, 1883.

• Cerdà et al. 2008

J. J. Cerdà, S. Kantorovich and C. Holm, J. Phys.: Condens. Matter, 2008, 20, 204125.

• Klapp and Schoen 2002

S. H. L. Klapp and M. Schoen, J. Chem. Phys., 2002, 117, 8050.

• Ouyang et al. 2011

W.-Z. Ouyang, S.-H. Xu and Z.-W. Sun, J. Chem. Phys., 2011, 134, 014901.

• Geiger and Klapp 2013

R. Geiger and S. H. L. Klapp, J. Mod. Phys., 2013, 04, 401.

• van Teeffelen and Löwen 2008

S. van Teeffelen and H. Löwen, Phys. Rev. E, 2008, 78, 020101.

• Liao and Klapp 2018

G.-J. Liao and S. H. L. Klapp, Soft Matter, 2018, 14, 7873.

• Kloeden and Platen 1992

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 $\tilde{\Phi}=N\pi\sigma^{2}/\left(4L^{2}\right)$. Therefore, $\Phi=\tilde{\Phi}\left(\sigma_{eff}/\sigma\right)^{2}$. In this work, we consider three different density values $\tilde{\Phi}=0.1$, $0.2$, and $0.5$. Consequently, the exact area fractions to the fifth decimal place are $\Phi=0.11632$, $0.23264$, and $0.58160$.

• Cates and Tailleur 2015

M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219.

• Blaschke et al. 2016

J. Blaschke, M. Maurer, K. Menon, A. Zöttl and H. Stark, Soft Matter, 2016, 12, 9821.

• Tailleur and Cates 2008

J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2008, 100, 218103.

• Bialké et al. 2013

J. Bialké, H. Löwen and T. Speck, EPL, 2013, 103, 30008.

• Digregorio et al. 2018

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 $\Delta r=0.0125\sigma$. The resulting exact values for $r_{L}$ are $r_{L}=1.15625\sigma$ for $\Phi=0.1$, $r_{L}=1.13125\sigma$ for $\Phi=0.3$, and $r_{L}=1.11875\sigma$ for $\Phi=0.5$.

• Weis and Levesque 2005

J.-J. Weis and D. Levesque, in Advanced Computer Simulation Approaches for Soft Matter Sciences II, Springer, Berlin, Heidelberg, 2005, pp. 163–225.

• Weis 2005

J.-J. Weis, J. Chem. Phys., 2005, 123, 044503.

• Trasca and Klapp 2008

R. A. Trasca and S. H. L. Klapp, J. Chem. Phys., 2008, 129, 084702.

• Weis 2002

J. J. Weis, Mol. Phys., 2002, 100, 579.

• Weis 2003

J.-J. Weis, J. Phys.: Condens. Matter, 2003, 15, S1471.

• Weis and Levesque 1993

J. J. Weis and D. Levesque, Phys. Rev. Lett., 1993, 71, 2729.

• Rovigatti et al. 2011

L. Rovigatti, J. Russo and F. Sciortino, Phys. Rev. Lett., 2011, 107, 237801.

• Rovigatti et al. 2012

L. Rovigatti, J. Russo and F. Sciortino, Soft Matter, 2012, 8, 6310.

• Rovigatti et al. 2013

L. Rovigatti, S. Kantorovich, A. O. Ivanov, J. M. Tavares and F. Sciortino, J. Chem. Phys., 2013, 139, 134901.

• Kantorovich et al. 2013

S. Kantorovich, A. O. Ivanov, L. Rovigatti, J. M. Tavares and F. Sciortino, Phys. Rev. Lett., 2013, 110, 148306.

• Kantorovich and Ivanov 2016

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.

• Ronti et al. 2017

M. Ronti, L. Rovigatti, J. M. Tavares, A. O. Ivanov, S. S. Kantorovich and F. Sciortino, Soft Matter, 2017, 13, 7870.

• Camp 2018

P. J. Camp, in Modern Problems of Molecular Physics, ed. L. A. Bulavin and A. V. Chalyi, Springer, Cham, 2018, pp. 185–204.

• Peroukidis and Klapp 2016

S. D. Peroukidis and S. H. L. Klapp, Soft Matter, 2016, 12, 6841.

• Schmidle et al. 2012

H. Schmidle, C. K. Hall, O. D. Velev and S. H. L. Klapp, Soft Matter, 2012, 8, 1521.

• Schmidle et al. 2013

H. Schmidle, S. Jäger, C. K. Hall, O. D. Velev and S. H. L. Klapp, Soft Matter, 2013, 9, 2518.

• Theurkauff et al. 2012

I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2012, 108, 268303.

• Kantorovich et al. 2015

S. S. Kantorovich, A. O. Ivanov, L. Rovigatti, J. M. Tavares and F. Sciortino, Phys. Chem. Chem. Phys., 2015, 17, 16601.

• Stenhammar et al. 2013

J. Stenhammar, A. Tiribocchi, R. J. Allen, D. Marenduzzo and M. E. Cates, Phys. Rev. Lett., 2013, 111, 145702.

• Stenhammar et al. 2014

J. Stenhammar, D. Marenduzzo, R. J. Allen and M. E. Cates, Soft Matter, 2014, 10, 1489.

• Luo and Klapp 2009

L. Luo and S. H. L. Klapp, J. Chem. Phys., 2009, 131, 034709.

• Baglietto et al. 2012

G. Baglietto, E. V. Albano and J. Candia, Interface Focus, 2012, 2, 708.

• Mahault et al. 2019

B. Mahault, F. Ginelli and H. Chaté, 2019, arXiv: 1908.03794v1 [cond-mat.stat-mech].

• Yoshinaga and Liverpool 2017

N. Yoshinaga and T. B. Liverpool, Phys. Rev. E, 2017, 96, 020603.

• Yoshinaga and Liverpool 2018

N. Yoshinaga and T. B. Liverpool, Eur. Phys. J. E, 2018, 41, 76.

• Matas Navarro et al. 2014

R. Matas Navarro, R. Golestanian, T. B. Liverpool and S. M. Fielding, Phys. Rev. E, 2014, 90, 032304.

• Matas Navarro and Fielding 2015

R. Matas Navarro and S. M. Fielding, Soft Matter, 2015, 11, 7525.

• Theers et al. 2018

M. Theers, E. Westphal, K. Qi, R. G. Winkler and G. Gompper, Soft Matter, 2018, 14, 8590–8603.

• Schwarzendahl and Mazza 2019

F. J. Schwarzendahl and M. G. Mazza, J. Chem. Phys., 2019, 150, 184902.

• Alarcón and Pagonabarraga 2013

F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56.

• Delmotte et al. 2015

B. Delmotte, E. E. Keaveny, F. Plouraboué and E. Climent, J. Comput. Phys., 2015, 302, 524.

• Hoell et al. 2018

C. Hoell, H. Löwen and A. M. Menzel, J. Chem. Phys., 2018, 149, 144902.

• Lauga et al. 2006

E. Lauga, W. R. DiLuzio, G. M. Whitesides and H. A. Stone, Biophys. J., 2006, 90, 400.

• Loose and Mitchison 2014

M. Loose and T. J. Mitchison, Nat. Cell Biol., 2014, 16, 38.

• Kümmel et al. 2013

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.

• Baraban et al. 2008

L. Baraban, D. Makarov, M. Albrecht, N. Rivier, P. Leiderer and A. Erbe, Phys. Rev. E, 2008, 77, 031407.

• Steinbach et al. 2016

G. Steinbach, S. Gemming and A. Erbe, Eur. Phys. J. E, 2016, 39, 69.

• Yener and Klapp 2016

A. B. Yener and S. H. L. Klapp, Soft Matter, 2016, 12, 2066.

• Klapp 2016

S. H. Klapp, Curr. Opin. Colloid Interface Sci., 2016, 21, 76.

• Mazars 2011

M. Mazars, Phys. Rep., 2011, 500, 43.

• Schoen and Klapp 2007

M. Schoen and S. H. L. Klapp, in Reviews in Computational Chemistry, John Wiley & Sons, Inc., 2007, pp. 301–340.

• Fily et al. 2014

Y. Fily, S. Henkes and M. Cristina Marchetti, Soft Matter, 2014, 10, 2132.