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 significant 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 flows 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 and four 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, electrostatic and steric interactions. The presence of a constraining environment exerts an influence on the behavior of self-propelled synthetic microswimmers, challenging the prediction and control of their individual and collective behaviour in realistic situations. Here, the authors use multiparticle collision dynamics to simulate self-propelled Janus toroidal particles near a wall and study how various contributions, such as thermal fluctuations, hydrodynamic and electrostatic interactions, chemical reactions, and gravity govern their collective behaviour.

I ndividual and collective behavior of self-propelled microobjects (also called agents or motors) is the main topic of a rapidly expanding field termed active matter [1][2][3][4][5][6][7][8][9][10] . 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 11,12 , artificial chemotaxis [13][14][15] , rheotaxis 16,17 , 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 18 and rods 11 , although more complex selfpropelled helices have been made using glancing angle deposition 19 . These factors significantly restrict the range of possible applications of self-propelled micromotors in robotics, drug delivery, and precision medicine [20][21][22] .
Advances in additive manufacturing and nanoscale 3D printing allow the design of particles with practically arbitrary shape, e.g., tori or propellers 23 . 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 Janus tori near walls. There have been a number of studies of the behavior of active particles near walls [24][25][26][27][28][29] in which various kinds of dynamical behavior have been observed. For instance, diffusiophoretically active Janus colloids near a hard wall have been shown to hover above the wall or slide along it at a constant height and a preferred tilt angle 27,29 . Also, there have been investigations of the dynamics of more complex active agents, such as rigid swimmers constructed from linked rotating beads 30 and swimmers with conformational dynamics 31 , near a wall. In addition, there is body of literature that considers the dynamics of a torus in bulk solution. There are studies using analytical and numerical continuum methods of the properties of inactive tori 32,33 , as well as investigations of active tori including a torus propelled by rotation about its centerline [34][35][36] , a torus made from linked rotating beads 36 and a torus propelled by a diffusiophoretic mechanism 37 .
Our investigations of the dynamics of Janus tori near a wall use multiparticle collision dynamics 38 , a particle-based mesoscale simulation technique for complex fluids that accurately incorporates thermal fluctuations, hydrodynamic interactions, chemical reactions, and solid inclusions 39,40 . We 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 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 these different interactions along with steric repulsion results in the formation of bound states of two and four tori.

Results and discussion
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 41 . 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 are 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 of reactants A and products B, a diffusiophoretic mechanism 8,[42][43][44] 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. 44 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 backward-moving torus. If it is a direction opposite toû, it is a forward-moving torus. In experiment, the backward vs forward mode of propulsion is controlled by the relative thicknesses of the Pt and Ni layers 23 . 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. 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 in a long-lived metastable state provided that F g > F d . If F g < F d such an initially oriented torus can easily escape from the wall to the bulk fluid phase. If the initial orientation is opposite so that its noncatalytic surface faces the wall, then both forces act in the same direction to confine the torus to the wall in what we call a trivial metastable state. These are metastable states since sufficiently large orientational fluctuations can allow the torus to change its orientation relative to the wall, thus changing the probability of residing near the wall. A corresponding set of statements apply for forward-moving tori; for example, for F g > F d the orientation of a forward-moving torus relative to the wall must be opposite of that of a backward-moving torus for longlived confinement to occur. Consequently, for F g > F d it is interesting to investigate the properties of the long-lived metastable states as a function of the strength of the confining gravitational force, and we present these results below.
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 R π 0 dθ sin θ pðθÞ ¼ 1, is shown in Fig. 1c. All distributions 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. From Fig. 1c we also see that large orientational fluctuations have very low probability because they are suppressed by interactions with the wall. Except for the largest gravitational strengths these wall interactions are diffusiophoretic in nature and do not depend strongly on direct wall repulsive forces. In the bulk fluid phase the orientational relaxation time is τ r~1 0 4 . In our simulations for times longer than τ r we have never observed escape from these long-lived metastable states, although such escape is possible.
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 backwardmoving torus with g = 0.15 are shown in Figs. 2a and b. The fields are displayed in a laboratory 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 Fig. 2c, a shows that the magnitude of the B-particle 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. In addition, 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 Fig. 2e-h. Comparison of the concentration fields reveals the distortions in these fields due to the phoretically-induced interaction with the wall. The distortion of the concentration gradient affects the diffusiophoretic force and, correspondingly, the propulsion. The flow fields have a different form for hovering tori and those in the bulk fluid where 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. While the near field velocity flow pattern for a torus bead in a cross-section of the torus is similar to that seen for a Janus colloid, the flow patterns in the vicinity of the entire torus differ from those for a spherical Janus colloid. The direction of the fluid velocity of the forward-moving torus is opposite to that for a backward-moving torus.
From these results we see that while the torus exhibits a hovering state above the wall in the presence of a gravitational field, the tilt angle distribution in Fig. 1c is centered at zero and the distribution becomes less sharply peaked as g decreases, consistent with weaker confinement that leads to enhanced thermal orientational fluctuations. Orientational fluctuations are suppressed by interactions with the wall since, if the torus tilts at fixed z, a portion of the torus will move closer to the wall and be repelled by it. Gravitational confinement prevents large changes in z as seen in Fig. 1b. The result is a dynamic hovering state where the torus is on average parallel to the wall but executes active diffusive motion due to thermally driven orientational fluctuations that produce fluctuating diffusiophoretic force components parallel to the wall.
The hovering and sliding dynamical states seen for Janus colloids near walls in the presence of a gravitational field have been shown to depend on whether the surface mobility is uniform or non-uniform 27,29 . We have assumed uniform mobility since the solvent species have the same interactions with the catalytic and noncatalytic torus surfaces and find hovering and sliding states. Also, spherical Janus colloids will not experience the same shape-dependent effect on the orientational fluctuations as torus colloids. Nevertheless, the essential factors are wall modulations of the concentration and velocity field effects arising from self diffusiophoresis in both cases.
Effect of electrostatic forces on the torus interaction with the wall. We now show that the experimentally observed states where the torus moves along the wall with a preferred tilt angle can be modeled in the simulations by including a charge-dipolar interaction between the torus and the wall. These results are summarized in Figs. 2a-l and 3.
Experimental results reported in Baker et al. 23 suggest that a torus acquires an electric dipolar moment resulting from the electrochemical reactions occurring on the Pt surface of a torus. where a wall is present, indicated by dashed horizontal lines, the ordinates refer to the z distance from the wall, while the abscissas gives the distance x along a plane parallel to the torus orientation vector that cuts the torus through its center. The black arrows indicated the direction of motion of moving tori in all panels. The blue colour corresponds to the non-catalytic surface, and the red colour corresponds to the catalytic surface, respectively. While direct experimental measurements of the electric dipole are not possible at this time, the supporting evidence for the dipole is two-fold: (i) 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, implying that the torus's platinum and polymer sides acquire opposite electric charges. (ii) Adding salt to the solution kills the effect of self-propulsion and, consequently, the ability to collect beads by the torus. This behavior is consistent with the electrohydrodynamic propulsion of the torus and the onset of the dipolar moment.
In addition, the bottom glass plate also acquires an electric charge in H 2 O 2 solution. Since the torus platinum cap is conducting, due to interaction with the absorbed charges in the glass plate, the electric dipole deviates from the axial direction and orients in-plane of the torus. This interaction, in turn, results in a constant tilt angle.
Our modeling of the torus dynamics near the wall does not take into account all of the effects described above that are postulated to give rise to wall charge-torus dipolar interactions. Instead, we assume the existence of such interactions and investigate if these interactions in the presence of a gravitational force can lead to the onset of torus tilt for some specific values of the strengths of gravitational and dipolar-wall interactions. Thus, we included electric point dipoles with the strength ϵ cd in the model: each bead i in the torus interacts with a uniformly charged wall through the interaction potential V ðbÞ cd ðz i Þ whose explicit form is given in the Methods section. The potential gives rise the force F cd ¼ Àð∂V ðbÞ cd =∂z i Þẑ F cd γðiÞẑ, where γðiÞ ¼ 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, for a given gravitational force, a sufficiently large 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 θ max rises sharply from values near zero in a narrow transition region to a domain where it has large θ max values that increase with further increases in ϵ cd . The plot is suggestive of a bifurcation to a tilted torus beyond a critical value ϵ Ã cd % 0:12 that is smoothed by thermal fluctuations. This observation is consistent with the plots in the bottom panels of Fig. 2i-l 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 that depends on the tilt angle θ. To examine this effect we consider the velocity component V u k defined as follows: the projections of the motor velocity V and orientation vectorû onto the surface plane are V ∥ = V ⋅ 1 ∥ and u k ¼û Á 1 k , respectively, where 1 k ¼ 1 Àẑẑ and 1 is the unit tensor. The mean velocity can be defined as V u k ¼ V k Áû k . Figure 3c shows that V u k increases as θ max increases from zero in the presence of thermal fluctuations since any tilt from horizontal will induce motion parallel to the wall. Such orientational fluctuations are responsible for the non-zero values of V u k , even at the values of the average θ max near zero seen in the figure. In the transition region V u k increases weakly since this region encompasses a small range of ϵ cd , followed by a sharper increase that is proportional to θ max .
Next, 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 iŝ u ¼ ðsin θ cos ϕ; sin θ sin ϕ; cos θÞ. In the Cartesian frame withẑ 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 hT fl ðtÞT fl ðt 0 Þi ¼ 2k B Tζ r δðt À t 0 Þ. 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 R π 0 dθ sin θ pðθÞ ¼ 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 dθ dt ¼ À γ 1 ζ r sin θð1 À γ 2 cos θÞ þ 1 ζ rT fl ðtÞÞ; ð5Þ and the corresponding Fokker-Planck equation is whose stationary solution is pðθÞ ¼ N e βγ 1 cos θð1Àðγ 2 =2Þ cos θÞ : 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  Dynamics of several tori near the wall. We now consider the dynamics of two and four interacting tori in the vicinity of the wall. The tori are subject to all the interactions discussed above for a single torus, as well as a short-range steric repulsion between beads on different tori. To accommodate several 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.
Two tori. The simulation results for two tori near the wall are summarized in Fig. 4a-l. Using the same parameter values, g = 0.15 and ϵ cd = 0.2, as in Fig. 2 for a single hovering backwardmoving torus near the wall, the results for a pair of such tori are given in Fig. 4a-d. This figure shows the probability distributions P(z) vs their c.o.m. distances from the wall, z, sample trajectories of the c.o.m. positions of the two tori, and the probability distributions of the tilt angles θ for each of the two tori, p(θ). We observe that for both tori the P(z) and distributions practically coincide. The p(θ) distributions also very closely correspond, and the average tilt angles are about the same as that for a single torus under the same conditions presented earlier. The trajectories show that the distances between the tori centers of mass, D, are concentrated between D = 18 and 28. These results, along with the instantaneous configuration in Fig. 4d and Supplementary Video 1, support the conclusion that the tori behave independently without forming a long-lived metastable bound state. This behavior is analogous to that for backward-moving spheredimer 45 and Janus 46 colloids which also show little or no tendency to cluster.
By contrast, two forward-moving tori near the wall can form long-lived bound states, as we now describe. We first consider systems with ϵ cd = 0.2 and two values of the gravitational force: g = 0.07 (again the same parameter values as in Fig. 2 for a single torus), and a larger value g = 0.09. The results for g = 0.07 in Fig. 4e-h indicate that the P(z) and p(θ) distributions of the individual tori are no longer the same, as was the case for backward-moving tori. 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. The plot of the probability density, P(D), in Fig. 5a, shows that the D values are concentrated between D = 7.0 and 10.0. Since a torus has a centerline radius a = 5, this indicates that the two tori form a stacked bound state. These features are evident in the instantaneous configuration in Fig. 4e-h, see also Supplementary Video 2.
For g = 0.09, the character of the metastable bound pair changes and adopts a predominantly side-by-side structure parallel to the wall. (See Supplementary Video 3) This observation is quantified in Fig. 4i-l where 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 forwardmoving torus, p(θ) is peaked at the origin, and the tilt angles differ from that for a single torus near the wall. The P(D) distribution in Fig. 5b is sharply peaked at D ≈ 14, characteristic of a well-defined bound pair. The trajectories indicate that the bound state translates and rotates. An instantaneous configuration is shown in Fig. 4i-l. Simulations have also been carried out for a lower value of ϵ cd = 0.1 and the same two values of g as above. The results show that for g = 0.09 the bound state again has a predominantly sideby-side structure, while for g = 0.07 the dynamics fluctuates between stacked and side-by-side bound states. A quantitative characterization of this behavior is given in Figs. 5c, d that plot the P(z) and p(θ) distributions. One can see that the P(z) distributions in Fig. 5c are intermediate between those in Figs. 4e and i, corresponding to stacked and side-by-side bound pairs: in particular the probability distribution of torus 1 has a smaller amplitude than that in Fig. 4e and a shoulder that extends into the high probability region of torus 2. The p(θ) distribution is similar to those in Fig. 4k. This behavior can be rationalized by noting that stronger gravitational confinement suppresses orientational fluctuations and causes the tori to adopt nearly side-by-side configurations parallel to the wall, while larger ϵ cd values tend to induce tilt. Thus, smaller g and larger ϵ cd will favor the stacked bound state.
The bound states of two forward-moving tori near the wall may be compared to the dynamics of two such tori in bulk solution in the absence of a gravitational force. Stacked bound long-lived metastable states are also observed, but they have rather different characters. Figure 6a, b show the two predominate configurations we observed. The tori that stack with their catalytic surfaces in the same direction and propagate along the common axis of the two tori are indicated by the arrow in Fig. 6a. By contrast, Fig. 6b illustrates the other configuration where the catalytic surfaces point towards each other with the noncatalytic surfaces on the outside of the bound pair. Since a forward-moving torus propagates in the direction of its catalytic surface, this bound state does not propagate but only executes active Brownian motion with rotation. Similar bound states were observed experimentally in Baker et al. 23 , Supplementary Video 6 from that paper.
Four tori. A greater variety of clustered metastable states is observed for four tori near the wall. We confine our attention to the same gravitational and charge-dipole interactions as for two tori described above, in particular, we consider ϵ cd = 0.2 with g = 0.07 and g = 0.09 for forward-moving tori. Figure 7a-c shows images of two cluster configurations drawn from different realizations of the dynamics. In one of these clusters in panels (a) and (b), side and top views, respectively, one sees two side-by-side tori, where one torus is interacting with two tilted tori. A distinctive feature of this configuration is that one of these two tiled tori has flipped its orientation and its catalytic surface is now oriented towards the wall. (The dynamics of this bound state is shown in Supplementary Video 4). The other cluster shown in top view in panel (c) comprises a pair of stacked tori similar to those in Fig. 4h. Quantitative characterization of these clusters is given in Fig. 7d, e that show the P(z) and p(θ) probability densities. The P(z) distribution has a double-peaked structure characteristic of a stacked configuration but it more closely resembles that for ϵ cd = 0.1 with g = 0.07 in Fig. 5c. The p(θ) distribution is also similar to that in Fig. 5d, except that it extends to large angles and has a small maximum near θ ≈ 150 o indicating that some tori have flipped their orientation to the wall. Thus, the statistically predominate cluster is a bound pair of two stacked tori whose structure is modified due to interactions among the four tori.
Turning next to systems of four tori with g = 0.09 and ϵ cd = 0.2, we see in Fig. 8a, b, side and top views, respectively, clusters comprising stacked and side-by-side torus pairs, and (c) two staggered interacting side-by-side pairs. (The dynamics giving rise to the bound state in Fig. 8c is shown in Supplementary Video 5.) The P(z) and p(θ) probability distributions shown in Fig. 8d, e resemble those in Fig. 4i and k, indicating that clusters with two staggered interacting side-byside pairs dominate the dynamics.  Thus, for these two and four tori systems, a number of different long-lived bound metastable states exist whose structure depends both on parameters that control their interactions with the wall and on steric forces that reflect the torus geometry. The geometry of the torus influences the concentration and flow fields in its vicinity and their interactions with the walls and other tori. This gives rise to multi-torus configurations that differ from those of simple spheres or dimer motors.

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. 4a-d). Another important prediction from our work is that, depending on the gravitational strength, a transition from stacked bound states (Fig. 4e-h) to flattened side-by-side bound states (Fig. 4i-l) 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 existence of metastable bound states of different types for four tori that depend on the magnitudes of the gravitational and charge-dipole interactions, suggests that  interesting regular and more complex collective tori states could be found in systems with large numbers of interacting tori near a wall. 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
Details of the MPCD algorithm. The simulations are carried out using a mesoscopic particle-based method. The dynamics of the torus motors and fluid particles are explicitly taken into account by combining molecular dynamics with multiparticle collision dynamics for the interactions among the fluid particles. Chemical reactions on the catalytic part of the torus surface and in the bulk of the fluid are also incorporated in this dynamical method. On long distance and time scales, the method has been shown to yield the Navier-Stokes and reaction-diffusion equations with analytical expressions for the transport coefficients [38][39][40]47 . Consequently, all relevant properties such as local velocity and concentration fields, as well as torus dynamical properties, can be extracted directly from the simulations from averages of the particle dynamics. Additional details of the simulation methods are given below.
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 more than a single torus) and contains N S = N A + N B solvent particles with mass m. A torus is constructed from linked spherical beads as described earlier 41 . Briefly, a torus is made from N b spherical beads with radius r b linked to form a ring. Although there are soft repulsive potential interactions between the beads and solvent particles, the solvent cannot penetrate into a bead closer than r b , which is used to define its radius. The tori with a centerline radius a = 5 considered in this work were made from N b = 39 beads, and the segmented structure that results from the overlapping linked beads with radius r b is visible in Fig. 1a. The neighboring beads in this ring are connected by harmonic springs with spring constant k s = 500. The equilibrium bond length between the ith and (i + j)th beads is r The torus has mass M 39 = 3890. While, the frictional and propulsion properties of a torus in bulk solution we described earlier, here we consider situations that are similar to the experimental studies where one or more tori experience gravitational and chargedipolar interactions with a wall, which lead to new phenomena. The interactions the tori experience are given below.
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 38 . 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 point 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 . 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ϵ. Although the reactants and products interact with the torus beads with different potentials, the potentials we consider do not take into account the more generic situation where the interactions may differ for the catalytic and noncatlaytic faces. Such more general interactions have been considered in particle-based simulations collective Janus colloid dynamics 46 and for a Janus colloid interacting with a wall 27 . 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, 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.
For the simulations involving several tori, we also included short-range repulsive interactions between the beads i, j of different tori, for r ij ≤ r c where r c ¼ 2 1 6 σ m and zero for r ij > r c , where r ij = r i − r j is the vector from the center of torus bead i to that of bead j. We take σ m = 2(r b + σ), ϵ m = ϵ in order to prevent binding due to solvent exclusion effects.
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 probability p when an A or B particle crosses a reaction radius R r = r c . We have taken p = 1 but the reaction rate may be varied by changing the value of p. 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 where reactions are carried out probabilistically in each multiparticle collision cell, and the dynamics yields the reaction-diffusion equations on long distance and time scales 48 .
Simulation parameters. 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 49 . 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 ≈ 1.44 × 10 6 . For τ = 0.5 the fluid dynamic viscosity is η = 1.43, the kinematic viscosity is ν = η/ρ = 0.143, and the common diffusion coefficient for the A and B species is D = 0.117. For a torus with velocity V = 0.006, the Reynolds number is Re = Va/ν ≈ 0.2, the Peclet number is Pe = Va/D ≈ 0.25, the Damköhler number of Da = κ + a/D ≈ 2 signals that the reaction on the torus surface with rate constant κ + per unit surface ares is in the diffusion-influenced regime, and the Schmidt number is Sc ≈ 1.2 so that momentum transport is greater than mass transport. From the scaling that enters mesoscopic dynamics 50 , simulations of systems with such dimensionless numbers, where fluid particles interact with solute particles through intermolecular potentials, behave as low Reynolds number fluids with liquid-like properties. This has been confirmed by simulations with a multiparticle collision time of τ = 0.1 where the Reynolds number is an order of magnitude smaller and the Schmidt number is an order of magnitude larger 41 .
Results are reported in dimensionless units with mass m, length σ, energy ϵ, and time t 0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffi mσ 2 =ϵ p .

Data availability
Raw data were generated at SciNet HPC Consortium computers. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.
Code availability