Forces that control self-organization of chemically-propelled Janus tori

Control of the individual and collective behavior of self-propelled synthetic micro-objects has immediate application for nanotechnology, robotics, and precision medicine. Despite signiﬁcant progress in the synthesis and characterization of self-propelled Janus (two-faced) particles, predictive understanding of their behavior remains challenging, especially if the particles have anisotropic forms. Here, by using molecular simulation, we describe the interactions of chemically-propelled microtori near a wall. The results show that a torus hovers at a certain distance from the wall due to a combination of gravity and hydrodynamic ﬂows generated by the chemical activity. Moreover, electrostatic dipolar interactions between the torus and the wall result in a spontaneous tilt and horizontal translation, in a qualitative agreement with experiment. Simulations of the dynamics of two tori near a wall provide evidence for the formation of stable self-propelled bound states. Our results illustrate that self-organization at the microscale occurs due to a combination of multiple factors, including hydrodynamic, chemical, and electrostatic interactions.


INTRODUCTION
Individual and collective behavior of self-propelled micro-objects (also called agents or motors) is the main topic of a rapidly expanding field termed active matter [1,2]. Recent research on synthetic active particles has demonstrated that they possess diverse functionalities which are similar to those found in the living world, such as propulsion and energy conversion [3,4], artificial chemotaxis [5][6][7] and rheotaxis [8,9], etc. However, most of the synthetic realizations of self-propelled particles lack the fidelity and efficiency of biological organisms. The limitations come mostly from micro-fabrication, resulting in relatively simple shapes such as spheres [10], rods [3], or helices [11]. These factors significantly restrict the range of possible applications of self-propelled micromotors in robotics, drug delivery, and precision medicine [12][13][14].
Advances in additive manufacturing and nanoscale 3D printing allow the design of particles with practically arbitrary shape, e.g., tori or propellers [15]. Functionalization of the tori surfaces by platinum (Pt) and nickel (Ni) enables propulsion in hydrogen peroxide (H 2 O 2 ) solution and a response to an applied magnetic field. Without a magnetic field, the tori (with sizes between 7 and 10 µm) sediment and hover near the bottom wall at a distance of the order 1-2 µm. Due to a spontaneous symmetry breaking of the hovering state, the tori tilt and translate parallel to the bottom with the speed of the order of few µm/s. The tori also form multiple rotating and translating bound states. If a dc magnetic field is applied parallel to the bottom wall, the tori turn perpendicular to the bottom and swim with much higher speed (10-20 µm/s) in the direction controlled by the magnetic field.
In this work we use molecular simulation to uncover fundamental mechanisms governing the organization of chemically self-propelled anisotropic particles. In particular, we investigate the dynamics of Janus tori near a wall using multiparticle collision dynamics [16], a particlebased mesoscale simulation technique for complex fluids that accurately incorporates thermal fluctuations, hydrodynamic interactions, chemical reactions, and solid inclusions [17,18]. Our results show how various interaction potentials govern the behavior of several tori interacting with each other and the confining walls. We have found that all of gravitational, hydrodynamic, and chemical interactions are responsible for the hovering behavior of a torus above the bottom wall. When the electrostatic dipolar interactions between the torus and the bottom wall are included, the torus tilts spontaneously and translates parallel to the bottom, in qualitative agreement with experiment. We also show how the interplay among the different interactions results in the formation of bound states of two tori.

ACTIVE JANUS TORUS NEAR A WALL
The simulations of the dynamics of a torus near a wall are carried out in a system with a slab geometry: the top and bottom walls of the simulation volume are parallel to the (x, y) plane and separated by a distance L z . Periodic boundary conditions are applied in the x and y directions with L x = L y . The volume contains one or more tori with centerline radius a constructed from linked overlapping spherical beads that form a ring [19]. The torus beads are Janus spheres with catalytic C and noncatalytic N hemispherical faces. The solvent comprises A and B fluid particles. Catalytic reactions A → B occur on the catalytic faces of the beads. Figure 1a shows a torus oriented so that the catalytic side of the torus faces the wall. The symmetry axis of the torus is indicated by the unit vectorû pointing from the catalytic to noncatalytic sides of the torus. The simulation algorithm and its implementation is described in Methods.
A torus is heavier than the solvent and, due to the gravitational force F g , an inactive torus will sediment to the bottom and experience a short-range repulsive force from the wall. For an active torus where chemical reactions occur on its catalytic surface, because of the asymmetric distribution reactants A and products B a diffusiophoretic mechanism [20][21][22] operates to propel the torus in solution. The magnitude and direction of the diffusiophoretic force F d on the torus depends on the interactions of the reactive species with the torus beads and the displacement from equilibrium. [23] We adopt the following convention: if the diffusiophoresis force F d is in the same direction as theû-vector, i.e., towards the noncatalytic (passive) surface, it is called a backwardmoving torus. If it is a direction opposite toû, it is a forward-moving torus.
If a backward-moving torus is initially oriented so that its catalytic surface faces the wall as in Fig. 1a, then the gravitational and diffusiophoretic forces will have opposite signs. Consequently, the torus will be confined to the vicinity of the wall, provided that F g > F d , if there are no sufficiently large orientational fluctuations allowing the torus to change its orientation relative to the wall so that escape is possible. Correspondingly, the orientation of a forward-moving torus relative to the wall must be opposite of that of a backward-moving torus for confinement to occur. Under these conditions, we can investigate the torus dynamics in the vicinity of the wall versus the strength of the confining gravitational force.
The probability distribution P (z) of the distance of the center of mass (c.o.m.) of the backward-moving torus from the wall is shown in Fig. 1b for several values of the gravitational strength parameter g (the value of g can be adjusted by the mass density of the torus material). For the g values in the range (0.1 − 0.3) presented in the figure, the torus remains confined to the vicinity of the wall and, as expected, the mean distance from the wall increases as g decreases. For g 0.15, the torus beads do not experience direct repulsive interactions with the wall. This behavior indicates that the observed hovering is due to diffusiophoresis generated by chemical reactions rather than the short-range repulsive forces.
The tilt angle θ of the torus with respect to the bottom plane is determined from the scalar product of their normal vectors: cos(θ) =û ·ẑ. The probability distribution of θ = arccos(û ·ẑ), p(θ), normalized so that π 0 dθ sin θ p(θ) = 1, is shown in Fig. 1c. All distribu-tions are peaked at the origin and broaden as g decreases, indicating larger angular fluctuations as the confinement to the wall vicinity weakens. For g = 0.3, the torus interacts strongly with the wall through short-range repulsive forces and p(θ) (not shown) has a very narrow peak at zero.
The forms of the concentration and fluid velocity fields that accompany the diffusiophoretic mechanism provide additional insight into the origin of torus hovering. The c B concentration and v(r) velocity fields in the (x, z) plane near the hovering backward-moving torus with g = 0.15 are shown in Fig. 2a and b. The fields are displayed in a c.o.m. reference frame. In this frame, for the concentration, we measure the average number of product B particles in a thin slab in the (x, z) plane with thickness ∆ = 1 whose y = L/2 coordinate lies at the torus c.o.m. The z direction in this frame is the same as theû vector of the torus. Similarly, the velocity field was calculated by averaging the velocities of all particles in volume elements with thickness ∆ = 1 in the (x, z) plane. Consequently, the fields are averages over the vertical and angular displacements that the torus experiences in its motion along the wall.
As expected, Fig. 2a shows the increased product concentration in the torus hole and near the catalytic surface. In the velocity field plot (Fig. 2b) we observe that the solvent particles near the C − N interface of the torus flow from the N to C surface. We also see that there is a suction above the torus. The solvent particles flow out along the torus bottom after they flow into the hole because the wall limits the fluid flow. There are also re-circulation regions with vortices slightly above the torus. For a hovering forward-moving torus, comparison of Figs. 2c and 2a shows that the magnitude of the Bparticle concentration field is significantly different because the catalytic surface no longer faces the wall. Also, comparison of the corresponding fluid velocity field near a hovering forward-moving torus surface in Fig. 2d with Fig. 2b shows that the solvent particles near the C − N interface flow from the C to N surfaces. Also, there is a suction above the torus hole that pushes the fluid particles along the torus bottom surface because the wall limits the flow.
The corresponding concentration and velocity fields for backward and forward-moving tori in the bulk fluid are shown in the middle panels of Fig. 2. Comparison of the concentration fields reveals the distortions in these fields due to the wall. The distortion of the concentration gradient affects the diffusiophoretic force and, correspondingly, the propulsion. The flow fields have a completely different form for hovering tori and those in the bulk fluid: for a backward-moving torus in the bulk phase the fluid flows outward from the noncatalytic surface of the torus and, after circulation, inward to the catalytic surface. The flow field for a cross-section of the torus is similar to that seen for a Janus colloid, as might be expected. The direction of the fluid velocity of the forward-moving torus is opposite to that for a backward-moving torus.

EFFECT OF ELECTROSTATIC FORCES ON THE TORUS INTERACTION WITH THE WALL
Experimental results reported in Ref. [15] indicated that a torus acquires an electric dipolar moment resulting from the electrochemical reactions occurring on the Pt surface of a torus. This conclusion is supported by the following observation: a torus collected negatively charged 1 µm polystyrene spheres along its Pt surface and accumulated 2 µm positively charged latex spheres along its polymer surface. To incorporate this effect, we included electric dipoles with the strength ǫ cd in the model: each bead i in the torus interacts with a uniformly charged wall with through the interaction potential V cd /∂z i )ẑ ≡ F cd γẑ, where γ = cos(2πi/(N b − 1)) controls the sign and relative magnitude of the dipole of the N b beads in the torus.
As one sees from Fig. 3a, the electrostatic interaction results in a symmetry breaking of the horizontally hovering state and the shift of θ max , the tilt angle at the maximum in p(θ), from θ max = 0 to θ max > 0. Figure 3b plots θ max versus the charge-dipole strength ǫ cd where it is seen that beyond a critical value ǫ * cd ≈ 0.11 it takes nonzero values and increases with increasing ǫ cd . This effect is also seen in the bottom panels of Fig. 2i-f that display the concentration and velocity field plots for backward and forward moving tori with charge-dipole interactions. These fields show the asymmetries that are expected in view of tilted geometries that break the symmetry of the horizontal state.
Furthermore, a tilted torus translates parallel to the bottom wall with a velocity proportional to the tilt anglē θ. To examine this effect we consider the velocityV u defined as follows: the projections of the motor velocity V and orientation vector θ max onto the surface plane are V = V · 1 and u =û · 1 , respectively, where 1 = 1 −ẑẑ and 1 is the unit tensor. The mean velocity can be defined asV u = V ·û . Figure 3c shows that V u remains constant for tilt angles below the bifurcation where θ max = 0 and increases as θ max increases beyond the bifurcation point.

DYNAMICS OF TWO TORI NEAR THE WALL
Next, we consider the dynamics of two interacting tori in the vicinity of the wall. The tori are subject to all the interactions in Eq. (11) as well as a short-range steric repulsion. To accommodate both tori in configurations where they can separate sufficiently far from each other and behave independently, the size of the simulation box is increased to 60×60×40. Periodic boundary conditions are again used in directions parallel to the walls.
For a pair of hovering backward-moving tori, Fig. 4  (top) shows the probability distributions P (z) vs their c.o.m. distances from the wall, z, trajectories of the c.o.m. of two different tori, and the probability distributions of the tilt angles θ of two tori, p(θ). We observe that for both tori the P (z) distributions practically coincide. The θ angle distributions also very closely correspond and the average tilt angles are about the same as that for a single torus under the same conditions. The trajectories show that the distances between the tori c.o.m. are concentrated between D = 18 and 28 so that the tori do not form a long-lived bound state. The instantaneous configuration in Fig. 4 (top) shows that the tori have no tendency to form a bound state, see Supplementary Video 1.
Two forward-moving tori with g = 0.07, ǫ cd = 0.2 behave very differently as can be seen from Fig. 4 (middle). Now both the P (z) and p(θ) distributions are distinct, and on average, one torus lies at a distance farther from the wall than the other, also the tilt angle is larger for the upper torus. Examination of the trajectories shows that the distances between the c.o.m. of two tori are concentrated between D = 7.5 and 9. This indicates that the two tori are stacked in a bound state since a torus has a radius a = 5. These features are evident in the instantaneous configuration in Fig. 4 (middle), see also Supplementary Video 2.
If the constant that gauges the magnitude of the gravitational force is increased to g = 0.09, the behavior of the pair of forward-moving tori changes, see Supplementary Video 3. In Fig. 4 (bottom) we see that the P (z) and p(θ) distributions for the two tori coincide, so both tori lie at similar distances from the wall, and their tilt angles are close. In contrast to the corresponding tilt angle distribution for a single forward-moving torus, p(θ) is peaked at the origin, and the tilt angles differ from that for a single torus near the wall. The trajectories indicate the formation of a long-lived bound state that translates and rotates. An instantaneous configuration is shown in Fig. 4 (bottom). This behavior results from stronger gravitational confinement causing the tori to adapt nearly parallel configurations to the wall. Consequently, they cannot explore regions further from the wall to adopt other configurations, such as the stacked bound state for g = 0.07.

DISCUSSION
Here we develop a simple phenomenological theory for p(θ), the probability distribution of the torus tilt angle. First, we consider the system in the absence of  charge-dipole interactions where the torus hovers above the wall at some distance due to hydrodynamic and diffusiophoretic forces and is subject to a gravitational force. The orientation vector of the torus isû = (sin θ cos φ, sin θ sin φ, cos θ). In the Cartesian frame witĥ z normal to the wallẑ ·û =û z = cos θ that defines the tilt angle θ. The dynamics of the orientation vector is governed by the Langevin equation describing the dynamics of the torus, The external torque due to the diffusiophoretic forces and gravitational field is denoted by T g = γu ×ẑ, where γ depends on the strength of the gravitational field, T fl is the fluctuating torque and ζ is the rotational friction tensor. The Langevin equation for the dynamics of the tilt angle can be determined from that for u z and is given by where the fluctuating force isT fl = T fl,x sin φ − T fl,y cos φ that satisfies the fluctuation-dissipation relation T fl (t)T fl (t ′ ) = 2k B T ζ r δ(t − t ′ ). Here ζ r is the component of the rotational friction tensor for u z and the corresponding rotational diffusion coefficient is D r = k B T /ζ r . The Fokker-Planck equation reads whose stationary solution is where β = 1/k B T and the normalization is determined from π 0 dθ sin θ p s (θ) = 1. This result may also be derived from Fokker-Planck equation that accounts for the effective orientational potential energy corresponding to u z , U r (θ) = −γ cos θ. The expression for p(θ) in Eq. (4) is in accord with the simulation data in Fig. 1c and shows that the angle corresponding to the maximum in this distribution is θ max = 0.
When the electric dipolar interaction is included the torque on the torus changes and includes another contribution T cd , that contributes to the previous torque and contains and additional contribution that is proportional to cos θ. The Langevin equation is now given by and the corresponding Fokker-Planck equation is whose stationary solution is In this case the corresponding rotational potential energy is U r (θ) = −γ 1 cos θ(1 − (γ 2 /2) cos θ) and has a contribution proportional to cos 2 θ that is responsible for breaking symmetry. In contrast to the situation where there are no charge-dipole interactions, now the probability distributions indicate that the angle corresponding to the maximum in the distribution is at θ max > 0, provided the charge-dipole strength γ 2 is sufficiently large. Figure 3a compares the angle probability distributions p(θ) for different values of the electric dipolar strength ǫ cd obtained from the simulations with the predictions of Eq. (7). The simulations results are in good agreement with the analytical expression in Eq. (7) showing that this simple model captures the origin of the symmetrybreaking mechanism.

CONCLUSIONS
We have shown that a variety of forces control the behavior of chemically-propelled torus colloids near walls. They include gravitational, steric-repulsive, diffusiophoretic, hydrodynamic, and electrostatic forces. If we neglect any of the forces, then the experimentally observed behavior cannot be reproduced. For example, no persistent horizontal translation near the wall can be obtained without electrostatic interaction. Our simulations also make a prediction that backward moving tori do not form bound states (see Fig. 4(top)). In experiment, the backward vs forward mode of propulsion is controlled by the relative thicknesses of the Pt and Ni layers [15]. For example, for thin Pt layers (about 10 nm), the tori swim Pt (catalytic) side forward, whereas for the Pt layers thicker than 40 nm, the tori swim polymer (non-catalytic) side forward. Another important prediction from our work is that, depending on the gravitational strength, a transition from stacked bound states (Fig. 4(middle)) to flattened" bound states (Fig. 4(bottom))occurs. In turn, the effective mass density can be controlled by depositing more catalyst (PT), making the particles heavier. Alternatively, 3D printing hollow particles will make them lighter. The algorithms developed here can be applied to particles of different shapes (e.g., propellers, helices, etc.) or more complex chemical reactions. In addition to electrostatic forces, magnetic interactions can be included as well.

METHODS
The simulation box has dimensions L x = L y = L z = 40 (or L x = L y = 60 and L z = 40 for systems with two tori) and contains N S = N A + N B solvent particles with mass m. The construction of the torus motor from a collection of beads was described in our earlier study of an active torus in bulk solution. [19] Here we consider tori with a centerline radius a = 5 comprising N b = 39 beads with radius r b , where neighboring beads are connected by harmonic springs with k s = 500. The equilibrium bond length between ith and (i + j)th beads is r i,i+j 0 = 2a sin(jπ/N b ), j = 1, 2, . . . , N b /2. The torus has mass M 39 = 3890.
The A and B fluid particles undergo bounce-back collisions with the walls, and interact with each bead in a torus motor through short-range repulsive Lennard-Jones potentials. The particles interact among themselves through multiparticle collisions [16]. The torus motor interacts with the bottom wall through a short-range repulsive Lennard-Jones wall potential energy function. We also consider situations where the wall has a positive charge density and a torus has nonuniform distribution of dipoles. In this case we also have charge density-dipole interactions for each bead with the wall. To account for the effect of gravity on the motor and the possible inhomogeneous distribution of mass due to the different compositions of the catalytic and noncatalytic portions of the beads, we have a gravitational potential energy function.
More specifically, these potential energy functions have the following forms: The α = {A, B} particles interact with the torus beads through the potential V (b) α (r ki ) with strength ε α , where r ki = r k − r i is the vector from the center of torus bead i to fluid particle k of species α. The short range potential function has the form for r ki ≤ r c where r c = r b + 2 1 6 σ and zero for r ki > r c .

7
We take σ = 1.0, r b = 1.0. For backward-moving tori we take ǫ A = 0.1ǫ, ǫ B = ǫ, while for forward-moving tori these potential strengths are interchanged and ǫ A = ǫ, ǫ B = 0.1ǫ. Each torus bead i interacts with the bottom wall through a repulsive Lennard-Jones wall potential energy function, and zero otherwise. Here z i is the vertical distance of motor bead i from the lower wall at z = 0. We take σ w = 3.0 and ǫ w = 1.0. The potential is truncated at the minimum of the potential energy function at z min = (2/5) 1/6 σ w = 2.575 so that the wall force vanishes for z ≥ z min . A similar expression applies to the upper wall. The torus beads also interact with the bottom wall through charge density-dipole potential energy functions, cd (z i ) = ǫ cd 2 ln (L 2 z +z 2 ci )/z 2 ci cos(2πi/(N b −1)), (10) where each torus bead is labeled by i = 0, . . . , N b − 1 and z ci is the vertical distance of motor bead i from the center of mass of the catalytic hemisphere to the bottom wall. The gravitational potential acting on a bead is given by V (b) g (z ci ) = −gz ci where g accounts for the mass and the acceleration due to gravity.
The total interaction energy is then given by where θ α k is an indicator function that is one if particle k is species α and zero otherwise.
The time evolution of the entire system is governed by molecular-dynamics-multiparticle collision dynamics. [16][17][18]24] An irreversible chemical reaction, C + A → C + B, takes place on the C hemispherical surface of each motor bead that converts species A to product B. The reactive events occur with unit probability when an A or B particle crosses a reaction radius R r = r c . The fluid phase reaction that is used to maintain the system in a nonequilibrium state has rate constant k b = 0.001 and is carried out using the reactive version of multiparticle collision dynamics. [25] The molecular dynamics portion of the evolution uses the velocity-Verlet algorithm with a molecular dynamics time step, δt = 0.001, and the multiparticle collision dynamics portion uses a multiparticle collision time of τ = 0.5. The collision rotation angle is α = π/2. Grid shifting is used to insure Galilean invariance. [26] The system temperature is k B T = 0.2. The solvent density is c 0 ≈ 10, so that the approximate number of fluid particles is given by N s = c 0 L x L y L z . For τ = 0.5 the fluid dynamic viscosity is η = 1.43 and the common diffusion coefficient for the A and B species is D = 0.117.
Results are reported in dimensionless units with mass m, length σ, energy ǫ, and time t 0 = mσ 2 /ǫ.