Ultrasound directed self-assembly of three-dimensional user-specified patterns of particles in a fluid medium

Ultrasound directed self-assembly (DSA) is defined as the process by which discrete components organize through interactions between themselves and an external ultrasound wave field. This process is similar to other DSA techniques driven by external fields, including electric,¹ magnetic,² or flow fields.³ The external field acts as a tunable mask that enables organizing particles dispersed in a fluid medium into user-specified patterns by adjusting the operating parameters of the transducers, such as their amplitude and phase, which generate the external field. However, electric and magnetic field-based DSA requires particles with specific electric and magnetic properties, respectively, and typically demands an ultra-high field strength, which limits dimensional scalability.^4,5 In contrast, ultrasound DSA relies on the acoustic radiation force associated with a standing ultrasound wave field to organize particles dispersed in a fluid medium into a specific pattern, independent of their material properties. Additionally, dimensional scalability is only limited by attenuation of the ultrasound wave, which can be mitigated by selecting a low-viscosity (bulk and shear) fluid medium.⁶ Thus, ultrasound DSA enables organizing patterns of particles for a number of engineering applications, including biological cell manipulation,⁷ microbubble filtration,⁸ and fabrication of engineered materials with exotic properties based on specific patterns of particles embedded in a matrix.^9–11 

Employing ultrasound DSA to create user-specified patterns of particles requires relating the ultrasound transducer arrangement and operating parameters (amplitude and phase) to the resulting pattern of particles assembled with a specific standing ultrasound wave field. This reduces to a “forward problem” that entails calculating the pattern of particles that result from any user-specified ultrasound transducer arrangement and operating parameters and an “inverse problem” that involves calculating the ultrasound transducer operating parameters required to assemble a user-specified pattern of particles.¹² The forward problem has been analyzed extensively in the literature.^13–15 The inverse problem has been solved analytically to enable manipulating single particles or a limited number of simple patterns of particles in one-dimensional (1D),^16,17 two-dimensional (2D),^18,19 or three-dimensional (3D) fluid reservoirs.^20–22 However, each direct solution technique only applies to one specific pattern, fluid reservoir geometry, and ultrasound transducer arrangement and, thus, limits the applicability of these existing methods. Alternatively, solving the inverse problem using numerical optimization techniques has enabled manipulating single particles in 3D²³ or creating user-specified patterns of particles in 2D.¹² Although critical to employing ultrasound DSA in engineering applications, no universal direct inverse solution method exists that enables ultrasound DSA of user-specified 3D patterns of particles in a fluid reservoir with arbitrary geometry and ultrasound transducer arrangement.

Thus, the objective of this work is to theoretically demonstrate and experimentally validate a direct and universal solution of the inverse ultrasound DSA problem for a 3D fluid reservoir with arbitrary geometry and ultrasound transducer arrangement. We first formulate the ultrasound wave field within a 3D fluid reservoir as a function of the ultrasound transducer operating parameters using the boundary element method based on Green's third identity.²⁴ Then, we apply the acoustic radiation force theory to relate the ultrasound wave field to a user-specified pattern of particles.^13,14 Finally, we pose the inverse problem as a constrained optimization problem that reduces to eigendecomposition.¹² 

Figure 1 shows an arbitrary-shaped 3D fluid reservoir lined with N_t ultrasound transducers of acoustic impedance Z_t(ω₀) around the boundary surface S, where ω₀ is the operating frequency of the ultrasound transducers. We have removed an elliptical section of the boundary surface in Fig. 1 to show the interior 3D solution domain D, which contains a fluid medium of density ρ_f and sound speed c_f. Additionally, the inset image in Fig. 1 shows a differential element ds(q) and normal direction n(q) of the boundary surface at point q on S, with reference to the fluid reservoir origin o. We compute the ultrasound wave field in the fluid reservoir in terms of the complex scalar velocity potential φ using the following assumptions. (1) The ultrasound wave field φ satisfies the Helmholtz equation $∇2φ+k02φ=0$ throughout D and on S, where k₀ = ω₀/c_f is the wave number of the ultrasound wave field in the fluid medium. (2) φ satisfies the impedance boundary condition $∂φ/∂n+ik0Z̃φ=v$ at each point q on S to account for reflection and transmission of the ultrasound wave field at the ultrasound transducer surfaces or inactive boundary surfaces (v = 0). Here, $Z̃=ρfcf/Zt$ is the impedance ratio of the fluid and the ultrasound transducer and v is the complex harmonic velocity amplitude of ultrasound transducer surface as it vibrates in the normal direction n. We have assumed smooth, inelastic ultrasound transducers of which the impedance has no angular dependence. While not strictly valid in practice, it still allows accurately approximating the standing ultrasound wave field within many fluid reservoirs, evidenced by good agreement between experimental and simulated patterns of particles.^12,16–19 We calculate the ultrasound wave field φ at any location x in D or on S using Green's third identity accounting for the impedance boundary condition as²⁴ 

Ω(x) = 1 if x is in D or 2 if x is on S, $i=−1$⁠, and G(q,x) is the Green's function that represents the free-field ultrasound wave field emitted from a unit-magnitude point source at location q on S and measured at location x in D or on S. In the case of a 3D ultrasound wave field, we calculate G(q,x) as⁶ 

where |q - x| is the Euclidean distance between points q and x. We discretize the fluid reservoir into N_d domain points and its boundary surface into N_b > N_t boundary elements, and calculate the ultrasound wave field in D using the boundary element method, which approximates the integral over boundary surface S in Eq. (1) as a summation over the N_b boundary elements representing S.²⁴ Thus, using the boundary element method, we relate the ultrasound wave field φ to the operating parameters v of each of the N_t ultrasound transducers, arranged into a single vector $v=[v1,v2,…,vNt]$⁠. We relate the ultrasound wave field in the 3D fluid reservoir to the resulting pattern of particles by calculating the acoustic radiation force f acting on each particle of density ρ_p and sound speed c_p as $f=−∇U$⁠, where U is the acoustic radiation potential. For a spherical particle in an inviscid fluid medium, with particle radius r_p ≪ λ₀ and ultrasound wavelength λ₀ = 2π /k₀ in the fluid medium, we calculate U as^13,14

Here, $βp=(ρpcp2)−1$ and $βf=(ρfcf2)−1$ are the compressibility of the particle and fluid medium, respectively, and the notation ⟨|φ|²⟩ indicates the mean-square fluctuations in the velocity potential φ over one period T = 2π/ω₀ of the ultrasound wave field. We treat the particle as a fluid and, thus, neglect the effect of shear stiffness.^13,14 From Eq. (3), we compute the acoustic radiation force f and identify the pattern of particles as the stable fixed points x_s where f is zero at x_s and points to x_s in the surrounding region. Additionally, the stable fixed points correspond to region(s) where U is locally minimum with respect to the fluid reservoir coordinates (x, y, z). Hence, to solve the inverse 3D ultrasound DSA problem and organize a user-specified pattern of particles at the set of desired positions X_des, we minimize the average value of U over all points xϵ X_des. Additionally, we constrain the amplitude of the N_t ultrasound transducers |v| = α, where α is a real scalar value, because in practice, they are limited to a finite input power. We use the formulation of Greenhall et al. to pose the inverse problem as a constrained optimization problem and calculate the optimal transducer operating parameters v^* that minimize the average value of U over the set of desired particle locations X_des for a fixed ultrasound transducer amplitude |v| = α, using eigendecomposition.¹² The 3D inverse solution technique presents a computational challenge compared to existing 2D ultrasound DSA techniques.^12,18,19,25 Modeling a 3D fluid reservoir using the boundary element method requires expanding the number of domain points from N_b = O(m²) (2D) to N_b = O(m³) (3D), where m is the number of uniformly spaced domain points in the x, y, and z directions. Additionally, the number of boundary elements increases from N_t = O(m) (2D) to N_t = O(m²) (3D), which significantly increases the processor time and memory required to solve the inverse problem via eigendecomposition.

Figure 2 shows a schematic of the experimental set-up we have used to validate our theoretical solution to the inverse ultrasound DSA problem. A poly(methyl methacrylate) fluid reservoir with interior dimensions L_x = L_y = L_z = 12.75 mm = 12.75λ₀ contains water (ρ_f c_f = 1.5 MRayl) with dispersed 80 nm carbon particles. Ultrasound transducers (lead zirconate titanate type SM111, center frequency ω₀/2π = 1.5 MHz, acoustic impedance Z_t = 35 MRayl²⁵) line three vertical walls and the floor of the fluid reservoir, creating a hard acoustic boundary at each ultrasound transducer surface (⁠$Z̃$ = 0.043). We note that such boundaries limit the patterns of particles that can be assembled if the operating frequency is a resonant frequency of the reservoir. However, this limitation can be removed by impedance matching of the ultrasound transducers to the fluid medium (Z_t = ρ_f c_f)²⁵ or by selecting a non-resonant operating frequency or reservoir size.¹² The ultrasound transducers are driven at their center frequency by a function generator, with independent control of their operating parameters v (amplitude and phase). A glass sheet on top of the fluid reservoir creates a reflector for the ultrasound waves generated by the transducer on the fluid reservoir floor and allows imaging of the patterns of particles from the top. One fluid reservoir wall is transparent to allow imaging the patterns of particles from the side.

We prescribe a user-specified pattern of particles in the domain D of the fluid reservoir and solve the inverse problem to compute the optimal ultrasound transducer operating parameters v^* required to assemble this pattern. We then apply the optimal ultrasound transducer operating parameters v^* to both the model and experiment to determine the resulting pattern of particles and compare both the simulated and experimentally obtained patterns to the originally defined user-specified pattern of particles. We prescribe patterns of particles assembled in vertical sheets, horizontal columns, and 3D dot arrangements. These patterns can be assembled in the cubic fluid reservoir of Fig. 2 by superimposing plane standing waves in the x, y, and z directions.

Figure 3 shows the simulated and experimental results obtained for user-specified patterns consisting of vertical sheets (Fig. 3(a)), horizontal columns (Fig. 3(b)), and 3D dots (Fig. 3(c)) of 80 nm carbon particles dispersed in water contained within the fluid reservoir shown in Fig. 2. Figures 3(a)–3(c) depict a schematic of the user-specified pattern in the center of the figure. For each user-specified pattern, we solve the inverse ultrasound DSA problem to calculate the optimal ultrasound transducer operating parameters, which we subsequently apply to the ultrasound transducers in both simulations and experiments. Table I shows the optimal operating parameters v₁, v₂, v₃, and v₄ of ultrasound transducers 1, 2, 3, and 4, respectively, for the vertical sheet, horizontal column, and 3D dot patterns. Additionally, we note that the optimal ultrasound transducer parameters for the vertical sheets and horizontal columns are comparable to those required to create 1D sheets in the x-direction¹⁷ and 2D dots in the x-z plane.²⁵ 

Figures 3(a)–3(c) show the simulated pattern of particles (gray) overlaying the corresponding radiation potential U (green) on the left and the experimentally obtained pattern of particles (black) on the right. The simulated patterns are computed as the regions around the stable fixed points x_s, where |f| is below ten percent of the maximum acoustic radiation force amplitude to ensure that the thickness of the experimentally obtained and simulated pattern features is equivalent. The inset images in Figs. 3(a)–3(c) display a top and side view of the simulated (left) and experimentally obtained (right) patterns with the user-specified patterns superimposed in red. We qualitatively observe good agreement between the user-specified patterns and the simulated and experimentally obtained patterns. Additionally, we compute the pattern error E_pat as the mean distance between the centroids of the user-specified pattern features and the centroids of the simulated and experimentally obtained pattern features, normalized by the user-specified pattern spacing λ₀/2, to quantify the differences between the user-specified patterns and the simulated and experimentally obtained patterns of particles. We compute E_pat of 1.2%, 1.8%, and 5.8% for the simulated patterns and E_pat of 16.6%, 17.7%, and 17.9% for the experimentally obtained patterns, indicating good agreement between the user-specified patterns and the simulated and experimentally obtained patterns, respectively. Differences between the user-specified patterns and the simulated and experimentally obtained patterns are due to imperfect dispersion of the carbon particles in the fluid medium in the experiments, near-field effects,⁶ secondary scattering of the ultrasound wave field between neighboring particles, and acoustic streaming, which are not included in the boundary element model of the fluid reservoir, but can disrupt the experimentally obtained patterns of particles.²⁶ The pattern error of the experimentally obtained patterns is also due in part to manufacturing imperfections of the cubic fluid reservoir and imperfect alignment of the ultrasound transducers with respect to each other. Because the complexity of aligning and independently controlling ultrasound transducers increases drastically with increasing N_t, we limit our experiments to a cubic fluid reservoir with N_t = 4 ultrasound transducers, which constrains the patterns of particles that can be achieved to 3D sheets, columns, or dot patterns. We note that sheet and column patterns have been experimentally demonstrated in 2D, under the assumption that the patterns are uniform in the third dimension.^12,25 However, by creating 3D patterns of vertical sheets (Fig. 3(a)) and horizontal columns (Fig. 3(b)), we show that the patterns are not uniform in the third dimension due to near-field effects, reflections from the reservoir walls, and other effects not accounted for in the 2D model. As a result, the simple 3D patterns of vertical sheets, horizontal columns, and dots provide proof-of-concept of the 3D inverse problem solution method by demonstrating accurate replication of user-specified patterns, as indicated by small pattern errors. We have also used simulations to demonstrate the ability of the 3D inverse solution method to create complex user-specified patterns within fluid reservoirs of any geometry, lined with a large number (N_t > 4) of ultrasound transducers in any arrangement. Figure 4 shows the results for a user-specified pattern of particles consisting of eight hollow spheres in a cubic fluid reservoir with N_t = 24 ultrasound transducers, and Fig. 5 shows an electron d_z2 orbital pattern of particles in a spherical fluid reservoir with N_t = 225 ultrasound transducers. Figures 4 and 5 show (a) the simulated fluid reservoir with the 3D user-specified pattern of particles indicated in red and (b) the simulated pattern of particles in gray obtained by solving the inverse problem for the user-specified pattern of (a) to compute the optimal ultrasound transducer operating parameters, which are subsequently employed to achieve the simulated pattern of particles. The inset images show 2D cross sections of the user-specified (red) and simulated (gray) patterns of particles overlaying the radiation potential (green). We observe good qualitative agreement between the user-specified and simulated patterns of particles. The results in Fig. 5 deviate primarily at the center of the pattern, where it contains sharp corners. When the density and compressibility of the particle and fluid medium are not identical (ρ_p ≠ ρ_m and β_p ≠ β_m), ultrasound DSA is limited to creating pattern features with non-zero radii. For instance, if ρ_p > ρ_m, β_p < β_m, and the standing spherical ultrasound wave field is expressed in terms of the Green's function from Eq. (2) as $φ(x)=Im[G(x0,x)]$ for an arbitrary point x₀, then the corresponding acoustic radiation force drives particles into a pattern of concentric spheres around x₀, where the smallest sphere has radius r = πc_f/ω₀. Thus, the smallest possible pattern feature radius is inversely proportional to the operating frequency ω₀ but will always remain greater than zero. Similarly, the 3D inverse solution technique is limited to creating user-specified patterns with pattern features separated by |Δx| ≥ λ₀/2 due to the minimum spacing between nodes or antinodes of the ultrasound wave field with wavelength λ₀. Thus, the 3D inverse solution technique fails to create user-specified pattern features with sub half-wavelength spacing (|Δx| < λ₀/2). Additionally, the 3D inverse solution technique may result in simulated pattern features that were not specified by the user (light gray regions in Figs. 4(b) and 5(b)) due to inevitable interference between the ultrasound waves generated by each ultrasound transducer and reflections from the fluid reservoir walls.¹² We observe that the feasibility of achieving a user-specified pattern is a function of the number of ultrasound transducers and their geometric arrangement, where the feasibility increases for an increasing number of transducers. We have defined the range of feasible patterns for the simple rectangular ultrasound transducer arrangement shown in Fig. 2. However, in general, defining the relationship between the range of feasible patterns and the arrangement and operating frequency of ultrasound transducers remains an open problem.

In conclusion, we have theoretically derived a direct method of solving the 3D inverse ultrasound directed self-assembly (DSA) problem to calculate the ultrasound transducer operating parameters required to create a 3D user-specified pattern of particles in a fluid reservoir with arbitrary simple, closed geometry, and ultrasound transducer arrangement. We demonstrate good quantitative agreement between the user-specified, simulated, and experimentally obtained patterns, including sheets, columns, and 3D dots of particles. We demonstrate with simulations that the solution method also applies to complex user-specified patterns of particles. This knowledge enables using 3D ultrasound directed self-assembly in a myriad of engineering applications, including biomedical and fabrication processes.

We acknowledge support from Army Research Office Contract No. W911NF-16-1-0457. JG also acknowledges support of the NASA Space Technology Research Fellowship, Award No. NNX15AP30H.