Acoustical boundary hologram for macroscopic rigid-body levitation

Acoustic radiation pressure is a time-average pressure applied to an object in the path of a propagating or standing sound wave, which can exert a force on the remote object. Such a non-contact force theoretically enables the manipulation of objects of various shapes in midair. However, although successful levitation has been reported in previous studies, these achievements were limited to levitation of objects with sizes comparable to (Marzo et al., 2018) or smaller than (Brandt, 2001; Foresti et al., 2013; Ochiai et al., 2014a; Whymark, 1975) the sound-source wavelengths λ, or to objects located near the sound sources.

There are two strategies for stable levitation of objects larger than λ in the far-field. The first strategy involves dynamic control of the sound-source signal according to the position, orientation, velocity, and angular velocity of the object in every moment measured by the employed sensors; we refer to this technique as feedback-based levitation. One problem with this strategy is that measurement completeness is crucial. Furthermore, a delay of 3 ms for a distance of 1 m between the object and sound source in air is inevitable, even if the measurement is perfect.

The second strategy involves the formation of a three-dimensional (3D) acoustic field, in which the object is stably levitated without time-domain control; this technique is referred to as stable potential levitation. If it were possible to generate a sound field that could exert an appropriate restoring force to counteract some minute arbitrary displacement of the object, stable levitation could be performed without sensor feedback. This strategy is feasible even considering the large propagation delay between the object and sound source. Further, such a quasistatic technique could be combined with the feedback-based levitation technique to achieve more flexible and stable object control. This approach could also be applied to acoustical levitation by using acoustic metamaterials (Marzo et al., 2017; Melde et al., 2016; Memoli et al., 2017). However, it is not at all obvious that a stable levitation field can be formed for any objects in reality. At present, it is known that such a restoring force field is formed for a particle within a standing wave, provided the particle is smaller than λ (Brandt, 2001; Foresti et al., 2013; Ochiai et al., 2014a; Whymark, 1975). Another research trend is to investigate the characteristics of acoustic helicoidal wavefronts, i.e., acoustical vortices (Baresch et al., 2013a, 2016; Marzo et al., 2018). It has been proven that techniques using acoustic vortices are capable of levitating a Mie sphere having a maximum diameter of 1.88λ, to the best of our knowledge. However, no examples of a field that produces an appropriate restoring torque as well as a suitable force for non-spherical objects larger than λ have been reported.

In this paper, based on the premise of phased array application, we present a numerical approach termed the boundary hologram method, in which a static levitation field around a free-shaped impedance boundary is designed. Hence, we show that numerically optimized fields can provide the desired restoring force and torque. Although it is not guaranteed that the necessary sound field can be generated for all combinations of objects and phased array geometries, we experimentally demonstrate that at least two examples, namely, a sphere with a diameter of 3.5 λ and a regular octahedron with a diagonal length of 5.9 λ can be stably levitated at a distance of 24 λ from the source without dynamic feedback control. Furthermore, the experimentally proven weight payload is 4–6 orders of magnitude larger than those of previously reported far-field acoustical levitation, while the acoustic intensity flux is equivalent (Baresch et al., 2016; Marzo et al., 2018, Marzo et al., 2015; Seah et al., 2014).

A typical particle levitator uses standing waves generated by a pair of transducers or a set of transducers and a reflector (Brandt, 2001; Foresti et al., 2013; Whymark, 1975). In a recent study, the manipulation degrees of freedom were expanded to three dimensions (Ochiai et al., 2014a). In addition, single-sided traps have recently been realized by extending the technique used in optical tweezers, first for lateral trapping (Lee et al., 2009) and then 3D trapping (Baresch et al., 2013a, 2016; Marzo et al., 2015; Ochiai et al., 2014b).

Historically, a trade-off between target size and levitation height (i.e., the distance from the wave sources) has been necessary. For example, Marzo et al. (2015) levitated spheres with diameters of 0.36 λ at a distance of 22 λ. A near field on an ultrasonic vibrator can lift an object larger than λ; this is called the squeeze-film effect. In an experiment by Ueha et al. (2000), a height of 0.006λ was obtained for an object size of 4.2 λ. Another approach is to create a standing wave between the wave sources and the target object itself (Andrade et al., 2016; Andrade et al., 2017; Zhao and Wallaschek, 2011). Using this strategy, Andrade et al. previously demonstrated the levitation of a sphere with a diameter of 3.6 λ at a distance of 0.5 λ. Further, Foresti et al. (2013) demonstrated the trapping of a toothpick with a length of 6 λ between wave sources and a reflector at 0.5 λ from the wave sources.

Very recently, Marzo et al. (2018) introduced a virtual vortex for Mie particle levitation. In that study, the virtual vortex was defined as a periodic multiplex of a chiral combination of acoustic vortices, and independent control of the lift force and orbital angular momentum was achieved. It was proven that the virtual vortex could levitate spheres with diameters up to 1.88 λ at 12 λ from the wave sources. Furthermore, Cox et al. (2018) demonstrated non-spherical levitation up to 0.4 λ using the multiplex technique of a twin-trap and a standing wave.

The boundary element method (BEM) has historically been a powerful analytical tool; for example, for the analysis of exterior scattering problems. The essence of the concept presented herein is the active drawing of acoustical imagery on the boundary manifold by controlling the incident wave field. Then, an explicit formulation between the incident wave generator and resultant radiation force is derived. The boundary integral equation is introduced from the Helmholtz equation by restricting its region of interest to the surface manifold, on which we design exactly. As detailed in Sec. III A, the acoustic radiation force applied to the body is also expressed as an integral of the acoustic pressure and particle velocity on the body surface.

The force F that acts on a rigid-body surface RB due to acoustic radiation pressure which is expressed by the Lagrangian density L and the momentum density tensor $ρ0uu$

where ρ₀ is the density of the medium, $u$ is the particle velocity, and $⟨·⟩$ denotes the time domain average (Westervelt, 1951). Using the surface normal admittance β and assuming a harmonic time-domain term, we obtain the following relations for u, the sound pressure p, and unit surface normal n

where j denotes an imaginary unit, c₀ is the speed of sound in the medium, and k is the wavenumber of the ultrasonic wave. We can rewrite the radiation force acting on the rigid body as

On a sound-hard boundary for which β = 0, this is simply

where $∇∥p$ indicates the tangential derivative on the surface. Similarly, the torque T is given as

where $r$ indicates the surface element position and $c$ is the center of mass of the rigid body (Maidanik, 1958; Zhang and Marston, 2011). Note that T is always zero and therefore uncontrollable for a true sphere in this formulation.

A phased array transducer (PAT) is an array of acoustic transducers for which the phases are configurable. The pressure of an incident wave at a point on the surface element of the target object $r$ is obtained as the integral of the wavefields generated by the phased array elements

where $g(r,s)$ is Green's function for the Helmholtz equation, $q(s)$ is the normal particle velocity on the phased array surface, and PAT denotes the PAT surface.

A detailed derivation of the boundary integral formulation and the following Galerkin method is omitted here. Interested readers are referred to Sauter and Christoph (2011). Briefly, the acoustic pressure $p(r)$ on a sufficiently smooth RB exerted by the incident wave is expressed as follows:

Specifically, the assumption of a sound-hard boundary yields the following:

where $δ(r,r′)$ is the Dirac delta function.

To handle these integral equations for non-algebraic shapes, we insert an ansatz that the density functions $p(r)$ and $q(s)$ are on the finite-dimensional span of a set of basis functions $ϕ$ and $ψ$⁠; i.e., $p(r)=∑ipiϕi(r)$ and $q(s)=∑iqiψi(s)$⁠. Applying the Galerkin method, we obtain the following system matrices and matrix equation:

where $p$ and $q$ are vectors with coefficients p_i and q_i. This equation gives a link between the imagery on the PAT $q$ and that on the body surface $p$⁠. For basis functions, we employ piecewise linear continuous functions for a triangular mesh; in other words, continuous boundary elements of degree 1 for $ϕ$ and piecewise constant functions for $ψ$ (see Sec. 4.1.7 of Sauter and Christoph, 2011).

Next, we discretize the integrals in $F$ and $T$⁠.

First, the above-mentioned ansatz is inserted into the integral of the squared density

Note that $νi,j$ is a 3D vector that corresponds to the normal vector. Let N_x, N_y, and N_z be matrices with coefficients $νxi,j,νyi,j$⁠, and $νzi,j$⁠. Then each component of Eq. (12) is expressed in matrix quadratic form. For example, its x component is

Similarly, $∫|p|2[(r−c)×n]xdS=p*Mxp$ is derived where $Mxi,j=∫ϕj*(r)ϕi(r)[(r−c)×n]xdS(r)$⁠.

Second, the tangential gradient of the pressure can be discretized as

Therefore, its norm is

and the v ∈ {x, y, z} component of its integral is also written in the quadratic form

Similarly, we define

Hence, using $p=B−1Gq$⁠, Eqs. (4) and (5) are described in a closed form of the transducer's complex gain. Their v ∈ {x, y, z} components are

Equations (18) and (19) are matrix quadratic forms with respect to phased array gain vectors, which can be efficiently handled on modern parallel processors.

The conditions for rigid-body levitation are as follows: (1) the external forces, which are typically gravity and the acoustic radiation force, must be balanced; (2) a restoring force should be generated against positional fluctuation; (3) the net moment must be zero; and (4) a restoring torque should be generated against rotational fluctuation. In previous studies on particles, condition (1) was automatically satisfied by a slight position shift of the particle, but conditions (3) and (4) were outside the scopes of those works. Marzo et al. (2015) satisfied condition (2) by maximizing the Gor'kov Laplacian to render it positive for particles. However, the Gor'kov approach assumes that the target is a small sphere (typically, R < 0.1 λ) in nature, and is inappropriate for application to larger and/or aspherical objects. Further, a positive Laplacian is not a sufficient condition for stabilizing a rigid body.

Let small positional and rotational perturbations from the target levitation point be expressed as $Δx=[Δx,Δy,Δz]T$ and $ΔΩ=[ΔΩx,ΔΩy,ΔΩz]T$⁠, respectively. We combine those two statements as a generalized position $X=[Δx,Δy,Δz,ΔΩx,ΔΩy,ΔΩz]T$ and force $F=[Fx,Fy,Fz,Tx,Ty,Tz]T$⁠.

For sufficiently small $X, F$ is expressed as a linear approximation using a Jacobian matrix $∇F$⁠, such that

Specifically, the sufficiently small region here is the region in which the first-order Taylor expansion is valid. That region is not precisely determined before the optimization introduced below is solved. We evaluate the size of the region in Sec. V E.

When all the forces and torques acting upon the object balance each other and resistance from the medium is ignored, the equation of motion becomes

where $M=diag(mI;I)$⁠, using m and $I$ as the mass and inertia matrix of the rigid body, respectively. Here, $X(t)$ is bounded for a bounded initial position $X(0)$⁠, if and only if all eigenvalues of $M−1∇F(0)$ are non-positive real numbers. This is known as the “Lyapunov stability condition” in control theory (Strogatz, 2001).

The Lyapunov stability is a weaker concept than asymptotic stability. A particle is called asymptotically stable when it is infinitely close to an equilibrium point as time expands. On the other hand, the Lyapunov stability only requires that the particle moves inside a neighborhood of an equilibrium point as time expands. In microscopic particle levitation, a lossy factor generated by the fluid viscosity contributes to the asymptotic stability. However, in this macroscopic scenario with a large Reynolds number, the contribution is small and we omitted the term from Eq. (21). Its explicit consideration requires an additional hydrodynamics model.

Next, we formulate an optimization problem to synthesize a phased array output that minimizes the eigenvalues while maintaining the force-torque equilibrium condition. In other words, an acoustic hologram image on the surface manifold is optimized so that the field can properly respond to the fluctuations of the body position and rotation. In terms of the objective function to be minimized, we employ a linear combination of the generalized force mean-squared error and the real parts of the eigenvalues, which is expressed as follows:

where $w$ and $v$ are positive weight hyper-parameters, $Fext$ represents the external force (e.g., gravity) and torque applied to the rigid body at the levitation point, and $λ=eig(M−1∇F(0))$ represents the eigenvalues. The hyper-parameters $w$ and $v$ are employed to tune the importance of supporting and restoring force, respectively. Similar parameters have been employed in Marzo et al. (2015). Note that the second term is replaced by the trace of the Jacobian in the case of $w=1$⁠. Although there is no explicit term to render the imaginary parts of the eigenvalues zero in this formulation, we expect that poles are placed on the real axis spontaneously.

This relaxed optimization problem is solvable using gradient methods such as the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method. The existence of a solution is not always guaranteed. That is, a solution satisfying the constraints cannot be found when the conditions are physically unfeasible, i.e., if the target is too heavy or its shape is too complex for the given phased array. Further, as this is a non-convex problem, we cannot always expect to obtain the global minimum.

For the experiments conducted in this study, we assumed that the body boundary was sound-hard. Note that this assumption is valid for most natural rigid materials in gaseous media. We implemented phased arrays comprising 996 air-coupled transducers (Nippon Ceramic T4010; hereafter “996 array”). Each transducer was driven by a 40 kHz rectangular wave of V_pk–pk = 24 V and $⟨V⟩=0 V$⁠. See Inoue et al. (2016) for details of the array system and Ito et al. (2016) for its frequency variations. In the experiments, the phased arrays were placed on the horizontal plane. For levitation targets, we fabricated a sphere with a diameter of 30 mm (mass: 6 × 10²mg) and a regular octahedron with a diagonal length of 50 mm (mass: 5 × 10²mg) made of expanded polystyrene. The geometry sketch in the case of the octahedron is shown in Fig. 1.

The sphere was levitated at 200 mm (24 λ) from the center of the array by both a single-beam wave (generated by a single-sided 996 array placed at the bottom only) and a standing wave (generated by bottom and top phased arrays, which were 400 mm apart and consisted of 1992 transducers in total). The latter was more stable. In the case of the octahedron, a stable solution was achieved by the optimizer for the double-sided phased array, but not for the single-sided phased array. In actuality, in our experiments, the double-sided array levitated the octahedron and the single-sided array did not.

Figures 2 and 3 show the levitation of the sphere and the regular octahedron by a 40 kHz (λ = 8.5 mm) ultrasound in air at a position 200 mm above the bottom phased array. In the figures, the three green laser lines form an orthogonal basis crossing at the levitation point (see also MM. 1 and MM. 2 for full videos of the levitations).¹

We did not succeed in levitating octahedrons with the single-sided 996 array regardless of the object size. One possible reason for this difficulty is that the acoustic complexity may not be sufficient to control the six-degrees-of-freedom (force and torque) stability. This hypothesis is supported by the fact that the nominal rank of matrix B⁻¹G for the double-sided array is more than twice that of the single-sided array.

In Fig. 4, optimized phase patterns of the sphere and the octahedron levitators are shown. The patterns are similar to that of a helical vortex field, which has been proven to stabilize particles (Baresch et al., 2013a,b, 2016; Marzo et al., 2018, Marzo et al., 2015; Shvedov et al., 2014). It is interesting that the present optimization technique yields a similar pattern without heuristics. On the other hand, the pattern is elliptic, unlike that for helical vortices. This is a unique characteristic of our pattern and suggests that symmetry is not always necessary for levitation.

The hyper parameters were $w=(1,1,1,0,0,0),v=(1,1,1,0,0,0)$ for the sphere and $w=(1,1,10,1,1,1),v=(1,1,1,1,1,1)$ for the regular octahedron, which were empirically sought from the all-ones vector. Note that the latter half of the components of the parameter vectors for sphere levitation are fixed to zero because they are for torque control.

Figure 5 shows a finite element method (FEM) simulation of the absolute acoustic pressure of the sphere and the octahedron levitation fields on the vertical plane. Both the acoustic field in the absence and presence of a target rigid body are shown. The acoustic fields were significantly changed by the rigid bodies and showed that reflections occurred on the rigid body surface.

Figure 6 shows a comparison of the simulated and measured absolute acoustic pressures on the horizontal and vertical planes of the sphere levitation field. It is apparent that the optimized acoustic field surrounded the target object. Note that the measurement and simulation were conducted in the absence of the sphere.

These FEM simulations were implemented on COMSOL Multiphysics 5.1b. The acoustic field was measured using a calibrated microphone (TYPE 4138-A-015, Brüel & Kær) attached to a 3D robot arm (M-710iC, FANUC). The measurement was performed under 20% voltage and rescaled because of the saturation of the microphone.

A preliminary experiment indicated that the radiation force predicted by Eq. (18) for the sphere and the octahedron almost converged to a certain value with an error rate ΔF/F < 3% when a mean mesh edge size of less than 2 mm was used. Figure 7 shows the acoustic radiation force exerted by a single focal point targeting the origin. The 996 array was used and the height from the array was 200 mm. The sphere and octahedron were moved along the horizontal x axis while the focal point remained fixed. The mean mesh edge size was 2 mm. The force acting on the sphere was overestimated, but that on the octahedron was relatively well matched with the prediction. A possible reason for this is the completeness of the curvature approximation. As octahedrons consist of planes, meshing of a octahedron perfectly models its normal vector, whereas meshing of a sphere does not. Meshing of spheres with shorter edge sizes improved the prediction slightly, but at the expense of high memory and time costs to construct matrices B and G. We concluded that a mean mesh edge size of 2 mm is a practical compromise. Hereafter, triangular meshes with a mean mesh edge size of 1 mm and the Gauss quadrature rule of order four were employed to ensure sufficient precision.

Figures 8 and 9 show the estimated and measured restoring force against the displacement around the equilibrium point. In detail, Figs. 8(a), 8(b), and 9(b) show the horizontal forces for horizontal displacement applied to the sphere and octahedron. The restoring force acted appropriately in a direction opposite to the displacement and converged at the origin. Figure 9(a) shows the net resultant force, which includes the restoring force and the gravitational force, applied to the sphere in the vertical plane, including the origin. Again, the force converged at the origin, and the gradient of the force was observed along the vertical line. Errors between the estimation and the measurement were observed here. A possible reason is that the objects made of expanded polystyrene are not perfectly sound-hard. The boundary between the air and expanded polystyrene has some finite impedance but our estimation model considers it infinite. Figure 8(c) shows the torque on the octahedron in response to the rotational displacement. It is apparent that an appropriate restoring torque is applied for rotational noise. Unfortunately, as we cannot build a mN·mm-order torque sensor without disturbing the acoustic field, only a numerical analysis is presented here.

The stable region of the rigid body is the neighborhood of the levitation point. Here, the restoring force acts like a linear spring, as our optimizer constrains the gradient of the restoring force at the equilibrium point (see Fig. 8). To observe this, we tracked the position of the sphere from various initial points. Figure 10 shows the time to ejection with respect to the initial object position. The initial positions were distributed three-dimensionally, and their variance-covariance matrix was

The spheres that dropped within 4 s were initially at positions greater than 5 mm from the origin. The stable region of radius 5 mm corresponded to the linear region in Fig. 8(a).

The field generated a force on the order of millinewtons, while the peak acoustic intensity was on the order of 10 W/cm². The generated forces were 4–6 orders of magnitude larger than those of the previously reported far-field acoustical levitation, while the acoustic intensity flux is equivalent to those in previous reports (Baresch et al., 2016; Marzo et al., 2018, Marzo et al., 2015; Seah et al., 2014).

The force field was measured by using a load cell (LTS-50GA Kyowa Electronic Instruments, rated capacity= 500 mN, nonlineality <0.1 mN, hysteresis <0.65 mN) and an amplifier (WGA-900A, Kyowa Electronic Instruments) mounted on a 3D robot arm (M-710iC, FANUC). The load cell and objects were connected via a stainless steel rod (diameter = 2.0 mm, length = 170 mm) to reduce the intrusiveness of the field.

The optimization given in Eq. (22) is a non-convex problem; thus, global minimization is not always guaranteed. Successful use of the limited-memory (L)-BFGS algorithm for large-scale applications including acoustical levitation has been reported (Liu and Nocedal, 1989; Marzo et al., 2015). In addition, the basin-hopping algorithm, which is a randomized technique designed to search for the global minimum, has also been used in large-scale applications (Wales and Doye, 1997).

We compared the BFGS and its combination with the basin-hopping algorithm for this optimization. Figure 11 shows the sequences of the largest eigenvalue in Eq. (22) with respect to the step count. Convergence to slightly different solutions was observed. Note that the value does not decrease monotonically as the vertical axis does not show the value of the objective function, but rather the largest eigenvalue. However, the levitation quality was not improved in this case. The wall times required for a single step of our SciPy-based implementation were 1.92 and 94.8 s for BFGS and basin-hopping, respectively, for the 996 array. We conclude that the BFGS converges to local minima but is practically useful in this application.

In the experiments throughout this work, we used a focus at the target center as the initial PAT phase pattern of the optimization. We experimented with various initial focuses aimed at other positions and obtained the same solution, provided the focal point was inside the target object.

The dimension of this optimization problem is the number of transducers in the array. Therefore, the time and space complexity is not affected by the number of meshes on the target object surface. Further, if we consider manipulating an object by using this technique, the matrix G requires recalculation once the object moves, whereas $B,N,N$⁠, and $M$ do not. The size of G corresponds to the (number of mesh vertices) × (number of transducers), whereas the other matrices have sizes corresponding to the (number of mesh vertices)². Thus, greater efficiency is achieved if the number of transducers is as small as possible, although this reduces the field diversity.

We demonstrated acoustical levitation for macroscopic rigid bodies. To achieve this, we employed a BEM formulation for acoustic radiation pressure applied to an arbitrary smooth surface, and formulated an optimization problem that maximizes the stability while balancing all forces and torques at a point sufficiently distant from the acoustic elements.

Even maintaining the acoustic intensity flux, the weight payload is four orders of magnitude larger than that of previously reported for far-field acoustical levitation. Two main factors contribute to this achievement. First, our method simultaneously models the supporting force and the restoring force that act on a large-surface object. Because this approach explicitly describes the object surfaces, it can target the exact surfaces that effectively contribute to levitation. Second, the model can calculate the radiation force for a wide class of acoustic fields; that is, any acoustic fields generated by discrete phased array transducers. In contrast, previous analytical studies have considered a relatively narrow class of acoustic fields, such as Bessel beams or acoustic vortices. Thus, our large-scale optimization is capable of finding a better solution for levitation.

Perfect asymptotic stability was not achieved. In microscopic particle levitation, the lossy factor is generated by the fluid viscosity and renders the trap tolerant to wind or other disturbances. However, in a macroscopic scenario with a large Reynolds number, the lossy factor from the viscosity is small. Considering some other damping factor in the system equation [Eq. (21)] will contribute to additional stability.

Throughout this paper, we assumed rigid bodies as targets. As the BEM formulation is concerned with the boundary only, deformation was not modeled. In further work, trapping of deformable objects, such as droplets, will expand the range of important applications (Ding et al., 2012; Laurell et al., 2007). BEM-FEM coupling is one possible approach here. Another is to add a penalty term to limit the local surface acoustic intensity in order to prevent object deformation.

The method can be applied without assuming any prior special conditions regarding the object shape or the configuration of the phased array. Even if a stable field cannot be obtained from numerical optimization, it is expected that a stable field can be achieved by expanding the geometry of the phased array as necessary.

The developed levitation technique is applicable to the 3D manipulation of multiple objects to make them follow complex 3D paths in air or liquid, which would be impossible to achieve using mechanical arms. While demonstrations of dynamic manipulation are beyond the scope of this study, the main advantage of this technique is the explicit controllability of the force and torque instead of indirect control by pointing to the destination location of the next step (Baresch et al., 2016; Marzo et al., 2015). Thus, this approach can exploit the systematic dynamical model. Maximizing push or pull force in an arbitrary direction is another interesting topic.

Another aspect of the presented method is that it enables the use of higher-frequency (shorter-wavelength) waves than those previously employed for microparticle manipulation, which improves the temporal response and strengthens the holding force generated by the steeper pressure gradient. This boundary-holographic-optimization approach is also applicable to electromagnetic fields, allowing performance improvement for optical tweezers (Grier, 2003; Shvedov et al., 2014).

We thank Professor S. Hara, Chuo University for the suggestions. This work was supported in part by a JSPS Grant-in-Aid for Scientific Research (S) 16H06303, JST CREST JPMJCR18A2, and a JSPS Grant-in-Aid for JSPS Fellows 15J09604.