Optical Calibration of Holographic Acoustic Tweezers

Recently, acoustic tweezers based on an array of ultrasonic transducers have been reported taking inspiration from holographic optical tweezers. In the latter technique, the calibration of the optical trap is an essential procedure to obtain the trap stiffnesses. On the contrary, in the case of acoustic tweezers the calibration of the acoustic forces is seldom carried out. To cover this gap, in this work, we adapt the calibration protocols employed in optical tweezers to acoustic tweezers based on arrays of ultrasonic transducers. We measure trap stiffnesses in the mN/m range that are consistent with theoretical estimates obtained by calculations of the acoustic radiation forces based on the Gor’kov potential. This work gives a common framework to the optical and acoustic manipulation communities, paving the way to a consistent calibration of hybrid acoustooptical setups.


I. INTRODUCTION
T HERE are different ways to trap and manipulate particles: for instance, it is possible to use light [1], but also magnetic [2] and electric fields [3] or a sound wave [4].In optical trapping, dielectric particles with size varying from tens of nanometers to few micrometers can be studied; however, its applicability is limited in case of light absorbing materials, due to increased scattering force or particle heating [5].Magnetic and electrostatic trapping can be used only with magnetic or charged particles, respectively.On the other hand, acoustic trapping requires less stringent conditions on both the size and the type of particle that can be trapped, allowing interesting applications in biomedical research [6], life sciences [7], [8], physics of liquids [9] and soft-matter [10], [11].
S. Marrara and D. Bronte Ciriza contributed equally to this work.
In acoustic trapping sound waves in the ultrasonic range are used to confine and manipulate millimeter and submillimeter particles in air [12]- [14] or a fluid [15]- [17].The simplest way to obtain particle trapping is the realization of an acoustic standing wave.To this purpose, either an ultrasound source and a reflector [18], or two sources facing each other [19], can be used.With these experimental configurations, both liquid and solid particles have been trapped, even with very large densities (Hg and Ir, respectively) [20].However, these apparatuses have the disadvantage of being expensive and bulky, requiring potentially dangerous high voltages and loosing efficiency after prolonged operation.To overcome these drawbacks, it has been recently reported the possibility of using arrays of commercial ultrasonic transducers both in the standing wave configuration [21] and in a single-sided array configuration [22], [23].In the latter case the shaping of the pressure field is achieved by controlling the phase of the signal emitted by each transducer.This method gets inspiration from holographic optical tweezers, in which an optical diffractive element (tipically a Spatial Light Modulator) is used to shape an incoming light beam [5].With the holographic acoustic tweezers [22], [24] different acoustic pressure field distributions have been obtained, giving twin traps (two close finger-shaped regions among which the particles are trapped), vortex traps and bottle traps [22], which enable the full 3D manipulation of single [22] or multiple [24] millimetric particles.The advantage of this configuration with respect to other techniques of acoustic beam shaping [25] is that the acoustic trap is computercontrolled [23], [26], allowing to change the trap coordinates on demand and in real-time.
In optical tweezers, a routine procedure that enables quantitative force measurements is the calibration of the trap, namely the measurement of its stiffness when a standard size spherical particle is trapped in a medium of known viscosity, and the rescaling in length units of the particle dynamics by means of appropriate calibration factors [5].Many different protocols have been proposed [27]- [32], and a comparison between them can be done [33], aiming at individuating the best method to calibrate the particular system under study.On the contrary, in the case of acoustic tweezers, few examples are found in literature [34]- [36], treating quite small (tens of µm) particles in a microfluidic environment.However, in these cases the full spatial localization of the particles is obtained close to a surface, which may affect the estimate of the trap stiffness.
In this work, we tailor two well assessed calibration methodologies developed in optical tweezers to acoustic tweezers Fig. 1. a) Scheme of the experimental set-up; b) CCD shadow image of a trapped polystyrene sphere; c) 8x8 transducers array (coloured elements) and map of the produced pressure field (in Pa).The trap reference frame is also represented.The trap is focused at (x f , y f , z f )= (0, 0, 2 cm) setups.After a brief introduction on the methods used to trap and manipulate the particles, we show the results of the calibration of holographic acoustic tweezers [22] obtained by both static and dynamic measurements.Finally, we compare these results with acoustic trapping forces calculated on the basis of the Gor'kov potential [37].

II. METHODS
In this section, after some details on the hardware used for the trap calibration (section II-A), we discuss the model used for the acoustic trap (section II-B) and the calculation of the acoustic force (section II-C).

A. Experimental setup
The calibration methodologies employ both video microscopy and a position-sensitive 4-quadrant photodetector to track the fluctuations from the equilibrium position of acoustically trapped particles.The experimental setup is sketched in Fig. 1.The space above an 8×8 array of ultrasonic transducers [22] is crossed by two perpendicular optical paths.In the first one, a diode laser (785 nm, Thorlabs, DL7140-201S), whose beam is made circular by a couple of anamorphic prisms, is directed towards a levitated millimeter sphere.The shadow of the acoustically trapped particle is imaged, through a mirror (Thorlabs, BB1-E03) and two lenses (Thorlabs, N-BK7, -B coated) in a telescope 3:1 configuration, onto a 4-quadrant photodetector (QPD, RS with home-made analog circuit) that records the axial (z axis) or transversal (x or y, depending on the orientation of the trap with respect to the detector) particle position fluctuations.The voltage signal acquired by the QPD is analyzed by custom-made LabView codes.The second optical path is perpendicular to the first one and consists of two lenses in a 2:1 telescope configuration and a CCD camera (Thorlabs Zelux CS165MU/M).Two lamps provide front and back illumination.The shadow of the levitated particle is imaged by the CCD, recording movies of the trapped particle with a sampling rate of approximately 100 Hz, after cropping the image.It is worth noting that (see Supplementary Figure 1) whereas in optical tweezers the laser beam allows both particle illumination and trapping (through a high numerical aperture objective), in our case the laser beam (and lamps in the perpendicular detection path) are used only for the particle illumination, aiming at detecting its fluctuations on QPD (or CCD).The trapping of the particle is carried out by the acoustic pressure field produced by the transducers array.Tracking of the particle is obtained with Python-based homemade codes, giving the trajectory of the particle center and the mean particle radius R. The array of transducers has been assembled following the instructions provided in A. Marzo et al. [23] and in the corresponding GitHub repository [26].The transducers (Murata MA40S4S), emitting at f =40 kHz, are controlled by a suitably programmed Arduino Mega 2560 Rev3 board generating 64 digital periodic signals and a printed circuit mounting 32 TC4427 MOSFET [23] to amplify these signals.The phase of each transducer is changed, according to the chosen shape and focus position of the trap, by means of a Java software [26].The circuits are powered by a DC power supply, with voltages ranging from 16 V to a maximum of 20 V.

B. Modeling of the acoustic trap
To realize the acoustic trap, each j-th transducer is modeled as a circular piston, whose far-field acoustic pressure p j at a target point A(x, y, z) is given by with ϕ j the initial phase of each transducer and Here, d j is the generic distance of the j-th transducer from the A point, κ = 2πf c0 is the wave vector, f is the frequency of the signal emitted by the transducer (40 kHz), c 0 =346 m/s is the velocity of sound in air, a is the transducer radius, θ j is the angle between d j and the normal to the array, J 1 is the Bessel function of the first kind, P 0 =0.17 is a power constant typical of the transducers and V pp is the voltage sent by the power supply.
The initial phase angle ϕ j is chosen in such a way that all the pressure fields p j are in phase at the same point, the trap focus, at (x f , y f , z f )= (0, 0, 2 cm); additional phase contributions need to be taken in account for the generation of twin, vortex or bottle traps [22].In the case of a twin trap oriented along the x direction (see Fig. 1 c), an additional π shift between the initial phases of transducers having x j > x f and x j < x f , respectively, is used.Once defined all the phases giving a trapping point at the focus, the total pressure field p is calculated by summing the p j from each transducer, that is, p = j p j .

C. Calculation of the acoustic force
In case of spherical particles with radius R smaller than the acoustic wavelength λ (in our case, λ =c 0 /f =8.65 mm), a good approximation of the acoustic potential is given by the Gor'kov potential U [37], [38] where p and ⃗ u are the incident acoustic pressure and the velocity fields at the location of the particle.The factors b 1 and b 2 are given by with c 0 and ρ 0 the sound velocity and the density in the medium, c p and ρ p the sound velocity and the density in the particle, and In case of harmonic p and ⃗ u fields and taking in account the relation (see Eq.25 in [38]) the Gor'kov potential can be expressed only in terms of the pressure field and its spatial derivatives: with and with V = 4 3 πR 3 the particle volume and ω the angular frequency of the pressure wave.
Finally, the acoustic force is obtained as

III. RESULTS AND DISCUSSION
An optical trap is characterized by studying the Brownian motion of the trapped particle in the confining optical potential [5], which for small displacements from equilibrium can be considered harmonic.
The dynamics of the trapped particle are described by the Langevin equation where m is the mass of the particle, k x is the trap spring constant, γ is the drag coefficient and χ is a white noise.For clarity, we have considered only the x direction, but analogous equations can be written for y and z directions, and consequently different trap spring constants k y and k z must be considered.At low Reynolds number, that is, when viscous forces on the particles are larger than inertial forces, the inertial term in (11) can be neglected, giving the overdamped Langevin equation where D = k B T γ is the diffusion coefficient and W x (t) is the white noise term.
The calibration of the trap is obtained by tracking the trapped particle fluctuations around the trap equilibrium position by means of a CCD camera or by means of a 4-quadrant photodiode.With both devices, the trap spring constants k i (i = x, y, z) can be estimated by fitting the signal Power Spectra (PS) [5], [27], [28].The fitting of PS gives also the pixels to meters (in case of the CCD detector) or the Volts to meters (in case of the QPD detector) calibration factors.
Inspired by the optical tweezers calibration protocols, we study the dynamics of a levitated spherical particle by using both video microscopy and the QPD device, and then analyze the obtained trajectories by means of the PS approach.In all measurements, we used a twin trap.The trap reference frame (see Fig. 1c), is chosen in such a way that the x direction is the one connecting the two pressure maxima of the twin trap.The y direction is perpendicular to x and parallel to the array of transducers, while z is the direction perpendicular to the array.As already outlined in the section II-A of the Methods, our CCD and QPD detectors are perpendicular to each other, allowing us to study, for the same twin trap, both the x and the y direction.If we want to choose which direction is observed on each detector, we can rotate the orientation of the trap by 90 • .Fluctuations in the z direction can be observed with both devices.
Calibration measurements consist in the study of the position fluctuations of the particle while it is levitating in the trap [39].We study how the particle position fluctuates around the trap equilibrium point when subject to a random perturbation.Moreover, we can also displace in a controlled manner the particle from the trap focus, thanks to the 3D manipulation capability of the setup, and study the damped oscillations of the particle position.In such a way, we are able to measure both the trap spring constants and provide also an estimation of the air drag coefficient γ.
For our measurements, we use styrofoam millimeter spheres.These particles are not sold as laboratory standards, and their physical properties are not provided by the producer.However, to calibrate the acoustic trap we need the density of the material.To this aim, we weighed twenty similar beads on an analytical balance with a sensitivity of 0.1 mg.By averaging on all the particles, we estimate a density of ρ = 36±6 kg/m 3 , which is compatible with the values reported by other authors [40]- [42].

A. Random perturbations
In Fig. 2, the PS calculated by tracking the trajectory of a trapped bead by means of the QPD detector (Fig. 2 a) and the CCD camera (Fig. 2 b) are shown.PS approach is a very useful tool for the analysis of the tracking signals because a periodic motion of the particle appears as a peak in the spectrum [28].
Here, we show the PS along x (yellow curve), y (red curve), and z (black curve) directions obtained by tracking a levitated particle in a twin trap with focus at (x f , y f , z f = 0, 0, 2 cm).To obtain the PS along x or y on the same detector, the trap is rotated by changing the signal phases on the transducers.
Peaks at approximately 8 Hz (f z ), 18.5 Hz (f y ) and 32.4 Hz (f x ) are clearly observed on both devices.In the case of the QPD detector, in x and y power spectra a small peak at f z is also observed, probably due to a small asymmetry of the bead, inducing a tilt on the particle fluctuations and therefore, a cross-talk [5] between QPD channels.A similar effect is not observed in the signals obtained by the CCD device because the tracking software fits a circle to the particle image and, thus, any particle asymmetry cannot affect the reconstruction of the particle trajectory.Due to the limited sampling rate of CCD (100 Hz) with respect to the QPD (5kHz), the PS calculated with CCD signals are in a shorter frequency range, but retain the most important information, i.e., the frequency peaks f x , f y and f z .From these peaks and considering a simple harmonic oscillator model, k i = 4π 2 f 2 i m (i = x, y, z), the spring constants are obtained as k x =16±2 mN/m, k y =5.2±0.7 mN/m and k z =0.97±0.13mN/m, as the particle mass is m=0.39±0.05mg.Finally, the PS obtained with the QPD signals have peaks at frequencies higher than 50 Hz.As these features are observed also without the particle in the trap, or even without laser illumination on the QPD, it is understood that they are periodic spurious signals likely due to electronic noise.
In Fig. 3a the power spectra obtained for particle fluctuations along Y direction are shown at increasing input voltages.As outlined by the black arrow, an increase of the peak frequency f y with increasing V pp is clearly observed.A similar behaviour has been found also for f x and f z .The corresponding trap spring constants k x , k y and k z as a function of V pp are shown in Figure 3b.In the voltage range used, a linear dependence of the k i on V pp is found, even if a quadratic dependence of the trap stiffness on V pp should be observed [37], due to the p 2 term in the Gor'kov potential.It is worth noting that the voltage range accessible in these measurements is limited by the maximum working voltage of the transducers (20 V) and by the minimum voltage required to stably trap these millimetric styrofoam particles (approximately 15-16 V).Due to this short voltage range, the quadratic dependence of the trap spring constants cannot be observed, whereas an "effective" linear dependence is found.
As in optical trapping [43], [44], the symmetry properties of the the trap can be studied by using parameters connected to the trap spring constants.The parameter k p = 1− kx ky can be used to estimate the symmetry on the xy plane.If the acoustic trap was completely isotropic in this plane, k p should be zero.We find that this parameter is close to -2, indicating a strong transversal anisotropy of the trap, due to a corresponding anisotropy in the confining potential (see Fig. 5).The trap aspect ratio can be studied by means of the k as = kx+ky 2kz parameter.In our case, we find a quite large value, indicating a prolate shape of the confining potential and of the trapping region in the axial direction, which further become more prolate at increasing V pp (see Fig. 3c and d).Note that at increasing V pp the acoustic force increases while the particle weight remains constant.Thus, the trap equilibrium point changes and, in particular, is located at higher z (Fig. 3d).The Gor'kov potential (7)(8)(9) shows a dependence on the third power of the radius.Thus, a correspondent dependence is also expected on the trap spring constants.To verify this, we trapped four different particles in the same experimental conditions (V pp =16 V, trap focus at x f , y f , z f = 0, 0, 2 cm) and measured their fluctuations in the trap.The analysis of the PS of the particle fluctuations in the trap allowed the estimation of the trap spring constants, which are shown in Fig. 4 as a function of R in log-log scale.The dependence of the k i on the third power of the particle radius is recognized.

B. Induced perturbations
We induced particle position perturbations by periodically changing the position of the trap focus.We shifted the focus by 0.5 mm along x, y and z directions at 5 s time intervals and we acquired the consequent damped oscillations of the particle position.
In the insets of Figure 5 the tracking signals registered on the QPD (x direction) and CCD (y and z directions) are shown.
Note that the faster oscillation is sent on the high acquisition rate detector (the QPD).The signals are characterized by steps, whose initial part (approximately 1 s) is an oscillation with a damped amplitude.After this, the signal becomes similar to what is observed in the random perturbation case.Thus, we decided to focus on the first part of the signal and to fit it with the damped harmonic oscillator model: where A is the amplitude, Γ 0 = γ m the damping coefficient, γ the viscous drag coefficient, m the mass of the particle, ω the angular frequency of the oscillator, t 0 the initial time, φ a phase constant and y 0 a signal offset.Here, ω is related to the oscillator natural angular frequency Ω 0 = k m by the following relation: By fitting the particle damped oscillations in x, y and z directions we can calculate both the trap spring constants k i = Ω 2 0,i m and also estimate the medium (air) viscous drag coefficient γ.
In the insets of Fig. 5 the whole signals recorded along the three x (a), y (b) and z (c) directions are shown.The step structure of the signals is easily recognized.One of the steps along with its fit (red curve) is also shown for each direction.The frequency of oscillation, and thus the trap spring constant, decreases from x to z, in agreement with the results obtained with random perturbations analysis.By averaging on the Γ 0 and ω values obtained from the fit of many different steps, the calculated trap spring constants are k x = 19 ± 3 mN/m, k y = 6 ± 1 mN/m, k z = 1.2 ± 0.2 mN/m, which are consistent with results obtained by analyzing only the random perturbations.Moreover, the values of Γ 0 obtained for each direction can be used to estimate also the average drag coefficient γ = (2.9 ± 0.8) • 10 −6 N•s/m, which has the same order of magnitude and is 6 times larger than the theoretical value γ 0 =0.47•10 −6 N•s/m expected for air at 25 • C. Note that, by averaging γ on the three directions, we model the particle as a perfect sphere, even if we know that some small asymmetries are present.This could explain the discrepancy between our result and the expected value for a perfect sphere.

C. Simulation of the acoustic trap.
To verify the agreement between experimental and theoretical results, we calculated the total pressure field, the Gor'kov  and c) the data are acquired with the CCD camera.The steps are fitted with a damped harmonic oscillator model (see Eq. 13) to estimate the oscillator frequency ω and the damping term Γ 0 .An example of the fitting is given (red curve) for each direction.Note that the time scale is reduced to t − t 0 aiming at a comparison between different step traces.Note also that, going from x to z, ω decreases, pointing out, respectively, a lower k.
potential, the force fields, the components of the acoustic forces along the three spatial directions and the trap spring constants for our 8×8 array of transducers.
For the calculations, we consider a particle with radius R=1.36±0.02mm, as in our measurements, and we use ρ 0 =1.18 kg/m 3 as the density of the medium, ρ p =36 kg/m 3 as the density of the particle material, and c p =900 m/s as the sound velocity in the particle.Finally, we used a voltage supply of V pp =18 V.
In Fig. 6a, a xy map of the Gor'kov potential in the z f =2 cm plane is shown.The potential shows an elongated shape which explains the anisotropy found in the transversal trap spring constants.The most stable point is at the focus.d) correspondent maps of the acoustic force for a R =1.36±0.02mm bead.The twin trap has been considered in the calculation.Focus has been set at x f =0, y f =0, z f =2 cm above the 8×8 array.e-g) Calculation of the x (e), y (f) and z (g) components of the acoustic force for the same bead.The trap spring constants kx, ky and kz have been calculated by a linear fit of the force in the trapping region.Note that in g) also the weight of the particle and the buoyancy have been taken in account.For this reason, the trapping point height is below z f and located at approximately 1.5 cm.
By calculating the gradient of the Gor'kov potential, it is possible to obtain the maps (Fig. 6b-d) of the acoustic force in the x, y and z direction.In case of the z direction, we have calculated the total axial force, which takes in account also the gravity and the buoyancy of the particle.From the maps, it is clear the anisotropic character of the force field, which has a different structure both in the transversal and axial plane.
The origin of the anisotropic trapping in holographic acoustic tweezers is in the anisotropic shape of the pressure field [22].Forces along x direction are due to gradients of the pressure field, whereas forces along y and z direction are due to gradients of the velocity field [22], with velocity itself obtained as the gradient of the pressure field (6).Moreover, along z also gravity and buoyancy control the trap equilibrium position.Thus, anisotropy in the trapping forces along the three spatial directions is expected and confirmed by the anisotropy of the trapping potential (see Fig. 6 a).
To individuate the trap equilibrium point and to calculate the trap spring constants, we can plot the F x , F y and F z components of the acoustic force along x, y and z directions, respectively.The point where the force vanishes with a negative slope is a stable equilibrium point and corresponds to the trap position.By linearly fitting the curves at the trap position we can determine the trap spring constants.In Figs.6(e)-6(g) the forces (red curves), the linear fits (dashed black lines) and the trap spring constants calculated (insets) are shown.The values obtained agree fairly well with the measured values, confirming that the Gor'kov potential approximation can be used to model the acoustic trapping of millimetric particles in air with this setup.

IV. CONCLUSIONS
In this work, we have presented a calibration procedure for an acoustic tweezers set-up based on a flat array of transducers.The procedure is inspired by protocols used for optical tweezers.The trap spring constants and the medium drag coefficient have been measured by using both random perturbations and controlled displacements from the equilibrium position.We have shown that the Gor'kov potential and the acoustic force calculated from it are a good approximation to acoustic forces that come into play in real setups.We think that these protocols, sharing the same tools used for optical tweezers, may bridge the gap between the two optical and acoustic trapping communities, paving the way to the realization of hybrid setups where both light and sound can be used to trap particles in air.

Fig. 2 .
Fig. 2. Power spectra of spontaneous oscillations of a levitated particle along the trap reference frame xyz.a) PS of the signal registered with QPD; b) PS of the signal registered with the video camera.Data are shown both in original units (V 2 s in the QPD case and pixel 2 s in the CCD case, left axes) and in calibrated units (m 2 s, rigth axes).Dashed lines point out correspondence between similar peaks.From data and the harmonic oscillator model, the trap stiffnesses are calculated as kx=16±2 mN/m, ky=5.2±0.7 mN/m and kz=0.97±0.13mN/m.

Fig. 3 .
Fig. 3. a) Power spectra measured at input voltage ranging from 16 V to 20 V for the same trap and particle shown in Fig. 2. The fy oscillation shifts towards higher values at increasing Vpp.Similar behaviour is observed for fx and fz (not shown).The power spectra have been displaced vertically for clarity.The black arrow is a guide for the eye.b) Trap stiffnesses kx, ky and kz as a function of Vpp.A linear dependence on the input voltage is observed in the voltage range used for measurements.c) Trap parameters kp and kas as a function of Vpp.While the transversal asymmetry pointed out by kp is approximately constant at increasing Vpp (orange circles), the trap aspect ratio, pointed out by kas (red squares), increases at increasing Vpp.d) Distributions of the trapped particle position.For clarity, they are shown only at the two extreme voltages 16 V and 20 V. The red dashed lines highlight the shift (approx. 1 mm) in the trap equilibrium position at increasing voltage.

Fig. 4 .
Fig. 4. Dependence of trap spring constants on particle radius.kx, ky and kz are shown as a function of R for four different particles.To highlight the k i functional dependence on the radius R, different models (black solid line for R 3 , dashed gray lines for R 2 and R 4 ) are plotted as guides to the eye.The R 3 dependence of the spring constants is recognized.

Fig. 5 .
Fig. 5. Damped harmonic oscillations of an acoustically trapped bead along x (a), y (b) and z (c) directions.The oscillations are induced by periodically changing the particle position along the same directions.The whole traces obtained are shown in the insets.In a), the QPD detector is used, while in b)and c) the data are acquired with the CCD camera.The steps are fitted with a damped harmonic oscillator model (see Eq. 13) to estimate the oscillator frequency ω and the damping term Γ 0 .An example of the fitting is given (red curve) for each direction.Note that the time scale is reduced to t − t 0 aiming at a comparison between different step traces.Note also that, going from x to z, ω decreases, pointing out, respectively, a lower k.

Fig. 6 .
Fig.6.Calculation of the a) Gor'kov potential and b-d) correspondent maps of the acoustic force for a R =1.36±0.02mm bead.The twin trap has been considered in the calculation.Focus has been set at x f =0, y f =0, z f =2 cm above the 8×8 array.e-g) Calculation of the x (e), y (f) and z (g) components of the acoustic force for the same bead.The trap spring constants kx, ky and kz have been calculated by a linear fit of the force in the trapping region.Note that in g) also the weight of the particle and the buoyancy have been taken in account.For this reason, the trapping point height is below z f and located at approximately 1.5 cm.