A first-principle mechanism for particulate aggregation and self-assembly in stratified fluids

An extremely broad and important class of phenomena in nature involves the settling and aggregation of matter under gravitation in fluid systems. Here, we observe and model mathematically an unexpected fundamental mechanism by which particles suspended within stratification may self-assemble and form large aggregates without adhesion. This phenomenon arises through a complex interplay involving solute diffusion, impermeable boundaries, and aggregate geometry, which produces toroidal flows. We show that these flows yield attractive horizontal forces between particles at the same heights. We observe that many particles demonstrate a collective motion revealing a system which appears to solve jigsaw-like puzzles on its way to organizing into a large-scale disc-like shape, with the effective force increasing as the collective disc radius grows. Control experiments isolate the individual dynamics, which are quantitatively predicted by simulations. Numerical force calculations with two spheres are used to build many-body simulations which capture observed features of self-assembly.


Equations of motion
Here we provide details regarding how we employ COMSOL to solve the incompressible Stokes equations coupled to a steady advection-diffusion equation for concentration. Making the steady assumption, we solve the following equations using the COMSOL finite element package, currently in an axisymmetric domain of height 2H and radius R: With boundary conditions on velocity field: and boundary conditions on density field: ∇ρ ·n| r=a = 0 (7) ∇ρ ·r| r=R = 0 (8) ρ| z=±H = ρ o ∓ σH (9) As described in the main text, these are the resulting equations after assuming that the Reynolds number is sufficiently small and that sufficient time has passed so that the transient density

Supplementary Note 2 Domain and meshing
We chose H and R large enough and mesh element size small enough for each simulation such that the peak speeds (both global maximum speed and maximum horizontal velocity in the equatorial plane) have converged to 3 significant digits. There is a refined region near the surface of the sphere (in the boundary layer) that transitions smoothly with a specified element growth rate to the far-field, where the largest element size is constrained (see Supplementary Figure 1

3D unsteady simulations
For the 3D simulations, a solution of the coupled Stokes flow and unsteady advectiondiffusion equations for density is sought by running to sufficiently long time for transients to decay. Specifically, we numerically solve the coupled equations For the case shown in Supplementary Figure 2, a rectangular domain size of X ×Y ×Z is used, with mirror symmetry conditions about the y − z and z − x planes. Parameters of the domain and mesh are provided in Supplementary Table 2., and the attractive force shows good convergence for both the domain size and mesh resolution used, as evidenced in Supplementary Figure 2(b). The forces are computed by numerically integrating the projected stress tensor which includes both viscous and pressure contributions. We monitor the evolution of this force in time to assess the timescale for decay of transients as a function of the separation distance. This is documented in Supplementary Figure 3 which shows the evolution of the force normalized by its long time limit for 4 different spacings. For particles which are close, the transient decays extremely rapidly, on the order of one minute. Clearly, the diffusion induced flow from one particle will take longer to affect the other particle if their separation distance is large, but the velocities induced when the particles are well separated is very weak, so it may well be the case that the steady assumption is still valid even at large separations. Even for separation distances of 10 diameters, the transient decays in less than two minutes for experiments which take 8 hours to collapse.
The force data shown has a reasonable fit to an expression of the form where F 0 = 2.0865 × 10 −6 dynes, which is then used for our modified Stokesian dynamics simulations.

Supplementary Note 4 Modified Stokesian dynamics
We have modified a Stokesian dynamics software package available on GitHub developed by Townsend (4) which implements the formalism pioneered by (5). This solver allows for a user defined external force field. In our implementation, we use the above force field expression (13) with a correction that accounts for the relative velocity of the moving particle with respect to the diffusion induced flow. This correction is a non-trivial modification of the algorithm, as it requires introducing semi-implicit time-stepping in the computation of velocities to avoid numerical instabilities. The algorithm is applied to a set of two small experiments as a benchmark. Supplementary Figure 4 depicts the separation distance as a function of time constructed with modified Stokesian dynamics compared to two experiments with spheres of radius 0.045cm in two different linear background stratifications (σ = 0.03, and σ = 0.025g/cm 4 ), with viscosity and diffusivity as in Supplementary Figure 3. With the strong agreement observed in this benchmark, the algorithm may be applied to hundreds of spheres, as we have demonstrated in the main body. We note that improvements can be made by running the two sphere calculation in a moving frame to represent more faithfully the coupled nature of the advection-diffusion dynamics. In future work we plan to report on these more refined approaches in detail.
S t o k e s i a n d y n a mi c s t r i a l 1 S t o k e s i a n d y n a mi c s t r i a l 2 E x p e r i me n t t r i a l 1 E x p e r i me n t t r i a l 2 Supplementary Figure 4: Stokesian dynamics (SD) compared with two-sphere experiments. SD modified by adding the two-body force computed via COMSOL given in force field expression (13). Here, a = 0.045 cm, with two linear background density fields, σ = 0.03g/cm 4 and σ = 0.025g/cm 4 in the first and second trial respectively, other experimental parameters matching Supplementary Figure 3.

Low Peclet number limit
We next discuss a semi-analytical approach (valid at low Peclet number) which exactly computes the density perturbation from the background for a single body and approximately for multiple bodies using the method of images. The flow is then computed using convolution with a special Green's function augmented again by the method of images. In the low Peclet number limit for the non-dimensionalization presented in the main text, the advection diffusion equation for the concentration field, ρ, reduces to a simple heat equation, which once transients decay (see prior timescale discussion), further reduces to Laplace's equation, with flow given by solution to the Stokes equations: with boundary conditions being impermeable for the concentration, and no-slip for the velocity.
For the case of the density perturbation itself, with a linear background density profile, explicit formulae are available: Interestingly, in this limit, the concentration field problem is actually equivalent to computing the induced potential field for a body held in a vertical constant flow, under the assumption of potential flow. For a single sphere (1): where a 1 is the radius of the sphere. In the presence of another sphere of radius a 2 , a correction to the density anomaly is given by (2) where c is the distance between centers (the companion sphere is obtained through summing identical translated formulae). The density anomaly for the oblate spheroid is given by (streamfunction in (1), potential through our own calculation) where a is the semi-major axis, b is the semi-minor axis, and (ξ, η, θ) are the oblate spheroidal coordinates that define the cartesian mapping x = √ a 2 − b 2 sin(η) cosh(ξ) cos(θ), y = √ a 2 − b 2 sin(η) cosh(ξ) sin(θ), and z = √ a 2 − b 2 cos(η) sinh(ξ). In Supplementary Figure 5, we present the density anomalies for these cases: on the left is a single sphere, in the middle is an oblate spheroid, and in the right is the case with two spheres. A horizontal wall (lid) may be implemented exactly by adding the solution arising from an imaging sphere positioned at twice the distance between the real sphere's center and the plane. A second horizontal wall (the bottom) may be approximated by repeating this procedure.

Low Peclet number limit
The fluid flow problem is then solved using an exact Green's function derived by Oseen (3) for a delta function of momentum inserted into a fluid in the presence of a fixed, no-slip sphere. With this explicit Green's function (given below), we may express the velocity field as a convolution with the exact concentration field: where x is the point at which the flow will be evaluated, y is the integration variable which represents the location of the point force, Ω is all of R 3 exterior to the sphere, and We note that the slow decay of the kernel generally will lead to divergences in free space. Imaging the Green's function G i with respect to two horizontal planes allows for the convolution domain to be limited to an infinitely wide channel of finite height, which is sufficient to restore convergence because of the far field decay rate of the density anomaly.
To represent a vertical no-flux wall at horizontal distance h from the center, which is equivalent to the two-sphere problem with a second sphere located at a position 2h, we introduce a second image of this type where the domain of integration is exterior to the sphere and truncated at the wall, x vert = (x 1 − 2h, x 2 , x 3 ), and y vert = (−y 1 , y 2 , y 3 ). To represent a horizontal no-flux wall at distance h, we introduce a second image as where the domain of integration is exterior to the sphere and truncated at the wall, x horiz = (x 1 , x 2 , x 3 − 2h), and y horiz = (y 1 , y 2 , −y 3 ). For our computation, we include an image representing a second sphere shifted in the horizontal direction, and then image both of these above and below to model a top and bottom wall, for a total of six images. The subsequent calculations are performed using a midpoint Riemann sum. Supplementary Figure 6 displays the flows induced by two fixed spheres in the low Peclet limit for four different separation distances. The trend is that, for well separated spheres, interactive effects are minimal, as the flow structures resemble the case for an individual from main article Figure 3. As the spheres are moved closer together, the sphere-sphere interaction modifies the flow structures, with the flow immediately between the spheres weakening as they approach each other (note the method of images induces some mild errors near the boundaries, but is faithful away from boundaries). Ultimately, as the spheres touch, the flow qualitatively resembles that of a single prolate spheroid and produces the strongest flow of the four spacings.