Many-body dynamics of chemically propelled nanomotors

TL;DR

This study explores the collective dynamics of chemically propelled nanomotors in confined environments, revealing phase segregation and fluctuations driven by chemical, hydrodynamic, and thermal interactions.

Contribution

It introduces a microscopic dynamical model for many-body nanomotor systems, highlighting the roles of chemical reactions, hydrodynamics, and interactions in their collective behavior.

Findings

Segregation into high and low density phases observed.

Strong fluctuations in homogeneous states identified.

Chemical and hydrodynamic interactions significantly influence dynamics.

Abstract

The collective behavior of chemically propelled sphere-dimer motors made from linked catalytic and noncatalytic spheres in a quasi-two-dimensional confined geometry is studied using a coarse-grained microscopic dynamical model. Chemical reactions at the catalytic spheres that convert fuel to product generate forces that couple to solvent degrees of freedom as a consequence of momentum conservation in the microscopic dynamics. The collective behavior of the many-body system is influenced by direct intermolecular interactions among the motors, chemotactic effects due to chemical gradients, hydrodynamic coupling, and thermal noise. Segregation into high and low density phases and globally homogeneous states with strong fluctuations are investigated as functions of the motor characteristics. Factors contributing to this behavior are discussed in the context of active Brownian models.

Tables

Table 1. Table 1: Properties of forward-moving (f-m) and backward-moving (b-m) single sphere-dimer motors for two different solvent conditions characterized by the values of τ MPC subscript 𝜏 MPC \tau_{\text{MPC}} .

	${\tau }_{\text{MPC}}$	$V_z$	$D_0$	${\tau }_{\text{D}}$	${\tau }_{\text{R}}$	${\tau }_{\text{B}}$	Pe	${\text{Pe}}^{\prime }$
		$\times {10}^3$	$\times {10}^3$
f-m	0.5	26	2	8241	3305	160	21	52
f-m	0.1	8.1	0.78	22187	13677	514	27	43
b-m	0.5	-11	2	8241	7362	378	19	22
b-m	0.1	-2.9	0.78	22187	16790	1434	12	15

Table 2. Table 2: Mean and standard deviation of scaled propulsion velocity V z ′ = V z / k B ​ T / M m superscript subscript 𝑉 𝑧 ′ subscript 𝑉 𝑧 subscript 𝑘 B 𝑇 subscript 𝑀 m V_{z}^{\prime}=V_{z}/\sqrt{k_{\text{B}}T/M_{\text{m}}} , effective diffusion coefficient D m subscript 𝐷 m D_{\text{m}} , and reorientation time τ R subscript 𝜏 R \tau_{\text{R}} for backward-moving motors, and diffusion coefficient D 0 subscript 𝐷 0 D_{0} and reorientation time τ R 0 superscript subscript 𝜏 R 0 \tau_{\text{R}}^{0} for inactive dimers, for different area fractions ϕ italic-ϕ \phi .

$\varphi$	$⟨V_z^{\prime }⟩$	${\sigma }_{V_z^{\prime }}$	$D_{\text{m}}$	$D_0$	${\tau }_{\text{R}}$	${\tau }_{\text{R}}^0$
				$\times {10}^3$
0	-1.42	0.99	0.54	2.0	7362	5691
0.05	-1.33	1.02	0.30	1.8	5542	5792
0.1	-1.23	1.05	0.22	1.7	4759	6020
0.15	-1.11	1.07	0.17	1.5	4393	6312
0.2	-0.99	1.09	0.13	1.3	4264	6761
0.25	-0.86	1.09	0.10	1.2	4313	7465
0.3	-0.71	1.09	0.07	1.0	4558	8593

Table 3. Table 3: Properties of forward-moving (f-m) and backward-moving (b-m) single sphere-dimer motors for full microscopic model (top) and for two modified microscopic models (bottom) that neglect hydrodynamics using a global Andersen thermostat.

	$d_{\text{C}}$	$d_{\text{N}}$	${\tau }_{\text{MPC}}$	${\tau }_{\text{AT}}$	${ϵ}_{\text{NA}}$	${ϵ}_{\text{NB}}$	$F_{\text{m}}$	$V_z$	$D_0$	${\tau }_{\text{D}}$	${\tau }_{\text{R}}$	${\tau }_{\text{B}}$	Pe	${\text{Pe}}^{\prime }$
							$\times {10}^3$	$\times {10}^3$
f-m	4	8	0.5	–	1	0.1	–	26	2	8639	3305	159	20	54
b-m	4	8	0.5	–	0.1	1	–	-21	2	8639	8020	194	41	44
f-m	8	8	0.5	–	1	0.1	–	25	1.5	16817	15517	203	76	82
b-m	8	8	0.5	–	0.1	1	–	-19	1.5	16817	17073	266	63	62
f-m	4	8	–	0.5	1	0.1	–	0.63	0.028	618576	165801	6649	24	93
b-m	4	8	–	0.5	0.1	1	–	-0.78	0.028	618576	135084	5357	25	115
f-m	4	8	–	0.5	1	–	4	0.7	0.028	618576	167372	5933	28	104
b-m	4	8	–	0.5	1	–	-4	-0.7	0.028	618576	152142	5961	25	103

Equations

$$V_z=\frac{k_{\text{B}}T}{\eta }\Lambda Q,$$

$$\Lambda ={\int }_0^{\infty }r[e^{-\beta V_{\text{NB}}(r)}-e^{-\beta V_{\text{NA}}(r)}]dr.$$

$$g(r)=\frac{1}{2\pi rn_{\text{N}}}\Big{\langle}\sum_{j<i=1}^{N_{\text{N}}}\delta(|{\bf r}_{\text{N}ij}|-r)\Big{\rangle}\;,$$

$$f(r,\theta)=\frac{1}{N_{\text{N}}}\Big{\langle}\sum_{j<i=1}^{N_{\text{N}}}\delta(\theta_{ij}-\theta)\delta(|{\bf r}_{\text{N}ij}|-r)\Big{\rangle}\;,$$

$$\mathcal{K}=1+2\pi {\rho }_{\text{m}}{\int }_0^{\infty }drr\left(g(r)-1\right),$$

$${\mathcal{K}}_{\text{HS}}={\left(1+\frac{8\varphi -2{\varphi }^2}{(1-\varphi {)}^4}\right)}^{-1}.$$

$$\Delta R^2(t)$$

$$\frac{d}{dt}{\mathbf{R}}_i$$

$$\frac{d}{dt}{\hat{\mathbf{z}}}_i$$

$$\frac{d}{dt}{\mathbf{R}}_{i\nu }=\mu ({\mathbf{F}}_{i\nu }^{\text{p}}+{\mathbf{F}}_{i\nu })+{\mathbf{f}}_{i\nu }^{\text{t}},$$

$$D_{\text{m}}(\varphi )=D_0(\varphi )+\frac{1}{2}|V_z(\varphi ){|}^2{\tau }_{\text{R}}(\varphi ).$$

$$\frac{d}{dt}{\mathbf{R}}_i$$

$$\frac{d}{dt}{\hat{\mathbf{z}}}_i$$

Full text

I Introduction

Self-organization and pattern formation are processes that occur in systems under both equilibrium and nonequilibrium conditions. The principles underlying the dynamics leading to the final structures are well understood, although in specific applications the evolution equations and the final inhomogeneous states may be very complicated. Cross and Hohenberg (1993); Hoyle (2006); Pismen (2006); Desai and Kapral (2009) In equilibrium systems, the phase ordering or segregation processes are often described in terms of order parameter fields, satisfying certain symmetries, whose evolution is controlled by equations of motion based on free energy functionals. In nonequilibrium systems, in general, no such free energy functionals exist, and the order parameter equations are typically derived from a nonlinear analysis. In this case, the variety of possible inhomogeneous states is large and their classification is more difficult.

Active systems, where the constituent objects either are forced by external means or move autonomously, are another class of nonequilibrium systems that self-organize in ways that differ from equilibrium and nonequilibrium systems with inactive elements. Systems of this type take many forms including shaken granular materials, a wide rage of biological systems, and synthetic active media. Although active systems have been intensively studied only recently, a considerable literature exists, documented in reviews Lauga and Powers (2009); Marchetti et al. (2013); Elgeti et al. (2015); Zöttl and Stark (2016); Bechinger et al. (2016), which describes the new phenomena that these systems exhibit.

The models that are used to describe the collective motions of active systems range from simple particle-based models Vicsek et al. (1995); Chaté et al. (2008); Peruani et al. (2006) to models based on continuum order-parameter equations Toner and Tu (1995); Simha and Ramaswamy (2002); Toner et al. (2005). Some of the simplest models deal with active Brownian particles that have a prescribed velocity, interact through hard or soft potentials, and experience translational and orientational Brownian motion. Fily and Marchetti (2012); Redner et al. (2013); Cates and Tailleur (2013); Stenhammar et al. (2013); Bialké et al. (2013); Wysocki et al. (2014); Fily et al. (2014); Speck et al. (2014); Stenhammar et al. (2014); Takatori and Brady (2015); Speck et al. (2015) Even with these minimal ingredients, the systems separate into high and low density phases for high enough volume or area fractions, and such dynamics has been mapped onto continuum order-parameter equations. Stenhammar et al. (2013); Speck et al. (2014)

In many systems where active particles move in a fluid environment, particle motion induces fluid flows in the surrounding medium that lead to hydrodynamic coupling among the particles. Investigations have shown that such hydrodynamic coupling can give rise to coherent collective motion. Saintillan and Shelley (2012); Zöttl and Stark (2014); Oyama et al. (2016)

The collective dynamics of self-propelled motors whose motion is governed by diffusiophoretic mechanisms Anderson (1989); Golestanian et al. (2005); Kapral (2013) is considered in this paper. Such chemically powered motors operate under nonequilibrium conditions, and catalytic chemical reactions on the motor generate concentration gradients of the reactive species that are responsible for motor motion. The concentration and fluid velocity fields that arise from the activity of all motors in the system influence the nature of the collective dynamics. Experimental studies of the collective behavior of chemically powered active colloidal particles have shown that these systems display clustering, schooling, and other collective states. Ibele et al. (2009); Theurkauff et al. (2012); Palacci et al. (2013); Buttinoni et al. (2013); Wang et al. (2015) Langevin models without hydrodynamic interactions have been developed and used to construct phase diagrams that show the possible types of collective behavior that such systems can exhibit. Pohl and Stark (2014); Saha et al. (2014); Pohl and Stark (2015)

We investigate the collective dynamics of sphere-dimer motors using a coarse-grained microscopic dynamical method that accounts for coupling through many-body hydrodynamic interactions, concentration gradients, and direct potential interactions among motors, as well as thermal fluctuations. Most studies of the collective dynamics of chemically powered motors are carried out on spherical Janus particles with catalytic and noncatalytic hemispherical faces. The sphere-dimer motors Rückner and Kapral (2007); Tao and Kapral (2008); Valadares et al. (2010); Yang et al. (2014); Reigh and Kapral (2015) we study here have a different shape, are made from linked catalytic and noncatalytic spheres, and are propelled by a diffusiophoretic mechanism. Simulations of small numbers of sphere-dimer motors showed that they self-assemble into transient clusters and display schooling behavior. Thakur and Kapral (2012); Colberg et al. (2014) Here we consider systems with thousands of motors and quantitatively characterize the collective dynamics. Because our particle-based dynamics describes all interactions and chemical reactions, and conserves mass, momentum, and energy, it captures the full many-body structures of the hydrodynamic velocity and concentration fields without resorting to simplifying approximations, such as the neglect of correlations in mean field models. It also accounts for the modification of the propulsion properties of individual motors due to chemical gradients arising from other motors in the system. As a result, we are able to probe aspects of the collective behavior that are not accessible using other models.

Sec. II contains a description of the coarse-grained microscopic model for the chemically powered sphere-dimer motor and the solvent in which it resides. The properties of a single sphere-dimer motor in solution are given in Sec. III. This section also summarizes the continuum description of this motor. Simulation results on the collective dynamics of large ensembles of motors are presented and analyzed in Sec. IV. The results are analyzed further in Sec. V where connections with active Brownian particle models are made. The conclusions of the paper are in Sec. VI.

II Microscopic model for collective motor dynamics

The coarse-grained particle-based model for the system includes both the motors and the multi-component fluid environment in which they move. The evolution of the entire system is carried out using a hybrid scheme that combines molecular dynamics (MD) and multiparticle collision dynamics (MPCD). Malevanets and Kapral (1999); *malevanets2000 The system is contained in a slab whose height is fixed at $L_{Z}=30$ and whose edge lengths $L_{X}=L_{Y}=L$ vary as specified below. The solvent consists of point-like fuel (A) and product (B) molecules with common mass $m_{\text{s}}=1$ and total number density $\varrho_{\text{s}}=9$. The solvent molecules undergo bounce-back collisions with the walls at $Z=0$ and $L_{Z}$ that reverse their velocities. Periodic boundary conditions are applied in the $X$ and $Y$ directions.

The catalytic (C) and noncatalytic (N) spheres in the dimer motors have effective diameters $d_{\text{C}}=4$ and $d_{\text{N}}=8$, while the dimer mass is $M_{\text{m}}=\frac{\pi}{6}\varrho_{\text{s}}m_{\text{s}}(d_{\text{C}}^{3}+d_{\text{N}}^{3})\approx 2714$, which makes it neutrally buoyant. In order to better visualize the collective dimer dynamics, their motion in the $Z$ direction was suppressed by using two long-range 9–3 Lennard-Jones wall potentials, V_{\text{w}S}(\zeta)=(3\sqrt{3}/2)\epsilon_{\text{w}S}\bigl{[}(\sigma_{\text{w}S}/\zeta)^{9}-(\sigma_{\text{w}S}/\zeta)^{3}\bigr{]}, where $\zeta$ is the distance of a sphere to the wall, and the interaction parameters are $\epsilon_{\text{w}S}=5$ and $\sigma_{\text{w}S}=L_{Z}/2$. These potentials largely confine the dimer dynamics to the $XY$ midplane of the slab, although the fluid flow and concentration fields are three dimensional. For a given number of motors $N_{\text{m}}$ and area fraction $\phi=N_{\text{m}}A_{\text{m}}/A$, where $A_{\text{m}}=\frac{\pi}{4}(d_{\text{C}}^{2}+d_{\text{N}}^{2})$ is the area of the dimer projected on the $XY$ plane, and $A=L^{2}$, the edge lengths of the midplane are chosen as L=\big{\lfloor}\sqrt{A_{\text{m}}N_{\text{m}}/\phi}\big{\rfloor}. (The projected area could be modeled as that of a cylindrical representation of the dimer and this would yield somewhat higher area fractions.)

Dimer spheres and solvent molecules, as well as spheres of different dimers, interact via shifted, truncated 12–6 Lennard-Jones potentials of the form V_{ij}(r)=\epsilon_{ij}\bigl{\{}4\bigl{[}(\sigma_{ij}/r)^{12}-(\sigma_{ij}/r)^{6}\bigr{]}+1\bigr{\}} for $r<\sqrt[6]{2}\,\sigma_{ij}$ and zero otherwise. The separation distances are $\sigma_{\text{CA}}=\sigma_{\text{CB}}=d_{\text{C}}/2$ and $\sigma_{\text{NA}}=\sigma_{\text{NB}}=d_{\text{N}}/2$ for sphere–solvent pairs, and $\sigma_{\text{CC}}=d_{\text{C}}+1$, $\sigma_{\text{CN}}=(d_{\text{C}}+d_{\text{N}})/2+1$, and $\sigma_{\text{NN}}=d_{\text{N}}+1$ for sphere–sphere pairs. The interaction energies are $\epsilon_{\text{CA}}=\epsilon_{\text{CB}}=\epsilon_{\text{NA}}=1$ and $\epsilon_{\text{NB}}=0.1$ or 10 for sphere–solvent pairs and $\epsilon_{\text{CC}}=\epsilon_{\text{CN}}=\epsilon_{\text{NN}}=10$ for sphere–sphere pairs. The dimer bond length is $d_{\text{CN}}=\sqrt[6]{2}\,(d_{\text{C}}+d_{\text{N}})/2$, which yields the maximum possible propulsion velocity Reigh and Kapral (2015) at the minimum possible distance that still ensures conservation of energy in the presence of chemical reactions at the catalytic sphere. Likewise, the sphere–sphere pair separation distances and interaction energies given above are large enough to ensure conservation of energy in the presence of chemical reactions and avoid the occurrence of solvent depletion forces between dimers. 111The simulation code explicitly asserts at every MD step that any solvent molecule interacts with at most one dimer sphere, which, besides satisfying the sphere–sphere separation criteria, allows for a performance optimization in the parallel computation of sphere–solvent forces. There are no intermolecular potentials among solvent molecules; these interactions are instead described by multiparticle collision dynamics.

In hybrid MD–MPCD, the solvent particles undergo multiparticle collisions at discrete time intervals $\tau_{\text{MPC}}$. As described in detail elsewhere Malevanets and Kapral (1999); *malevanets2000; Kapral (2008); Gompper et al. (2009), to carry out multiparticle collisions, at discrete times the system is partitioned into collision cells, and rotation operators are assigned to the cells. Post-collision particle velocities are then obtained in each cell by rotating the particle velocity relative to the cell center of mass velocity and adding back the center of mass velocity. [MPCDcollisionsinacellwerecarriedoutusingrotationsby$π/2$aboutarandomlychosenaxis.Galileaninvariancewasguaranteedbyimplementinggridshiftingwhencarryingoutthecollisions;]ihle2001; *ihle2003 Between these multiparticle collisions, all particles in the system evolve by Newton’s equations of motion in the potentials given above. Since MPCD conserves mass, momentum, and energy, one may derive the Navier-Stokes equations for the solvent on long distance and time scales with known values of the transport coefficients. Malevanets and Kapral (1999); *malevanets2000; Ihle and Kroll (2001); *ihle2003; Kapral (2008); Gompper et al. (2009) The MPCD collision cell size is $a=1$, and the solvent number density is $n=9$. The velocity-Verlet time step is $\tau_{\text{MD}}=0.01$, and the thermal energy of the system is $k_{\text{B}}T=1/6$. Given these parameters, two values of the multiparticle collision time $\tau_{\text{MPC}}$ are considered. For $\tau_{\text{MPC}}=0.5$, the solvent viscosity is $\eta=1.22$, the A and B common diffusion coefficient is $D_{\text{A}}=0.099$ and the Schmidt number is $\text{Sc}=1.4$. For $\tau_{\text{MPC}}=0.1$, the solvent viscosity is $\eta=4.51$, the A and B common diffusion coefficient is $D_{\text{A}}=0.020$ and the Schmidt number is $\text{Sc}=25$. We note that for our collision model, the solvent molecules interact with the motors through intermolecular potentials so that the fluid particle diffusion does not enter in the Navier-Stokes equation. In this circumstance, the fluid will still exhibit liquid-like properties in spite of the relatively low values of the Schmidt numbers, which are a consequence of coarse-graining to achieve computational efficiency. Padding and Louis (2006)

The catalytic chemical reactions that drive the directed motion are carried out by converting fuel to product in irreversible [Whileirreversiblereactionsareconsideredhere; thereactivedynamicsmaybegeneralizedtotreatreversiblereactions.See; forexample; ]thakur2011 reactions $\text{A}+\text{C}\rightarrow\text{B}+\text{C}$ with unit probability in the vicinity of a catalytic sphere, $r_{\text{CA}}<\sqrt[6]{2}\,\sigma_{\text{CA}}$. Given this reactive dynamics, the rate constant for reactions on a single catalytic sphere, $k$, may be approximately determined Tucci and Kapral (2004) and is given by $k^{-1}=k_{0}^{-1}+k_{D}^{-1}$, where $k_{0}=p_{0}r_{\text{C}}^{2}(8\pi k_{\text{B}}T/m)^{1/2}$ is the intrinsic reaction rate coefficient and $k_{D}=4\pi D_{\text{A}}r_{\text{C}}$ is the Smoluchowski diffusion-controlled rate coefficient, where the reaction probability is $p_{0}=1$, and $r_{\text{C}}$ is the radius of the catalytic sphere. For our parameter values, $k_{0}=8.19$, and for $\tau_{\text{MPC}}=0.5$, we have $k_{D}=2.48$ corresponding to diffusion-influenced kinetics, while for $\tau_{\text{MPC}}=0.1$, we have $k_{D}=0.50$ so that the reaction is diffusion-controlled.

In order to maintain the system in a steady state, irreversible reactions $\text{B}\stackrel{{\scriptstyle k_{2}}}{{\rightarrow}}\text{A}$ occur locally in the bulk of the solution with rate coefficient $k_{2}=0.01$ outside the sphere–solvent interaction zones $r_{S\text{B}}>\sqrt[6]{2}\,\sigma_{S\text{B}}$ for all spheres $S$. For this purpose, the reactive version of MPCD was employed. Rohlf et al. (2008) On long distance and time scales, this reactive dynamics yields the reaction-diffusion equations.

All parameter values listed above, as well as the results that follow, are reported in dimensionless units with mass in units of $m$, energy in units of $\epsilon$, distances in units of $\sigma$, and time in units of $t_{0}=\sqrt{m\sigma^{2}/\epsilon}$.

III Single sphere-dimer motor

It is useful to summarize the properties of a single motor in solution before studying the collective dynamics of many motors. 222The simulations of single motors were carried out using systems of size $L_{X}=L_{Y}=60$ and $L_{Z}=30$, and averaging was performed over 60 realizations of $10^{8}$ MD steps each. The majority of the results presented below are for two values of $\epsilon_{\text{NB}}$. For $\epsilon_{\text{NB}}=0.1$, the sphere dimer moves in a direction with the catalytic sphere at its head, which for convenience we shall refer to as a forward-moving motor (Fig. 1a), while for $\epsilon_{\text{NB}}=10$, the motor has the noncatalytic sphere at its head and will be called a backward-moving motor (Fig. 1b).

Letting $\hat{\bf z}$ be a unit vector directed along the dimer bond from the N to C spheres, the projection of the motor velocity along $\hat{\bf z}$ will be denoted by $V_{z}$. The motor reorientation time, $\tau_{\text{R}}$, is determined from the decay of the orientation correlation function $\langle\hat{\bf z}(t)\cdot\hat{\bf z}\rangle$. The time it takes a motor to move a distance equal to its effective radius, $R_{\text{m}}=\sqrt[3]{d_{\text{C}}^{3}+d_{\text{N}}^{3}}/2=4.16$, is the ballistic time, $\tau_{\text{B}}=R_{\text{m}}/|V_{z}|$, while the time it takes the motor to diffuse this distance is the diffusion time, $\tau_{\text{D}}=R_{\text{m}}^{2}/D_{0}$, where $D_{0}$ is the diffusion coefficient of an inactive dimer. Péclet numbers can be defined as the ratio of the diffusion and ballistic times, $\text{Pe}^{\prime}=\tau_{\text{D}}/\tau_{\text{B}}$, and the ratio of the orientational and ballistic times, $\text{Pe}=\tau_{\text{R}}/\tau_{\text{B}}$. The values of these quantities can be found in Table 1.

The continuum theory of the sphere-dimer motor is most conveniently formulated in a bispherical coordinate system. Stimson and Jeffery (1926); Morse and Feshbach (1953); Happel and Brenner (1973); Reigh and Kapral (2015); Michelin and Lauga (2015); Popescu et al. (2011) The concentration fields can be found by applying a radiation boundary condition involving $k_{0}$ on the catalytic sphere and a reflecting boundary condition on the noncatalytic sphere. The motor velocity has a functional form typical for phoretic propulsion,

[TABLE]

where

[TABLE]

Here $\beta=1/(k_{\text{B}}T)$, $V_{\text{NB}}$ and $V_{\text{NA}}$ denote the interaction potentials of the B and A species with the noncatalytic sphere, respectively, while the factor $Q$ depends on the system parameters and is determined from the solutions of the reaction-diffusion and Stokes equations in bispherical coordinates. The explicit form for the motor velocity is given in Eq. (8) of Ref. Reigh and Kapral, 2015 from which the factor $Q$ can be found. The continuum solutions have been compared with microscopic simulations. In particular, experiment Valadares et al. (2010), continuum theory, and simulation Reigh and Kapral (2015); Tao and Kapral (2008) indicate that the maximum motor velocity occurs near a 1:2 ratio of catalytic to noncatalytic sphere diameters, and this has motivated the radius ratio chosen for most of the work in this study, although results for sphere-dimer motors made from spheres with equal sizes will also be given. The continuum results for the fluid velocity fields may also be obtained analytically and depend on the dimer bond length. For the bond length chosen in this study, the forward-motor far-field flow decays as $1/r^{2}$, characteristic of a force dipole, and the flow pattern is that of a “puller” where the fluid flow is inward from the front and rear of the motor and outward on its sides (cf. Fig. 4 of Ref. Reigh and Kapral, 2015). The near-field flow pattern is more complex with fluid circulation that changes the puller character. The flow pattern is reversed for the backward motor with a far-field flow characteristic of a “pusher”.

IV Collective dynamics

In this section, we describe the dynamics at selected points in the system parameter space. The nature of the collective behavior depends on the values of the microscopic parameters that define the system: intermolecular potentials, MPCD parameters for the solvent, and motor density. Once these parameters and system conditions are specified, all properties of the system may be determined from the solution of the evolution equations. [ThesimulationswereperformedusingacodewritteninOpenCL~CandLuathatispublishedunderafreesoftwarelicense;]colberg2013 Unless stated otherwise, results will be reported for $\tau_{\text{MPC}}=0.5$. It is useful to organize the discussion of the collective dynamics based on whether self-propulsion leads to forward or backward motor motion.

Forward-moving motors

Consider interaction energies $(\epsilon_{\text{NB}}=0.1,\epsilon_{\text{NA}}=1)$ where the parameter $\Lambda=1.18>0$ (see Eq. (2)) and the motor moves in the forward direction with the C sphere at its head. Since the catalytic reaction produces a high concentration of species B and depletes that of A near the portion of the N sphere surface closest to the C sphere, for these potential parameters, a net force will act to move the motor in the forward direction. It also follows from this that a motor will respond to the concentration gradients of other motors by moving in the direction of high B concentration, i.e., forward-moving motors are chemotactically attracted to each other. An individual motor will respond to both its self-generated concentration gradient and that due to other motors, and forward-moving motors tend to align with their catalytic heads pointing in the same direction. Geometric factors related to the motor shape also play a role in such alignment processes since they influence how the motors interact and assemble when they are close to each other.

Hydrodynamic interactions also contribute to collective dynamics and these effects have been investigated for active swimmers. The far-field flow of pullers is attractive in the swimming direction and repulsive in the perpendicular directions, giving rise to a tendency for two swimmers to move away from each other. Lauga and Powers (2009) Simulations of suspensions of pullers do not show large-scale correlated motions Saintillan and Shelley (2012), but when near-field flows and the finite volume of swimmers are taken into account, large density fluctuations have been observed for pullers. Oyama et al. (2016); Alarcón and Pagonabarraga (2013) Recall that on the basis of the far-field flow, for our dimer bond length, the forward-moving motor can be classified as a puller, but the near-field flow is more complex.

The evolution of systems with $N_{\text{m}}=1000$ and 6000 sphere-dimer motors for several area fractions $\phi$ is shown in Fig. 2. The number of solvent particles $N_{s}$ in the system ranges from $N_{s}\approx 10^{7}$ to $10^{8}$ depending on the area fraction. The initial state of the system was prepared by randomly placing the dimers with random orientation and without overlap such that $V_{SS^{\prime}}(r)=0$ for spheres $S$ and $S^{\prime}$ and with the solvent comprising only A molecules. Segregation into low-density gas-like and high-density disordered solid-like phases was observed, with the time scale for coarsening dependent on the area fraction. A video showing the evolution of clusters can be found in the supplementary material. For the area fractions considered, dimer motors propagate with a strong ballistic component between encounters. Cluster growth occurs through single dimer attachment to a cluster and by cluster–cluster aggregation events involving propagating clusters. Due to the presence of thermal noise, clusters evaporate, fragment, and combine, although in the late stages where large clusters exist, the gas phase is very dilute and the time scales for cluster breakup and formation are very long.

The average number of dimers in the largest cluster, $N_{\xi}$, is plotted as a function of time for different values of $\phi$ in Fig. 3a. For each area fraction, we compare the results for two different values of the total number of motors in the simulation volume to gauge the magnitude of finite-size effects. The log–log scale plots of $N_{\xi}$ indicate power-law scaling of the cluster growth that is dependent on $\phi$. In the initial stage of the growth, for $t\leq 10^{4}$, $N_{\xi}(t)\propto t^{\alpha}$. For $\phi=0.1$ and 0.2, the initial regime does not exhibit finite-size effects; however, for $\phi=0.3$, one can see that the initial regime for $N_{\text{m}}=1000$ is suppressed compared to that for $N_{\text{m}}=6000$. The cluster structure in this regime can be seen in the panels in the left-most column in Fig. 2. The exponent $\alpha(\phi)$ was fitted for systems with $N_{\text{m}}=1000$ and increases with $\phi$ (see Fig. 3b). As above, finite-size effects may alter the exponent values for $\phi\geq 0.2$. The simulations results for 4000 and 6000 dimers suggest the presence of another power-law regime with smaller exponents at longer times, but more extensive simulations of even larger systems are needed to verify its presence and quantitatively characterize its behavior.

The temporal variation of the average number of motors in the gas phase, $N_{\text{g}}$, also displays power-law behavior, $N_{\text{g}}(t)\propto t^{\nu}$, in the initial time regime. The value of the exponent $\nu(\phi)$ decreases with increasing $\phi$ as shown in Fig. 3b. This slow decay in the short-time regime is followed by a rapid decay to very small values since the gas phase is very dilute.

Since solid-like clusters are quickly formed, the longtime dynamics of individual motors is determined by diffusive motions of the clusters containing these motors. This diffusive cluster dynamics occurs on very long time scales and requires significantly longer simulation times than those considered in this study in order to determine the diffusion coefficients accurately.

In the steady state, large metastable clusters persist for very long periods of time. The structure of these disordered solid-like clusters is determined by the interactions among motors, dependent on the sphere-dimer geometry, and the chemotactic interactions among motors. Fig. 4 shows expanded views of the interior and periphery of a cluster. There is a high degree of orientational order on the periphery with motors pointing with their catalytic heads towards the bulk of the cluster. This is a consequence of the attraction of these forward-moving motors to regions of high product (B) concentration, which is concentrated in the cluster vicinity, along with geometric effects that arise from the packing of the asymmetric dimers into clusters. The interior is globally disordered but has a high degree of local structural order. In the peripheral region, and elsewhere in the interior of the cluster, there exist square lattice arrangements of N spheres with C spheres occupying the interstitial positions. Since the C and N spheres are linked, there is a high degree of dimer orientational order. Other prominent types of local ordering can be seen in the expanded views. In particular, there are configurations in the cluster interior where dimers are oriented so that their C heads are aligned in order to partially surround an N sphere.

These qualitative observations are reflected in the forms of the radial and orientational correlations. The NN radial distribution function is defined as

[TABLE]

with ${\bf r}_{\text{N}ij}={\bf r}_{\text{N}i}-{\bf r}_{\text{N}j}$, where ${\bf r}_{\text{N}i}$ is the position of the noncatalytic monomer $i$, $N_{\text{N}}$ is the number of noncatalytic monomers, and $n_{\text{N}}=N_{\text{N}}/V$. The angle brackets denote an average over time and realizations. This function is plotted in Fig. 5a.

There is a very prominent peak at $r=9.9$ that corresponds to the NN packing discussed above, with smaller-amplitude weakly damped oscillations that persist for long distances reflecting the approximate long-range order. Fine structure is seen in some of the smaller-amplitude next-neighbor peaks. Fig. 5b shows the joint probability distribution $f(r,\theta)$ of N–N sphere pairs,

[TABLE]

where $\hat{{\bf z}_{i}}\cdot\hat{{\bf z}_{j}}=\cos\theta_{ij}$ and $\theta_{ij}$ is the angle between the unit bond vectors centered on noncatalytic spheres $i$ and $j$ for values of $r$ corresponding to some of the prominent peaks in $g(r)$, namely, those at $r=9.9$, 13.8, and 14.3. In particular, note that strong orientational order is associated with the peak at $f(r=9.9,\theta=40)$. This peak arises from configurations in the cluster where nearest-neighbor dimers are oriented with their heads pointing towards a central noncatalytic sphere.

Backward-moving motors

For interaction energies $(\epsilon_{\text{NB}}=10,\epsilon_{\text{NA}}=1)$, where $\Lambda=-0.60<0$, the motor will move in the backward direction with the N sphere at its head. Backward motors tend to move in a direction of lower product concentration and will chemotactically respond to other motors by moving away from them, although geometric factors play a role and aligned pairs of motors exist and move in the backward direction as a unit for short periods of time. As a result, one might expect that clustering will be less pronounced for such motors. This is what is found, but transient clustering has important effects on the collective dynamics of these motors.

Fig. 6 shows instantaneous configurations of backward-moving motors at several times. No large-scale cluster formation takes place as for forward-moving motors; however, substantial regions of inhomogeneous density are evident in the figure. (See also the video in the supplementary material.) The tendency to form small transient clusters can be quantified by considering the steady state radial distribution functions, $g(r)$, with $r$ being the magnitude of the distance between the N spheres. The distribution functions for inactive and backward-moving motors are compared in Fig. 7a. There is a strong peak in the NN distribution function at $r=10$ for backward-moving motors, while only very weak structural ordering is seen for inactive dimers.

The tendency to form small transient clusters is reflected in the $\phi$ dependence of ${\mathcal{K}}$, defined as

[TABLE]

which is the nonequilibrium analog of ${\mathcal{K}}_{\text{T}}=k_{\text{B}}T\rho_{\text{m}}\,\kappa_{\text{T}}$, where $\kappa_{\text{T}}$ is the isothermal compressibility. However, ${\mathcal{K}}$ depends on the nonequilibrium steady state radial distribution function instead of its equilibrium analog. Fig. 7b plots ${\mathcal{K}}$ as a function of $\phi$ for active and inactive dimers. For inactive dimers, ${\mathcal{K}}$ decays quickly with increasing area fraction. Its behavior is similar to that of a fluid of rigid hard spheres. Using the Carnahan-Starling equation of state Carnahan and Starling (1969) for a hard sphere fluid, ${\mathcal{K}}_{\text{HS}}$ is given by

[TABLE]

This quantity is also plotted in Fig. 7b and matches the inactive sphere dimer result. In contrast to inactive dimers, ${\mathcal{K}}$ for backward-moving motors is larger and has a convex rather than a concave shape. For equilibrium systems, $\kappa_{\text{T}}$ is related to the number fluctuations by $k_{\text{B}}T\rho_{\text{m}}\,\kappa_{\text{T}}=(\langle N^{2}\rangle-\langle N\rangle^{2})/\langle N\rangle$, and the markedly different $\phi$ dependence of this function for active dimers reflects the presence of strong density fluctuations. For low area fractions, all of these expressions tend to the ideal gas value of unity.

Interactions among motors reduce the mean propulsion velocity of a motor as seen in Fig. 8, which plots the probability distributions of the velocity projected along the propagation direction for various values of $\phi$. The distributions are Gaussian with mean values corresponding to the propulsion velocity and widths close to those for a thermal distribution of velocities. (The velocity is scaled by its thermal value, $\sqrt{k_{\text{B}}T/M_{\text{m}}}$.) The magnitude of the mean motor velocity steadily decreases and its width increases as $\phi$ increases. (See data in Table 2.)

The effects of crowding on dimer orientational dynamics are very different for active motors and inactive dimers. Fig. 9a plots the reorientation time, $\tau_{\text{R}}$, determined from the decay of the orientation autocorrelation function, for both of these cases as function of the dimer area fraction. For inactive dimers that execute thermal orientational Brownian motion, $\tau_{\text{R}}$ is an increasing function of $\phi$. This is the expected effect due to crowding that will tend to hinder reorientation. Active dimers undergo ballistic motion in a direction along the dimer bond vector. Collisions with other motors at low dimer densities will perturb the ballistic dynamics and shorten the reorientation time. However, as the dimer density increases, dimers tend to form oriented transient clusters, and this increases the reorientation time. As a result, $\tau_{\text{R}}$ is a nonmonotonic function of $\phi$. For small ensembles of forward-moving motors, $\tau_{\text{R}}$ was observed to vary in a similar nonmonotonic fashion with dimer number. Thakur and Kapral (2012)

The existence of transient clustering is also reflected in the self-diffusion coefficients of backward-moving motors. For a single motor in solution, the mean square displacement $\Delta R^{2}(t)$ takes the form

[TABLE]

where $\tau_{\text{v}}$ is the velocity relaxation time. In the ballistic regime, $t\ll\tau_{\text{v}}$, $\Delta R^{2}(t)\approx(2k_{\text{B}}T/M_{\text{m}}+V_{z}^{2})\,t^{2}$, while for long times, $\Delta R^{2}(t)\approx 4D_{\text{m}}t$ with an effective sphere-dimer diffusion coefficient of $D_{\text{m}}=D_{0}+\frac{1}{2}V_{z}^{2}\tau_{\text{R}}$, where the bare diffusion coefficient is $D_{0}=(k_{\text{B}}T/M_{\text{m}})\,\tau_{\text{v}}$. For a single backward-moving motor, the effective diffusion coefficient has a value $D_{\text{m}}=0.299$, while the bare diffusion coefficient is $D_{0}=0.002$.

Fig. 10 compares $\Delta R^{2}(t)$ versus time for inactive dimers and backward-moving motors, with the ballistic and diffusive regimes indicated by the straight line segments that intersect at the crossover time separating these regimes. For a single dimer where Eq. (7) applies, the crossover time is approximately $t_{c}\approx 2\tau_{\text{R}}$. The inertial regime for inactive dimers is given by $\Delta R^{2}(t)\approx 2k_{\text{B}}T/M_{\text{m}}\,t^{2}$, and at $\phi=0.3$, the crossover time is $t_{c}\approx 30$, which is much shorter than $2\tau_{\text{R}}(\phi)$. In addition, sub-diffusive dynamics is observed before the diffusive regime is reached. For backward-moving motors at $\phi=0.3$, the sub-diffusive regime is absent and the crossover time is $t_{c}\approx 2000$, which is still considerably smaller than $2\tau_{\text{R}}(\phi)\approx 8000$.

The effective motor diffusion coefficient as a function of the area fraction $\phi$, $D_{\text{m}}(\phi)$, was determined from the long-time behavior of $\Delta R^{2}(t)$, and the ratio $D_{\text{m}}(\phi)/D_{\text{m}}(\phi=0)$ is plotted versus $\phi$ in Fig. 9b. One observes that this ratio is smaller and decreases more rapidly for backward-moving motors than for inactive dimers. Although the diffusion coefficients of the active backward-moving motors are much larger than the corresponding diffusion coefficients of inactive dimers, the effects of crowding are more significant since the motor velocity and reorientation time, which determine the effective diffusion coefficients of the motors, are strong functions of $\phi$ (Table 2).

Other motor and system parameters

The phenomena observed for asymmetric sphere-dimer motors with $(d_{\text{N}}=8,d_{\text{C}}=4)$, for the specific solvent conditions and interaction potentials used above, are present in systems with other motor and solvent parameters. While we have not carried out a systematic study as a function of all parameters that determine the dynamics, a few examples given below will serve to illustrate the results.

The interaction parameters $(\epsilon_{\text{NB}},\epsilon_{\text{NA}})$ determining forward and backward motions were chosen to be $(0.1,1)$ and $(10,1)$, respectively, corresponding to $\Lambda=1.18$ and $\Lambda=-0.60$ in the continuum theory expression for the motor velocity in Eq. (1). We may also choose parameter pairs $(0.1,1)$ and $(1,0.1)$ so that $\Lambda$ takes the values $\Lambda=\pm 1.18$. The results for this more symmetric case again show clustering for forward motors and only strong fluctuations for backward motors. In addition, similar results, shown in Fig. 11, are found for symmetric sphere-dimer motors with $d_{\text{N}}=d_{\text{C}}=8$ (see videos in the supplementary material). The change to a symmetric motor geometry increases the tendency of backward-moving motors to align, and transient strings of motors exist and play a role in enhancing density fluctuations due to transient cluster formation. For forward-moving motors, this tendency to form transient strings is less significant since noncatalytic spheres are chemotactically attracted to catalytic spheres. The large clusters that form for forward-moving motors have hexagonal ordering with defects. Thakur and Kapral (2012)

Additional information concerning these cases can be found in Table 3.

Finally, for both forward- and backward-moving motors, the qualitative features of the collective dynamics for the higher viscosity fluid with $\tau_{\text{MPC}}=0.1$ are similar to those described above for $\tau_{\text{MPC}}=0.5$. In particular, the forward-moving motors form large clusters, although the time scale on which this cluster formation occurs is longer for $\tau_{\text{MPC}}=0.1$ due to the smaller single-motor velocity (see Table 1). Similarly, for the backward-moving motors, only transient cluster formation is observed.

V Factors contributing to collective behavior

Segregation into domains of low- and high-density phases occurs even in active systems with only repulsive interactions among the active particles. The Langevin and Fokker-Planck models for such active Brownian systems include terms to account for particle propulsion in a given direction, forces arising from direct repulsive interactions, and thermal noise. Hydrodynamic effects are not considered. The Langevin equations take the form

[TABLE]

Here ${\bf R}_{i}$ and $\hat{\bf z}_{i}$ are the position and orientation of motor $i$, the propulsion force is ${\bf F}^{\text{p}}_{i}=\zeta V_{z}\hat{\bf z}_{i}$, and the mobility is $\mu=1/\zeta$, with $\zeta$ being the friction coefficient. The random functions ${\bf f}^{\text{t}}_{i}$ and ${\bf f}^{\text{r}}_{i}$ satisfy fluctuation-dissipation relations, $\langle{\bf f}^{\text{t}}_{i}(t){\bf f}^{\text{t}}_{j}(t^{\prime})\rangle=2D_{0}\delta_{ij}{\bf 1}\delta(t-t^{\prime})$ and $\langle{\bf f}^{\text{r}}_{i}(t){\bf f}^{\text{r}}_{j}(t^{\prime})\rangle=2D_{0}^{\text{R}}\delta_{ij}{\bf 1}\delta(t-t^{\prime})$, where $D_{0}=k_{\text{B}}T/\zeta$ and $2D_{0}^{\text{R}}=1/\tau_{\text{R}}$ are the translational and rotational diffusion coefficients, respectively, for the inactive motor.

Numerical and analytical studies of such models have shown that phase separation is favored when the speed $|V_{z}|$ is a decreasing function of the active Brownian particle density. Fily and Marchetti (2012); Redner et al. (2013); Cates and Tailleur (2013); Stenhammar et al. (2013, 2014); Bialké et al. (2013); Wysocki et al. (2014); Fily et al. (2014); Speck et al. (2014, 2015) Important parameters that control whether the system will segregate into distinct phases are the Péclet number, $\text{Pe}=|V_{z}|\tau_{\text{R}}/R_{\text{m}}$, and the volume or area fraction, $\phi$. (We consider two-dimensional systems where $\hat{\bf z}_{i}=(\cos\theta_{i},\sin\theta_{i})$.) Phase diagrams have been constructed that show regions in the Pe–$\phi$ parameter plane where segregation into dense and gas-like phases occurs. Speck et al. (2014); Stenhammar et al. (2014); Speck et al. (2015) The effective diffusion coefficient for the active Brownian particles is predicted to depend on the area fraction as $D_{\text{m}}(\phi)=D_{0}+\frac{1}{2}V_{z}(\phi)^{2}\tau_{\text{R}}$. This form is consistent with active Brownian model simulation results, which also show that $|V_{z}(\phi)|$ decreases linearly with $\phi$. Stenhammar et al. (2014)

Two-dimensional simulations of a Brownian dynamics model for symmetric sphere-dimer motors that includes reorientation effects due to conservative forces have also been carried out. Gonnella et al. (2014) This model is very similar to that in Eq. (V), except applied to the monomers in the dimer, and its overdamped limit takes the form

[TABLE]

where $\nu$ labels the monomers in the dimer. The propulsion force ${\bf F}^{\text{p}}_{i\nu}$ again acts along the dimer bond. Phase separation dynamics is observed for certain values of Pe and $\phi$.

Some features of the collective motion of the backward-moving sphere-dimer motors considered in this paper are captured by such Brownian dynamics models. In particular, the dependence of the effective diffusion coefficient on area fraction is qualitatively similar to that of active Brownian particles although there are differences. The motor speed is a decreasing function of $\phi$ but there are deviations from the linear behavior seen in active Brownian models. The main difference in the effective diffusion coefficient arises from the dependence of $\tau_{\text{R}}$ on $\phi$, in contrast to the active Brownian models which assume that $\tau_{\text{R}}$ is independent of $\phi$. Taking this effect into account, $D_{\text{m}}(\phi)$ of the sphere-dimer can be written again as a simple density-dependent generalization of that for a single motor,

[TABLE]

This theoretical estimate is plotted in Fig. 9b, and one can see that the results agree well with the simulations. In fact, the main discrepancy occurs for $\phi=0$ where the statistical uncertainty is largest.

However, these simple Brownian particle models cannot fully describe the collective dynamics of chemically powered sphere-dimer motors. For example, for the symmetric sphere-dimer motors with beads of identical size, the Brownian particle models cannot distinguish between forward and backward motions, yet as Fig. 11 shows, forward-moving motors yield large cluster formation, but no such cluster formation occurs for backward-moving motors. These differences arise from interactions due to chemical gradients which are not included in many active Brownian particle models.

Phenomenological models, analogous to those of active Brownian particles, but incorporating the effects of chemical concentration gradients and again neglecting hydrodynamic interactions, have been used to describe the collective behavior of diffusiophoretic motors. In particular, for the $\text{A}\to\text{B}$ motor reactive dynamics considered here, the equations describing the position and orientation of motor $i$ are supplemented with terms involving concentration gradients,

[TABLE]

where $c_{\text{B}}({\bf r})$ is the local concentration of the product B and the parameters $a_{1,2}$ characterize the strength of the concentration gradient coupling. Anderson (1989); Anderson and Prieve (1991) The motor velocity $V_{z}(c_{\text{B}})$ is a function of the concentration. Depending on the values of the parameters $a_{1,2}\propto\frac{k_{\text{B}}T}{\eta}\Lambda$, different types of collective behavior, including gas-like, dynamic cluster states and other time-dependent and collapsed cluster states can exist. Pohl and Stark (2014); Saha et al. (2014); Pohl and Stark (2015)

Modifications of microscopic dynamics

Given the above summary of these active Brownian models that neglect hydrodynamic interactions, it is interesting to examine the consequences of this neglect in the context of our microscopic model. Here we consider two modifications of the microscopic dynamics where momentum is no longer locally conserved so that hydrodynamic interactions are absent, and the diffusiophoretic propulsion force that depends on concentration gradients is replaced by a constant propulsion force.

Coupling to the hydrodynamic fluid velocity fields can be eliminated by replacing the collision step in MPCD with a global Andersen thermostat. Instead of multiparticle collisions, the post-collision solvent velocities are assigned from a Boltzmann distribution at temperature $T$, which locally destroys conservation of momentum. The motors experience the same forces by self-produced chemical gradients as in the full dynamics; however, in the absence of local momentum conservation, these forces on the motor are not balanced by solvent flow and the diffusiophoretic mechanism no longer operates. Coupling to fluid hydrodynamic modes also changes other dynamical properties of the motor, such as its reorientation time, diffusion coefficient, and corresponding friction coefficient (see Table 3). Thus, neglect of coupling to fluid hydrodynamic modes in the microscopic theory alters both the dynamics of a single motor and its interactions with other motors. This is reflected in the significantly reduced values of the single-dimer propulsion velocities given in Table 3, which are approximately 30–40 times slower than their values when momentum is locally conserved. In phenomenological models, single motor properties are effectively rescaled while neglecting hydrodynamic coupling among motors.

Within the context of this microscopic model without hydrodynamic interactions, one may explore the effects of chemical coupling on the collective dynamics. Effects due to chemical gradients can be eliminated by removing the reactive collisions that produce product in the $\text{A}\to\text{B}$ reaction and replacing them with an external, constant force $F_{\text{m}}$ applied along the bond of each sphere dimer. The magnitude of the force is chosen such that a single dimer at area fraction $\phi\to 0$ has the same propulsion velocity as a single dimer propelled by a self-produced chemical gradient. The addition of external forces implies viscous heating of the solvent, which may be avoided by using a global Anderson thermostat that also eliminates hydrodynamic interactions.

Fig. 12 compares simulation results using the two modified dynamical schemes described above for asymmetric sphere-dimer motors with $d_{\text{C}}=4$ and $d_{\text{N}}=8$.

When chemical gradients are taken into account but hydrodynamic interactions are neglected, for forward-moving motors, one observes the formation of large clusters that persist for long times (Fig. 12a). Similarly, for backward-moving motors, a statistically stationary regime involving small transient clusters is observed (Fig. 12c). These behaviors are qualitatively similar to those observed in the phenomenological Langevin models that account for chemical gradients but neglect hydrodynamic interactions; however, the time scales of the dynamics are different since neglect of hydrodynamic interactions at the microscopic level destroys the diffusiophoretic mechanism as discussed earlier.

When replacing chemical coupling by an applied constant force for propulsion, in addition to neglect of hydrodynamic interaction, for forward-moving motors instead of the formation of large clusters one observes transient clusters with a significant amount of translational and rotational motion (Fig. 12b). The collective behavior of backward-moving motors is qualitatively similar to that of the full microscopic model (Fig. 12d).

Recall that for symmetric sphere-dimer motors, there is no difference between forward- and backward-moving motors when chemical gradients are neglected. Thus, the difference in the results for forward- (Fig. 12b) and backward-moving (Fig. 12d) asymmetric motors obtained using an applied force to move the motor instead of a chemical gradient is a geometric effect arising from the motor asymmetry. For full microscopic dynamics or models that include chemical gradients but no hydrodynamics, chemical gradients play a major role in determining the dynamical structure, but geometrical effects arising from direct intermolecular interactions may also be important. The relative importance of these two effects merits further study.

VI Conclusion

Through the study of microscopic models of active systems that include both the active particles and solvent species, insight into the various factors that determine the character of the collective behavior can be obtained. For the diffusiophoretic sphere-dimer motors considered in this paper, these factors include the explicit microscopic description of the reaction kinetics on the motor surface, direct motor–motor anisotropic interactions, hydrodynamic coupling arising from solvent fluid flows, and thermal noise. Of course, there is a price to be paid for this level of description: there are limitations on the size of the system that may be studied conveniently. Although the simulation results presented in this paper are for very large systems comprising thousands of motors and up to $10^{8}$ solvent molecules, the number of motors is several orders of magnitude smaller than the largest active Brownian particle simulations. As a consequence, our study was restricted to a limited number of points in parameter space, while Brownian dynamics simulations have been used to construct phase diagrams in large parameter regions and characterize the scaling behavior of the phase segregation process. Since our dynamics follows directly from the intermolecular interactions, it can be used to test the validity of the assumptions that underlie the phenomenological models and aid in the construction of more accurate reduced descriptions. The research described in this paper can be extended to other motors, such as simpler Janus motors, to investigate further aspects of the collective dynamics. Huang et al. (2016)

The phase behavior of active systems using Brownian dynamics models is usually studied in parameter spaces that include quantities such as the motor velocity, Péclet numbers, volume fraction, and diffusiophoretic coupling coefficients. In microscopic models, these parameters are functions of the intermolecular interactions that define the system. For example, changes in the solvent dynamics through the MPCD collision rule in our model will change the solvent viscosity and diffusion transport properties and will also simultaneously change the motor velocity and reorientation time, as well as the diffusion-influenced reaction kinetics on the motor. Changes in the intermolecular interactions between the reactive species and the motor (described by the $\Lambda$ factor in the continuum theories) will simultaneously change the self-propulsion properties of the motor and the coupling through the concentration gradients due to other motors. If hydrodynamic interactions are suppressed in the microscopic model by altering the collision dynamics, the basic diffusiophoretic mechanism that causes a single motor to be self-propelled is also turned off. Rescaling of transport properties and other alterations to the dynamics can restore some of the effects but at the expense of using a dynamical description that is less firmly grounded.

The other factor that plays a role in determining the character of the collective dynamics is the manner in which the system is driven from equilibrium. Fuel must be supplied to the motors and product must be removed to maintain the system out of equilibrium, and these species may be introduced at the boundaries or globally through bulk nonequilibrium reactions as in this study. If species are supplied or removed at the boundaries, geometry and dimensionality will play important roles; for example, correlations arising from many-body concentration fields can lead to non-analytic dependence of reaction rates on the motor volume fraction. Tucci and Kapral (2004); Lebenhaft and Kapral (1979)

Active systems present challenges because they operate out of equilibrium and display diverse phenomena. Theoretical models at different levels of description can be used to unravel and understand the origins of these phenomena.

This work was supported in part by grants from the Natural Sciences and Engineering Research Council of Canada and Compute Canada.