Acoustic levitation of a Mie sphere using a 2D transducer array

Most acoustic levitation techniques are limited to objects smaller than half the wavelength. To overcome this limit, different strategies have been proposed for suspending macroscopic objects in mid-air. Two approaches to levitate spherical and non-spherical macroscopic objects have been recently presented: the acoustical virtual vortices and the boundary hologram method. However, the former approach places high demands on the available hardware due to the mandatory high switching rate while the latter uses a computationally expensive model that prevents future real-time manipulation. In the present work, we demonstrate the single-beam levitation of a Mie sphere using a 2D transducer array. To achieve this, we employ a computationally fast sound field model based on spherical harmonics expansion. To obtain a suitable array output, we formulate an optimization problem that maximizes the stability of the sphere while keeping the net force balanced. In addition, we prove the local asymptotic stability for the equilibrium position and determine a domain of attraction using Lyapunov-based methods. In experiments, we show that the macroscopic sphere is stably levitated in a twin tuning forks trap, which results from a superposition of two twin trap signatures and a bottle trap signature. This result could open up the possibility of a computationally fast and convenient non-contact manipulation of macroscopic objects by a superposition of holographic elements in future applications.


I. INTRODUCTION
When a sound wave impinges on the surface of an object, it exerts an acoustic radiation pressure on it. 1In general, the radiation pressure is weak, but if the wave intensity is increased considerably, the resulting radiation force can be strong enough to counteract the gravitational force and levitate small particles in mid-air.3][4][5] One of the most interesting features of acoustic levitation is that it can be used to trap a wide variety of materials, such as solids, 6 liquids, [7][8][9] soap bubbles, 10 and small living creatures. 11,12coustic levitation has also numerous potential applications in biology, [13][14][15] chemistry, 16 pharmacy, 17,18 and microassembly. 19,20ifferent strategies have been proposed for suspending objects.The most common approach traps objects much smaller than the acoustic wavelength at the pressure nodes of a standing wave field, which can be generated either by a device consisting of a transducer and opposing reflector, [21][22][23] between two opposing transducers 24 or arrays of transducers. 25,26Another approach uses ultrasonic beams generated by a single-sided array of transducers to trap small objects in mid-air. 27In contrast to devices based on standing waves, acoustic levitation based on single beams does not require the object to be confined between opposing acoustic elements.
However, most acoustic levitation techniques are limited to objects that are smaller than half the wavelength.In order to overcome this limit, different strategies have been proposed for levitating objects larger than the acoustic wavelength.For instance, the near-field acoustic levitation method 28,29 was used to levitate large planar objects closely above a surface that is vibrating at an ultrasonic frequency.However, the maximum levitation height for this technique is on the order of tens of micrometers.1][32] Another approach for levitating large spherical objects uses an array of transducers for generating virtual vortices of high aperture. 33Recently, Inoue et al. 34 have proposed the boundary hologram method to levitate large non-spherical objects.They combined an optimization algorithm with the boundary element method to find the optimal phase angles of an array of transducers.This method provides not only a restoring force to trap the object in position but also a restoring torque to trap the object, e.g., a regular octahedron made of polystyrene, in orientation.However, the boundary hologram method uses a computationally expensive model that prevents future real-time manipulation of macroscopic objects while the method of acoustic virtual vortices places high demands on the available hardware due to the mandatory high switching rate.
In the present work, we demonstrate the acoustic levitation of a Mie sphere, i.e., the sphere size and the wavelength have the same order, using a 2D array of ultrasonic transducers.A numerical approach based on spherical harmonics expansion 35 simulates the acoustic radiation force acting on the sphere, whereas an optimization-based algorithm is used to find the emission phase of each transducer.The method offers a remarkable and potentially real-time alternative for levitating and manipulating macroscopic objects in mid-air.This article is structured in six sections.After an introduction to the sound field model (Sec.II A) and the calculation of the acoustic radiation force (Sec.II B), we present a dynamic model to describe the translational motion of the sphere (Sec.III A).Subsequently, we state an optimization-based approach to determine a suitable array output for a static levitation (Sec.III B), followed by a detailed analysis of the stability of the equilibrium and its domain of attraction (DoA) using Lyapunov-based methods (Sec.III C).After a presentation of our hardware used for the experiments (Sec.IV), we discuss our results (Sec.V).Finally, we give an outlook on further research that is linked with the present work (Sec.VI).

II. ACOUSTIC RADIATION FORCE ON A SPHERE GENERATED BY A TWO-DIMENSIONAL ARRAY
For the calculation of the acoustic radiation force we employ, with minor modifications for the transducer model, the approach of Andersson and Ahrens. 35This approach based on the work of Sapozhnikov and Bailey 36 is tailored to phased arrays and offers notable advantages for numerical optimization, see also Secs.III B and V.

A. Sound field generated by the phased array
The incident sound field p i that affects an object exposed to it can be calculated by a spherical harmonics expansion as 37 where j n (kr) are the nth order spherical Bessel functions of the first kind and the wave number is given as k ¼ ω=c, where c denotes the speed of sound and ω is the angular frequency.The position vector r ¼ xe x þ ye y þ ze z from the center of the levitated sphere to a point in Cartesian space can be expressed in spherical coordinates (r, θ, w) with the relations where where the associated Legendre polynomials P m n (x) are related to the Legendre polynomials P n (x) by the formula 37 Similar to Eq. ( 1), the scattered sound field is given as where h n (kr) ¼ j n (kr) þ iy n (kr) are the nth order spherical Hankel functions of the first kind and y n (kr) are the spherical Bessel functions of the second kind.Since the incident and scattered waves are calculated in the frequency domain, the model does not consider time delays when the sphere moves with a finite velocity such as that considered by Rudnick and Barmatz. 38For compressible spheres, the coefficients Ŝm n ¼ c n S m n are 39 where a is the radius of the sphere, e Z ¼ (ρ p c p )= ρ 0 c 0 ð Þ is the relative impedance, and the prime 0 indicates the derivative of a function with respect to its argument.In addition, the subscripts 0 and p of the speed of sound c and the density ρ denote the medium and the material of the sphere, respectively.For soundhard or sound-soft surfaces, the corresponding formulas for the coefficients Ŝm n are given in Gumerov and Duraiswami (p.146). 40ndersson and Ahrens 35 have shown that the expansion coefficients S m n can be written as sum of coefficients j S m n from each transducer as if the sound pressure field is created by a discrete transducer array.In Eq. ( 7), r 0 is the radius of the spherical integration surface Γ and

Journal of Applied Physics
the superscript * indicates the complex conjugated element.For the elements of the transducer array, we employ a single-frequency piston source model in the far-field that is given as where r is the position in which the pressure is calculated and r j denotes the position of transducer j.Its source strength q j ¼ A j e if j is determined by the amplitude A j and the phase angle f j .The directivity function D f (ξ j ) of the transducer is taken as where J 1 is the first-order Bessel function of the first kind, r p is the piston radius, and ξ j denotes the angle between the transducer normal n j and the vector r À r j .Since the expansion coefficients j S m n are known for point sources, we adapt the solution by Williams, Eq. (8.22) 37 to Eqs. ( 8) and ( 9) and employ for the region r , r j as the solution of Eq. ( 7) for the computation of the expansion coefficients j S m n of the sound field model for each transducer.

B. Acoustic radiation force on a sphere
To calculate the acoustic radiation force F ac (r) ¼ F x e x þF y e y þ F z e z for the model of the sound pressure field in Sec.II A, we utilize the equations obtained by Andersson and Ahrens 35 that are given as where < Á f g and = Á f g denote the real and imaginary parts of the argument, respectively.In addition, the coefficients n , and H m n are given as For the calculation of the spatial derivatives, one can use either the equations given in Gumerov and Duraiswami (Chap.2) 40 or employ a numerical approximation as for small values h, where e a denotes the unit vector of the Cartesian axis a [ x, y, z f g.Regarding the optimization in Sec.III B, the derivatives of the radiation force with respect to the phase angles f can be obtained by using For implementation purposes, the infinite sums in Eqs. ( 11)- (13)  have to be truncated to a certain order N. Its choice involves a trade-off between the accuracy of the calculation and the computational effort, which particularly plays a role in future real-time applications.To ensure an upper limit of the absolute error ϵ of the acoustic radiation force F ac (r), we employ the formula by Gumerov and Duraiswami (p.432) 40 for the approximation of order N as where σ ¼ d=a, d is the radial distance from the center of the sphere, and the threshold value ka ð Þ * is given by Additional formulas for the approximation of N are given in studies of Andersson and Ahrens 35 and Xu et al. 41

III. DYNAMIC MODEL, OPTIMIZATION APPROACH, AND STABILITY ANALYSIS A. Dynamic model
The translational movement of a sphere in an acoustic field can be described by the nonlinear state space model

Journal of Applied Physics
where the state vector x ¼ r `v1 `contains the position r and velocity v ¼ v x e x þ v y e y þ v z e z of the sphere, F ac (r) is the exerted acoustic radiation force, and F g ¼ Àmge z the gravitational force that acts on the sphere, whose mass is m ¼ 4 3 πρ p a 3 and the mass matrix we employ the formula used by Fushimi et al. 42 as where , where μ .0 denotes the kinematic viscosity of the fluid.Since typical Reynolds numbers for our experiments ranges from 0 to 100, a nonlinear damping model was chosen instead of a linear model based on Stokes law.The system in Eq. ( 23) has an equilibrium at which has an equilibrium point at the origin b x ¼ 0.

B. Optimization approach
For a successful levitation of the sphere at a given position r e , the net force as well as the net moment must be zero.In addition, restoring forces and torques should be generated to deal with translational and rotational perturbations of the macroscopic object. 34nfortunately, the fast method of holographic acoustic elements by Marzo et al. 27 cannot be used to calculate a suitable array output because it is only applicable to Rayleigh spheres, i.e., the particle size is much smaller than the acoustic wavelength.Moreover, a singular value decomposition (SVD) recently applied by Helander et al. 43 is also not feasible since there is no linear relation between the source strengths q j and the acoustic radiation force.Hence, the array output has to be determined by solving a nonlinear optimization problem, where usually the amplitude A j of each transducer is kept constant and the phase angles f j are solely used as optimization variables. 27,41n order to create more complex sound fields, the simultaneous optimization of the amplitudes A j and phase angle f j as well as the transducer frequency f j was investigated in the studies of Andersson and Ahrens 44 and Puranen et al. 45 However, the choice of optimization variables always depends on the available hardware since for many arrays one can only adjust the phase of the transducers, see also Table 1 in the work of Zehnter and Ament 46 for a brief overview.Furthermore, additional variables increase the complexity of the optimization problem which has a negative impact on real-time applications.Since an exclusive optimization of the phase angle vector f is sufficient in our case, we formulate an optimization problem similar to Inoue et al., 34 as we employ the objective function The function consists of two terms: the weighted mean squared error between the acoustic radiation force F(f) and an external force F e for each Cartesian axis C ¼ x, y, z f gensures that the latter acts on the sphere at a given position r e .For static levitation, typically F e ¼ ÀF g ¼ mge z is chosen.In the context of dynamic applications, F e can also be provided by a higher-level open or closed-loop controller to implement a motion request along a desired trajectory.Examples of this hierarchical controller structure can be found in the studies of Zemánek et al. 47 and Matouš et al. 48Regarding the second term, similar to the studies of Inoue et al. 34 and Xu et al., 41 we suppose that the equilibrium position r e (b r ¼ 0) is located in a small region B # R 3 where the firstorder expansion of the Taylor series is valid.Consequently, the translational motion of the sphere for b x [ B can be approximated by the linear system that results from the nonlinear model in Eq. ( 25) by a linearization at the equilibrium point b x ¼ 0 near which the assumption b v % 0 is valid.Thus, the linear dynamics in Eq. ( 27) is mainly determined by its subsystem By minimizing the real parts of λ(f) ¼ eig A(f) f g, a placement of the eigenvalues λ(f) on the negative real axis is enforced to achieve local asymptotic stability at the position r e .In addition to stability, the term < λ i (f) f gin Eq. ( 26) plays an important role for the generation of appropriate restoring forces to deal with translational perturbations of the sphere.
Similar to the study of Inoue et al., 34 we employ positive hyper-parameters w i and v i in the objective function in Eq. ( 26) to tune the importance of supporting and restoring forces, respectively.A careful adjustment of these weights is essential since the acoustic radiation force F(f) is typically in the range of several μN to mN, whereas the eigenvalues λ(f) are scaled by the term (1=m), see Eq. (28).It is also noteworthy that these two terms represent contradictory objectives.Thus, a high overweight of one term often leads to the case that the other term is weakly regarded during optimization to achieve an overall better minimum of the objective function J(f).Hence, we have observed in experiments that these results often lead to deviating or even unstable equilibria.Therefore, we presume that for each Cartesian axis i [ x, y, z f g, a balanced weighting of the corresponding terms in Eq. ( 26), ensured by each pair (w i , v i ) of hyper-parameters, strongly favors the quality of the obtained optimization results.Besides the absolute difference of both terms, a good indication for the choice of each pair (w i , v i ) is also provided by the selected upper limit of the absolute error ϵ of F(f) in Eq. (21).
To solve the unconstrained optimization problem in Eq. ( 26), we employ the gradient-based Broyden-Fletcher-Goldfarb-Shanno (BFGS) method because it has been successfully applied to similar problems in literature studies. 27,34,44Since the problem is nonconvex, it cannot be guaranteed that the method always finds the global minimum.It is also possible that no feasible solutions exists.This may be the case if the sphere is too heavy for the array or should be stabilized at a certain position for which no suitable array output can be found. 34

C. Stability analysis
In Sec.III B, the local asymptotic stability at the position r e was achieved by minimizing the real parts of the eigenvalues λ(f) ¼ eig A(f) f gin Eq. ( 26).However, assuming < λ Ai f g , 0 8i, the eigenvalues λ As ¼ eig A s f g of the matrix A s in Eq. ( 27) are not hyperbolic since < λ As,i f g¼ 0 8i.Therefore, unlike the position r e , no conclusions about the local stability of the origin b x ¼ 0, i.e., x e ¼ r è 0 1Â3 Â Ã `, of the dynamic model in Eq. ( 25) can be made. 49This result is caused by the fact that, compared to the linear viscous damping assumed in the studies of da Silva et al., 50 we use a nonlinear friction model in Eq. ( 24) whose terms vanish during linearization at b x ¼ 0, and, therefore, they do not contribute to the stability of the equilibrium in the linearized model.Thus, the successful levitation that was observed in the experiments (see Sec. V) has to be proven in a different way.To achieve this, we use Lyapunov-based methods, which are common in nonlinear control theory.In particular, we utilize the invariance principle of Krasovskii-LaSalle. 49heorem 1 (Invariance principle).Let X be a compact, positive invariant set for b _ x ¼ g(b and V : Let Y be the set of all points in X, where _ V(b x) ¼ 0. Let Z be the largest invariant set in Y.Then, the following applies and suppose that no solution can stay identically in S, other than the trivial solution b x(t) ; 0. Then, the origin is asymptotically stable.
Based on the theorem and the corollary, the following steps are necessary to prove the local asymptotic stability of b x ¼ 0: first, we need to define a domain Ω with b x ¼ 0 [ Ω and then find a function V : Ω !R that is continuously differentiable and satisfies V(0) ¼ 0, V(b x) .0 for b x = 0, and _ V(b x) 0 in Ω.Finally, we have to show that in S no solution except the trivial solution b x(t) ; 0 can remain.As it can be seen from Figs. 4(a)-4(c), the components of the acoustic radiation force F ac (b r) are Lipschitz and fulfill where the two other Cartesian axes b, c [ C with b, c = a in b r a kept constant at the equilibrium position.The components of the drag force F d (b v) are also Lipschitz and fulfill we conclude that _ V 1 (b x) is negative semi-definite.Since the implication chain for the model in Eq. ( 25) for _ where e F(b r) ¼ F ac (b r) þ F g , it becomes clear that only the trivial solution b x(t) ; 0 remains in the set : By applying the invariance principle, we conclude that the equilibrium b x ¼ 0 is locally asymptotically stable.
In addition to stability, the domain of attraction (DoA) is of particular interest since it can be regarded as a quality measure for the stability of an equilibrium x e .The larger Ω D , the more likely it is that the system will return to its equilibrium after a deflection that may be caused by an external disturbance.In terms of our application, this means that high imprinted restoring forces that act like a spring (see also Fig. 4) ensure a robust static levitation and allow, for example, a higher initial misplacement of the sphere inside the acoustic trap.Apart from well-researched system classes, e.g., robotics, or systems of low complexity, e.g., a pendulum, the DoA for an equilibrium of a nonlinear system can usually only be approximated by a set Ω E , Ω D .To calculate Ω E from a given dynamic model, several simulation-based methods have been proposed, which can be classified into Lyaponov-based methods, e.g., sum of squares programming, and non-Lyapunov-based methods, e.g., trajectory reversing. 51Among these methods, we decided to use Lyapunovbased sampling.Although this method allows only a conservative estimation of Ω E , 52 the existence of a Lyapunov function can be guaranteed, 49 whose search is often a time-consuming process.
Furthermore, it is usually significantly less computationally intensive than methods based on techniques like trajectory reversing.Since < λ A,i f g , 0 8i applies (see Sec. III B), we can state a Lyapunov function for the nonlinear subsystem b _ where P is determined by the Lyapunov equation for a positive definite matrix Q [ R 3Â3 , e.g., Q ¼ I 3Â3 , where P has to be positive definite, i.e., P .0. It shall be noted that it is always possible to find a valid solution for P .0 when < λ A,i f g , 0 8i is satisfied. 49Since the additional differentiation has no influence regarding the DoA Ω E , we assume that _ V 2 (b r) can be calculated by Using Eqs. ( 35) and ( 37), the DoA Ω E can be expressed as the largest compact, positive invariant set, where the constant c [ R can be determined via numerical methods. 52The application of this technique for the determination of Ω E is shown in Sec.V.

IV. EXPERIMENTAL SETUP
To carry out the experiment, we utilize a 16 Â 16 array of 40 kHz ultrasonic transducers (Manorshi MSO-P1040H07T) of 9:8 mm in diameter to levitate a polystyrene sphere with a diameter of d ¼ 8:9 mm and a density of ρ p ¼ 15 kg m À1 .The sphere can be regarded as Mie sphere because its size and the wavelength λ ¼ c f ¼ 343 ms À1 40 kHz % 8:6 mm are in the same order.The transducers are driven by square wave signals generated by a Field-programmable gate array (FPGA) (Altera Cyclone IV-EP4CE6).Shift registers (Texas Instruments 74HC595) are used for converting 32 outputs of the FPGA into 256 independent signals, which are amplified up to 20 Vpp by MOSFET drivers (Microchip MIC4127).Although the transducers are driven by square wave signals, the emitted acoustic wave is sinusoidal due to the narrow bandwidth of the transducers. 53or the experiments, a phase resolution of π 16 rad is employed, and the phase of each transducer is transferred from MATLAB to the FPGA via a serial interface using a data transfer rate of 256 kbit=s.The setup of Fig. 1 is used for measuring the acoustic field emitted by the array without the presence of the sphere.The ultrasound signals are captured by a calibrated microphone (Brüel & Kjaer, type 4138-A-015) which is moved by a 3-axis XYZ translation stage (NRT150/M Thorlabs stages and BSC203 Thorlabs motor controller).The microphone signals are amplified by a conditioning amplifier (Brüel & Kjaer, Nexus 2690-A-0S2) and captured by an oscilloscope (Keysight DSOX2014A), which communicates with the PC via USB.In the experiments, the sound field was measured over the XY, XZ, and YZ planes with a spatial resolution of 0.5 mm.To reduce the signal distortion caused by nonlinear propagation and microphone saturation, the sound fields were measured with the array operating at a low voltage amplitude of 3 Vpp.The setup was also used for measuring the individual responses (amplitudes and phases) of each transducer for a better prediction of the acoustic field generated by the array.In addition to the measurements made with the setup of Fig. 1, a high-speed camera (FASTCAM Mini UX50, Photron) was used for recording the sphere oscillation in y and z directions.To record the oscillation in the z direction, the trapping position was switched from z ¼ 49 mm to z ¼ 50 mm.Similarly, the sphere oscillation in the y direction was recorded by keeping z constant at z ¼ 50 mm and switching the horizontal trapping position from y ¼ 1 mm to y ¼ 0 mm.The videos were recorded with a spatial resolution of 28 pixels=mm.A tracking algorithm was used to extract the sphere position from the videos.

V. RESULTS
In preparation of the experiments, we stored the positions r j and normals n j of each transducer as well as their individual manufacturing-related tolerances ΔA j and Δf j in a calibration file.Subsequently, we implemented the model in Sec.II B for the calculation of the sound field and the acoustic radiation force as well as the force gradients with respect to (x, y, z) and phase f j of each transducer in MATLAB.
To solve the optimization problem in Eq. ( 26), we used the C-based solver L-BFGS-B-C 54 for enhanced execution speed.Since the coefficients j S m n , Ψ n , A m n , B m n , G m n , and H m n in Eq. ( 10) and Eqs. ( 14)-( 18) can be pre-calculated for a given position r, the model of Andersson and Ahrens 35 offers considerable advantages since no computationally expensive terms like Bessel functions have to be evaluated during optimization.These two factors contributed significantly to being able to achieve wall times of % 40 ms for one optimization step.Considering the smaller amount of transducers, the required time should still not be as high as the values given in the literature, e.g., in the work of Inoue et al. 34 Although the BFGS method only converges to a local minimum, the quality of the obtained results in our tests was usually very

Journal of Applied Physics
good.Thus, we omitted other methods such as genetic algorithms or basin hopping, which have also been applied to similar optimization problems. 34,41n addition, we have measured the sound field of the created acoustic trap without the sphere, starting from r e and altering each spatial direction by +20 mm.The resulting XY, XZ, and YZ planes are depicted in Figs.2(b), 2(d), and 2(f ).Comparing them with simulation results in Figs.2(a), 2(c), and 2(e), one can see that the contours of the field coincide to a large extent and there is only a maximum deviation of the sound pressure of 60 Pa.This difference between simulated and experimental results can be attributed to the measurement uncertainty.According to our estimates, the measured pressure can differ up to 22 Pa from the true value due to the microphone uncertainty itself (0.2 dB) and the uncertainty caused by the angle of incidence of the incoming wave (up to 1:56 dB).In addition to the microphone uncertainty, one has to take into account that the resulting sound pressure field differs from the measured free sound field due to the presence of the sphere.Further simulative studies lead to the conclusion that the contribution of the scattered field to the total sound pressure field for the experiment in Fig. 2(g) at 12 Vpp is approximately up to 17%.Based on these results, we can conclude that the measured sound pressure field is overall sufficiently well reproduced by the model.Moreover, Figs.2(a)-2(f ) lead to the conclusion that the sphere is levitated in a twin tuning forks trap (TTFT).This name was chosen because in the isosurfaces of Fig. 2(h), the sphere seems to be trapped by a pair of tuning forks.A video showing a rotating 3D view of the acoustic trap is available online (Multimedia view).To retain with the terminology of Marzo et al., 27 its visual shape could be regarded as superposition of two twin trap signatures and a bottle trap signature, see also Figs. 3 and 5 in the mentioned study. 27Following this idea, we subtracted a focusing element from the obtained phase angles related to Fig. 2 to reveal the signature of the trap, see Fig. 3.The signature can be divided into eight different parts.Comparing them with the signatures in Fig. 6 in Marzo et al., 27 there is a noteworthy correlation.Apart from small deviations, it can be clearly seen that there is a phase difference of at least Δf % π between the areas "(1)-( 4)" and "( 5)- (8)," which directly corresponds to a bottle signature.In addition, there is also a notable phase difference of Δf % π between the areas "(5) and ( 6)" as well as "( 7) and (8)."It is quite obvious that this corresponds to the signatures of two twin traps that are arranged crosswise.This thought is underlined by a comparison of the pressure fields in Figs.2(a) and 2(b) with Fig. 3(d) in Marzo et al. 27 While the bottle trap provides the necessary acoustic radiation force in the vertical direction, the two twin traps increase the horizontal restoring forces to enable stable trapping of the sphere.This can also be seen in Fig. 4, where the simulated acoustic radiation force related to the experiment in Fig. 2(g) is shown.Based on this simulation results, the array at 12 Vpp still has sufficient capacity to levitate a bigger sphere than a Mie particle.However, the remaining phase angles f % π=4 (green) and f % Àπ=4 (light blue) could not be clearly assigned to any specific signature.We, therefore, assume that these deviations result from the choice of the abort criterion of the BFGS method in the optimization.
To investigate the DoA for the equilibrium r e , we calculated in Fig. 5, the phase portrait of the reduced dynamic system g 2 (b x) related to the experiment in Fig. 2(g) for the XY (a), XZ (b), and YZ planes (c).Based on these three figures, we assume that the DoA Ω D [ R 3 for the equilibrium r e can be roughly approximated by the set It is also noteworthy that this set can possibly be derived from the smallest common set of a pair of areas between two inflection points which define the limits of one-dimensional stability.To give an example, the course of F x in Fig. 4(a) has inflection points at x in ¼ +4:0 mm, the same applies for F y at y in ¼ +4:0 mm in Fig. 4(b).This observation directly corresponds to the area inside the dotted black circle in Fig. 5(a).In our simulative study, we applied the method of Lyapunov-based sampling to calculate the set Ω E , Ω D for which local asymptotic stability for r e can be guaranteed and have determined this set as where V 2 (b r) ¼ b r `Pb r, c ¼ 2:4543 Â 10 À10 m, and the matrix P [ R 3Â3 is given as 0:2007 0:0125 0:0047 0:0125 0:2168 À0:0105 0:0047 À0:0105 0:9490 from Eq. (36), where A at b r ¼ 0 was computed using Eq. ( 28).
The set Ω D is shown in Fig. 5 in red for the XY (a), XZ (b), and YZ planes (c).If we compare Ω D and Ω E , we can conclude that Ω E approximates Ω D quite well, although there is still a notable difference in Figs.5(b) and 5(c) due to the chosen conservative approach.Furthermore, the small difference for Ω E regarding the XZ and YZ planes in Figs.5(b) and 5(c) can be explained by the small deviation between the values 0:2007 and 0:2168 for X and Y on the main diagonal of P in Eq. (41).In order to correct these deviations and to enlarge Ω E even further, methods like trajectory reversing 52 can be used where Ω E is set as initial value.
To evaluate the quality of levitation in Fig. 2(g), we recorded the position of the sphere for approximately 150 s at 125 frames per second (FPS).A video showing the sphere levitation corresponding to Fig. 6 is available online (Multimedia view).The horizontal and vertical sphere position over time are shown, respectively, in Figs.6(a) and 6(c).The corresponding Fast Fourier Transform (FFT) of the sphere oscillation is also shown in Figs.6(b) and 6(d).The levitating sphere had a total deflection of approximately 0:35 mm along y and 1 mm along the z direction.The low-frequency oscillation of the sphere up to 2 Hz is possibly caused by acoustic streaming and external air currents, which are not considered in our model.Future studies are required to understand and minimize these low-frequency oscillations since it can hinder applications requiring high precision manipulation of the levitated object.The peak at 16 Hz in Fig. 6(b) corresponds to the natural frequency along the y direction, whereas the first peak is

Journal of Applied Physics
attributed to the natural frequency in the z direction.We do not have a definite explanation for the third peak, located at 26:6 Hz, but our hypothesis is that the sphere behavior can be described by a three-dimensional nonlinear oscillator, and the third peak arises due to a nonlinear coupling between different modes of the oscillator.
The levitation stability was also analyzed by switching the trapping position and recording the sphere position over time As described in previous articles, [55][56][57][58] the dynamic behavior of the sphere can be described by a simple model based on a springmass-damper-system, in which the oscillation frequency f a is given as where k a is the linear trapping stiffness along direction a and m is the sphere mass.Similar to what occurs for small particles in a standing wave field, 56,57 the trapping stiffness and, thus, the oscillation frequency of the sphere depend on pressure gradients, which cause the sphere to move to the positions of low acoustic pressure amplitude.The higher the pressure gradient, the higher is the oscillation frequency.The numerical values for k y and k z are found by calculating the angular coefficient from the curves of Fig.

Journal of Applied Physics
maximum deflection in y and z would be k y =k z ¼ 2:87, which is consistent with the experimental ratio of Δz=Δy % 2:86.
Finally, we have performed further simulative studies to investigate the generalizability and the limitations of our presented approach.Based on these results, we presume that our approach can be used not only for Mie spheres but also feasible array configurations can be determined for bigger spheres whose diameter is much larger than the acoustic wavelength λ.Apart from a correspondingly longer computation time to solve the optimization problem in Eq. ( 26), the authors are not aware of any restrictions on the total number of transducers nor on the complexity of the necessary calculations.In addition, the approach can potentially be applied to spheres with other densities or in other media, such as water, possibly using transducers with a higher frequency in the Megahertz range.However, one limitation remains: since we use a spherical harmonics expansion to calculate the acoustic radiation force, the object to be levitated must be a sphere or at least have a sphere-like geometrical shape.Therefore, for levitation of macroscopic objects with arbitrary geometry, one still has to rely on approaches such as the boundary hologram method of Inoue et al. 34

VI. CONCLUSIONS
In the present work, we demonstrated the single-beam acoustic levitation of a Mie sphere using a 2D transducer array.In order to achieve this, we employed the sound field model from Andersson and Ahrens 35 that is based on spherical harmonics expansions to calculate the acoustic radiation force on the sphere.Moreover, we formulated an optimization problem similar to Inoue et al. 34 that maximizes the stability of the sphere while keeping the net force balanced at the equilibrium position.
By simulation, we have proven the local asymptotic stability for the equilibrium with Lyapunov-based methods and determined a sufficient, yet conservative DoA Ω E utilizing our dynamic model of the sphere.In experiments, we demonstrated the asymptotic stable levitation of the sphere inside the acoustic trap over long periods of time.We presume that the low-frequency motion of the sphere up to 2 Hz with a maximum deflection of Δz ¼ 1 mm and Δy ¼ 0:35 mm can be explained by acoustic streaming, which is currently not considered in our dynamic model.The ratio between the maximum deflection in z and y (Δz=Δy % 2:86) is consistent with the ratio between the trapping stiffness along the y and z directions (k y =k z ¼ 2:87), which supports the acoustic streaming hypothesis.
A high-speed camera was also used for recording the sphere oscillation after the trapping position was switched vertically and horizontally.The natural frequencies obtained experimentally presented a reasonable agreement with the oscillation frequencies predicted by the simulation.The differences between simulated and experimental results are attributed to the uncertainty of the microphone response and by nonlinear effects that are not considered in the model, such as the nonlinear propagation of the sound wave and acoustic streaming.
The choice of both model and the L-BFGS-C 54 solver contributed mainly to the fact that we were able to achieve wall times of % 40 ms for one optimization step.However, since usually at least approximately 500 steps are necessary to determine a suitable array output, a practical real-time manipulation via optimization-based methods can currently only be achieved with cumbersome trajectoryshaping via pre-calculation.A promising alternative to this approach can be found in the observed visual shape of the twin tuning forks trap (TTFT).We have shown that this trap according to Marzo et al. 27 is composed of two twin trap signatures and a bottle trap signature.Since only a static pattern of phase angles is necessary to stably levitate the sphere, the approach is suitable for most of the transducer arrays that are mentioned in the literature 46 and places fewer demands on the hardware than the method of acoustic virtual vortices that was proposed by Marzo et al. 33 Consequently, the observed trap signature may facilitate the real-time manipulation of macroscopic objects in mid-air since we presume that it might be possible to levitate macroscopic objects by creating appropriate holographic elements instead of using computationally expensive optimization-based methods.Furthermore, it could also be possible to rapidly translate and rotate macroscopic objects by gradually shifting and rotating the holographic signature of the acoustic trap toward the next desired position or orientation, similar to what was done by Marzo et al. 27 for Rayleigh particles.Finally, it will be very interesting to see which other specific trap signatures can be created by a superposition of basic holographic elements, i.e., a vortex, a twin, and a bottle signature, that can keep macroscopic objects in stable suspension.
polar angle and w ¼ atan2(y, x), and w [ Àπ, π ð the azimuth angle.Moreover, S m n are the expansion coefficients and the time dependency e Àiωt is implicit, where i denotes the imaginary unit.The spherical harmonic bases Y m n are defined as

FIG. 2 .
FIG. 2. In (a), (c), and (e), the simulated incident pressure fields for the XY, XZ, and YZ planes at 3 Vpp are shown.In these diagrams, the center and the circumference of the sphere are indicated by a white dot and a dashed black circle, respectively.The corresponding measurements of the sound field without the sphere are depicted in (b), (d), and (f ).In (h), the simulated amplitude isosurface of 275 Pa (red) corresponding to (a)-(f ) as well as the isosurface of the sphere (green) are shown.A video showing a 3D view of the acoustic trap is available online (Multimedia view).A photography of the experiment with a polystyrene sphere which is levitated at [0, 0, 50] mm is pictured in (g) using 12 Vpp for all 256 transducers.Multimedia view: https://doi.org/10.1063/5.0037344.1

FIG. 3 .FIG. 4 .FIG. 7 . 3 FIG. 6 .
FIG. 3. Holographic signature of the acoustic trap.The phase modulation for the 256 transducers for generating the trap (left) as well as its decomposition into the focusing element (center) and the holographic signature (right) are shown.