# Patterning and micromanipulation of miscible fluids using vortex-based single-beam acoustic tweezers

Samir Almohamad, Gustav Modler, Ravinder Chutani, Udita Ghosh, Sarah

Cleve, Henrik Bruus, Michael Baudoin

► To cite this version:

Samir Almohamad, Gustav Modler, Ravinder Chutani, Udita Ghosh, Sarah Cleve, et al., Patterning and micromanipulation of miscible fluids using vortex-based single-beam acoustic tweezers. Physical Review Applied, 2025, 23 (5), pp.054094. 10.1103/PhysRevApplied.23.054094. hal-04692599v2

HAL Id: hal-04692599

https://hal.science/hal-04692599v2

Submitted on 13 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Patterning and micromanipulation of miscible fluids using vortex-based single-beam
acoustic tweezers

Samir Almohamad,¹,* Gustav K. Modler,²,* Ravinder Chutani,¹ Udita

U. Ghosh,¹ Sarah Cleve,¹, † Henrik Bruus,², ‡ and Michael Baudoin¹,³, §

¹Univ. Lille, CNRS, Centrale Lille, Univ. Polytechnique Hauts-de-France,
UMR 8520 - IEMN - Institut d'Electronique de Microelectronique et de Nanotechnologie, F-59000 Lille, France

²Department of Physics, Technical University of Denmark, DTU Physics Building 309,
2800 Kongens Lyngby, Denmark

³Institut Universitaire de France, 1 rue Descartes, 75005 Paris
(Dated: November 13, 2025)

Vortex-based single-beam tweezers have the ability to precisely and selectively move a wide range of objects, including particles, bubbles, droplets, and cells with sizes ranging from the millimeter to micrometer scale. In 2017, Karlsen and Bruus [Phys. Rev. Appl. 7, 034017 (2017)] theoretically suggested that these tweezers could also address one of the most challenging issues: the patterning and manipulation of miscible fluids. In this paper, we experimentally demonstrate this ability using acoustic vortex beams generated by interdigital transducer-based active holograms. The experimental results are supported by a numerical model based on acoustic body force simulations. This work paves the way for the precise shaping of chemical concentration fields, a crucial factor in numerous chemical and biological processes.

# I. INTRODUCTION

In both optics and acoustics, particle trapping and dis-
placement using radiation forces were initially demon-
strated with standing waves created by interfering
counter-propagating waves. In 1986, Nobel Prize Laure-
ate Arthur Ashkin demonstrated that 3D selective trap-
ping of dielectric particles could be achieved with a single
focused laser beam [1], ushering in a new era of optical
manipulation. Inspired by this seminal work, Junru Wu
[2] began exploring the trapping possibilities offered by
focused beams in acoustics. He demonstrated the abil-
ity to trap particles using two opposing focused beams
interfering near the focus to create a localized standing
wave trap. Later, K.K. Shung and his team explored the
potential of single focused beams for particle trapping
[3-6]. However, as demonstrated later on by Gong and
Baudoin [7] a wide range of objects, including particles
that are stiffer and denser than the surrounding fluids,
can only be trapped laterally with a focused beam (no
axial trap) and only in the Mie regime at specific fre-
quencies near the particle resonances. Therefore, other
methods are required to achieve 3D and robust trapping
of such objects with a single beam.

Along these lines, Baresch, Marchiano, and Thomas
first proposed theoretically [8] and later demonstrated
experimentally [9] that 3D trapping with single beams
can be achieved using specific wavefields known as fo-
cused acoustical vortices, i.e. Bessel beams of topological
order *l* greater than one [10]. Marzo et al. [11] demon-
strated in the Long Wavelength Regime (LWR), through

an optimization algorithm, that these beams are opti-
mal for single beam 3D trapping with the most even
trapping capabilities along the three main axes. They
also investigated the trapping capabilities of these beams
(along with other types, such as twins or bottle beams)
in air. Since then, the trapping capabilities of acoustic
tweezers based on acoustical vortices have been exten-
sively explored [10, 12]. These tweezers have been shown
to precisely and selectively trap and manipulate various
objects, including bubbles [13, 14], droplets [15], kidney
stones [16], microparticles [17–19] and cells [20].

While the manipulation of objects with a clear inter-
face has been widely demonstrated, the ability of single-
beam tweezers to trap and manipulate miscible liquids
with a diffuse interface, has not yet been demonstrated
experimentally. Such manipulation presents a significant
challenge due to (i) the generally weak contrast in den-
sity and compressibility between miscible liquids, (ii) the
absence of an interface between the two fluids, which in
the case of immiscible liquids contributes to the stiffness
contrast of droplets through surface tension, and (iii) the
presence of diffusion, which severely limits the tweezers
manipulation time. In other words, rapid manipulation
is required before diffusion occurs, merges the two liq-
uids, and reduces the contrast at the core of the acoustic
force.

The potential to manipulate inhomogeneous fluids
with acoustic waves was explored experimentally by
Deshmukh et al. [21]. They demonstrated that acous-
tic forces generated by standing waves can modify mis-
cible fluid distributions and stabilize density profiles
against hydrostatic pressure gradients due to the higher
impedance liquid being attracted to pressure nodes.
They subsequently leveraged this discovery to introduce
a new method, termed 'isoacoustics', for sorting cells
based solely on their acoustic contrast, independently of
their size [22–24]. The response of inhomogeneous flu-

<footnote>* These authors contributed equally.</footnote>

<footnote>† sarah.cleve@univ-lille.fr</footnote>

<footnote>‡ bruus@fysik.dtu.dk</footnote>

<footnote>§ michael.baudoin@univ-lille.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. (A) Schematics of the experiments setup: the microchannel is placed on top of an acoustic tweezer, and in the lightpath of an optical microscope. The resulting image is recorded using a sCMOS camera. (B) Design of the tweezer, where the two electrodes are colored in blue and green respectively. (C) Design of the co-flow system for creating parallel water and Ficoll lines; L1, L2, and L3 denote the lengths of the channel sections, and θ represents the angle of the inclined sections relative to the section of length L3. The channel of interest (L3)has a depth of 40 µm and width of 1500 µm. (D) Three example snapshots of the patterning process, which sets in upon turning on the acoustic tweezer.</figcaption></figure>

ids to acoustic fields was later rationalized by Karlsen
et al. [25], who introduced an acoustic force density
model based on density and compressibility gradients
to explain the forces acting on the fluid, as summa-
rized below in Sec. III. Karlsen and Bruus further sug-
gested, through numerical simulations [26], the potential
of acoustic tweezers based on focused beams or acous-
tical vortices to pattern and manipulate miscible fluids,
although this was not demonstrated experimentally.

In this paper, we demonstrate the ability to pattern
and manipulate miscible fluids using vortex-based acous-
tic tweezers synthesized with InterDigital Transducers
(IDTs)-based active holograms. The experimental re-
sults are carefully compared with 2D numerical simu-
lations based on the acoustic body force, showing good
quantitative agreement.

## II. METHOD: EXPERIMENTS

A. Setup description

The experiments consist of two parts: (i) the cre-
ation of miscible fluid concentration patterns imposed by
the vortex wavefield structure and (ii) the trapping and
translation of a fluid blob trapped at the center of the vor-
tex. These experiments were performed inside a 40 µm
deep and 1.5 mm large microfluidic channel using an IDT-
based tweezer (Fig. 1.A). The tweezer (Fig. 1.B) is con-
nected to the bottom of the microfluidic chamber (made
of a glass coverslip) using a drop of silicone oil. The result
is visualized with a sCMOS high sensitivity Prime-BSI
photometric camera operating at a steady frame rate of

12.5 fps and imaging through a Nikon Ti2E optical mi-
croscope equipped with a module for fluorescent imaging.
A typical manipulation sequence consists of the following
steps. (i) A centered narrow band of fluorescently labeled
Ficoll solution (~100 µm in width) surrounded by water
is created using the co-flow system depicted on Fig. 1.C,
with two syringe pumps precisely controlling the fluids
flow rates. (ii) The flow is stopped and the tweezer is
activated. The tweezer creates a vortex beam, which is
focused through a 6.5 mm thick glass wafer glued on top
of the piezoelectric substrate, and then transmitted to
the microfluidic chamber. This produces a pattern rem-
iniscent of the vortex ring structure through respective
motion of the two fluids (Fig. 1.D). (iii) In the manipu-
lation experiments, the blob of fluid trapped at the core
of the vortex is displaced by moving the tweezer hori-
zontally with a highly accurate XY Thorlabs motorized
stage with a precision of 100 nm. (iv) Eventually diffu-
sion results in the merging of the two fluids and hence
the disappearance of the concentration gradients.

B. Acoustic tweezers design, fabrication and characterization

Design

The acoustic tweezers were designed using the method
described in Refs. [18, 20]. In short, the phase map
resulting from the intersection of a propagating (Han-
kel) acoustic vortex [10] of topological order 1 with the
source plane is discretized over 4 levels (0, π/2, π, 3π/2).
The two levels corresponding to opposite phases (0, π)
are materialized into electrodes sputtered on the sur-
face of a 0.5 mm thick, 3-inch Y-36 lithium niobate
(LiNbO₃) piezoelectric substrate, resulting in two inter-
twined spiralling lines of decreasing distance and width
(see Fig. 1.B). The two other phase levels are not materi-
alized into electrodes, which enables to leave some space
between the electrodes.

In the present work, tweezers were designed to cre-
ate a focused vortex beam operating at the frequency
f = 18 MHz. A 6.5 mm thick D263 glass wafer was
glued on top of the electrodes, to ensure focusing (and
hence better localization [18] of the trap) before the beam
reaches the microfluidic chamber. The design of the elec-
trodes was calculated to synthesize transverse waves in-
side the glass, as their speed of sound is closer to the one
of water compared to longitudinal waves, hence ensuring
better wave transmission to the fluid contained in the
microfluidic chamber. The transverse wave speed c<sub>t</sub> used
to compute the propagation through the glass was cho-
sen equal to 3500 ms<sup>-1</sup>, which is an approximate value
of the theoretical value calculated from the glass prop-
erties provided by the glass manufacturer using the for-
mula c<sub>t</sub> = √(E/(2ρ<sub>g</sub>(1+ν))) = 3467 ms<sup>-1</sup>, with E = 72.9 GPa
the Young's modulus, ν = 0.208 the Poisson's ratio and
ρ<sub>g</sub> = 2510 kg m<sup>-3</sup> the glass density. The spiralling struc-
ture was made of 16 turns, with a void central region with
no electrodes to enable visualization through the tweez-
ers. This resulted in a lateral radius of the patterned
region of ~9 mm and hence an aperture of ~54°.

Fabrication

Classical photolithography technique was used to pat-
tern spiraling metallic electrodes from metallic layers de-
posited on the substrate with Plasma Enhanced Vapor
Deposition (PECVD): (i) The LiNbO₃ substrate under-
goes ultrasonication in acetone and propanol for 3 min,
and is then dried with nitrogen gas. (ii) An adhesion pro-
moter called HexaMethylDiSilazan (HMDS) is spread on
the substrate using a spin coater, followed by a negative
AZnLoF2020 photoresist with a thickness of 3 µm. The
resist is cured on a hot plate at 110°C for 90 s, after which
the patterns of the tweezers are transferred into the resist
using an optical mask and a MA6/BA6 SUSS Microtec
UV optical aligner. (iv) The substrate is then heated on a
hot plate at 110°C for 120 s to complete the cross-linking
process, and immersed in AZ326 developer for 30 s before
being rinsed with deionized water. (v) The substrate is
coated with a layer of titanium (30 nm) and gold (30 nm)
with PECVD. (vi) The lift-off process is achieved using
a SVC14 solution at room temperature for one day, fol-
lowed by ultrasonication at 35 kHz and 15% power to
improve the efficiency of the lift-off operation.

A glass wafer of borosilicate D263 T with a diameter of
56.8 mm and a thickness of 6.5mm is then glued on top
of the piezoelectric substrate using optically transpar-

ent epoxy glue (EPOTEK 301-2) to ensure good wave
transmission. Prior to the glueing step, a 15 nm thick
chromium layer is deposited on the upper face of the
glass wafer without etching to serve as markers for local-
izing the vortex center. The two substrates are cleaned
with acetone and propanol and treated with O₂ plasma to
make the surfaces hydrophilic and improve glue spread-
ing. The glue is degassed and mixed using an ARV 310
vacuum mixer to prevent bubble formation, and a drop
of 3.45 µL glue is deposited on the center of the piezoelec-
tric substrate. The glass wafer is then positioned atop the
piezoelectric substrate and left on a horizontal plate un-
til the glue covers the whole surface between the lithium
niobate and glass wafer. After ensuring complete glue
coverage, the structure is left to cure on the plate for two
days at room temperature. The final step involves etch-
ing the markers on the glass wafer, which will enable to
identify the vortex center in the experiments. The whole
structure is cleaned with acetone and propanol and dried
with nitrogen gas to remove dust, and the glass wafer
is coated with AZ1505 resist (thickness about 1 µm) us-
ing a spin coater and placed to dry for 1 day at room
temperature. Once the patterns have been transferred
onto the resist using the optical mask and MA6/BA6
SUSS aligner, the structure is immersed in MIF726 de-
veloper and subsequently washed with deionized water.
The chromium layer is removed using a Cr etchant, and
the tweezers are placed in acetone to remove resist traces,
followed by propanol and drying with nitrogen gas.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Vibrational field created by the tweezer driven at 18 MHz. Left: Normalized amplitude (upper row) and phase (lower row) of the normal displacement measured experimentally with a UHF-120 Laser Doppler Vibrometer. Right: Normalized amplitude (upper row) and phase (lower row) of the field predicted by angular spectrum propagation of the electrodes binary source. Note that the maximum vibration amplitude is of the order of a few nm.</figcaption></figure>
## Characterization

The vibrational acoustic field produced by the tweezers was characterized using a Polytech UHF-120 Laser Doppler Vibrometer (LDV). The LDV measured the normal displacement (amplitude and phase) at the bottom of the microfluidic chamber, i.e. at the top of the glass slide used as the support of the microfluidic PolyDiMethylSiloxane (PDMS) chamber. The wavefield measured experimentally represented on Fig. 2 is compared to simulations performed with Angular Spectrum (AS) code which takes the electrode design as a binary source plane. Both the amplitude and phase are faithful to the numerical predictions. The differences in the amplitude color map stems from slight under-sampling of the acoustic field which does not enable to capture precisely the maxima and minima (also clearly visible on the figure in Sec. III A). The more pronounced experimental anisotropy might be related to the non-isotropic coupling coefficients along different directions owing to the crystal nature of the $LiNbO_3$ resulting in non isotropic amplitude along different directions as well as the presence of a set of two large electrodes (see Fig. 1.B) intersecting the spiraling electrodes to efficiently deliver the current with minimal Joule losses.

Note that the two sets of electrodes of the tweezers are driven by a signal generator IFR 2023A providing a signal at 18 MHz and 5 dBm amplified with an AR50A250 150 W power amplifier, resulting in vibrations at the top of the glass slide with an amplitude of the order of 1 nm. The precise value of the amplitude can vary according to the thickness of the coupling layer.

## C. Channel design and fabrication

### Design

The microchannels are made of a PDMS slab fixed on top of a 150 µm thick glass slide. The same glass (Borosilicate D263) as the one used for the glass wafer was chosen as the bottom part of the channel to ensure good transmission of the wave between the glass wafer and the microchannel. PDMS was chosen to fabricate the microchannel due to its weak impedance mismatch with water and good absorbing properties, resulting in reduced resonance effects in the z direction. The channel design is represented on Fig. 1.C and is similar to the one used in Ref. [21]. It is made of two inlets and a co-flow design which enables the creation of a central Ficoll line surrounded with water. The dimensions of the channel are given in Fig. 1.C.

### Fabrication

The PDMS slab is fabricated with the following procedure. (i) A mold is fabricated with lithography tech-

nique. To that end, a silicon wafer is first cleaned for 15 min using piranha solution (mixture of sulfuric acid ($H_2SO_4$) and hydrogen peroxide ($H_2O_2$) in a 3:1 ratio). The wafer is then rinsed twice with de-ionized water and dried using nitrogen gas. After that, a negative photoresist, $SU8-2035$, is spin-coated onto the wafer to a thickness of 40 µm, followed by baking at 65°C for 3 min, and then at 95°C for 6 min. To transfer the desired shape of the channel from an optical mask to the channel, the mask and wafer are aligned and exposed to UV light, after which the wafer is baked first at 65°C for 2 min, and then at 95°C for 6 min. The wafer is then treated with $SU8$ developer for 5 min, rinsed with propanol, and dried using nitrogen gas. This ends up the fabrication of the mold required for PDMS patterning.

(ii) The mold is then placed in a 3-inch glass Petri dish and subsequently PDMS, mixed with a ratio 10:1 with the curing agent and blended to remove bubbles, is poured onto the mold. It is placed in a preheated oven at 110°C and left to cure for at least 10 min until it solidifies. After cooling down, the PDMS is separated from the mold and the in- and outlets are punched using a biopsy puncher.

(iii) Finally, the PDMS slab is bonded to the glass slide. To achieve this, the glass slide is cleaned with acetone and propanol, while the PDMS is rinsed with propanol and treated with oxygen plasma. Next, the PDMS is gently pressed onto the glass slide, followed by heating the assembly at 70°C for a minimum of 10 min to ensure that the bond solidifies. The PDMS-glass assembly that results from this process is now ready to be used for microfluidic experiments.

## D. Ficoll preparation

The Ficoll solution is made of a mixture of 1.9 g of Ficoll (Cytiva Ficoll PM400), 0.1 g fluorescent-labelled Ficoll (Sigma-Aldrich Polysucrose 400-fluorescein isothiocyanate conjugate) combined with 11.33 g of deionized water to obtain a Ficoll solution with a mass concentration of 15%. While we were able to manipulate Ficoll concentrations ranging from 5% to 30%, a concentration of 15% was selected for the experiments presented here, as it represents a good compromise between acoustic contrast, moderate viscosity, and effective fluorescence imaging. Indeed, for concentrations of 5%, 15%, and 30%, the contrasts between the Ficoll solution and water are 1.02, 1.05, and 1.10 for density, and 0.97, 0.90, and 0.82 for compressibility, respectively, see Ref. [27]. Meanwhile, the viscosities of these Ficoll solutions are approximately 2.2, 11.4, and 129 times higher than that of water. As a result, lower concentrations yield weaker acoustic body forces but with weaker flow resistance, whereas higher concentrations enhance acoustic forces but significantly increase flow resistance. The chemicals were mixed using a magnetic stirrer bar inside a sealed Beaker for several hours until the Ficoll is fully dissolved.
<figure><img src="image_4.png" /><figcaption>FIG. 3. Typical trends (shown here for three representative pixels within the region of interest) used for the calibration of the fluorescent imaging procedure. (A) Evolution of the grey scale pixel value captured by the camera as a function of the Ficoll concentration. (B) Variation of this grey value for 15% Ficoll as function of time featuring a decrease of the fluorescent intensity with time.</figcaption></figure>

### E. Post-processing

To follow the miscible fluids dynamics induced by the activation of the tweezers, it is necessary to track the evolution of the Ficoll concentration. This concentration is tracked by monitoring light emission of fluorescent Ficoll. This fluorescent Ficoll has an excitation peak at 496 nm and an emission peak at 521 nm. Light emitted by the Nikon Intensilight C-HGFI Mercury Fluorescent light source and by the fluorescent Ficoll respectively were filtered by a GFP-4050B bandpass filter cube with excitation and emission bandwidth of 466 ± 25 nm and 525 ± 25 nm respectively and a dichroic mirror at 495 nm. The filtered light emitted by the fluorescent molecule was then recorded using a high-sensitivity sCMOS Prime BSI camera.

The fluorescent emission captured by the camera was converted into concentration with the following calibration procedure: (i) Prior to the experiments, the micro-channel was filled successively with different solutions of Ficoll of controlled mass concentration ranging between 0% and 15%, obtained from the dilution of the initial solution described in the previous paragraph. (ii) Pictures of the channel filled with these different solutions were taken with the exact same settings (relative position of the channel and tweezers, light intensity, microscope and camera settings) as the one used for the final experiments. (iii) Due to a slight shadowing effect induced by the presence of the electrodes resulting in non uniform illumination, a calibration between the grey scale of the image taken by the camera and the concentration was performed for each pixel of the camera in the region of interest. Typical calibration curves are represented on Fig. 3.A, highlighting the linear relation between Ficoll concentration and the grey scale pixel value recorded with the camera. (iv) Finally, as

we observed photobleaching (resulting in a decrease of fluorescence emission over time), the bleaching effect in presence of a quiescent liquid was evaluated for every pixel prior to experiments (see Fig. 3.B for an example), approximated by a linear function and subsequently corrected in the post-processing procedure. To sum up, the post-processing procedure consists thus of a custom code translating grey-color values into fluorescent concentration taking into account the linear dependence of fluorescence on Ficoll concentration and correcting shading effects of the optical path as well as photobleaching over time. This calibration procedure was repeated for each set of experiments.

### III. METHOD: NUMERICS

A numerical model of the acoustic patterning is formulated based on the theory by Karlsen and Bruus [26], in which acoustic fields acting on a short acoustic time scale $t$ give rise to acoustic body forces $f_{ac}$ on inhomogeneous fluid solutions on a slow time scale $\tau$, due to density and compressibility gradients produced by spatially-varying solute concentration fields. By exploiting this disparity between $t$ and $\tau$ within the usual perturbative framework of theoretical acoustics, $f_{ac}$ can be expressed as [25],

$$f_{ac} = -\frac{1}{4}|p_1|^2 \nabla \kappa_0 - \frac{1}{4}|v_1|^2 \nabla \rho_0 \quad (1a)$$

$$= \frac{1}{4} (\kappa_0 |p_1|^2 - \rho_0 |v_1|^2) \frac{\nabla \rho_0}{\rho_0} + \frac{1}{2} \kappa_0 |p_1|^2 \frac{\nabla c_0}{c_0}, \quad (1b)$$

where $p_1$ and $v_1$ are the acoustic pressure and velocity fields, and $\kappa_0$, $\rho_0$, and $c_0$ are the concentration-dependent compressibility, density, and speed of sound of an undisturbed fluid. In this work, the gradients in the material parameters are generated by an inhomogeneous solution of Ficoll PM400 in water characterized by its mass fraction $s$. Once the solution dependencies $\kappa_0(s)$, $\rho_0(s)$, and $c_0(s)$ as well as the acoustic fields $p_1$ and $v_1$ are specified as detailed in the following Sec. III A, the acoustic body force $f_{ac}$ is fully characterized by Eq. (1).

#### A. Description of the simulated configuration

In accordance with the experiments, we consider a shallow microfluidic channel of length $L = 8000$ µm along $x$, width $W = 1500$ µm along $y$, and height $H = 40$ µm along $z$, featuring rigid walls and containing an inhomogeneous solution of Ficoll PM400 in water, as sketched in Fig. 4.A. The local concentration of Ficoll is represented by a mass fraction field $s(r, \tau)$, which initially is a straight band along the $x$-axis, constant along $z$, and with a Gaussian smearing across the width in $y$,

$$s(r, \tau = 0) = s_{\max} \exp \left[ -\frac{1}{2} \left( \frac{y}{\sigma_0} \right)^2 \right], \quad (2)$$
<figure><img src="image_5.png" /><figcaption>FIG. 4. (A) Schematics (not to scale) of the simulated configuration. The dashed region indicates the computational domain of the 2D simulations. (B) Ficoll concentration along the channel width y. A Gaussian g(y) is fitted to the experimental data and used for the numerical simulations. (C) Normalized absolute vibration amplitude in the channel. The experimental data is averaged over the angular coordinate θ. The Bessel function j₁(kₑr) is fitted to the experimental values and used (with respective phase) for the numerical simulations, using a fitting parameter kₑ = 24 063 m⁻¹.</figcaption></figure>

where s_max and σ_0 are determined from best fit to experimental data (Fig. 4.B), which, following the procedure presented in the supplementary material of Ref. [22], leads to the constant diffusivity D_0. Furthermore, the solution is characterized by the s-dependencies of the material parameters given in Ref. [27],

$$ \rho_0(s) = \rho_w(1 + a_1 s), \quad a_1 = 0.349, \quad (3a) $$

$$ c_0(s) = c_w(1 + c_1 s), \quad c_1 = 0.167, \quad (3b) $$

$$ \eta_0(s) = \eta_w \exp(b s), \quad b = 16.2, \quad (3c) $$

$$ \kappa_0 = [\rho_0(s) c_0^2(s)]^{-1}. \quad (3d) $$

An acoustic vortex travels up through the glass and is then transmitted to the fluid of the channel. Owing to (i) the weak acoustic contrast between water and PDMS and strong absorbing properties of the PDMS resulting in weak reflection, (ii) the small depth of the channel compared to the attenuation length (~14 cm in water at 18 MHz), the field intensity variations along the z axis are neglected. Hence, the acoustic wave pressure lateral profile p₁(r) in the fluid is assumed to be the same profile as the one measured with LDV at the bottom of the channel well approximated by the function (see Fig. 4.C),

$$ p_1 = p_a \frac{j_1(k_e r)}{j_1^{\max}} e^{i(\theta - \omega t)}, \quad (4) $$

where (r, θ, z) are the cylindrical polar coordinates with its origin at the center point of the channel, p_a is the acoustic pressure amplitude, k_e an effective wavenumber, ω = 2πf the angular frequency, j_1 the spherical Bessel function of order one, and j_1^max = 0.436 the maximum of j_1. The consistency of such approximation has been verified in ref. [20] through angular spectrum computation of the transmission of the acoustic field from the

glass slide to the fluid. The corresponding acoustic velocity v_1 is then calculated using the first-order momentum equation,

$$ v_1 = - \frac{i}{\omega \rho_0} \nabla p_1. \quad (5) $$

Note that if a spherical vortex was generated with transducers positioned all around the focal point, the effective wavenumber kₑ would correspond to the actual wavenumber k = ω/c_t. But the finite aperture (here 54°) limits lateral focusing resulting in larger rings. Hence, this parameter was optimized to obtain the best match with experiments (Fig. 4.C), leading to the value kₑ = 24 063 m⁻¹.

The equations governing the slow time-scale hydrodynamic velocity v(r, τ), pressure p(r, τ), and Ficoll concentration s(r, τ) consist of the Navier-Stokes equation with f ac as a body force, the continuity equation, and the advection-diffusion equation,

$$ \rho_0 [\partial_{\tau} v + (v \cdot \nabla) v] = \nabla \cdot \sigma + f_{ac}, \quad (6a) $$

$$ \partial_{\tau} \rho_0 + (v \cdot \nabla) \rho_0 = -\rho_0 \nabla \cdot v, \quad (6b) $$

$$ \partial_{\tau} s + (v \cdot \nabla) s = \nabla \cdot (D_0 \nabla s). \quad (6c) $$

Here, σ is the fluid stress tensor, given by

$$ \sigma = -p \boldsymbol{I} + \eta_0 \left[ \nabla \boldsymbol{v} + (\nabla \boldsymbol{v})^{\mathrm{T}} \right] + (\eta_0^{\mathrm{b}} - \frac{2}{3} \eta_0) (\nabla \cdot \boldsymbol{v}) \boldsymbol{I}, \quad (6d) $$

where I is the identity tensor, (·)T indicates the transpose, and η₀ and η₀^b ≈ 2.79η₀ (valid for pure water) are the shear and bulk viscosities, respectively. These governing equations are supplemented with no-slip and no-flux conditions on the rigid walls,

$$ v = 0, \quad (\hat{\boldsymbol{n}} \cdot \nabla) s = 0, \quad (7) $$

with $\hat{\boldsymbol{n}}$ being the outwards pointing unit normal to the boundary surface.
When using the explicit expressions (3) - (5), we find that the acoustic body force $f_{ac}$ given in Eq. (1a) is dominated by the pressure term, $f_{ac} \approx -\frac{1}{4}|p_1|^2 \nabla \kappa_0 \approx \frac{a_1+2c_1}{4}|p_1|^2 \kappa_w \nabla s$. Consequently, $f_{ac}$ in Eq. (6) tends to maintain the high-density Ficoll regions as it is proportional to $\nabla s$, but relocating these regions from high to low pressure areas as it is proportional to $|p_1|^2$. This observation provides a simple interpretation of the evolution of the experimentally observed Ficoll patterns in this work. It is also important to note that in the limit of low density and compressibility contrasts, as well as low acoustic pressure levels, the diffusion process is not affected by the acoustic body force

We have numerically simulated this 3D model; however, memory requirements prevent us from simulating it at the actual size of the system. We hence proceeded by creating a reduced 2D model, obtained by averaging the 3D model over the height dimension along $z$, which approximates the 3D model well in the limit of long acoustic wavelengths compared to the channel height, as quantified by the dimensionless parameter $\epsilon = (k_e H/2)^2 \ll 1$. In experiments, we find that $\epsilon = 0.23$.

The 2D height-averaged model is formulated in terms of the height-averaged fields $\langle v \rangle$, $\langle p \rangle$, and $\langle s \rangle$, defined by

$$ \langle f \rangle(x, y, \tau) \equiv \frac{1}{H} \int_{-H/2}^{H/2} f(x, y, z, \tau) \, dz. \quad (8) $$

When $\epsilon \ll 1$, the acoustic-induced spatial variations in the lateral dimensions, such as the concentration field pattern and underlying flow fields, evolve significantly slower than those in the height dimension. As a result, the velocity field adopts a parabolic profile in the height dimension, $v = \frac{3}{2} \langle v \rangle [1 - (2z/H)^2]$, and the concentration remains largely uniform along $z$, with only a small variation caused by the parabolic velocity profile, which leads to classical Taylor dispersion. This enables us to height-average the governing equations (6), resulting in the following 2D equations,

$$ \eta_0 \nabla_{\parallel}^2 \langle v \rangle - \frac{12\eta_0}{H^2} \langle v \rangle = \nabla_{\parallel} \langle p \rangle - f_{ac} [\nabla_{\parallel} \langle s \rangle], \quad (9a) $$

$$ 0 = \nabla_{\parallel} \cdot \langle v \rangle, \quad (9b) $$

$$ \partial_{\tau} \langle s \rangle + (\langle v \rangle \cdot \nabla_{\parallel}) \langle s \rangle = \nabla_{\parallel} \cdot (D_0 \cdot \nabla_{\parallel} \langle s \rangle), \quad (9c) $$

where $\nabla_{\parallel} = \hat{e}_x \partial_x + \hat{e}_y \partial_y$ is the lateral part of the gradient $\nabla$, and $D_0 = D_0 I + \frac{1}{210} \frac{H^2}{D_0} \langle v \rangle \langle v \rangle$ is the dispersion tensor, valid in the long-wavelength limit (where $\langle v \rangle$ varies slowly) and for sufficiently high Péclet numbers, $Pe = \frac{|\langle v \rangle| H}{2 D_0} \gg 1$. The term $-\frac{12\eta_0}{H^2} \langle v \rangle$ represents the viscous damping exerted by the top and bottom walls. The height-averaged versions of Eq. (7) are

$$ \langle v \rangle = 0, \quad (\hat{n} \cdot \nabla_{\parallel}) \langle s \rangle = 0. \quad (10) $$

This 2D height-averaged model has been simulated and successfully compared with the full 3D model for a down-sized version of the system, where the 3D model is tractable, see Appendix A.

## B. Description of the numerical method

Numerical simulations are carried out in the finite element software COMSOL Multiphysics [28], following a similar procedure to that of Ref. [26], including a standard mesh-convergence analysis [29] (Appendix B). The acoustic fields in Eq. (4) are implemented analytically through $f_{ac}$, and the concentration $s$ is represented by a logarithmic field, $\ln(s)$, to improve numerical stability. The computational domain $\Omega$ has been reduced to the square subdomain indicated by the dashed lines in Fig. 4.A by exploiting the rapid radial decay of $f_{ac}$, which ensures that the associated error of introducing the artificial boundaries is sufficiently small. Symmetries have been exploited to reduce the computational domain further when appropriate. The domain is meshed with free triangular elements, whose sizes increase with distance from the acoustic center. The parameters used for simulations are summarized in Table I.

The time-dependent equations are solved using a backward-differentiation time-stepping method with outputs every 0.04 s and a strict adaptive stepping scheme. The cell Péclet number $Pe_{cell}$ and the total amount of solute, $S = \int_{\Omega} s \, dV$, are monitored at all times during simulations to ensure that $Pe_{cell}$ remains reasonably small and that $S$ is conserved.

# IV. RESULTS

## A. Miscible fluids patterning

A typical sequence of the dynamic evolution of the Ficoll concentration field upon activation of the tweezer is presented in Fig. 5 (top row) and compared to numerical simulations (bottom row). As expected from the density and compressibility contrasts between Ficoll and water, Ficoll is moved toward the minima of the pressure field and water to the maxima (A to D) by the acoustic body force, ending up in the formation of ring structures reminiscent of the vortex amplitude map (Fig. 2). However,

<table><caption>TABLE I. Parameter values used in the numerical simulations. The frequency $f$ is given in Sec. II B, the effective wavenumber $k_e$ and the diffusivity $D_0$ are obtained from fitting to experiments (see Fig. 4), and the remaining parameter values are from Table II in Ref. [27].</caption><thead><tr><th colspan="3">Acoustics</th><th colspan="3">Ficoll solution</th></tr><tr><th>Parameter</th><th>Value</th><th>Units</th><th>Parameter</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td><i>f</i></td><td>18.0 MHz</td><td><i>D</i><sub>0</sub></td><td></td><td>9.19 &times; 10<sup>-11</sup> m<sup>2</sup> s<sup>-1</sup></td><td></td></tr><tr><td><i>k</i><sub>e</sub></td><td>24063 m<sup>-1</sup></td><td><i>&rho;</i><sub>w</sub></td><td></td><td>996.85 kg m<sup>-3</sup></td><td></td></tr><tr><td><i>&epsilon;</i> = (<sup><i>k</i><sub>e</sub><i>H</i></sup>&frasl;<sub>2</sub>)<sup>2</sup></td><td>0.23</td><td><i>c</i><sub>w</sub></td><td></td><td>1496.30 m s<sup>-1</sup></td><td></td></tr><tr><td></td><td></td><td><i>&eta;</i><sub>w</sub></td><td></td><td>0.893 Pa s</td><td></td></tr><tr><td></td><td></td><td><i>a</i><sub>1</sub></td><td></td><td>0.349</td><td></td></tr><tr><td></td><td></td><td><i>c</i><sub>1</sub></td><td></td><td>0.167</td><td></td></tr><tr><td></td><td></td><td><i>b</i></td><td></td><td>16.20</td><td></td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Snapshots of the experimental (top row) and simulated (bottom row) evolution of Ficoll concentration after activation of the tweezer. The reader may refer to the videos available in the supplementary material (Movie S1). The grey triangles in the upper, experimental row correspond to the markers of the tweezer, which are aligned around its center and black in the original snapshots. The circle (blue) and point (green) used to define the normalized concentration difference $\Delta\tilde{s}$ of Eq. 11 between the first ring and the center of the vortex are shown. The colour bar defines the concentration range in the snapshots.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Plot of the time evolution of the normalized concentration difference $\Delta\tilde{s}$ of Eq. (11). Solid lines and dashed lines indicate experimental and numerical results respectively. The labeled dots correspond to the snapshots represented in fig. 5. The inset shows the time evolution of the corresponding center concentration $s_{cent}$ (green), the concentration $\langle s \rangle_{circ}$ averaged over the first maximum intensity ring (blue), and the difference between the two (black). The locations of the center and ring circle chosen for the evaluation of the concentration are represented in green and blue respectively on Fig. 5.</figcaption></figure>

as diffusion becomes dominant, it leads to progressive homogenization of the two fluids and blurred rings (E to F). Indeed, as the acoustic body force is proportional to the gradient of the fluid acoustic properties (density, compressibility) it cannot prevent diffusion. The main features of this time sequence are well captured by the numerical simulations (Fig. 5, bottom row).

For a more quantitative comparison of the evolution of the concentration gradients between regions of maximum and minimum acoustic intensity, we plotted on Fig. 6 (insert) the evolution of the Ficoll concentration $\langle s \rangle_{circ}$ averaged over the maximum intensity ring (represented by a blue circle on Fig. 5) and at the vortex center $\langle s \rangle_{cent}$ (green point on Fig. 5). Since the absolute values of the concentration are very sensitive to the exact position of the circle and center, and the relative position of the tweezers compared to the central Ficoll line, we further compared the normalized concentration difference $\Delta\tilde{s}$ between the circle and center concentration,

$$ \Delta\tilde{s} = \frac{\langle s \rangle_{circ} - \langle s \rangle_{cent}}{\max(\langle s \rangle_{circ} - \langle s \rangle_{cent})} \quad (11) $$

We have plotted $\Delta\tilde{s}$, $\langle s \rangle_{cent}$, and $\langle s \rangle_{circ}$ in Fig. 6, where the labels A to F corresponds to the image sequence of Fig. 5. The curves exhibit a rapid increase in the concentration gradient upon activation of the tweezers (A to D) and then a slow decrease due to the diffusion process (D to F). While the first part of the sequence is very well captured by the predictions of the numerical simulations, some deviations appear once the concentration difference has reached a maximum. A comparison between the two curves indicates that the diffusion process is undervalued by the simulations in this second part of the sequence. These deviations might stem from the exclusion in the numerical models of (i) thermal processes, which could affect the fluid dynamics [30] as well as photobleaching, and (ii) acoustic streaming [31], which could create some additional recirculating flows accelerating the fluid mixing process.
<figure><img src="image_8.png" /></figure>

FIG. 7. (A) Experimental and (B) numerical normalized concentration difference $\Delta\bar{s}$ as a function of time. The curves are normalized by the maximum overall concentration for the experiments and numerical simulations respectively. Each line corresponds to a given vibration amplitude (see scale-bar). Dashed lines indicate amplitudes that are only available numerically (low amplitudes) or experimentally (high amplitudes). The markers correspond to the maxima of these curves, which we define as the characteristic time $t^*$. (C) Inverse of the characteristic times $t^*$ as a function of the vibration amplitude. The inlet shows the same data but with $t^*$ and in a logarithmic scale.

## B. Characteristic time of pattern formation

As it is a crucial parameter for applications, we further investigated the evolution of the time $t^*$ required to reach the maximum concentration difference between the vortex center and the first minimum intensity as a function of the wave amplitude (Fig. 7). Both the experiments and the simulations exhibit a relatively linear decrease (exponent -1.09 in the log-log plot for best fit with the combined experimental and numerical data) of the time $t^*$ with the applied vibration amplitude.

Note that the experimental wave normal amplitude at the beam center was characterized with the LDV and showed a linear dependence with the electrical power applied by our signal generator. Hence, the relative evolu-

<figure><img src="image_9.png" /></figure>

FIG. 8. Translation of a blob of Ficoll solution (here with a concentration of 20% in mass) trapped at the center of an acoustic vortex. The channel is displaced with respect to the tweezer 1.1 s after turning on the tweezer. The channel walls of the PDMS channel are indicated as reference of the displacement undergone.

tion of the vibration amplitude in a set of experiments is known. Nevertheless, the absolute value of the vibration amplitude shall depend on the thickness of the coupling silicone oil layer, which might differ between the LDV measurements and a set of experiments. Hence the absolute vibration amplitude was obtained in Fig. 7 by determining the best match between experimental and numerical values of $t^*$, and assuming a relation $p_1 \approx \rho_w c_w \omega d_1^n$ between the normal vibration amplitude $d_1^n$ and the pressure field $p_1$. The range of vibration amplitudes obtained in this way is of the same order of magnitude as the vibration amplitude measured with the LDV.

## C. Experimental manipulation of a Ficoll blob

To conclude this work, we further explored the possibility to move the blob of fluid trapped at the center of acoustic vortex after formation of the patterns. Fig. 8 and movie S2 show translation along the *y* axis and movie S3 shows the successive displacement along the two orthogonal directions *x* and *y*. These results constitute the first experimental evidence of a trapped fluid blob translation with vortex-based tweezers as theoretically predicted by Karlsen and Bruus [26]. The translational capability could be further improved in the future, by employing fluid combinations with stronger acoustic contrast and tweezers with weaker secondary rings. Note
that in this section, the results have been obtained with a 20% mass concentration Ficoll solution instead of the 15 % used in the patterning experiments due to a better visualisation in the displacement case.

V. CONCLUSION, PERSPECTIVES

In this paper we demonstrated the ability to pat-
tern and manipulate inhomogeneous miscible fluids us-
ing acoustic tweezers based on interdigital transducers.
The experimental trends are recovered through numer-
ical simulations based on the formulation proposed by
Karlsen et al. [25], which relies on the acoustic force
density. Besides the precise manipulation of fluid blobs
with selective acoustical tweezers, the development of a
wealth of holographic techniques – passive [32–40], active
[17–20, 41, 42] and dynamic [9, 43–47]) – paves the way
toward complex patterning of chemical and biomolecule
gradient concentrations within microfluidic chambers.
This opens a wide range of perspectives ranging from
controlled studies of the fundamental roles played by
biomolecular gradients in various biological processes [48]
to the design of reconfigurable microlens arrays [49].

Appendix A: Numerical validation of the
height-averaged 2D model

To validate the height-averaged 2D model, we have
compared it to a full 3D model. However, due to com-
puter memory requirements, this numerical comparison
has been carried out on a reduced size of the system,
namely a cylindrical subdomain with the correct height
H, but with a small radius R = 444 µm. The radius
approximately coincides with a node in the fluid velocity
field to minimize the error from the artificial rigid bound-
ary. The axis of the cylinder goes through the center of
the acoustic wave, and the initial Ficoll band coincides
with the center. All of these measures were necessary to
enable the 3D simulations, which ultimately took several
days to complete.

Figure 9 shows the normalized $L^2$-norm difference $\mathcal{E}_{2D-3D}$ between the concentration difference curves $\Delta \tilde{s}_{2D}$ and $\Delta \tilde{s}_{3D}$ obtained by the 2D height-averaged model (both with and without Taylor dispersion) and the 3D model over the slow time variable $\tau$, as a function of $\epsilon = (k_e H/2)^2$,

$$
\mathcal{E}_{2D-3D} = \sqrt{\frac{\int |\Delta s_{2D} - \Delta s_{3D}|^2 d\tau}{\int |\Delta s_{3D}|^2 d\tau}} . \quad (A1)
$$

The error of the height-averaged 2D model decreases with
decreasing $\epsilon$. Furthermore, the 2D model with Taylor
dispersion consistently produces smaller errors than the
2D model without Taylor dispersion. The parameter $\epsilon$
was varied by changing $H$, and to compensate for the
increased viscous damping at small $\epsilon$ and maintain the

<figure><img src="image_10.png" /><figcaption>FIG. 9. Convergence of the 2D model versus $\epsilon = (k_e H/2)^2$. The error $\mathcal{E}_{2D-3D}$ is calculated as the normalized $L^2$-norm difference between the concentration difference curves obtained by the 2D model (with and without Taylor dispersion) and the 3D model. As $\epsilon$ decreases, better agreement between 2D and 3D models is attained.</figcaption></figure>

Péclet number, the acoustic body was increased by the
square of the relative change in *H* from the default value.

Appendix B: Mesh convergence study

A standard mesh-convergence study was conducted to
ensure that simulations results are independent of the
mesh [29]. Solutions were computed on increasingly re-
fined meshes, each identified by a mesh refinement factor
CONV, which reduces the maximum mesh element size
for increasing values. The normalized L²-norm difference
E was calculated for the solution u for a given value of
CONV, using the solution on the most refined mesh as
the reference u<sub>ref</sub>,

$$
\mathcal{E}_{\text{mesh}} = \sqrt{\frac{\int_{\Omega} |u - u_{\text{ref}}|^2 \, dV}{\int_{\Omega} |u_{\text{ref}}|^2 \, dV}} . \quad (\text{B1})
$$

Figure 10 shows $\mathcal{E}$ vs. CONV, time-averaged over 0 to 5 s for the 2D height-averaged model for $d_1^n = 0.5$ nm. For the final simulations, we used the mesh corresponding to CONV = 2, at which numerical errors become smaller than 0.1%.

ACKNOWLEDGEMENT

The authors acknowledge the support of the Interna-
tional Chair Will ICMINA, financed by Agence Nationale
de la Recherche "ANR-21-IDES-0006".
<figure><img src="image_11.png" /><figcaption>FIG. 10. Time-averaged normalized $L^2$-norm difference between solutions of the 2D height-averaged model on increasingly refined meshes, which are identified by increasing values of CONV, and the solution on the most refined mesh.</figcaption></figure>

[1] A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, Observation of a single-beam gradient force optical trap for dielectric particles, Optics letters **11**, 288 (1986).

[2] J. Wu, Acoustical tweezers, The Journal of the Acoustical Society of America **89**, 2140 (1991).

[3] J. Lee, S.-Y. Teh, A. Lee, H. H. Kim, C. Lee, and K. K. Shung, Single beam acoustic trapping, Applied physics letters **95** (2009).

[4] F. Zheng, Y. Li, H.-S. Hsu, C. Liu, C. Tat Chiu, C. Lee, H. Ham Kim, and K. K. Shung, Acoustic trapping with a high frequency linear phased array, Applied physics letters **101** (2012).

[5] X. Chen, K. H. Lam, R. Chen, Z. Chen, P. Yu, Z. Chen, K. K. Shung, and Q. Zhou, An adjustable multi-scale single beam acoustic tweezers based on ultrahigh frequency ultrasonic transducer, Biotechnology and bioengineering **114**, 2637 (2017).

[6] H.-C. Liu, Y. Li, R. Chen, H. Jung, and K. K. Shung, Single-beam acoustic trapping of red blood cells and polystyrene microspheres in flowing red blood cell saline and plasma suspensions, Ultrasound in medicine & biology **43**, 852 (2017).

[7] Z. Gong and M. Baudoin, Single beam acoustical tweezers based on focused beams: A numerical analysis of two-dimensional and three-dimensional trapping capabilities, Physical Review Applied **18**, 044033 (2022).

[8] D. Baresch, J.-L. Thomas, and R. Marchiano, Three-dimensional acoustic radiation force on an arbitrarily located elastic sphere, The Journal of the Acoustical Society of America **133**, 25 (2013).

[9] D. Baresch, J.-L. Thomas, and R. Marchiano, Observation of a single-beam gradient force acoustical trap for elastic particles: acoustical tweezers, Physical review letters **116**, 024301 (2016).

[10] M. Baudoin and J.-L. Thomas, Acoustic tweezers for particle and fluid micromanipulation, Annual Review of Fluid Mechanics **52**, 205 (2020).

[11] A. Marzo, S. A. Seah, B. W. Drinkwater, D. R. Sahoo, B. Long, and S. Subramanian, Holographic acoustic elements for manipulation of levitated objects, Nature communications **6**, 8661 (2015).

[12] S. Guo, Z. Ya, P. Wu, and M. Wan, A review on acoustic vortices: Generation, characterization, applications and perspectives, Journal of Applied Physics **132** (2022).

[13] D. Baresch and V. Garbin, Acoustic trapping of microbubbles in complex environments and controlled payload release, Proceedings of the National Academy of Sciences **117**, 15490 (2020).

[14] W.-C. Lo, C.-H. Fan, Y.-J. Ho, C.-W. Lin, and C.-K. Yeh, Tornado-inspired acoustic vortex tweezer for trapping and manipulating microbubbles, Proceedings of the National Academy of Sciences **118**, e2023188118 (2021).

[15] S. Lin, A. Riaud, and J. Zhou, Selective sparse sampling of water droplets in oil with acoustic tweezers, ACS sensors **9**, 2066 (2024).

[16] A. C. Barnes, A sound option for the removal of kidney stones, Proceedings of the National Academy of Sciences **117**, 17473 (2020).

[17] A. Riaud, M. Baudoin, O. Bou Matar, L. Becerra, and J.-L. Thomas, Selective manipulation of microscopic particles with precursor swirling rayleigh waves, Physical Review Applied **7**, 024007 (2017).

[18] M. Baudoin, J.-C. Gerbedoen, A. Riaud, O. B. Matar, N. Smagin, and J.-L. Thomas, Folding a focalized acoustical vortex on a flat holographic transducer: Miniaturized selective acoustical tweezers, Science advances **5**, eaav1967 (2019).

[19] R. A. Sahely, J.-C. Gerbedoen, N. Smagin, R. Chutani, O. B. Matar, and M. Baudoin, Ultra-high frequency
vortex-based tweezers for microparticles manipulation with high spatial selectivity and nanonewton forces, arXiv preprint arXiv:2203.05214 (2022).

[20] M. Baudoin, J.-L. Thomas, R. A. Sahely, J.-C. Gerbe-
doen, Z. Gong, A. Sivery, O. B. Matar, N. Smagin,
P. Favreau, and A. Vlandas, Spatially selective manipula-
tion of cells with single-beam acoustical tweezers, Nature
communications **11**, 4244 (2020).

[21] S. Deshmukh, Z. Brzozka, T. Laurell, and P. Augustsson, Acoustic radiation forces at liquid interfaces impact the performance of acoustophoresis, Lab on a Chip **14**, 3394 (2014).

[22] P. Augustsson, J. T. Karlsen, H.-W. Su, H. Bruus, and J. Voldman, Iso-acoustic focusing of cells for size-insensitive acousto-mechanical phenotyping, Nature communications **7**, 11556 (2016).

[23] M. Rezayati Charan and P. Augustsson, Acoustophoretic characterization and separation of blood cells in acoustic impedance gradients, Physical review applied **20**, 024066 (2023).

[24] M. Rezayati Charan, F. Berg, and P. Augustsson, Acoustofluidic three-dimensional motion of suspended cells at near-zero acoustic contrast in homogeneous media, Physical Review Applied **19**, 014046 (2023).

[25] J. T. Karlsen, P. Augustsson, and H. Bruus, Acoustic force density acting on inhomogeneous fluids in acoustic fields, Physical review letters **117**, 114504 (2016).

[26] J. T. Karlsen and H. Bruus, Acoustic tweezing and pat-
terning of concentration fields in microfluidics, Physical
Review Applied **7**, 034017 (2017).

[27] W. Qiu, J. T. Karlsen, H. Bruus, and P. Augustsson, Experimental characterization of acoustic streaming in gradients of density and compressibility, Phys. Rev. Appl. **11**, 024018 (2019).

[28] Comsol Multiphysics 6.1 (2024), https://www.comsol.com/release/6.1.

[29] P. B. Muller, R. Barnkob, M. J. H. Jensen, and H. Bruus, A numerical study of microparticle acoustophoresis driven by acoustic radiation forces and streaming-induced drag forces, Lab Chip **12**, 4617 (2012).

[30] W. Qiu, J. H. Joergensen, E. Corato, H. Bruus, and P. Augustsson, Fast microscale acoustic streaming driven by a temperature-gradient-induced non-dissipative acoustic body force, Phys. Rev. Lett. **127**, 064501 (2021).

[31] J. T. Karlsen, W. Qiu, P. Augustsson, and H. Bruus, Acoustic streaming and its suppression in inhomogeneous fluids, Phys. Rev. Lett. **120**, 054501 (2018).

[32] K. Melde, A. G. Mark, T. Qiu, and P. Fischer, Holograms for acoustics, Nature **537**, 518 (2016).

[33] X. Jiang, Y. Li, B. Liang, J.-c. Cheng, and L. Zhang, Convert acoustic resonances to orbital angular momentum, Physical review letters **117**, 034301 (2016).

[34] N. Jiménez, Formation of high-order acoustic bessel beams by spiral diffraction gratings, Physical Review E **94**, 053004 (2016).

[35] N. Jiménez, V. Romero-García, L. M. García-Raffi, F. Camarena, and K. Staliunas, Sharp acoustic vortex

focusing by fresnel-spiral zone plates, Applied Physics Letters **112** (2018).

[36] S. Jiménez-Gambín, N. Jiménez, J. M. Benlloch, and F. Camarena, Holograms to focus arbitrary ultrasonic fields through the skull, Physical Review Applied **12**, 014016 (2019).

[37] Z. Ma, A. W. Holle, K. Melde, T. Qiu, K. Poeppel, V. M. Kadiri, and P. Fischer, Acoustic holographic cell patter-
ing in a biocompatible hydrogel, Advanced Materials **32**,
1904181 (2020).

[38] J. Li, A. Crivoi, X. Peng, L. Shen, Y. Pu, Z. Fan, and S. A. Cummer, Three dimensional acoustic tweezers with vortex streaming, Communications Physics **4**, 113 (2021).

[39] M. Xu, W. S. Harley, Z. Ma, P. V. Lee, and D. J. Collins, Sound-speed modifying acoustic metasurfaces for acous- tic holography, Advanced Materials **35**, 2208002 (2023).

[40] M. Xu, C. Vidler, J. Wang, X. Chen, Z. Pan, W. S. Harley, P. V. Lee, and D. J. Collins, Micro-acoustic holograms for detachable microfluidic devices, Small , 2307529 (2024).

[41] R. D. Muelas-Hurtado, J. L. Ealo, J. F. Pazos-Ospina, and K. Volke-Sepúlveda, Generation of multiple vortex beam by means of active diffraction gratings, Applied Physics Letters **112** (2018).

[42] R. D. Muelas-Hurtado, J. L. Ealo, and K. Volke-Sepúlveda, Active-spiral fresnel zone plate with tunable focal length for airborne generation of focused acoustic vortices, Applied Physics Letters **116** (2020).

[43] A. Riaud, J.-L. Thomas, E. Charron, A. Bussonnière, O. Bou Matar, and M. Baudoin, Anisotropic swirling surface acoustic waves from inverse filtering for on-chip generation of acoustic vortices, Physical Review Applied **4**, 034004 (2015).

[44] A. Riaud, J.-L. Thomas, M. Baudoin, and O. Bou Matar, Taming the degeneration of bessel beams at an anisotropic-isotropic interface: Toward three-dimensional control of confined vortical waves, Physical Review E **92**, 063201 (2015).

[45] A. Riaud, M. Baudoin, J.-L. Thomas, and O. B. Matar, Saw synthesis with idts array and the inverse filter: toward a versatile saw toolbox for microfluidics and biological applications, IEEE transactions on ultrasonics, ferroelectrics, and frequency control **63**, 1601 (2016).

[46] A. Marzo Pérez and B. W. Drinkwater, Holographic acoustic tweezers, Proceedings of The National Academy of Sciences, 116 (1), 84-89 (2019).

[47] Z. Ma, K. Melde, A. G. Athanassiadis, M. Schau, H. Richter, T. Qiu, and P. Fischer, Spatial ultrasound modulation by digitally controlling microbubble arrays, Nature communications **11**, 4537 (2020).

[48] T. M. Keenan and A. Folch, Biomolecular gradients in cell culture systems, Lab on a Chip **8**, 34 (2008).

[49] X. Gao, X. Hu, J. Zheng, Q. Hu, S. Zhao, L. Chen, and Y. Yang, On-demand liquid microlens arrays by non-contact relocation of inhomogeneous fluids in acoustic fields, Lab on a Chip **22**, 3942 (2022).