Dynamics of Acoustically Bound Particles

It is well known that acoustic fields can produce forces on single particles, however they can also induce inter-particle forces due to multiple scattering events. This multi-particle force -- here referred to as acoustic binding -- is comparable to other acoustic forces when the particles are of order wavelength in diameter. In principle, this force could be used as a tunable method for directing the assembly of particles of mm-scales, but has not been extensively explored in previous work. Here, we use a novel numerical method to compute binding interactions between strongly scattering bodies and find that they can produce stable clusters of particles with approximately wavelength separation. Moreover, we also observe that -- depending on the level of damping -- these structures can produce driven linear, rotational, or vibrational motion. These effects are a result of the non-conservative and non-pairwise nature of the acoustic binding force, and represent novel contactless manipulation and transport methods with potential application to metamaterial synthesis and drug delivery.


I. INTRODUCTION
While matter can alter the path of an external field through scattering and absorption, the field can also induce forces on matter.The resulting forces are a topic of great interest within active and condensed matter communities, as well as the field of material synthesis [1][2][3].For micron-scale objects, optical trapping has found numerous applications in physics, biology, and medicine [4][5][6][7][8][9][10].Optical forces can be placed on single objects, however second-order effects, known as optical binding can introduce inter-particle forces that has the potential to self-organize structures with wavelength-scale features [11,12].These binding forces arise when scattered and incident fields interfere to produce a new field gradient between two or more scattering bodies.
Could optical forces be used to manipulate larger particles?The field carried momentum goes as the intensity over the velocity of the wave p = I/c, which restricts the size of bodies that can be manipulated by an external field.Current optical trapping methods can levitate particles on the order of 10 µm in size (10 ng in mass) [13], limited by the maximum practical power density.
To manipulate larger particles we can employ a wave with slower speed, such as sound, which produces a force which is roughly 10 6 times higher for the same level of input power.Moreover, sound waves have a longer wavelength, so the analogous 'acoustic binding' effects will also self-organize structure on larger scales.
Ultrasound has previously been demonstrated to produce forces on single particles or even between many par-ticles [14][15][16][17][18][19][20][21][22][23].Acoustic levitation and trapping have been widely used to manipulate single particles on millimeter scales [24][25][26][27][28][29][30][31].(Commonly used 40 kHz sound waves have a wavelength of about 8.6 mm.) Far less known is the corresponding multiple particle force; in previous literature this has gone by several names -including acoustically induced mutual force [32] and acoustic interaction force [33], but here we will refer to this force as acoustic binding.The acoustic binding force arises from interference between the scattered field and the incident, resulting a long range oscillatory force.As discussed below, this force is distinct from the secondary Bjerknes force, which is a short range interparticle force acting on deformable bodies like bubbles in acoustic fields [34][35][36].
Unlike its optical counterpart, the acoustic binding force has only been studied in a small number of specific cases [19,32,33,37,38], despite the fact it represents a potentially powerful tool for self-organizing structures on mm-scales.Several different methods have been used to model acoustic forces, including multipole expansion of the acoustic scattering from weak scatterers [22,37,39], finite element methods [38,40], and a simplification of the problem to pairwise interactions [41].Here, we explore the acoustic binding forces in more detail using a novel numerical method, known as the method of fundamental solutions (MFS) [42].This method is both fast and accurate, allowing for one to conduct molecular dynamics simulations of small particle clusters.
For particles in the Mie scattering regime (ka > ∼ 1, where k is the wavenumber of the ultrasound field and a is the particle radius) we find that binding forces become comparable to trapping forces (fig.1).Consequently, Relative strengths of various acoustic forces as a function of size of the scattering body, computed using MFS numerical model.For ka < ∼ 1 the gradient force (blue) dominates; for ka > ∼ 1 the secondary forces (red, green, orange) are of comparable magnitude.The binding force is shown for two particles at fixed separation distances of 1.25λ (green) and 7.25λ (red) while the gradient force is for a single particle located between a pressure node and antinode.For reference the particle volume (dashed blue) and square of the volume (dashed orange) are included.b,c) Two particle acoustic binding force (black), 2nd order pressure field ( p2 ) (red), and gradient of the 2nd order pressure field integrated over particle surface (dashed gray) vs particle separation for (b) ka = 1 and (c) ka = 10.
these forces can be used to self-assemble structure on wavelength scales, and it should be possible to alter this structure by changing the properties of the particles or acoustic field.Furthermore, we observe active clusters of particles which drift, oscillate, and/or rotate depending on the cluster configuration (for example, see Supplemen-tary Movies 1, 2, and 3).As a result, acoustically bound clusters are not limited to passive self-assembly, opening up the possibility of driven particle assembly analogous to recently developed active matter systems [43,44].Previous examples of acoustically driven systems have relied on body asymmetry [45], complex sound fields [26,46], or acoustohydrodynamic interactions [47][48][49] as the driving mechanism.Conversely the behavior demonstrated here is not a result of any complex particle or beam characteristics, but rather the result of the nonlinearity of acoustic binding interactions between the constituents of this simple system.

A. Acoustic Forces
Acoustic fields can induce many types of forces on individual particles or between them.Broadly speaking, these can be separated into scattering forces [14,22,32], those due to streaming induced by acoustic waves [48,49], and those due to the deformation of the particles (e.g.Bjerknes forces) [15,35,36].As detailed below, for wavelength-sized (Mie) solid particles in a gaseous medium, the scattering forces will dominate over these other types.The scattering forces can be further divided into gradient forces, radiation pressure forces, and binding forces; each of these is directly analogous to the equivalent optical forces [4,[7][8][9]11].For small particles (ka 1), the acoustic gradient force is given by [14]: where ρ represents the mass density, κ = (ρc 2 ) −1 is the compressibility, and c denotes the sound speed; subscripts p and 0 refer to properties of the scattering particle and background medium, respectively.Angle brackets denote a time average over one oscillation period (i.e.g = 1 T T 0 g dt).In general, this force will tend to pull particles into pressure maxima or minima, depending on the particle properties.
In this manuscript we consider 'sound-hard' particles, for which the particle density is much greater than the background medium and the particle compressibilty is much less than the background medium.In this approximation -which is quite good for solid particles in a gaseous medium -f 1 = f 2 = 1.Conceptually, the gradient force can be regarded as being caused by a single scattering of the acoustic wave from the particle, and consequently it scales like the particle volume.For larger particles (ka > ∼ 1), this force must be computed numeri- Both maps are created by keeping one particle (gray quarter-sphere) fixed at the origin, which here corresponds to a location in the incident field halfway between node and antinode, and computing the scattered fields (and hence total force) on a second particle centered at a given location on the map.The single-particle force is subtracted from the computed total force at each location so that these maps show only the radial binding force between two identical particles.For the acoustic maps the MFS code was implemented, see supplemental materials section E for a discussion of the simulation of optical forces.The forces in these plots are all scaled to the corresponding acoustical/optical reference force.The hashed yellow regions correspond to exclusion zones within which the two particles would overlap.
cally from the scattered field.
Other forces arise due to higher order scattering events.We will separate these into the radiation pressure force, which acts on a single particle, and the acoustic binding force, which acts on clusters of two or more particles.For small particles (ka 1), both of these effects scale like the square of the particle volume (fig.1), and as a result are only comparable in magnitude to the gradient force when the particle is of order wavelength in size.The radiation pressure force pushes the particle in the propagation direction of the acoustic wave; this is often ignored because it cancels in counter-propagating beams (e.g. as used in a conventional acoustic levitation setup).
The acoustic binding force arises due to the interference of the external field and the field scattered from a neighboring particle.Because interference modifies the total field intensity in the vicinity of the particle, nearby particles feel an additional force due to this modulation.For small particles, this can be approximated as the gradient of the total field (fig.1b), but this approximation breaks down in the Mie regime (fig.1c).The binding force can be particularly difficult to model when particles are close together due to strong multiple scattering events.
Other acoustic forces can arise due to acoustic streaming or deformation of the particles; each of these effects can also induce interparticle forces under the right circumstances.However, both of these effects are negligible for wavelength sized solid particles in air.
Acoustic streaming is a non-zero mean flow induced by an acoustic field [48], which subsequently induces a drag force on any particles embedded in this flow [48][49][50].One can estimate a critical particle radius below which streaming dominates over other second order acoustic effects [51]: This threshold depends on the kinematic viscosity of the fluid (ν), the sound frequency (ω), acoustophoretic contrast factor (Φ), and a factor that depends on the geometry of the stream generating surface (Ψ) which is of order unity.For 40 kHz ultrasonic fields in air this critical radius is around 10µm, and so the scattering forces dominate for the Mie sized (a ≈ mm) bodies investigated here.
Bjerknes forces result from deformation of 'particles' in the acoustic field, giving rise to a non-zero time averaged pressure on the particle surface, and are generally significant only for situations where the particle compressibilty is the same or greater than the background medium (e.g.bubbles in water) [35,36].For the soundhard approximation made here, this force is identically zero, however for solid particles in a gaseous medium it will be orders of magnitude smaller than scattering induced forces.
Finally, we note that in this manuscript we do not consider acoustically mediated torques which may cause individual bodies to spin.The modelled incident field (planar standing wave) carries no angular momentum, as is generally the case for the fields in conventional acoustic trapping devices [29,32,33,37,41].Scattering between three or more bodies can yield asymmetries in the total sound field at the surface of these bodies, however unless they are sound absorbing an interaction torque will not be imparted [33,46,52].
For the remainder of this manuscript then we will consider only forces generated due to sound scattering events.The scaling of these three forces (gradient, scattering, and binding) can be computed analytically for particles in the Rayleigh regime [14,39], but this quickly becomes intractable for larger particles.This is especially problematic as this is precisely the regime in which the binding force becomes comparable to the gradient force.Although approximate forms exist for very weakly scattering or widely separated particles [41], such approximations are not valid for realistic experimental conditions.This is especially true if one wishes to use acoustic forces to guide the self-organization of large particle assemblies, in which case we expect many particles to be in physical contact.

II. NUMERICAL METHODS
We model the acoustic field using a complex oscillating potential field, φ(r) ∝ e −iωt , which is related to firstorder velocity, pressure, and density fluctuations in the following way [53]: where ω is the sound frequency, p is the pressure, ρ is the density, and v is the local fluid velocity.In all cases, the physical fields are given by the real part of the complex representation.This linear wave field produces zero net force on an embedded particle if averaged over a full cycle, however, such a formulation does not satisfy the Navier-Stokes equation without including higher order terms.A perturbation approximation gives rise to second order corrections -p 2 , ρ 2 , and v 2 -which can be formulated in terms of the first order terms.These second order terms invoke a nonzero net force when averaged over a cycle of the wave [14]: where the integral in the force calculation is taken over the surface of each particle, and φ represents the total potential field, including both external and scattering contributions.Here we have assumed that p 1 p 0 , in which case higher order terms (p 3 , etc.) are negligible [53].In this work we are interested in the case of solid particles in a gaseous medium.Using the sound-hard approximation, this is equivalent to v 1 • n = 0, removing the second term in eqn.(9).(See Supplemental Materials section D for a discussion of the validity of this approximation.) We solve for the total first order field, φ, including an external driving input and multiple scattering effects from more than one particle, using the method of fundamental solutions [42].In this model, the scattered field outside the particles is computed from a lattice of virtual scatterers placed inside each particle.We solve for the amplitude of these virtual scatterers by enforcing the sound-hard boundary condition on another lattice of evaluation points located on the surface of each particle.In practice, this is done by casting the problem as a linear matrix equation and solved using standard numerical techniques.The number of points in each lattice affects model accuracy at the cost of computation time; the larger the scattering body, the greater number of source and evaluation points are required to fully resolve the local field oscillations on the surface of that body.Since the method of fundamental solutions is a surface-based scattering approach, it requires less computational resources than volume-based scattering approaches such as FDTD, DDA, and FEM [40,48,54,55].Once the total field has been solved, it is numerically integrated over the surface of each particle using a Gaussian quadrature rule [56] to obtain a per-particle force.More details on the numerical method, and an analysis of its accuracy, are provided in Supplemental Materials section B.

A. Acoustic Force Scaling
The relative strength of binding, scattering, and gradient forces can be computed for particles of arbitrary size using MFS (fig.1).In particular, we note that gradient forces dominate for ka 1, while the other forces become comparable in magnitude only for ka > ∼ 1.Although the optical binding force between a pair of particles is ∼1 order of magnitude weaker than the scattering force (and has approximately the same scaling with size), FIG. 3. a, b) First five stable separation distances for two particles side-by-side (a) and for a hexagonal cluster of seven particles (b).The size of the encompassing blue marker at each distance encodes the strength of the relative spring constant associated with that equilibrium location.The effective spring constants are de-dimensionalized according to the acoustic wavelength and corresponding reference force for a given particle size.Bottom: Acoustic binding force between two (c) and seven (d) particles of size parameter ka =1 (black), 2 (red), and 6 (green), shown as vertical lines of the corresponding color on a, b.For the studies involving seven particles, the cluster was arranged in a regular hexagonal pattern and the force was computed on the right-most particle in the cluster.
this is mitigated by two factors: 1) for clusters of many particles the binding forces will be additive, while the scattering force is not, and 2) the scattering forces can be cancelled with appropriate field design.
In the results that follow, particles are suspended in a pair of counter-propagating acoustic plane waves travelling in the ±z-direction.This is consistent with the design of common acoustic levitation experiments [26,28,32], apart from the fact that we do not assume the beam is strongly focused in the transverse direction.The interference between these beams confines the particles to a single plane in z, but allows them to move freely within this plane without experiencing gradient or scattering forces.Unless otherwise specified, particles are simulated at a pressure node of the incident field, which has field amplitude of 200 Pa and frequency 40 kHz.
We find it useful to introduce here a reference force, to which one can scale all other relevant effects.Here, we use the force experienced by a perfectly absorbing sphere within an external acoustic field, F 0 : where p 1 is the amplitude of the oscillating external pressure field.

B. Pairwise Interactions
We first consider the force between pairs of particles in our acoustic standing wave (fig.2).In general, we observe an oscillating force -with a period given by the sound wavelength -which falls off like 1/r and is primarily in the radial direction.This arises because the main contribution to the acoustic binding force is a gradient force in the combined incoming field and the field scattered from the neighboring particle, whose amplitude falls off like 1/r.(Although the binding force can in part be explained as a 'gradient force', we do not treat it as part of the gradient force because it is only present for multiple particles.)Note that the force is generally quite weak if the particles are exactly on the pressure node (z = 0); this is because the incoming field has no amplitude here.In practice, gravity will displace particles from the anti-node, in which case this transverse force will be stronger.
For larger particles (ka > ∼ 2), the scattering becomes more intense along the z-axis, and so the binding force does as well.Interestingly, the qualitative features of the acoustic binding force are quite similar to optical binding (fig.2d-f ).This is explained by the fact that both are second order radiation scattering forces which arise for the same general reasons; the differences can be at-tributed to the fact that optical radiation is vector while acoustic fields are scalar.We note also that the optical forces shown in fig. 2 are for a relative refractive index of 1.5 -approximately equivalent to glass particles in air -while the acoustic forces are computed for an infinite impedance contrast.
The features of either acoustic or optical binding have many characteristics which make them interesting candidates for self-or directed-assembly of particle clusters.First, we note that the oscillatory nature of the force means that the particles have multiple stable separations whose distance can be tuned by changing the field wavelength.We characterize this with an effective spring constant, k ef f = −∂F r /∂r, at each of the stable particle separations (fig.3).For closely spaced particles, there are additional near-field effects which are strongly dependent on size parameter, ka.Indeed, the binding force can be either attractive or repulsive for contacting spheres (fig.3c), and as a result this is tunable via sound frequency.Although not considered here, we also note that a sound wave composed of several different frequencies could superimpose the interference patterns and produce an even more complicated -and tunable -interparticle force.
A pair of particles in a viscous fluid is effectively a damped harmonic oscillator, and so can be characterized by a dimensionless quality factor, Q.If we assume a Stokes drag law for the particles, we can compute the quality factor as a function of particle size and density (see Supplemental Materials section C).For experimental parameters relevant for solid particles in a gaseous medium, this quality factor scales like Q ∼ 6 (ka) 7/2 ρp 1kg/m 3 1/2 for ka < ∼ 3, and has a somewhat more complicated structure for ka > ∼ 3. (Note that this scaling assumes the incoming wave pressure is increased relative to the gravitational force so that F 0 = 4.25mg.)As a result, we expect oscillations of particle pairs (or many particle clusters) to be under-damped for realistic experimental parameters, in contrast to optically bound systems which -assuming wavelength sized particles and visible optical fields -should be in the over damped regime unless the particles are suspended in vacuum.As we shall see in section 'Evolution of Many Particle Clusters', this has significant implications for the formation of acoustic clusters with many particles.We also note that a dependence of cluster stability on damping level has previously been predicted for optically bound systems [12].

C. Acoustic Binding of Many Particles
How does the acoustic binding force scale when we have more than two particles in the field?The total force on each particle is determined by a combination of the incoming and scattered fields, which manifest in a net force through the p 2 term.As p 2 relies on the square of the field values, there are three terms involved in a force calculation: φ 2 in , φ 2 sc , and φ in φ sc .For weakly scattering particles, or for bodies interacting in the (scattered) far-field one can neglect φ 2 sc as it is very small in comparison to the other terms.In this case, the force is (approximately) linear in φ sc , in which case one would expect the force to be nearly pairwise; indeed this simplifying assumption has been made in some previous research [41].However, for strongly scattering particles which are close together the scattered field can become quite strong (φ sc ∼ φ inc ), and so this approximation should break down.
Thus, in the ka > ∼ 1 regime, we can no longer assume that the stable particle separation distances for many particle clusters will be the same as for pairs.To explore this effect, we compute an effective spring constant for a hexagonal cluster of 7 particles, using the radial force on one of the outer particles (fig.3b, d).For small particles (ka ≤ 1.8) we observe that the stable separations shift from r ∼ 1λ to r ∼ 1.18λ.For larger particles, we observe a more complicated shift in the the stable positions, consistent with the observation that multiplescattering effects should be stronger for larger particles.From these observations, we can conclude that a pairwise force approximation is not appropriate for closely spaced wavelength sized particles, and that the full N -body force calculations must be performed for accurate results.
A pair of particles at the same z coordinate -but displaced in x/y -must have the equal and opposite forces by symmetry.For more than two particles, this is no longer the case, especially once multiple scattering effects are included.To probe for these effects, we consider the forces between a stable three-particle triangular cluster and a fourth particle placed nearby (fig.4a).If we compare the net forces on the cluster and the lone particle, we find they are not equal and opposite, and this effect becomes stronger as the lone particle gets closer to the cluster.(fig.4b).Conversely, if we compute the forces by summing over pairs of particles, we find they are equal and opposite, as expected.Although this nonconservative force appears at first to violate Newton's third law, this is not the case: the net momentum is carried away by the scattered acoustic field.This simple result demonstrates that acoustic binding can produce non-conservative driving forces, but only when the full N -body force is considered.As we shall show in the next section, it is also possible to produce stable clusters of particles with a non-zero force on the entire cluster.

D. Evolution of Many Particle Clusters
The results described up to this point have been computed only for fixed particle locations.To understand how many particles freely moving in an acoustic binding force evolve in time, we have coupled our force model to a simple molecular dynamics simulation.To prevent particle overlap (which would cause the MFS model to produce non-physical results), we also include a short The force on an extra particle placed near the cluster, as well as the additional force on the cluster itself.The non-zero net force implies the presence of non-conservative effects which can give rise to drifting of stable clusters in the right configuration.The forces computed from a pairwise model are plotted as dashed lines and show clear discrepancy compared to the actual forces, especially when the lone particle is in close proximity to the cluster.range repulsive force of the form:

Model Parameters
This force is designed to turn on slightly before the particles are in actual contact, but to do so at such a distance that it does not significantly affect the results (see Supplemental Materials section C).We also include a gravitational force (F g = −9.8m/s 2 ẑ), and a Stokes drag (F d = 6πµav), where we assume the background medium is air at STP.Typical particle velocities are of order 30 mm/s or less, and so the Reynolds number of the flow around a particle is Re < ∼ 7; as a result a Stokes drag law is a reasonably good approximation.We do not consider hydrodynamic coupling between the particles.Time stepping in this model is handled using an adaptive timestep Runge-Kutta integrator (the Dormand-Prince method), specified to have an absolute velocity error tolerance of 10 −6 m/s.No thermal fluctuations are considered, as they should have negligible impact at these scales.
The simulations were performed for sets of 100 trials at various damping parameters, with initial particle locations assigned according to a Gaussian distribution with σ = 3a (see supplementary materials section B).All sets of trials were performed for ka = 1.466, which represents particles of size a = 2 mm and a driving frequency of 40 kHz in air.The particle density and incident field amplitude were both varied between trial sets in order to achieve identical binding strengths with varying levels of damping; the density and amplitude values used can be found in Table I and correspond roughly to Aerogel, expanded polystyrene, and urethane foam, as well as the driving pressures necessary to levitate bodies of these densities.(As before, the pressure is modulated so that F 0 = 4.25mg irrespective of density.)The quality factors for pairs of particles with these properties is Q ∼ 11, 34, 107.In these time-evolution studies we allow for motion in the z-direction so that the gravitational force pulls particles to a natural levitation plane which sits slightly below a pressure node.The simulations run for a predetermined time (10 s) or until a terminal event is reached, which here corresponds to acceleration of all particles falling below a minimum threshold of |a| ≤ 10 −6 m/s 2 .
In general, we observe that randomly placed configurations of 5 particles will arrange themselves into stable clusters, with a strong preference for near-integer wavelength particle spacing (fig.6).Interestingly the first preferred interparticle separation distance within all sets of trials converges to 1.14λ, which is slightly greater than the value which would be expected from pairwise interactions.This separation can be reproduced by instead considering the forces on a 7-particle hexagonal cluster, which has a nearly identical equilibrium separation (fig.3b).
Damping plays a critical role in the formation of the clusters, as seen in fig.5a-c.As expected, the relaxation time is much shorter for the trials with higher damping.Although the level of damping should not change the equilibrium configurations, weakly damped systems have difficulty shedding their kinetic energy so they can  and c) low-damped systems, all beginning from the same state.The high and mid-damped trials both resulted in the common 'boat' configuration, while for lower levels of damping one of the particles is ejected from the cluster.The colormap indicates the second order pressure field on the x-y plane surrounding the particles, and the particle tracks are colored sequentially from white to dark red to represent time progression.d) A simulation of a quasi-stable 'manta-ray' configuration, showing the instantaneous states at various progressive times.Scale bar is 1cm in all panels.See supplementary movies 1 and 2 to see 'boat' and 'manta ray' formations respectively.settle into a stable cluster.On average, this produces clusters with larger particle separations as the damping is reduced, or, equivalently, the density is increased.The low-damped trials produced virtually no compact fiveparticle configurations; at least one particle was generally accelerated through a local minima without enough damping to slow it down, and was effectively ejected away from the rest of the particles in a cluster (e.g.fig.5c and supplementary movie 4).This behavior seems to become prominent when the pairwise quality factor is Q > ∼ 50, and is consistent with results previously seen in simulations of optically bound clusters [12], which in some cases require a minimum level of damping to form stable structures [57].
For all levels of damping, we find that most of the clusters are not stationary, but rather have some nonzero drift velocity when the cluster is asymmetric.Thus, we characterize the final structures in terms of the moment of inertia and drift speed of the final state of the simulation (fig.7), with cluster moment of inertia given by: where the sum is over all particles and x/ȳ denote the center of mass of the configuration in the x/y direction.A number of cluster geometries were commonly produced throughout the mid and high damped trials, with the most common shape being what we term the 'boat' configuration (fig.5a,b).The boat -which is one of the few observed clusters with no drift velocity -comprised 12% of the mid-damped and 18% of the high damped final configurations (videos of this and other structures forming can be seen in the supplemental materials).We do not observe a correlation between moment of inertia and drift speed, suggesting that extended particle clusters are not required to produce drifting configurations.Additionally, the drift speed can be rescaled by a reference velocity which would be experienced by a single particle subject to a force of F 0 -this is equivalent to the speed a perfectly sound absorbing particle would move in a single plane wave of sound.Notably, the fastest moving clusters all contain a pair of particles in contact, which may effectively behave like a larger particle and so experience higher non-conservative driving forces.
In rarer cases we also observe quasi-stable clusters which oscillate or rotate in time; an example is shown in fig.5d and supplementary movie 2. In this configuration -which we term the 'manta ray' -the cluster itself undergoes oscillations while it drifts in space.This oscillatory motion continues indefinitely with no perceptual loss in amplitude, with the driving of the local oscillations perfectly balanced by the system damping.The oscillatory behavior of the manta ray configuration is only observed only for the mid-damped case; the high-damped case produces nearly identical configurations: these do not oscillate but do drift at a high velocity.
IV. DISCUSSION Acoustic binding produces complex, long range forces which make it an intriguing candidate for directed assembly of mm-scale systems.Although previous work has considered simplified models of this force [37,38,41], we find that non-pairwise and non-conservative effects significantly modify the resulting structures and their stability.This is especially true in the size regime where the multi-particle binding force is strongest compared to gradient and scattering forces, ka ∼ 1 − 10.If one wishes to use this force for on-demand assembly of complex structures, this has both positive and negative at-tributes.On one hand, the modifications of the force due to multiple-scattering effects could be exploited to assemble more complex structures without using a more complex incoming field.Moreover, the non-conservative driving forces make acoustically bound clusters analogous to active-matter systems which are currently being explored for a variety of applications [43,44,58].Unlike typical active matter systems, these driving forces are the result of specific particle configurations, and soin principle -they are both switchable and tunable.It is likely that non-spherical particles or even mixtures of differently sized particles could enhance these effects, as they arise due to symmetry breaking of the stable configurations.
Conversely, the observation that the many-body forces must be considered raises practical difficulties in computing the acoustic binding forces and predicting the structures they produce.This will be especially true as the number of particles is increased, since the computational cost of a fully coupled N -body system scales like N 3 .Although this could be somewhat mitigated using advanced numerical techniques (such as the fast multiple method [59] or GPU-based simulations), it may also be possible to produce empirical models of the many-body force which are more computationally efficient.In either case, with further numerical optimization these models could allow acoustic binding to be used for 'inverse-design' approaches where an incoming field (or combination of several) is used to produce a desired structure on demand [60][61][62].
Finally, we note that -to our knowledge -non-contact acoustic binding effects have only been observed in a single experiment, and the largest clusters observed had only 3 particles [32].There are likely two reasons for this: 1) most acoustic levitation experiments used focused fields which do not have room for particles that are transversely displaced by one or more wavelengths, and 2) driving forces make configurations of solid particles (ρ ∼ 1000 kg/m 3 ) unstable in air.Both of these limitations could be overcome with appropriate experimental design.Indeed, there is no fundamental requirement to use focused sound fields for acoustic levitation.Moreover, the relative amount of damping can be modified either by using a different background medium (e.g. a liquid rather than a gas, as used in the aforementioned experiment), or using higher frequency sound waves and smaller particles.Such a system would have both fundamental and practical applications.Understanding the self-organization of actively driven systems has attracted considerable interest in the last decade; acoustic binding provides a tunable platform to study these effects on easily accessible length scales.With continued research, acoustic binding could ultimately provide a practical method for controlling the assembly of particle clusters or even meta-materials composed of large numbers of particles in active or passive configurations.

Dynamics of Acoustically Bound Particles: Supplemental Information
A. WEAK SCATTERING THEORY Bruus et. al.[? ] showed that for small spherical scatterers (ka 1), the trapping force can be formulated in terms of the gradient of an acoustic potential, which itself is defined in terms of incident field values at the location of the particle as well as relevant particle/fluid parameters (eqn.2).This arises in weak field scattering theory by invoking a multipole expansion to model the scattered field.In this case one can easily express the monopole and dipole terms of the scattered field in terms of the incident field.For sound hard bodies the acoustical contrast terms (eqns.3-4) are both unity, simplifying much of the underlying theory.
where subscripts p and 0 denote particle and medium properties respectively.
We compare our MFS model to this small particle formulation in order to validate our model.Aside from this we also compare our model to a solution given by a truncated harmonic expansion (which is exact for the case of a single particle in an axially symmetric incident field).To compare the three models a particle is held fixed in between a node and an antinode of the external pressure wave (planar standing wave, amplitude 200Pa, frequency 40kHz), and we compute the forces according to each model as the size of the particle is varied.We see that for ka 1 all solutions agree, whereas for larger sized particles the weak scattering theory described by Bruus will overestimate the radiation force, while the force from the MFS model exactly matches the exact solution (fig.1).
FIG. 1. Percent error between acoustic trapping force computed from small particle theory and either our MFS code (dashed orange), or an exact solution (solid blue).Note how all three models give the same result as particle size approaches the Rayleigh limit.

arXiv:2111.08479v5 [physics.class-ph] 8 Mar 2022
The exact solution can be expressed as a truncated harmonic expansion: In the above h n is the spherical Hankel function of the first kind, and P n is the n th Legendre Polynomial.The weights C n are chosen in order to solve the boundary condition on the surface of the sphere.For the sound hard case that boundary condition is ∂ r (φ in + φ sc ) = 0 on the surface of the sphere.

B. NUMERICAL METHODS AND MODEL ACCURACY
We model the scattered sound field from a sound-hard spherical body using the Method of Fundamental Solutions [? ].In this method N image sources are placed outside the region of interest, in this case on a concentric lattice inside the spherical body, and the scattered field is given by a weighted superposition of fundamental solutions (Greene's functions) emanating from these image sources: where k is the angular wavenumber, ω is the angular frequency, and R = |r−r j | is the distance from the j th source point to the field location r.The weights are chosen by enforcing the sound-hard boundary condition (∂ r (φ in + φ sc ) = 0) at a second lattice of evaluation points (r j ) on the surface of the sphere through casting the problem as a set of linear equations: Note that the number of source and evaluation points must be the same so that the system is exactly solvable.Various lattice types for the source and evaluation points were tested using this model, including a fibonacci lattice, and lattices involving points subdivided onto a cube or an icosahedron and projected onto a unit sphere.Ultimately we found that using the icosahedral projection lattice for both the source and evaluation points gave the best results due to its high degree of symmetry.To integrate forces over the surface of the particle, we employ the Gauss-Legendre quadrature method [? ]: This method involves summing weighted values of the integrand over the surface of the sphere at various nodal points (x i ), in this case points spaced evenly at different (also evenly spaced) polar angles of the sphere.The interval at each polar location is rescaled so it fits on the interval -1 to 1, and the weights (w i ) depend on the location of each node, as well as the legendre polynomial associated with each location on the interval from -1 to 1: To determine the accuracy of our model, we again used the same simple scenario described in the section A, this time varying the number and of source points, and the number of quadrature points (upon which the Legendre-Gauss quadrature rule is applied).As number of source and evaluation points affects computation time, we are interested in the minimum number of points within each lattice necessary to obtain a satisfactory force computation.To answer this question we compare the force computed using MFS to the exact solution according to the Harmonic expansion described above.Note that the exact number of quadrature points is given by 2q 2 where q is the quadrature number  3. Quality factor vs particle size for acoustic binding between a pair of particles within a planar standing wave incident pressure field of amplitude p1 = 200 Pa and frequency ω = 2π × 40 kHz.The particles are fixed on the plane halfway between a node and an antinode and the separation distance is varied between the two particles in order to find the locations of equilibrium.Effective spring constants are then taken as the slope of the binding force at these locations.This plot was created using the spring constants associated with the first separation distance for the various particle densities simulated in our time-evolution studies.
Python code to compute these forces, as well as notebooks used to compute errors and other quantities can be found in a Github Repository for the project [? ].

C. QUALITY FACTOR
Since we have a measure of an effective spring constant (k ef f ) associated with various equilibrium locations of particles of different sizes and these bodies experience viscous linear damping, we can define the quality factor associated with these damped oscillatory systems.
In the small particle limit the acoustic binding force (and hence associated effective spring constant) scales as the particle volume squared, using this fact together with (eqn.12) one concludes that the quality factor scales as ka 7/2 in the Rayleigh limit (fig.3).

D. SOUND HARD APPROXIMATION
The sound hard boundary condition describes a scenario where no bulk transmission of sound occurs in the scattering body.We chose this approximation as most experimental work focuses on solid bodies levitating in air, for which there is a high density and compressibility contrast between the particles and medium.Making this approximation also simplifies the MFS solution, as we can ignore the sound field inside the particles.Relaxing this approximationand allowing for sound to penetrate the scattering particles -is possible using MFS, but requires a more complicated model which also includes an extra set of virtual scatterers.
Just how sound penetrable a body is depends on the relative measures of the density and sound speeds within a scattering body and the host medium.For normal incidence of a planar standing wave on an infinite 2D interface, the reflection coefficient is given by [? ]: where subscripts u and l denote the media above and below the refracting interface respectively.We note that for solid bodies in a gaseous medium the speed of sound and density of the scattering body will both always be much greater than those of the gas, so that m > 1 and n < 1.Even for a very light solid material -such as expanded polystyrene (EPS) -the reflection coefficient is R ∼ 1.The density and speed of sound in air are ρ air = 1.225 kg/m 3 and c air = 343m/s.The density of EPS varies, but for estimation purposed we will assume ρ EP S = 10 kg/m 3 .The speed of sound in EPS is dependent on the porosity, and ranges from 2000 -3000 m/s as the porosity is decreased [? ].Thus we expect a reflection coefficient of R ∼ 0.96, implying that the sound hard approximation is appropriate in this case.

E. OPTICAL BINDING CALCULATIONS
Optical binding forces are computed using the Discrete Dipole Approximation (DDA) [? ?].To implement DDA, each particle is subdivided into many individual pieces, each with a size of ∼ λ/10 which is then treated like a polarizable point source.(In practice, spherical particles are divided an exact integer number of times, so that the size of each chunk is as close to λ/10 as possible.)The multiple scattering problem in the presence of an incoming field can then be expressed as a matrix equation, similar to MFS.This matrix can be solved using a variety of methods [? ].In our case, we use a hybrid solution where the internal scattering in a single particles is solved using an explicitly inverted internal scattering matrix, and scattering between particles is handled using an iterative approach.In practice 3 iterations are used for ka = 1 -effectively considering up to four scattering events -and 4 iterations are used for ka = 2 − 4. In all cases the convergence of the solution is 10 −2 or better (measured as the relative root mean squared change in the dipole strength per iteration).
Forces are computed per dipole using standard techniques [? ] and summed to compute per-particle forces.We have compared the results of this solution to existing DDA simulations [? ] and found nearly identical results.

FIG. 1 .
FIG.1.Acoustic force scaling.a) Relative strengths of various acoustic forces as a function of size of the scattering body, computed using MFS numerical model.For ka < ∼ 1 the gradient force (blue) dominates; for ka > ∼ 1 the secondary forces (red, green, orange) are of comparable magnitude.The binding force is shown for two particles at fixed separation distances of 1.25λ (green) and 7.25λ (red) while the gradient force is for a single particle located between a pressure node and antinode.For reference the particle volume (dashed blue) and square of the volume (dashed orange) are included.b,c) Two particle acoustic binding force (black), 2nd order pressure field ( p2 ) (red), and gradient of the 2nd order pressure field integrated over particle surface (dashed gray) vs particle separation for (b) ka = 1 and (c) ka = 10.

FIG. 2 .
FIG.2.Acoustic vs. optical binding.Force maps for acoustic binding (a-c) and optical binding (d-f) for different sizes of particles, ka = 1 (c, f), 2 (b, e), and 4 (a, d).Both maps are created by keeping one particle (gray quarter-sphere) fixed at the origin, which here corresponds to a location in the incident field halfway between node and antinode, and computing the scattered fields (and hence total force) on a second particle centered at a given location on the map.The single-particle force is subtracted from the computed total force at each location so that these maps show only the radial binding force between two identical particles.For the acoustic maps the MFS code was implemented, see supplemental materials section E for a discussion of the simulation of optical forces.The forces in these plots are all scaled to the corresponding acoustical/optical reference force.The hashed yellow regions correspond to exclusion zones within which the two particles would overlap.

FIG. 4 .
FIG.4.Non-conservative effects in acoustically bound systems.a) The second order pressure field, p2, near a stable triangular cluster with ka = 1.466.b) The force on an extra particle placed near the cluster, as well as the additional force on the cluster itself.The non-zero net force implies the presence of non-conservative effects which can give rise to drifting of stable clusters in the right configuration.The forces computed from a pairwise model are plotted as dashed lines and show clear discrepancy compared to the actual forces, especially when the lone particle is in close proximity to the cluster.

FIG. 5 .
FIG.5.Acoustically bound cluster evolution.Initial (gray transparent) and final (white) states produced using simulation for a) highly damped, b) mid-damped, and c) low-damped systems, all beginning from the same state.The high and mid-damped trials both resulted in the common 'boat' configuration, while for lower levels of damping one of the particles is ejected from the cluster.The colormap indicates the second order pressure field on the x-y plane surrounding the particles, and the particle tracks are colored sequentially from white to dark red to represent time progression.d) A simulation of a quasi-stable 'manta-ray' configuration, showing the instantaneous states at various progressive times.Scale bar is 1cm in all panels.See supplementary movies 1 and 2 to see 'boat' and 'manta ray' formations respectively.

FIG. 6 .
FIG. 6. Final interparticle separations.Histograms of the final interparticle separation distances for three sets of 100 trial MD simulations.(a) high damping, (b) mid damping,(c) low damping.In all cases the first preferred separation distance is 1.14 λ, indicating the importance of multibody effects.

FIG. 7 .
FIG. 7. Cluster drift speeds and moments of inertia.Histograms of the drift speed (a) and moment of inertia (b) of some of the common clusters found in the trial simulations.Blue correspond to mid damping, while orange figures correspond to high damping.Data from trials involving low damping were omitted as the damping was not sufficient to lock particles into a local minimum, i.e. configurations tended to fly apart.

FIG. 2 .
FIG.2.Discrepancy between MFS computed force and that from the Harmonic expansion.The heatmap represents the number of decimal points of error.This was done for a single particle of size ka = 1 located between the node and antinode of the external field.
FIG.3.Quality factor vs particle size for acoustic binding between a pair of particles within a planar standing wave incident pressure field of amplitude p1 = 200 Pa and frequency ω = 2π × 40 kHz.The particles are fixed on the plane halfway between a node and an antinode and the separation distance is varied between the two particles in order to find the locations of equilibrium.Effective spring constants are then taken as the slope of the binding force at these locations.This plot was created using the spring constants associated with the first separation distance for the various particle densities simulated in our time-evolution studies.