---
abstract: |
    The lack of detailed balance in active colloidal suspensions allows dissipation to determine stationary states. Here we show that slow viscous flow produced by polar or apolar active colloids near plane walls mediates attractive hydrodynamic forces that drive crystallization. Hydrodynamically mediated torques tend to destabilize the crystal but stability can be regained through critical amounts of bottom-heaviness or chiral activity. Numerical simulations show that crystallization is not nucleational, as in equilibrium, but is preceded by a spinodal-like instability. Harmonic excitations of the active crystal relax diffusively but the normal modes are distinct from an equilibrium colloidal crystal. The hydrodynamic mechanisms presented here are universal and rationalize recent experiments on the crystallization of active colloids.\
    \
    DOI: [10.1103/PhysRevLett.117.228002](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.117.228002)\
author:
- Rajesh Singh
- 'R. Adhikari'
title: Universal hydrodynamic mechanisms for crystallization in active colloidal suspensions
---

In active colloidal suspensions [@palacci2013living; @petroff2015fast], energy is continuously dissipated into the ambient viscous fluid. The balance between dissipation and fluctuation that prevails in equilibrium colloidal suspensions [@einstein1905theory; @kubo1966fluctuation] is, therefore, absent. Nonequilibrium stationary states in active suspensions, then, are determined by both dissipative and conservative forces, quite unlike passive suspensions where detailed balance prevents dissipative forces from determining phases of thermodynamic equilibrium. In this context, it is of great interest to enquire how thermodynamic phase transitions driven by changes in free energy are modified in the presence of sustained dissipation.

In two recent experiments disordered suspensions of active colloids have been observed to spontaneously order into two-dimensional hexagonal crystals when confined at a plane wall. Bottom-heavy synthetic active colloids which catalyze hydrogen peroxide when optically illuminated are used in the first experiment [@palacci2013living] while chiral fast-swimming bacteria of the species *Thiovulum majus* are used in the second experiment [@petroff2015fast]. Given this remarkably similar crystallization in two disparate active suspensions it is natural to ask if the phenomenon is universal and to search for mechanisms, necessarily involving dissipation, that drive it.

Our current understanding of phase separation in particulate active systems is derived from the coarse-grained theory of motility-induced phase separation (MIPS) where active particles are advected by a density-dependent velocity [@tailleur2008statistical; @cates2010arrested; @cates2013active; @cates2015]. Microscopic models with kinematics consistent with MIPS also show phase separation and crystallization of hard active disks have been reported in two dimensions [@henkes2011active; @fily2012athermal; @bialke2012crystallization; @redner2013structure]. However, these models ignore exchange of the locally conserved momentum of the ambient fluid with that of the active particles and are, thus, best applied to systems where such exchanges can be ignored. Fluid flow is an integral part of the physics in [@palacci2013living; @petroff2015fast] and a momentum-conserving theory, currently lacking, is essential to identify the dissipative forces and torques that drive crystallization.

In this Letter we present a microscopic theory of active crystallization that connects directly to the experiments described above. Specifically, we account for the *three-dimensional* active flow in the fluid and the effect of a plane wall on this flow. Representing activity by slip in a thin boundary layer at the colloid surface [@ghose2014irreducible; @singh2014many; @singh2016traction] we rigorously compute the long-ranged many-body hydrodynamic forces and torques on the colloids. Thus we estimate Brownian forces and torques to be smaller than their active counterparts by factors of order $10^{2}$ (for synthetic colloids in [@palacci2013living]) to $10^{4}$ (for bacteria in [@petroff2015fast]) making them largely irrelevant for active crystallization. We integrate the resulting deterministic balance equations numerically to obtain dynamical trajectories.

Our main numerical results are summarized in Fig. (\[fig:Dynamics-of-crystallization\]). Panels (a)-(c) show the spontaneous destabilization of the uniform state by attractive active hydrodynamic forces, the formation of multiple crystallites, and their coalescence into a single hexagonal crystal at late times. Panels (d)-(f) show the structure factor at corresponding times. The route to crystallization is not through activated processes that produce critical nuclei, but through a spinodal-like instability produced by the unbalanced long-ranged active attraction. The uniform state is, therefore, always unstable and crystallization occurs for all values of density, in contrast to the finite density necessary for crystallization in MIPS models [@cates2015]. Active hydrodynamic torques tend to destabilize the ordered state but stability is regained when these are balanced by external torques (from bottom-heaviness in [@palacci2013living]) or by chiral activity (from bacterial spin in [@petroff2015fast]). Crystallites of chiral colloids rotate at an angular velocity that is inversely proportional to the number of colloids contained in them, as shown in panel (g). This is in excellent agreement with the experiment [@petroff2015fast]. The critical values of bottom-heaviness and chirality above which orientational stability, and, hence, positional order, is ensured is shown in panel (h). We now present our model and detail the derivation of our results.

![Panels (a)-(c) are instantaneous configurations during the crystallization of $1024$ active colloids of radius $b$ at a plane wall. The colloids are colored by their initial positions. Panels (d)-(f) show the structure factor $S(\mathbf{k})$ at corresponding instants. Wavenumbers are scaled by the modulus of the reciprocal lattice vector $k_{0}$ and the contribution from $\mathbf{k}=0$ is discarded. Panel (g) shows the variation of the angular velocity $\mathbf{\Omega}_{c}$ of a crystallite with the number $N$ of colloids in it. A typical configuration is shown in the inset. Panel (h) is the state diagram for orientational stability in terms of the measure of chirality $V_{0}^{(3a)}$ and bottom-heaviness $T_{0}$ (see text). Each dot represents one simulation. Here $v_{s}$ is the self-propulsion speed of an isolated colloid, $\tau=b/v_{s}$, and $\epsilon$ is the scale of the repulsive steric potential.\[fig:Dynamics-of-crystallization\]](Figure1){width="45.00000%"}

*Model:* We consider $N$ spherical active colloids of radius $b$ near a plane wall with center-of-mass coordinates $\mathbf{R}_{i},$ orientation $\mathbf{p}_{i}$, linear velocity $\mathbf{V}_{i}$, and angular velocity $\boldsymbol{\Omega}_{i}$, where $i=1\ldots N$. Activity is imposed through a slip velocity $\mathbf{v}_{i}^{\mathcal{A}}$ which is a general vector field on the surface $S_{i}$ of the $i$-th colloid satisfying $\int{\bm{\hat{\rho}}_{i}\cdot}\mathbf{v}_{i}^{\mathcal{A}}\,d\text{S}_{i}=0$ [@slipConstraint], where $\boldsymbol{\rho}_{i}$ is the vector from the center of the colloid to a point on its surface. The fluid velocity $\mathbf{v}$ is subject to slip boundary conditions $$\mathbf{v}(\mathbf{R}_{i}+\boldsymbol{\rho}_{i})=\mathbf{V}_{i}+\bm{\Omega}_{i}\times\bm{\rho}_{i}+\mathbf{v}_{i}^{\mathcal{A}}(\bm{\rho}_{i}).\label{eq:slip-RBM-BC}$$ on the colloid surfaces, to a no-slip boundary condition $\mathbf{v}=0$ at the plane wall located at $z=0$, and to a quiescent boundary condition at large distances from the wall. The slip is conveniently parametrized by an expansion $\mathbf{v}^{\mathcal{A}}(\mathbf{R}_{i}+\boldsymbol{\rho}_{i})=\sum_{l=1}^{\infty}\tfrac{1}{(l-1)!(2l-3)!!}\,\mathbf{V}_{i}^{(l)}\cdot\mathbf{Y}^{(l-1)}(\bm{\hat{\rho}}_{i})$ in irreducible tensorial spherical harmonics $\mathbf{Y}^{(l)}(\bm{\hat{\rho}})=(-1)^{l}\rho^{l+1}\bm{\nabla}{}^{(l)}\rho^{-1}$, where $\bm{\nabla}^{(l)}=\bm{\nabla}_{\alpha_{1}}\dots\bm{\nabla}_{\alpha_{l}}$. The expansion coefficients $\mathbf{V}_{i}^{(l)}$ are $l$-th rank reducible Cartesian tensors with three irreducible parts of ranks $l,$ $l-1,$ and $l-2$, corresponding to symmetric traceless, antisymmetric and pure trace combinations of the reducible indices. We denote these by $\mathbf{V}_{i}^{(ls)}$, $\mathbf{V}_{i}^{(la)}$ and $\mathbf{V}_{i}^{(lt)}$ respectively. The leading contributions from the slip, $$\begin{aligned}
\mathbf{v}_{i}^{\mathcal{A}}(\bm{\rho}_{i}) & = & \underbrace{-\mathbf{V}_{i}^{\mathcal{A}}+\tfrac{1}{15}\mathbf{V}_{i}^{(3t)}\cdot\mathbf{Y}^{(2)}}_{\mathrm{achiral,\,polar}}-\underbrace{\tfrac{1}{9}\boldsymbol{\varepsilon}\cdot\mathbf{V}_{i}^{(3a)}\cdot\mathbf{Y}^{(2)}}_{\mathrm{chiral,\thinspace apolar}}\nonumber \\
 & + & \underbrace{\mathbf{V}_{i}^{(2s)}\hspace{-0.1cm}\cdot\mathbf{Y}^{(1)}}_{\mathrm{achiral,\,apolar}}\hspace{-0.12cm}-\underbrace{\mathbf{\Omega}_{i}^{\mathcal{A}}\hspace{-0.05cm}\times\hspace{-0.05cm}\bm{\rho}_{i}-\tfrac{1}{60}\boldsymbol{\varepsilon}\cdot\mathbf{V}_{i}^{(4a)}\hspace{-0.12cm}\cdot\mathbf{Y}^{(3)}}_{\mathrm{chiral,\,polar}}\label{eq:slip-truncation}\end{aligned}$$ have coefficients of polar, apolar and chiral symmetry. Here $\boldsymbol{\varepsilon}$ is the Levi-Civita tensor. The retained modes have physical interpretations: for a single colloid in an unbounded fluid, $\mathbf{V}^{\mathcal{A}}$ ($l\sigma=1s)$ and $\bm{\Omega}^{\mathcal{A}}$ $(l\sigma=2a)$ are the linear and angular velocities in the absence of external forces and torques, $\mathbf{V}^{(2s)}$ is the active contribution to the stresslet, while $\mathbf{V}^{(3a)},\mathbf{V}^{(3t)}$, and $\mathbf{V}^{(4a)}$ are strengths of the chiral torque dipole, polar vector quadrupole, and chiral octupole respectively. The tensors are parametrized uniaxially, $\mathbf{V}_{i}^{\mathcal{A}}=v_{s}\mathbf{p}_{i}$, $\mathbf{\boldsymbol{\Omega}}_{i}^{\mathcal{A}}=\omega_{s}\mathbf{p}_{i}$, $\mbox{\ensuremath{{\bf V}_{i}^{(2s)}=V_{0}^{(2s)}(\mathbf{p}_{i}\mathbf{p}_{i}-\tfrac{\mathbf{I}}{3})}}$ and so on, where $v_{s}$ and $\omega_{s}$ are the speeds of active translation and rotation and $V_{0}^{(2s)}$ positive (negative) corresponds to a pusher (puller). [@supplementalM].

The synthetic active colloids in [@palacci2013living] are polar and achiral (they self-propel but do not spin) while the bacteria in [@petroff2015fast] are polar and chiral (they self-propel and spin). Both these cases are included in the leading contributions. In [@ghose2014irreducible] a procedure is outlined for estimating the leading coefficients from experimentally measured flows and it is shown that the active flow produced by flagellates and green algae can be modeled by slip. Our model is of sufficient generality, then, to include both synthetic and biological active colloids, and situations where swirling and time-dependent slip may be necessary [@drescher2010; @drescher2011; @guasto2010; @goldstein2015green].

![image](Figure2){width="90.00000%"}

*Active forces and torques:* Newton’s equations of motion for the colloids reduce, when inertia is negligible, to instantaneous balance of forces and torques $$\mathbf{F}_{i}^{H}+\mathbf{F}_{i}^{P}+\boldsymbol{\xi}_{i}^{T}=0,\quad\mathbf{T}_{i}^{H}+\mathbf{T}_{i}^{P}+\boldsymbol{\xi}_{i}^{R}=0.\label{eq:Newton}$$ Here $\mathbf{F}_{i}^{H}=\int\mathbf{f}\,d\text{S}_{i}$, $\mathbf{F}^{P}$ and $\mathbf{\boldsymbol{\xi}^{T}}$ are respectively the hydrodynamic, body and Brownian forces while, $\mathbf{T}_{i}^{H}=\int\boldsymbol{\rho}_{i}\mathbf{\times}\mathbf{f}\,d\text{S}_{i}$, $\mathbf{T}_{i}^{P}$ and $\boldsymbol{\xi}_{i}^{R}$ are, corresponding torques, $\boldsymbol{\sigma}$ is the Cauchy stress in the fluid and $\mathbf{f}=\bm{\hat{\rho}_{i}}\cdot\boldsymbol{\sigma}$ is the traction. The linearity of the Stokes equation implies that these must be of the form

\[force-formulation\] $$\begin{aligned}
{1}
\mathbf{F}_{i}^{H}= & -\boldsymbol{\gamma}_{ij}^{TT}\mathbf{\cdot V}_{j}-\boldsymbol{\gamma}_{ij}^{TR}\mathbf{\cdot\boldsymbol{\Omega}}_{j}-\sum_{l\sigma=1s}^{\infty}\boldsymbol{\gamma}_{ij}^{(T,\,l\sigma)}\negmedspace\cdot\mathbf{V}_{j}^{(l\sigma)},\label{eq:linear-force-torque}\\
\mathbf{T}_{i}^{H}= & -\boldsymbol{\gamma}_{ij}^{RT}\mathbf{\cdot V}_{j}-\boldsymbol{\gamma}_{ij}^{RR}\mathbf{\cdot\boldsymbol{\Omega}}_{j}-\sum_{l\sigma=1s}^{\infty}\boldsymbol{\gamma}_{ij}^{(R,\,l\sigma)}\negmedspace\cdot\mathbf{V}_{j}^{(l\sigma)},\label{eq:torque-expression}\end{aligned}$$

where repeated particle indices are summed over. The $\boldsymbol{\gamma}_{ij}^{\alpha\beta}$, with $\alpha,\beta=T,R$, are the usual friction matrices associated with rigid body motion and $\mathcal{\boldsymbol{\gamma}}_{ij}^{(\alpha,\,l\sigma)}$ are friction tensors associated with the irreducible modes of the active slip. They are of rank $l+1$, $l$, and $l-1$, respectively, for $\sigma=s,a,t$. The forces and torques depend on relative position (through the $\mathcal{\boldsymbol{\gamma}}_{ij}^{(\alpha,\,l\sigma)}$) and on relative orientation (through the $\mathbf{V}_{j}^{(l\sigma)}$). Their signature under time-reversal shows that the active contributions are dissipative.

We calculate the friction tensors using a Galerkin discretization of the boundary integral equation [@singh2014many; @singh2016traction] with the Lorentz-Blake Green’s function [@blake1971c] which, by construction, vanishes at the plane wall. The $\boldsymbol{\gamma}_{ij}^{(T,\,l\sigma)}$ decay as $r_{ij}^{-(l+1)}$ and $r_{ij}^{-(l+2)}$ in the directions parallel and perpendicular to the wall. The $\boldsymbol{\gamma}_{ij}^{(R,\,l\sigma)}$ decay one power of $r_{ij}$ more rapidly. While the force and torque so obtained are sufficient to study colloidal motion, additional insight is obtained from studying the flow, which we compute from its boundary integral representation. Further details are given in [@supplementalM].

The modes $l\sigma=1s$ and $l\sigma=2a$ contribute most dominantly to forces and torques and they attain their lower bounds far away from the wall, where their magnitudes are $\textcolor{black}{\ensuremath{F=6\pi\eta bv_{s}}}$ and $T=8\pi\eta b^{3}\omega_{s}$. The bacteria in [@petroff2015fast] have radius $b\sim4\,\mu\text{m}$, swimming speed $v_{s}\sim500\,\mu\text{m/s}$ and angular speed $\omega_{s}\sim50\,s^{-1}$ in a fluid of viscosity $\eta=10^{-3}\,\text{kg/ms}$. This gives an estimate of $F\sim40\times10^{-12}\,\text{N}$ and$T\sim10^{-16}$ Nm. For the synthetic colloids in [@palacci2013living], $b\sim2\,\mu\text{m}$, $v_{s}\sim10\,\mu\text{m/s}$, which corresponds to $F_{\mathcal{}}\sim10^{-13}\,\text{N}$. Typical Brownian forces and torques are of order $\mathcal{O}\left(k_{B}\text{T/b}\right)\sim10^{-15}\,\text{N}$, and $\mathcal{O}\left(k_{B}\text{T}\right)\sim10^{-21}\,\text{Nm}$ respectively. Thus active forces and torques overwhelm Brownian contributions by factors of 100 or more in these experiments and, henceforth, we neglect their effects. Trajectories are obtained by integrating the kinematic equations $\mathbf{\dot{R}}_{i}=\mathbf{V}_{i}$ and $\dot{\mathbf{p}}_{i}=\mathbf{\Omega}_{i}\times\mathbf{p}_{i}$, where $\mathbf{V}_{i}$ and $\mathbf{\boldsymbol{\Omega}}_{i}$ satisfy Eq. (\[eq:Newton\]) with Brownian contributions removed. Integration methods and parameter choices are detailed in [@supplementalM].

*Crystallization kinetics:* The kinetics of crystallization obtained from numerical solutions is shown in Movie 1 [@supplementalM], together with the evolution of the structure factor $S(\mathbf{k})$. The uniform state is destabilized, most notably for any initial density, by attractive active hydrodynamic forces. Steric repulsion between particles balances these to produce crystallites with hexagonal positional order. Rings in the structure factor first appear at wavenumbers that correspond to Bragg vectors of the lattice, reminiscent of a spinodal instability, representing the averaged scattering from randomly oriented crystallites. These sharpen into Bragg peaks as the crystallites coalesce and orientational order grows. Finally particles assemble into a single crystallite which continues to rotate, while the structure factor shows a clear sixfold symmetry. In Movie 2 [@supplementalM] we show the formation of a hexagonal unit cell from the simulation of seven polar and chiral active colloids. The crystallite rotates with an angular velocity parallel to the chiral axis of the colloids.

*Universal mechanisms:* To better understand the mechanisms behind active crystallization we show, in Fig. (\[fig:distortion-of-active-flow\]) , the active flow near a wall and the dominant contributions to the flow-mediated forces and torques. The top three panels show the increasing distortion of the flow produced by the leading polar $(l\sigma=1s)$ and apolar $(l\sigma=2s)$ modes for $\mathbf{p}_{i}$ normal to the wall and $V_{0}^{(2s)}<0.$ The flow develops a monopolar character as the colloid is brought to rest at a height $h$ by the balance of hydrodynamic attraction, Fig. (\[fig:distortion-of-active-flow\]d), and steric repulsion from the wall. The induced monopole on the colloids leads to attractive forces between them below a critical height $h$ from the wall as shown in Fig. (\[fig:distortion-of-active-flow\]f). Nearby colloids entrained in this flow are attracted towards the central colloid as shown in the rightmost panel and in Movie 3 [@supplementalM]. The balance of the hydrodynamic attraction and steric repulsion determines the lattice spacing $d$. We note that even an apolar colloid is attracted to the wall, Fig. (\[fig:distortion-of-active-flow\]d), and induces hydrodynamic attractive forces. Thus, unlike MIPS [@cates2015], polarity is not necessary for crystallization. The induced monopole also tends to reorient the colloids, by generating a torque in the plane of wall, as shown by the curved red arrows in Fig. (\[fig:distortion-of-active-flow\]c) and quantified in Fig. (\[fig:distortion-of-active-flow\]g). **[@supplementalM].** Activity *and* body forces pointing away from the wall are necessary for positional order while bottom-heaviness *or* chirality is necessary for orientational stability.

*Harmonic excitations:* We now study harmonic excitations $\mathbf{u}_{i}$ of a perfect hexagonal crystal by expanding the positions as $\mathbf{R}_{i}=\mathbf{R}_{i}^{0}+\mathbf{u}_{i}$ around the stationary state $\mathbf{R}_{i}^{0}=(X_{i}^{0},\,Y_{i}^{0},\,h)$ and ignoring orientational fluctuations. Force balance to leading order gives $$\begin{aligned}
-\boldsymbol{\gamma}_{ij}^{TT}\cdot\dot{\mathbf{u}}_{j} &  & +\left(\mathbf{\boldsymbol{\nabla}}_{j}\boldsymbol{\gamma}_{ij}^{TT}\cdot\mathbf{V}^{\mathcal{A}}-\mathbf{D}_{ij}\right)\cdot\mathbf{u}_{j}=0,\label{eq:stability-eq-ac}\end{aligned}$$ where $\mathbf{D}_{ij}=-\boldsymbol{\nabla}_{j}\boldsymbol{\nabla}{}_{i}U\big|_{0}$ and $U$ is the sum of all steric potentials. **This shows that relaxation is determined by both activity and elasticity, unlike in an equilibrium colloidal crystal where elasticity alone relaxes strains. The normal modes of relaxation can be obtained by Fourier transforming in the plane and in time. The dispersion is found from the solutions of** $$\det\Big|-i\omega\boldsymbol{\gamma}_{\mathbf{k}}^{TT}+i\mathbf{k}\,\boldsymbol{\gamma}_{\mathbf{k}}^{TT}\cdot\mathbf{V}^{\mathcal{A}}-\mathbf{D}_{\mathbf{k}}\Big|=0.\label{eq:det-stability}$$ **Here $\mathbf{k}=(k_{1},\,k_{2})$ is the wavevector restricted to the first Brillouin zone [@supplementalM], $\omega$ is the frequency and $\mathbf{D}_{\mathbf{k}}$** is the dynamical matrix**. The pair of dispersion relations for motion parallel to the wall are shown in Fig. (\[fig:Normal-modes\]). The dispersion for $k\ll k_{0}$, where $k_{0}$ is the magnitude of the reciprocal lattice vector, is quadratic in wavenumber $$\omega_{\pm}=-i\tfrac{\gamma_{\perp}^{T}hv_{s}}{2\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}}f_{\pm}(\theta)\,k^{2},\label{eq:small-k}$$ where $f_{\pm}(\theta)$ are angular factors, $\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}$ and $\gamma_{\perp}^{T}$ are one-body frictions parallel and perpendicular to the wall, and $\tan\theta=\tfrac{k_{2}}{k_{1}}$. The small-$k$ approximation is compared with the numerical solution in Fig. (\[fig:Normal-modes\]) and it is found to hold for $k\lesssim0.1k_{0}$. These can be interpreted as overdamped phonon modes of the active crystal [@joannyDiscuss]. The presence of the active term** $i\mathbf{k}\,\boldsymbol{\gamma}_{\mathbf{k}}^{TT}\cdot\mathbf{V}^{\mathcal{A}}$ in Eq. (\[eq:det-stability\]) makes them differ from phonon modes of a colloidal crystal.

![Branches of the dispersion relation for the two planar normal modes of relaxation of a hexagonal active crystal. The curves in upper panel show the dispersion along high symmetry directions of the Brillouin zone (first inset). The surfaces in the second and third insets show the dispersion over the entire Brillouin zone. Polar plots in the lower panel, have comparisons of full numerical solution of Eq. (\[eq:det-stability\]) with the approximate solution at small $k$ of Eq. (\[eq:small-k\]) for $k=0.01k_{0}$ (left panel) and $k=0.3k_{0}$ (right panel).\[fig:Normal-modes\]](Figure3){width="47.00000%"}

*Discussion*: In this work, we have considered only hydrodynamic forces and torques, unlike the case of MIPS [@tailleur2008statistical; @cates2010arrested; @cates2013active; @cates2015] where Brownian torques drives reorientations [@henkes2011active; @fily2012athermal; @bialke2012crystallization; @redner2013structure]. We have shown that the latter are at least two orders of magnitude weaker than the former for experiments in the class of [@palacci2013living; @petroff2015fast]. However, it is conceivable that thermal fluctuations will play a more significant role when the activity is comparatively weak, modifying both the nature of crystallization transition and the stability of the crystalline phase. The spinodal-like instability appears due to the uncompensated long-ranged attractive active forces. These can be compensated by entropic forces to stabilize the disordered phase at finite temperatures. A nucleational route to crystallization, with activity-enhanced rates, is then possible in the regime where the active forces reduce the nucleation barrier without driving it to zero. In the crystalline phase, thermal fluctuations will excite both phonon and topological modes. Phonon fluctuations will destroy long-range translational order [@peierls1935quelques; @landau1937theorie], but due to the activity-enhanced stiffness of these modes, large system sizes (compared to equilibrium) will be needed to observe the power-law decay of correlations. Topological defects will be excited at higher temperatures and a defect unbinding transition [@kosterlitz1973ordering; @halperin1978theory; @nelson1979theory; @young1979; @chaikin2000principles], modified by activity, may destroy translational order entirely, producing instead an active hexatic phase. These present exciting avenues for future research. We remark that wall-bounded clustering phenomena in algae [@drescher2009] and charged colloids [@squires2001effective] are mediated by specific forms of the universal hydrodynamic mechanisms presented here.

Finally, we suggest that the flow-induced phase separation found here may provide a paradigm, complementary to MIPS, in which theoretical and experimental studies of momentum-conserving driven [@yeo2015collective] and active matter [@trau1996field; @solomentsev1997particle; @matas2014hydrodynamic; @pandey2014flow; @wang2015one; @wykes2016dynamic] may be situated.

We thank M. E. Cates, P. Chaikin, D. Frenkel, D. J. Pine, A. Laskar and T. V. Ramakrishnan for helpful discussions and IMSc for computing resources on the Nandadevi clusters.

[54]{}ifxundefined \[1\][ ifx[\#1]{} ]{}ifnum \[1\][ \#1firstoftwo secondoftwo ]{}ifx \[1\][ \#1firstoftwo secondoftwo ]{}““\#1””@noop \[0\][secondoftwo]{}sanitize@url \[0\][‘\
12‘\$12 ‘&12‘\#12‘12‘\_12‘%12]{}@startlink\[1\]@endlink\[0\]@bib@innerbibempty [****,  ()](\doibase
  10.1126/science.1230020) [****,  ()](\doibase 10.1103/PhysRevLett.114.158102) [****,  ()](\doibase http://dx.doi.org/10.1002/andp.19053220806) [****,  ()](http://stacks.iop.org/0034-4885/29/i=1/a=306) [****,  ()](\doibase 10.1103/PhysRevLett.100.218103) [****, ()](\doibase dx.doi.org/10.1073/pnas.1001994107) [****,  ()](http://iopscience.iop.org/article/10.1209/0295-5075/101/20010/meta) [****,  ()](\doibase 10.1146/annurev-conmatphys-031214-014710) [****,  ()](\doibase 10.1103/PhysRevE.84.040301) [****,  ()](\doibase 10.1103/PhysRevLett.108.235702) [****,  ()](\doibase 10.1103/PhysRevLett.108.168301) [****,  ()](\doibase 10.1103/PhysRevLett.110.055701) [****,  ()](\doibase 10.1103/PhysRevLett.112.118102) [****, ()](http://stacks.iop.org/1742-5468/2015/i=6/a=P06017) @noop [ ()]{} @noop [  ]{} @noop [  ]{} [****,  ()](\doibase
  10.1103/PhysRevLett.105.168101) [****,  ()](\doibase
  10.1073/pnas.1019079108) [****, ()](\doibase 10.1103/PhysRevLett.105.168102) @noop [****,  ()]{} [****,  ()](\doibase 10.1017/S0305004100049902) @noop [ ]{} @noop [ ****,  ()]{} @noop [****,  ()]{} [****,  ()](http://iopscience.iop.org/article/10.1088/0022-3719/6/7/010/meta) [****,  ()](\doibase 10.1103/PhysRevLett.41.121) [****,  ()](\doibase 10.1103/PhysRevB.19.2457) [****, ()](\doibase 10.1103/PhysRevB.19.1855) @noop [**]{}, Vol.  (, ) [****,  ()](\doibase
  10.1103/PhysRevLett.102.168101) [****,  ()](\doibase 10.1017/S0022112001005432) [****,  ()](\doibase 10.1103/PhysRevLett.114.188301) [****,  ()](\doibase 10.1126/science.272.5262.706) [****,  ()](\doibase 10.1021/la970294a) [****,  ()](\doibase 10.1103/PhysRevE.90.032304) [****, ()](\doibase 10.1039/C6SM02104B) [****,  ()](\doibase
  10.1021/acs.accounts.5b00025) [****,  ()](\doibase
  10.1039/C5SM03127C) [****,  ()](\doibase 10.1016/0021-9797(79)90336-9) [****,  ()](\doibase 10.1007/BF01194638) @noop [**]{}, Mathematics and its applications (, ) @noop [**]{} (, ) @noop [**]{} (, ) [****,  ()](\doibase 10.1016/0378-4371(82)90127-3) @noop [**]{} (, ) [****,  ()](\doibase 10.1063/1.454658) [****, ()](\doibase 10.1007/BF01517501) [****,  ()](\doibase 10.1016/0378-4371(80)90067-9) [****,  ()](\doibase 10.1063/1.1674820) @noop [**]{} (, ) [[,]{} ](https://github.com/rajeshrinet/pystokes) () [****,  ()](\doibase 10.1103/PhysRevE.74.041127) [ [,]{} ](http://hplgit.github.com/odespy/doc/tutorial/html/wrap_odespy.html#using-adaptive-methods) ()

****

Active force, torque and flow\[sec:propulsion-matrices-in-wall-bounded\]
========================================================================

We derive, in this section, the expressions for the active forces, torques, and exterior flow in a suspension of $N$ active colloids bounded by a plane wall. The system of coordinates is shown in Fig. (\[fig:coordinate-system\]). The spheres are centered at $\mathbf{R}_{i}$ and their velocities and angular velocities are $\mathbf{V}_{i}$ and $\mathbf{\Omega}_{i}$ respectively. ${\bf p}_{i}$ denotes the orientation of the particle while points on the boundaries of the spheres is given by $\mathbf{s}_{i}=\mathbf{R}_{i}+\boldsymbol{\rho}_{i}$, where $\boldsymbol{\rho}_{i}$ is the radius vector. To ensure no-slip on the wall, we associate am image centered at ${\bf R}_{i}^{*}$ with the $i$-th colloid [@blake1971c], and a similar correspondence for all other quantities of the colloid and its image.

We closely follows our previous work [@singh2014many; @singh2016traction] where a boundary integral formulation has been used to solve the Stokes equation with arbitrary boundary conditions. The principal difference here is in the choice of the Green’s function which satisfies the no-slip condition at the plane wall [@blake1971c]. In the interest of being self-contained, we repeat certain key steps en route to the solution. A clear expression of the linearity of Stokes flow is found in its integral representation, where the flow in the bulk is given in terms of integrals of the tractions and velocities at the boundaries [@fkg1930bandwertaufgaben; @ladyzhenskaya1969; @pozrikidis1992; @kim2005], $$\begin{aligned}
{1}
v_{\alpha}(\mathbf{r})= & -\int G_{\alpha\beta}^{W}(\mathbf{r},\,\mathbf{s}_{j})f_{\beta}(\mathbf{s}_{j})\,d\mathrm{S}_{i}\nonumber \\
+ & \int K_{\beta\alpha\gamma}^{W}(\mathbf{r},\,\mathbf{s}_{j})\hat{\rho}_{\gamma}v_{\beta}(\mathbf{s}_{j})\,d\mathrm{S}_{i},\label{eq:BIE}\end{aligned}$$ where repeated particle indices are summed over, $\mathbf{s}_{j}=\mathbf{R}_{j}+\boldsymbol{\rho}_{j}$ is a point on the surface of $j$-th particle and $G_{\alpha\beta}^{W}(\mathbf{r},\,\mathbf{s}_{j})$ is the Green’s function of the Stokes system satisfying no-slip condition, $\mathbf{v}=0$ on the wall at $z=0$. The stress tensor $K_{\alpha\beta\gamma}^{W}(\mathbf{r},\,\mathbf{s}_{j})$ and the pressure vector $P_{\alpha}^{W}$($\mathbf{r},\,\mathbf{s}_{j})$ satisfy $K_{\alpha\beta\gamma}^{W}(\mathbf{r},\,\mathbf{s}_{j})=-\delta_{\alpha\gamma}P_{\beta}^{W}+\eta\left(\nabla_{\gamma}G_{\alpha\beta}^{W}+\nabla_{\alpha}G_{\beta\gamma}^{W}\right)$ and $-\nabla_{\alpha}P_{\beta}^{W}(\mathbf{r},\mathbf{r'})+\eta\nabla^{2}G_{\alpha\beta}^{W}=-\delta\left(\mathbf{r}-\mathbf{r'}\right)\delta_{ij}$ respectively.

We solve the Fredholm integral equation of Eq. (\[eq:BIE\]) by expanding the boundary fields in irreducible tensorial spherical harmonics, $\mathbf{Y}^{(l)}$, which are orthogonal basis function on the surface of the sphere $\frac{1}{4\pi b^{2}}\int\mathbf{Y}^{(l)}(\widehat{\bm{\rho}})\,\mathbf{Y}^{(l')}(\widehat{\bm{\rho}})\,d\mathrm{S}=\delta_{ll'}\,\frac{l!\,(2l-1)!!}{(2l+1)}\mathbf{\Delta}^{(l)},$ where $\mathbf{\Delta}^{(l)}$ is tensor of rank $2l$, projecting any $l$-th order tensor to its symmetric irreducible form [@mazur1982; @hess1980formeln]. The boundary velocity including active slip and its expansion in this basis has been provided above. The orthogonality of the basis functions gives the expansion coefficients in terms of surface integrals of traction and velocity as [@ladd1988; @ghose2014irreducible], $$\begin{aligned}
{1}
\mathbf{F}_{i}^{(l)} & =\frac{1}{(l-1)!(2l-3)!!}\int\mathbf{f}(\mathbf{R}_{i}+\bm{\rho}_{i})\mathbf{Y}^{(l-1)}(\bm{\hat{\rho}}_{i})\,d\mathrm{S}_{i},\nonumber \\
\mathbf{V}_{i}^{(l)} & =\frac{2l-1}{4\pi b^{2}}\int\mathbf{v}^{\mathcal{A}}(\mathbf{R}_{i}+\bm{\rho}_{i})\mathbf{Y}^{(l-1)}(\bm{\hat{\rho}}_{i})\,d\mathrm{S}_{i}.\end{aligned}$$ The coefficients of the traction and velocity are tensors of rank $l$ and can be written as irreducible tensor of rank $l,$ $l-1$ and $l-2$ [@singh2014many]. The first term in the traction expansion is the force $\mathbf{F}_{i}^{(1)}=\mathbf{F}_{i}^{H}$, while the antisymmetric part of the second term is the torque $b\boldsymbol{\varepsilon}\cdot\mathbf{F}_{i}^{(2)}=\mathbf{T}_{i}^{H}$. The first term in the velocity expansion is $\mathbf{V}_{i}^{(1)}=-\mathbf{V}_{i}^{\mathcal{A}}$ and the antisymmetric part of the second term is $\tfrac{1}{2b}\boldsymbol{\boldsymbol{\varepsilon}}\cdot\mathbf{V}_{i}^{(2)}=-\mathbf{\Omega}_{i}^{\mathcal{A}}$. Here $\mathbf{V}_{i}^{\mathcal{A}}=-\tfrac{1}{4\pi b^{2}}\int\mathbf{v}^{\mathcal{A}}(\bm{\rho}_{i})\,d\mathrm{S}_{i}$ denotes the self-propulsion, while $\bm{\Omega}_{i}^{\mathcal{A}}=-\frac{3}{8\pi b^{4}}\int\bm{\rho}_{i}\times\mathbf{v}^{\mathcal{A}}(\bm{\rho}_{i})\,d\mathrm{S}_{i}$ denotes the self-rotation of an isolated active colloid in unbounded flow. The expression for fluid flow can be obtained in terms of coefficients of traction and velocity, $$\begin{aligned}
{1}
\mathbf{v}(\mathbf{r})=\sum_{l=1}^{\infty}\Big(-\boldsymbol{G}_{j}^{(l)}\cdot\mathbf{F}_{j}^{(l)}+ & \boldsymbol{K}_{j}^{(l)}\cdot\mathbf{V}_{j}^{(l)}\Big),\label{eq:fluid-flow}\end{aligned}$$ where the boundary integrals $\boldsymbol{G}_{j}^{(l)}$ and $\boldsymbol{K}_{j}^{(l)}$ can be written in terms of Green’s function and its derivatives (Appendix \[sec:Expression-for-boundary\]). We multiply the fluid velocity by the $l$-th tensorial harmonic and integrate over the $i$-th boundary. Using the orthogonality of these basis functions, we obtain an infinite-dimensional linear system of equations for the unknown traction coefficients [@singh2014many], $$\begin{aligned}
{1}
\tfrac{1}{2}\,\mathbf{V}_{i}^{(l)}=\sum_{l'=1}^{\infty}\Big(-\boldsymbol{G}_{ij}^{(l,\,l')}\cdot\mathbf{F}_{j}^{(l')} & +\boldsymbol{K}_{ij}^{(l,\,l')}\cdot\mathbf{V}_{j}^{(l')}\Big),\label{eq:linear-system}\end{aligned}$$ where the matrix elements $\boldsymbol{G}_{ij}^{(l,\,l')}$ and $\boldsymbol{K}_{ij}^{(l,\,l')}$ can be evaluated in terms of the Green’s function and its derivatives, as given in Appendix \[sec:Expression-for-boundary\].

The traction and velocity coefficients are reducible and their irreducible decomposition is given as [@brunn1976effect; @schmitz1980force], $$\begin{aligned}
\mathbf{F}_{i}^{(ls)}=\overbracket[0.7pt][2.0pt]{\mathbf{F}_{i}^{(l)}},\quad & \mathbf{F}_{i}^{(la)}=\overbracket[0.7pt][2.0pt]{\bm{\varepsilon}\cdot\mathbf{F}_{i}^{(l)}},\quad & \mathbf{F}_{i}^{(lt)}={\bm{\delta}\cdot\mathbf{F}_{i}^{(l)}},\\
\mathbf{V}_{i}^{(ls)}=\overbracket[0.7pt][2.0pt]{\mathbf{V}_{i}^{(l)}},\quad & \mathbf{V}_{i}^{(la)}=\overbracket[0.7pt][2.0pt]{\boldsymbol{\mathbf{\varepsilon}}\cdot\mathbf{V}_{i}^{(l)}},\quad & \mathbf{V}_{i}^{(lt)}={\bm{\delta}\cdot\mathbf{V}_{i}^{(l)}}.\end{aligned}$$ Here the operator $\overbracket[0.7pt][2.0pt]{(\dots)}=\boldsymbol{\Delta}^{(l)}(\dots)$ extracts the symmetric irreducible part of the tensor it acts on. We use these irreducible coefficients and the linear system of equations to solve for the unknown traction in terms of the known boundary velocity [@singh2016traction]. The relations between the irreducible coefficients of the traction and velocity, then, becomes [@singh2016traction], $$\begin{aligned}
{1}
\mathbf{F}_{i}^{(l\sigma)}= & -\boldsymbol{\gamma}_{ij}^{(l\sigma,\,1s)}\cdot\mathbf{V}_{j}-\boldsymbol{\gamma}_{ij}^{(l\sigma,\,R)}\cdot\mathbf{\Omega}_{j}\nonumber \\
- & \sum_{l'\sigma'=1s}^{\infty}\boldsymbol{\gamma}_{ij}^{(l\sigma,\,l'\sigma')}\cdot\mathbf{V}_{j}^{(l'\sigma')}.\label{eq:main-traction-l}\end{aligned}$$ This infinite set of equations, called the traction laws [@singh2016traction], manifestly shows the linear relation between the traction and velocity coefficients and defines the friction tensors. The expressions for the friction tensors can be obtained by an iterative scheme [@singh2016traction]. We use the one-body solution as the initial guess for the iteration, $$\mathbf{F}_{i}^{H}=-\gamma^{T}(\mathbf{V}_{i}-\mathbf{V}_{i}^{\mathcal{A}}),\qquad\mathbf{T}_{i}^{H}=-\gamma^{R}(\mathbf{\Omega}_{i}-\mathbf{\Omega}_{i}^{\mathcal{A}}),$$ where $\gamma^{T}$ and $\gamma^{R}$ are one particle friction corresponding to translation and rotation. Near a wall no-slip wall, they are, ${\gamma}^{T}=\gamma_{\perp}^{T}\hat{\boldsymbol{z}}+\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}(\hat{\boldsymbol{x}}+\hat{\boldsymbol{y}})$ and ${\gamma}^{R}=\gamma_{\perp}^{R}\hat{\boldsymbol{z}}+\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{R}(\hat{\boldsymbol{x}}+\hat{\boldsymbol{y}})$ [@kim2005]. Here ${{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}$ and $\perp$ subscripts indicating directions parallel and perpendicular to the wall. The expressions after one iteration, corresponding to the “first reflection” in Smoluchowski’s classical method, are shown in Appendix \[appendix:Evaluation-of-gamma\].

![Coordinate system used to describe active spherical particles and its images near a no-slip wall. The $i$-th particle and its image is shown. See text for description. \[fig:coordinate-system\]](SI-Figure1){width="47.00000%"}

Crystalline steady states \[sec:Crystalline-steady-states\]
===========================================================

In this section we work out the steady states of the active crystals using the leading terms of the force and torque equations. Using the leading order force balance for $i$-th particle, the steady state condition for position is given as $$\begin{aligned}
\bm{\gamma}_{ij}^{TT}\cdot\mathbf{V}_{j}^{\mathcal{A}}+\mathbf{F}_{i}^{P}=0.\label{eq:minimal-steady-state}\end{aligned}$$ Here $\mathbf{V}_{i}^{\mathcal{A}}=-v_{s}\hat{\mathbf{z}}$ is self-propulsion of the colloid at a speed $v_{s}$, assumed to be moving $\perp$ to the wall. The body force $\mathbf{F}_{i}^{P}=-\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}U$ is due to a short-ranged repulsive potential $U,$ which depends on displacement $\mathbf{r}_{ij}=\mathbf{R}_{i}-\mathbf{R}_{j}$ and is given as, $U(r_{ij})=\epsilon\left(\frac{r_{min}}{r_{ij}}\right)^{12}-2\epsilon\left(\frac{r_{min}}{r_{ij}}\right)^{6}+\epsilon,$ for $r_{ij}<r_{min}$ and zero otherwise [@weeks1971role], where $\epsilon$ is the potential strength. The same potential has been used to model colloid-colloid repulsion $\mathbf{F}_{i}^{PP}$ and the colloid-wall repulsive force $\mathbf{F}_{i}^{PW}$.

*One- and two-body dynamics:* To estimate the height at which the particle is brought to rest close to the wall, we use the $z-$component of the force balance, $-\gamma_{\perp}^{T}v_{s}=F_{3}^{PW}$. Here $F_{3}^{PW}$ is the repulsive force from the wall in $\mathbf{\hat{z}}$ direction, while $\gamma_{\perp}^{T}v_{s}$ is the attractive force of the colloid to the wall in same direction. The balance between the attraction and repulsion sets the height $h$ at which the colloid is brought to rest. We now consider force balance for a pair of particles in planar direction, $$\begin{aligned}
{1}
-v_{s}\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}\gamma_{\perp}^{T}\mbox{\ensuremath{\mathcal{F}}}_{i}^{0}\mbox{\ensuremath{\mathcal{F}}}_{j}^{0}G_{\alpha3}^{W}(\mathbf{R}_{i},\mathbf{R}_{j})+F_{\alpha}^{PP} & =0,\end{aligned}$$ where $$\mathcal{F}_{i}^{l}=\left(1+\tfrac{b^{2}}{4l+6}\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}^{2}\right),$$ is an operator encoding the finite size of the sphere and $\alpha$ may takes either of the values 1 or 2 corresponding to two equivalent directions parallel to wall. We have used results provided in Appendix \[appendix:Evaluation-of-gamma\], to write the expression for friction. The solution of this equation gives the lattice spacing $d$. For fixed particle-wall potential, increasing $v_{s}$ decreases the resting height $h$ and separation between *pairs, $d$,* as show in Fig. (\[fig:steady-state\]).

*Rotational dynamics:* In Fig. (\[fig:Dynamics-of-crystallization\]), we show the state diagram, obtained from simulation, which shows that the crystal is stable over a critical strengths of either bottom-heaviness or chirality. For an initially symmetric distribution, a crystal stabilized by external torque alone *does not rotate*, while the crystal *stabilized by chirality does rotate*. When the crystal is rotating at an angular velocity $\mathbf{\Omega}_{c}$ about its center of mass $\mathbf{R}_{c}$, the velocity the $i$-th colloid at position $\mathbf{R}_{i}$ can be then written as $\mathbf{\dot{R}}_{i}=\mathbf{\Omega}_{c}\times\mathbf{R}_{i}$. Force balance parallel to the wall is then $$\begin{aligned}
\bm{\gamma}_{ij}^{TT}\cdot\left[\mathbf{\Omega}_{c}\times(\mathbf{R}_{j}-\mathbf{R}_{0}^{c})\right]+\bm{\gamma}_{ij}^{TR}\cdot\mathbf{\Omega}_{j}=0.\label{eq:minimal-steady-state-1}\end{aligned}$$ The angular speed perpendicular to wall is $\Omega=\Omega_{i}^{\mathcal{A}}$. This implies that in absence of chiral self-rotation there is *no* rotation of the crystal. The angular velocity of the crystal can be obtained by power counting - $\bm{\gamma}_{ij}^{TT}$ scales as $r_{ij}^{-3}$ in direction parallel to wall while $\bm{\gamma}_{ij}^{TR}$ scales as $r_{ij}^{-4}$. The angular velocity of the crystal, then, scales as $\Omega_{c}\propto1/R_{c}^{2}$. In Fig. (\[fig:Dynamics-of-crystallization\]) we show that rotation period of a crystal scales inversely as number of particles $N$ in the crystal for an assembly of chiral particles, which is an excellent agreement with a recent experiment [@petroff2015fast].

![Steady states of active crystallization. Left panel has the plot of leading terms for the analytical solution of height $h$, shown in solid line, along with the full numerical result, shown as dotted curve. Right panel has similar set of plots for lattice spacing $d$. The leading order estimates are found to be in agreement with the numerical solution.\[fig:steady-state\]](SI-Figure2){width="47.00000%"}

Harmonic excitations\[sec:dispersion-relation\] 
================================================

In this section we study harmonic excitations $\mathbf{u}_{i}$ of the crystal about a stationary state $\mathbf{R}_{i}^{0}=(X_{i}^{0},\,Y_{i}^{0},\,h)$, such that $\dot{\mathbf{R}}_{i}^{0}=0$ and $\mathbf{\Omega}_{i}^{0}=0$ at this location. The force balance condition is then $\boldsymbol{\gamma}_{ij}^{TT}\cdot\mathbf{V}_{j}^{\mathcal{A}}\big|_{0}+\mathbf{F}_{i}^{P}\big|_{0}=0$. We, now, consider a small displacement about this state $\mathbf{R}_{i}=\mathbf{R}_{i}^{0}+\mathbf{u}_{i}$. Expanding the friction tensors about the stationarity point, we have, $$\boldsymbol{\gamma}_{ij}^{TT}=\boldsymbol{\gamma}_{ij}^{TT}\big|_{0}+(\mathbf{u}_{i}\cdot\mathbf{\boldsymbol{\nabla}}_{{\scriptscriptstyle \mathbf{R}_{i}}}-\mathbf{u}_{j}\cdot\mathbf{\boldsymbol{\nabla}}_{{\scriptscriptstyle \mathbf{R}_{j}}})\mathbf{}\boldsymbol{\gamma}_{ij}^{TT}\big|_{0}+\mathcal{O}(\mathbf{u}^{2}).$$ The force can be expanded in a similar way $\mathbf{F}_{i}^{P}=\mathbf{F}_{i}^{P}\big|_{0}-\mathbf{D}_{ij}\cdot\mathbf{u}_{j},$ where $\mathbf{D}_{ij}=-\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\mathbf{\boldsymbol{\nabla}}_{{\scriptscriptstyle \mathbf{R}_{i}}}U\big|_{0}$. Using the equations of motion and considering terms linear in the displacement, the equation becomes $$\begin{aligned}
-\boldsymbol{\gamma}_{ij}^{TT}\cdot\dot{\mathbf{u}}_{j} &  & +\left(\mathbf{\boldsymbol{\nabla}}_{{\scriptscriptstyle \mathbf{R}_{j}}}\boldsymbol{\gamma}_{ij}^{TT}\cdot\mathbf{V}^{\mathcal{A}}-\mathbf{D}_{ij}\right)\cdot\mathbf{u}_{j}=0.\end{aligned}$$ We seek a solution of the form $\mathbf{u}_{i}(t)=\mathbf{u}_{\mathbf{k}}(t)\,e^{i\mathbf{k}\cdot\mathbf{R}_{i}}.$ Using this, the force balance condition becomes, $$-\boldsymbol{\gamma}_{\mathbf{k}}^{TT}\cdot\dot{\mathbf{u}}_{k}+\left(i\mathbf{k}\,\boldsymbol{\gamma}_{\mathbf{k}}^{TT}\cdot\mathbf{V}^{\mathcal{A}}-\mathbf{D}_{\mathbf{k}}\right)\cdot\mathbf{u}_{\mathbf{k}}=0.\label{eq:stability-fourier}$$ Here $\mathbf{D}_{\mathbf{k}}$ is the Fourier transform of $\mathbf{D}_{ij}$ and $\boldsymbol{\gamma}_{\mathbf{k}}^{TT}$ is the Fourier transform of the friction tensor

\_ & =\_[i=1]{}\^[N]{}\_[i1]{}e\^[i(\_[i]{}-\_[1]{})]{},\
\_\^[TT]{} & =\_[i=1]{}\^[N]{}\_[i1]{}\^[TT]{}e\^[i(\_[i]{}-\_[1]{})]{}.

Here, $\mathbf{D}_{\mathbf{k}}$ is called the dynamical matrix [@born1954dynamical]. We now write $\boldsymbol{\gamma}_{i1}^{TT}$ in terms of its planar Fourier transform $$\boldsymbol{\gamma}_{i1}^{TT}=\int\hat{\boldsymbol{\gamma}}_{\mathbf{k}}^{TT}(\mathbf{k};\,h)e^{-i\mathbf{k}'\cdot(\mathbf{R}_{i}-\mathbf{R}_{1})}\,\frac{d^{2}k'}{(2\pi)^{2}},$$ to obtain an expression for $\boldsymbol{\gamma}_{\mathbf{k}}^{TT}$, $$\begin{aligned}
\boldsymbol{\gamma}_{\mathbf{k}}^{TT}(\mathbf{k};\,h) & = & \sum_{i}\int\hat{\boldsymbol{\gamma}}_{\mathbf{k}}^{TT}(\mathbf{k};\,h)e^{-i(\mathbf{k}'-k)\cdot(\mathbf{R}_{i}-\mathbf{R}_{1})}\,\frac{d^{2}k'}{(2\pi)^{2}},\nonumber \\
 & = & \frac{1}{A_{c}}\sum_{\lambda}\hat{\boldsymbol{\gamma}}_{\mathbf{k}}^{TT}(\mathbf{k}+\mathbf{q}_{\lambda};\,h).\end{aligned}$$ Here we have used the identity $$\sum_{i}e^{-i\mathbf{k}\cdot\mathbf{R}_{i}}=\frac{(2\pi)^{2}}{A_{c}}\sum\mathbf{\boldsymbol{\delta}}(\mathbf{k}-\mathbf{q}_{\lambda}),$$ where $A_{c}$ is area of the unit cell and $\mathbf{q}_{\lambda}$ are reciprocal lattice vectors. We now identify two parts of $\boldsymbol{\gamma}_{\mathbf{k}}^{TT},$ $$\begin{aligned}
{1}
 & \boldsymbol{\gamma}_{\mathbf{k}}^{TT}(\mathbf{k};\,h)=\hat{\boldsymbol{\gamma}}_{{\scriptscriptstyle \mathbf{k}_{0}}}^{TT}(\mathbf{k};\,h)+\sum_{\lambda'}\hat{\boldsymbol{\gamma}}_{{\scriptscriptstyle \mathbf{k}_{q}}}^{TT}(\mathbf{k}+\mathbf{q}_{\lambda};\,h)\end{aligned}$$ Here $\hat{\boldsymbol{\gamma}}_{{\scriptscriptstyle \mathbf{k}_{0}}}^{TT}(\mathbf{k};\,h)$ corresponds to the $\mathbf{q}_{\lambda}=0$ and terms at arbitrary non-zero $q$ are denoted by $\hat{\boldsymbol{\gamma}}_{{\scriptscriptstyle \mathbf{k}_{q}}}^{TT}(\mathbf{k};\,h)$. Their leading order forms can be written as $$\begin{aligned}
{1}
\hat{\boldsymbol{\gamma}}_{{\scriptscriptstyle \mathbf{k}_{0}}}^{TT}(\mathbf{k};\,h) & =\gamma^{T}\mathbf{I}+\frac{\gamma^{T}\gamma^{T}}{A_{c}}\mathcal{F}^{k}\,\hat{\mathbf{G}}^{W}(\mathbf{k};\,h),\\
\hat{\boldsymbol{\gamma}}_{{\scriptscriptstyle \mathbf{k}_{q}}}^{TT}(\mathbf{k};\,h) & =\frac{\gamma^{T}\gamma^{T}}{A_{c}}\mathcal{F}^{k}\,\sum_{\lambda'}\hat{\mathbf{G}}^{W}(\mathbf{k}+\mathbf{q}_{\lambda};\,h).\end{aligned}$$ Here $\mathcal{F}^{k}=1-b^{2}k^{2}/3$ and $\hat{\mathbf{G}}^{W}(\mathbf{k};\,h)$ is the two-dimensional Fourier transform of $\mathbf{G}^{W}$ (see Appendix \[appendix:FT-GW\]). The prime on the summation on the right indicates that $\lambda=0$ is excluded from the sum. We now turn to the calculation of the dynamical matrix,

\_ & =\_[i=1]{}\^[N]{}(U’+U”)\_[0]{}(1-e\^[i\_[i]{}]{}).

Here $U'=-12\epsilon\,\big[\left(r_{min}/d\right)^{12}-\left(r_{min}/d\right)^{6}\big]$ and $U''=12\epsilon\,\big[14\left(r_{min}/d\right)^{12}-8\left(r_{min}/d\right)^{6}\big]$. We evaluate the above in the nearest neighbor approximation in the direction parallel to the wall. The expression for $\mathbf{D}_{\mathbf{k}}$ and $\boldsymbol{\gamma}_{\mathbf{k}}^{TT}$ can be evaluated numerically by summing over the reciprocal lattice vectors. The sum is unconditionally and rapidly convergent as the Green’s function decays as $r_{ij}^{-3}$ in the direction parallel to the wall. The dispersion is obtained numerically from Eq. (\[eq:stability-fourier\]) and is shown in Fig. (\[fig:Normal-modes\]).

*Long-wavelength approximation:* Analytical expression for the normal modes can be obtained in the $k\rightarrow0$ limit. Keeping terms of the $\mathcal{O}(k^{2})$, Eq. (\[eq:stability-fourier\]) becomes $$\begin{aligned}
\,\left(\begin{array}{c}
\dot{u}_{k_{1}}\\
\dot{u}_{k_{2}}
\end{array}\right) & = & -\frac{h\gamma_{\perp}^{T}v_{s}\,k^{2}}{\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}}\boldsymbol{f}(\theta)\left(\begin{array}{c}
u_{k_{1}}\\
u_{k_{2}}
\end{array}\right),\end{aligned}$$ $$\boldsymbol{f}(\theta)=\left(\begin{array}{cc}
c_{1}\cos^{2}\theta+\mathcal{C}_{2}\sin^{2}\theta & c_{3}\sin\theta\cos\theta\\
c_{3}\sin\theta\cos\theta & c_{4}\sin^{2}\theta+c_{5}\cos^{2}\theta
\end{array}\right).$$ Here $k_{1}=k\cos\theta$, $k_{2}=k\sin\theta$ and $c_{i}$ are positive constants that can be determined in terms of the parameters of the steric potential and the friction tensors: $c_{1}=\gamma_{\perp}^{T}h/2\eta A_{c}+\left(\tfrac{3}{2}U'+\tfrac{9}{8}U''+\right)/h\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}v_{s},$ $c_{2}=\left(\tfrac{3}{2}U'+\tfrac{3}{8}U''+\right)/h\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}v_{s},\,$ $c_{3}=\gamma_{\perp}^{T}h/2\eta A_{c}+\tfrac{3}{4}U''/h\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}v_{s}$, $c_{4}=\gamma_{\perp}^{T}h/2\eta A_{c}+3c_{5}$ and $c_{5}=\left(\tfrac{1}{2}U'+\tfrac{3}{8}U''+\right)/h\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}v_{s}$. We can now diagonalize this matrix equation to obtain the relaxation of the overdamped modes after Fourier transforming in time. The eigenvalues of the resulting equations give the dispersion relation **$$\omega_{\pm}=-i\tfrac{\gamma_{\perp}^{T}hv_{s}}{2\gamma_{{{\mkern3mu\vphantom{\perp}\vrule depth 0pt\mkern2.5mu\vrule depth 0pt\mkern3mu}}}^{T}}f_{\pm}(\theta)\,k^{2},\label{eq:small-k-1}$$ in terms of an angular factor,** $f_{\pm}(\theta)=c_{15}\cos^{2}\theta+c_{24}\sin^{2}\theta\pm[(c_{15}\cos^{2}\theta+c_{24}\sin^{2}\theta)^{2}-4(c_{1}c_{5}\cos^{4}\theta+c_{2}c_{4}\sin^{4}\theta-c_{3}^{2}\cos^{2}\theta\sin^{2}\theta+(c_{1}c_{4}+c_{2}c_{5})\cos^{2}\theta\sin^{2}\theta)]{}^{1/2}$, with $c_{15}=c_{1}+c_{5}$ and $c_{24}=c_{2}+c_{4}$. The comparison of the long wavelength solution with the full numerical solution has been plotted in Fig. (\[fig:Normal-modes\]). For **$k\lesssim0.1k_{0}$,** the approximate solution shows excellent agreement.

Numerical method
================

In this section, we outline the method used to simulate the dynamics of active colloidal particles near a no-slip wall. We invert Eq. (\[eq:main-traction-l\]) to obtain rigid body motions in terms of the known slip modes, body forces and torques [@singh2016traction]. This gives the “mobility” formulation, $$\begin{aligned}
{1}
\mathsf{\mathbf{V}}_{i} & =\boldsymbol{\mu}_{ij}^{TT}\cdot\mathbf{F}_{j}^{P}+\boldsymbol{\mu}_{ij}^{TR}\cdot\mathbf{T}_{j}^{P}+\sum_{l\sigma=2s}^{\infty}\boldsymbol{\pi}_{ij}^{(T,l\sigma)}\cdot\mathbf{\mathsf{\mathbf{V}}}_{j}^{(l\sigma)}+\mathsf{\mathbf{V}}_{i}^{\mathcal{A}},\\
\mathsf{\mathbf{\Omega}}_{i} & =\boldsymbol{\mu}_{ij}^{RT}\cdot\mathbf{F}_{j}^{P}+\boldsymbol{\mu}_{ij}^{RR}\cdot\mathbf{T}_{j}^{P}+\sum_{l\sigma=2s}^{\infty}\boldsymbol{\pi}_{ij}^{(R,\,l\sigma)}\cdot\mathbf{\mathsf{\mathbf{V}}}_{j}^{(l\sigma)}+\mathsf{\mathbf{\Omega}}_{i}^{\mathcal{A}}.\end{aligned}$$ The mobility matrices $\boldsymbol{\mu}_{ij}^{\alpha\beta}$, with ($\alpha,\beta=T,R$), are inverses of the friction matrices $\boldsymbol{\gamma}_{ij}^{\alpha\beta}$ [@kim2005]. The propulsion tensors $\boldsymbol{\pi}_{ij}^{(\alpha,\,l\sigma)}$, first introduced in [@singh2014many], relate the rigid body motion to modes of the active velocity. They are related to the slip friction tensors by [@singh2016traction], $$\begin{aligned}
{1}
-\bm{\pi}_{ij}^{(\text{T},\,l\sigma)} & =\boldsymbol{\mu}_{ik}^{TT}\cdot\boldsymbol{\gamma}_{kj}^{(T,\,l\sigma)}+\boldsymbol{\mu}_{ik}^{TR}\cdot\boldsymbol{\gamma}_{kj}^{(R,\,l\sigma)},\\
-\bm{\pi}_{ij}^{(R,\,l\sigma)} & =\boldsymbol{\mu}_{ik}^{RT}\cdot\boldsymbol{\gamma}_{kj}^{(T,\,l\sigma)}+\boldsymbol{\mu}_{ik}^{RR}\cdot\boldsymbol{\gamma}_{kj}^{(R,\,l\sigma)}.\end{aligned}$$ We retain modes corresponding to $l\sigma=1s,\,2s,2a,\,3a,\,3t$ and $4a$ in the active slip. The role of these individual modes is summarized in Table (\[tab:slip-mode-verbose\]). The mobilities are calculated using the PyStokes [@pystokes] library. The initial distribution of particles is chosen to be the random packing of hard-spheres [@skoge2006packing]. We use an adaptive time step integrator using the backward differentiation formula (BDF) to integrate these equations of motion [@langtangen2012tutorial]. In Table (\[tab:Table-of-simulation-parameters\]) of Appendix, we present the parameters used to generate the figures.

**

[|&gt;b[1.4cm]{}|&gt;b[2cm]{}|&gt;b[2cm]{}|&gt;b[2cm]{}|]{} Slip mode & Positional clustering & Orientational stability & Cluster rotation [\
]{} $\mathbf{V}^{a}$ & Yes & No & No[\
]{} $\mathbf{\Omega}^{a}$ & No & Yes & No[\
]{} $\mathbf{V}^{(2s)}$ & Yes & No & No[\
]{} $\mathbf{V}^{(3t)}$ & Yes & No & No[\
]{} $\mathbf{V}^{(3a)}$ & No & Yes & Yes[\
]{} $\mathbf{V}^{(4a)}$ & No & Yes & Yes[\
]{}

[|&gt;b[2cm]{}|c|&gt;p[2cm]{}|&gt;p[2cm]{}|&gt;p[2cm]{}|c|c|]{} Figure & \# of colloids & $v_{s}$ & $V_{0}^{(3a)}/v_{s}$ & $T_{0}/\epsilon$ & Wall WCA & Inter-particle WCA[\
]{} 1 (1-f) & 1024 & 0.1 & 100 & 0 & $r_{min}=3.4b$,$\,$$\epsilon=0.083$ & $r_{min}=5b$,$\,$$\epsilon_{p}=0.004$[\
]{} 1 (g) & - & 0.01 & 100 & 0 & $r_{min}=3.4b$,$\,$$\epsilon=0.083$ & $r_{min}=5b$,$\,$$\epsilon_{p}=0.004$[\
]{} 1 (h) & 256 & 0.01 & - & - & $r_{min}=3.4b$,$\,$$\epsilon=0.083$ & $r_{min}=5b$,$\,$$\epsilon_{p}=0.004$[\
]{} 2 & 1 and 2 & 0.1 & 100 & - & - & -[\
]{} 3 & - & 0.1 & - & - & - & $r_{min}=5b$,$\,$$\epsilon_{p}=0.004$[\
]{}

Expression for boundary integrals and matrix elements\[sec:Expression-for-boundary\]
====================================================================================

The boundary integrals in fluid flow, Eq. (\[eq:fluid-flow\]), can be solved exactly. The resulting solution is given in terms of the Green’s function, Eq. (\[eq:wall-G\]), and its derivatives [@singh2014many], $$\begin{aligned}
{1}
\boldsymbol{G}_{j}^{(l)}(\mathbf{r},\mathbf{R}_{j}) & =\frac{2l-1}{4\pi b^{2}}\int\mathbf{G}^{W}(\mathbf{r},\mathbf{R}_{j}+\bm{\rho}_{j})\mathbf{Y}^{(l-1)}(\hat{\bm{\rho}}_{j})\,d\mathrm{S}_{i}=b^{l}\mathcal{F}^{l-1}\mathbf{\bm{\nabla}}_{{\scriptscriptstyle \mathbf{R}_{j}}}^{(l-1)}\mathbf{G}^{W}(\mathbf{r},\mathbf{R}_{j}),\\
\boldsymbol{K}_{j}^{(l)}(\mathbf{r},\mathbf{R}_{j}) & =\frac{1}{(l-1)!(2l-3)!!}\int\mathbf{K}^{W}(\mathbf{r},\mathbf{R}_{j}+\bm{\rho}_{j})\cdot\mathbf{n}\mathbf{Y}^{(l-1)}(\hat{\bm{\rho}}_{j})\,d\mathrm{S}_{i}=\frac{4\pi b^{l+1}}{(l-2)!(2l-1)!!}\mathcal{F}^{l-1}\mathbf{\bm{\nabla}}_{{\scriptscriptstyle \mathbf{R}_{j}}}^{(l-2)}\mathbf{K}^{W}(\mathbf{r},\mathbf{R}_{j}).\end{aligned}$$ The integrals appearing in the linear system of the equations, Eq. (\[eq:linear-system\]), are, $$\begin{aligned}
{1}
\boldsymbol{G}_{ij}^{(l,\,l')}(\mathbf{R}_{i},\mathbf{R}_{j}) & =\frac{(2l-1)(2l'-1)}{(4\pi b^{2})^{2}}\int\mathbf{Y}^{(l-1)}(\hat{\bm{\rho}}_{i})\mathbf{G}^{W}(\mathbf{R}_{i}+\bm{\rho}_{i},\,\mathbf{R}_{j}+\bm{\rho}_{j})\mathbf{Y}^{(l'-1)}(\hat{\bm{\rho}}_{j})\,d\mathrm{S}_{i}\,d\mathrm{S}_{j}\\
\boldsymbol{K}_{ij}^{(l,\,l')}(\mathbf{R}_{i},\mathbf{R}_{j}) & =\frac{2l-1}{4\pi b^{2}\,(l-1)!(2l-3)!!}\int\mathbf{Y}^{(l-1)}(\hat{\bm{\rho}}_{i})\mathbf{K}^{W}(\mathbf{R}_{i}+\bm{\rho}_{i},\,\mathbf{R}_{j}+\bm{\rho}_{j})\cdot\mathbf{n}\mathbf{Y}^{(l'-1)}(\hat{\bm{\rho}}_{j})\,d\mathrm{S}_{i}\,d\mathrm{S}_{j}.\end{aligned}$$ These integrals are solved exactly to give matrix elements in terms of the Green’s function and its derivatives [@singh2014many], $$\begin{aligned}
{1}
\boldsymbol{G}_{ij}^{(l,\,l')}(\mathbf{R}_{i},\mathbf{R}_{j}) & =\begin{cases}
{\displaystyle \mathcal{G}_{ii}^{(l,\,l')}+b^{l+l'-2}\mathcal{F}_{i}^{l-1}\mathcal{F}_{j}^{l'-1}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}^{(l-1)}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}^{(l'-1)}\mathbf{G}^{*}(\mathbf{R}_{i},\mathbf{R}_{j});} & \qquad\qquad\qquad\qquad{\displaystyle \qquad\quad}{\displaystyle j=i,}\\
{\displaystyle b^{l+l'-2}\mathcal{F}_{i}^{l-1}\mathcal{F}_{j}^{l'-1}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}^{(l-1)}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}^{(l'-1)}\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j});} & \qquad\qquad\qquad\qquad{\displaystyle \qquad\quad{\displaystyle j\neq i},}
\end{cases}\\
\boldsymbol{K}_{ij}^{(l,\,l')}(\mathbf{R}_{i},\mathbf{R}_{j}) & =\begin{cases}
-{\displaystyle \tfrac{1}{2}\delta_{ll'}\bm{\Delta}^{(l-1)}+\frac{4\pi b^{(l+l'-1)}}{(l'-2)!(2l'-1)!!}\mathcal{F}_{i}^{l-1}\mathcal{F}_{j}^{l'-1}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}^{(l-1)}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}^{(l'-2)}\mathbf{K}^{*}(\mathbf{R}_{i},\mathbf{R}_{j});}\quad\qquad\qquad & {\displaystyle j=i,}\\
{\displaystyle \frac{4\pi b^{(l+l'-1)}}{(l'-2)!(2l'-1)!!}\mathcal{F}_{i}^{l-1}\mathcal{F}_{j}^{l'-1}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}^{(l-1)}\bm{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}^{(l'-2)}\mathbf{K}^{W}(\mathbf{R}_{i},\mathbf{R}_{j});}\qquad\qquad\qquad\qquad\qquad & {\displaystyle j\neq i,}
\end{cases}\end{aligned}$$ $$\mathcal{G}_{ii}^{(l,\,l')}=\delta_{ll'}\frac{2l-1}{2\pi b}\int\mathbf{Y}^{(l-1)}(\hat{\bm{\rho}})\left(\mathbf{\bm{I}}-\hat{\bm{\rho}}\hat{\bm{\rho}}\right)\mathbf{Y}^{(l-1)}(\hat{\bm{\rho}})\,d\Omega.$$

First order off-diagonal approximation for friction tensors\[appendix:Evaluation-of-gamma\]
===========================================================================================

The expressions for the friction tensor can be calculated from the solution of the linear system, provided above, using the Jacobi method [@singh2016traction]. The first order approximation to friction tensors used in this work are provide below, $$\begin{aligned}
\Big(\boldsymbol{\gamma}_{ij}^{(TT)} & \Big)^{[1]}= & \gamma^{T}\gamma^{T}\,\mathcal{F}_{i}^{0}\mathcal{F}_{j}^{0}\,\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j}),\quad\Big(\boldsymbol{\gamma}_{ij}^{(RT)}\Big)^{[1]}=\tfrac{1}{2}\gamma^{T}\gamma^{R}\,\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}\times\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j}),\end{aligned}$$ $$\Big(\boldsymbol{\gamma}_{ij}^{(TR)}\Big)^{[1]}=\tfrac{1}{2}\gamma^{T}\gamma^{R}\,\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\times\left(\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j})\right),\quad\Big(\boldsymbol{\gamma}_{ij}^{(RR)}\Big)^{[1]}=\tfrac{1}{4}\gamma^{R}\gamma^{R}\,\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}\times\left(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\times\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j})\right),$$ $$\Big(\boldsymbol{\gamma}_{ij}^{(T,\,2s)}\Big)^{[1]}=\frac{28\pi\eta b^{2}}{3}\gamma^{T}\,\mathcal{F}_{i}^{0}\mathcal{F}_{j}^{1}\,\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j}),\quad\Big(\boldsymbol{\gamma}_{ij}^{(R,\,2s)}\Big)^{[1]}=\frac{28\pi\eta b^{2}}{6}\gamma^{R}\,\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}\times\left(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j})\right),$$ $$\Big(\boldsymbol{\gamma}_{ij}^{(T,\,3a)}\Big)^{[1]}=\frac{13\pi\eta b^{3}}{9}\gamma^{T}\,\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\times\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j})),\quad\Big(\boldsymbol{\gamma}_{ij}^{(R,\,3a)}\Big)^{[1]}=\frac{13\pi\eta b^{3}}{18}\gamma^{R}\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}\times\left(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\times\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j}))\right),$$ $$\Big(\boldsymbol{\gamma}_{ij}^{(T,\,3t)}\Big)^{[1]}=-\frac{4\pi\eta b^{3}}{5}\gamma^{T}\,\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}^{2}\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j}),\quad\Big(\boldsymbol{\gamma}_{ij}^{(T,\,4a)}\Big)^{[1]}=-\frac{121\pi\eta b^{4}}{10}\gamma^{T}\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\times\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j})),$$ $$\Big(\boldsymbol{\gamma}_{ij}^{(R,\,4a)}\Big)^{[1]}=-\frac{121\pi\eta b^{4}}{20}\gamma^{R}\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{i}}}\times\left(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}(\boldsymbol{\nabla}_{{\scriptscriptstyle \mathbf{R}_{j}}}\times\mathbf{G}^{W}(\mathbf{R}_{i},\mathbf{R}_{j}))\right),\quad\Big(\boldsymbol{\gamma}_{ij}^{(R,\,3t)}\Big)^{[1]}=0.$$

Fourier transform of the Lorentz-Blake Green’s function\[appendix:FT-GW\]
=========================================================================

In this section, we derive the Fourier transform of the Green’s function for a fluid flow bounded by a plane infinite wall. Blake [@blake1971c] has derived the Green function of the Stokes equation which satisfies no-slip condition on the wall, $$\begin{aligned}
{1}
G_{\alpha\beta}^{\text{w}}(\mathbf{R}_{i},\,\mathbf{R}_{j}) & =G_{\alpha\beta}(\mathbf{R}_{i},\,\mathbf{R}_{j})+G_{\alpha\beta}^{*}(\mathbf{R}_{i},\,\mathbf{R}_{j}^{*})=G_{\alpha\beta}(\mathbf{R}_{i},\,\mathbf{R}_{j})-G_{\alpha\beta}(\mathbf{R}_{i},\,\mathbf{R}_{j})+G'_{\alpha\beta}(\mathbf{R}_{i},\,\mathbf{R}_{j})\label{eq:wall-G}\\
 & =G_{\alpha\beta}(\mathbf{r}_{ij}^{*})-G_{\alpha\beta}(\mathbf{r}_{ij}^{*})-2h\nabla_{{\scriptscriptstyle \mathbf{r}_{ij}^{*}}}G_{\alpha3}(\mathbf{r}_{ij}^{*})\mathcal{M}_{\beta\gamma}+h^{2}\nabla_{{\scriptscriptstyle \mathbf{r}_{ij}^{*}}}^{2}G_{\alpha\gamma}(\mathbf{r}_{ij}^{*})\mathcal{M}_{\beta\gamma}.\nonumber \end{aligned}$$ Here $\mathbf{r}_{ij}=\mathbf{\mathbf{R}}_{i}-\mathbf{\mathbf{R}}_{j}$, $\mathbf{r}_{ij}^{*}=\mathbf{\mathbf{R}}_{i}-\mathbf{\mathbf{R}}_{j}^{*}$ and $\boldsymbol{\mathcal{M}}=\mathbf{I}-2\mathbf{\hat{z}}\mathbf{\hat{z}}$. $\mathbf{G}$ is the Green’s function in the unbounded fluid flow, $$\mathbf{G}(\mathbf{R}_{i},\,\mathbf{R}_{j})=\frac{1}{8\pi\eta}\left(\frac{\mathbf{I}}{r_{ij}}+\frac{\mathbf{r}_{ij}\mathbf{r}_{ij}}{r_{ij}^{3}}\right).$$ We define the Fourier transform in the plane of the wall as, $$\hat{\varphi}(k_{1},\,k_{2},\,r_{3})=\mathbb{\mathbb{F}}\left[\varphi\right]=\frac{1}{(2\pi)^{2}}\int\varphi(r_{1},\,r_{2},\,r_{3})e^{i(k_{1}r_{1}+k_{2}r_{2})}\,dr_{1}dr_{2}.$$ The Fourier transform of $G'_{ij}$, last term of Eq. (\[eq:wall-G\]), is then [@blake1971c], $$\hat{G}'_{\alpha\beta}(\mathbf{k};\,h)=\frac{h}{2\eta k}\left[ik_{\alpha_{1}}(\delta_{\alpha3}\delta_{j\alpha_{1}}+\delta_{\beta3}\delta_{\alpha\alpha_{1}})+h\left(ikk_{\alpha_{1}}\{\delta_{\alpha3}\delta_{j\alpha_{1}}-\delta_{\beta3}\delta_{\alpha\alpha_{1}}\}-k_{\alpha_{1}}k_{\alpha_{2}}\delta_{\alpha\alpha_{1}}\delta_{\beta\alpha_{2}}-\delta_{\alpha3}\delta_{\beta3}k^{2}\right)\right]e^{-2kh},$$ where $\alpha_{1}$ and $\alpha_{2}$ only take values 1 or 2 corresponding to directions parallel to wall. The rest of terms in Eq. (\[eq:wall-G\]), can be transformed using the relation $\mathbb{\mathbb{F}}\left[\frac{1}{r}\right]=\frac{2\pi\,e^{-kz}}{k}$. The two-dimensional Fourier transform of the wall Green’s function for a source at height $h$ from the wall is then, with $\mathcal{E}=1-e^{-2kh}$, $$\begin{aligned}
\hat{\mathbf{G}}^{W}(\mathbf{k};\,h) & =\frac{1}{4\eta k^{3}}\left(\begin{array}{ccc}
\mathcal{E}k_{2}^{2}+2hkk_{1}^{2}e^{-2kh}\quad & -\mathcal{E}k_{1}k_{2}-2hkk_{1}k_{2}e^{-2kh}\quad & -i2hk^{2}k_{1}e^{-2kh}\\
\\
-\mathcal{E}k_{1}k_{2}-2hkk_{1}k_{2}e^{-2kh}\quad & \mathcal{E}k_{1}^{2}+2hkk_{2}^{2}e^{-2kh}\quad & -i2hk^{2}k_{2}e^{-2kh}\\
\\
-i2hk^{2}k_{1}e^{-2kh}\quad & -i2hk^{2}k_{1}e^{-2kh}\quad & \mathcal{E}k^{2}-2hk^{3}e^{-2kh}
\end{array}\right)+\hat{\mathbf{G}}'(\mathbf{k};\,h).\label{eq:Gw-FT}\end{aligned}$$

![image](SI-Figure3){width="82.00000%"}
