# Numerical study of acoustic cell trapping above elastic membrane disks
driven in higher-harmonic modes by thin-film transducers with patterned
electrodes

André G. Steckel<EMAIL_ADDRESS>Department of Physics, Technical
University of Denmark, DTU Physics Building 309, DK-2800 Kongens Lyngby,
Denmark  Henrik Bruus<EMAIL_ADDRESS>Department of Physics, Technical
University of Denmark, DTU Physics Building 309, DK-2800 Kongens Lyngby,
Denmark

(23 December 2021)

###### Abstract

Excitations of MHz acoustic modes are studied numerically in 10-µm-thick
silicon disk membranes with a radius of 100 and 500 µm actuated by an attached
1-µm-thick (AlSc)N thin-film transducer. It is shown how higher-harmonic
membrane modes can be excited selectively and efficiently by appropriate
patterning of the transducer electrodes. When filling the half-space above the
membrane with a liquid, the higher-harmonic modes induce acoustic pressure
fields in the liquid with interference patterns that result in the formation
of a single, strong trapping region located 50 - $100~{}\textrm{\textmu{}m}$
above the membrane, where a single suspended cell can be trapped in all three
spatial directions. The trapping strength depends on the acoustic contrast
between the cell and the liquid, and as a specific example it is shown by
numerical simulation that by using a 60% iodixanol solution, a cancer cell can
be held in the trap.

##  I Introduction

Recently, the concept of thin-film-actuated devices has been introduced in the
field of microscale acoustofluidics at MHz ultrasound frequencies. In 2018,
Reichert et al. demonstrated experimentally and numerically, how to generate
useful acoustofluidic responses in microchannels with a thin silicon-membrane
lid driven by a lead-zirconate-titanate (PZT) thin-film transducer [1]. In
this case, the transducer makes up around 15% by volume (v/v) of the actuated
membrane, which is excited while leaving the bulk part of the device inert. In
2021, we successfully modeled and experimentally validated the excitation by
aluminum scandium nitrite (AlSc)N thin-film transducers of MHz-modes in
millimeter-sized glass-block devices without microchannels [2]. In this
system, the transducer is only 0.2% v/v. In a follow-up study [3], which
constitutes the theoretical foundation of our present work, we demonstrated by
numerical simulation, how thin-film transducers can induce acoustofluidic
responses in bulk microfluidic devices on par with those obtained by using
conventional bulk PZT transducers. In Ref. [3], we also studied the robustness
and pointed out several advantages of using thin-film transducers: the low
sensitivity of the thin-film device to the material, the thickness, and the
quality factor of the thin-film transducer, and that the microfabrication
techniques, by which the thin-film devices can be produced, also allow for
careful shaping of the transducer electrodes, a design freedom that may be
used to boost the acoustofluidic response of the device.

In this work, we follow up on the latter idea, and present a numerical study
of how proper shaping of the electrodes on circular thin-film-driven membranes
can increase the excitation amplitudes of the higher-harmonic vibration
eigenmodes in the membrane. When placing such membranes in the bottom wall of
a cavity containing a liquid, we show that specific higher harmonics of the
order $n\gtrsim 10$ in the membrane induce acoustic pressure fields in the
liquid with interference patterns that result in the formation of a single,
strong trapping region located 50 - $100~{}\textrm{\textmu{}m}$ above the
membrane, where a single suspended cell can be trapped in all three spatial
directions. The choice of this model system is motivated by the increasing use
of disk-shaped membranes in acoustofluidic applications, such as capacitive
micromachined ultrasonic transducers (CMUT) for imaging, inkjet printing, and
testing [4], thin-film resonators for mixing and biosensing [5, 6], and
silicon-membrane devices for particle manipulation [7].

In the field of acoustofluidics, electrode shaping is of course used
extensively when dealing with surface acoustic waves (SAW), as the electrodes
directly define these waves [8]. However, when dealing with bulk acoustic
waves (BAW), the topic of this work, electrode shaping is rarely used, and not
at all for the above-mentioned membrane devices. A simple split-electrode
configuration with an applied anti-symmetric driving voltage, has been used on
experiments on bulk piezoelectric (PZE) transducers [9, 10] and on thin-film
transducers [1] to obtain a strong excitation of anti-symmetric modes for
optimal particle focusing. Such systems has also been studied in numerical
simulations [11, 12, 13, 3]. Furthermore, in a combined experimental and
numerical study, Hammarström et al. used a square-shaped back and a full front
electrode on a bulk PZE transducer to create a dynamically-defined array of
particle traps in the liquid just above the transducer [14].

In the field of microelectromechanical system (MEMS), shaping of PZE-
transducer electrodes has been studied in much more detail. For example the
development of energy harvesting MEMS devices involving shaping of electrodes
is an active research field [15, 16, 17, 18, 19]. A particular thorough study
is the combined experimental and theoretical work on the excitation and
detection of more than 50 permitted arbitrary modes in disc-, plate-, ring-,
and beam-shaped PZT-on-silicon resonators by electrode shaping [20]. It is
this kind of selective excitation of chosen resonance modes that we in the
present work extend from solid-state MEMS to acoustofluidic devices for
particle handling.

The contents of the paper is as follows. In Section II we introduce our model
system, a disk-shaped silicon-membrane-based microscale acoustofluidic system
driven at MHz-frequencies by an (AlSc)N PZE thin-film transducer with
patterned electrodes for selective excitation of higher-harmonic resonance
modes, a system chosen for its compatibility with standard MEMS fabrication
techniques. A brief summary is given of the basic theory and modeling
developed in our previous work [3]. It includes time-harmonic perturbation
theory, the electromechanical theory of the elastic membrane and the linear
PZE transducer, acoustics and time-averaged acoustic streaming of the liquid,
boundary conditions, and the numerical implementation in the software COMSOL
Multiphyscis [21]. In Section III, we present the main result of the paper,
the selective enhanced excitation by electrode shaping of higher-harmonic
membrane resonance modes for acoustofluidic applications. In Section IV, we
present an application example by showing how specific higher-harmonics
membrane modes may induce the above-mentioned trapping of a single suspended
cell. In Section V, we discuss the possibility of size-dependent cell trapping
and the advantages and disadvantages of the presented method of trapping, and
finally in Section VI we present our main conclusions.

##  II Model system, theory, and numerical implementation

The model system, sketched in Fig. 1, consists of a thin disk-shaped silicon
membrane completely covered by a (AlSc)N thin-film PZE transducer on one of
its surfaces. The ground electrode of the transducer is always fully covering
the transducer, whereas the excitation electrode may be divided into several
pieces each with individual alternate-current (AC) excitation voltages. Below
the membrane is air (treated as vacuum), and above it a liquid. The system is
axisymmetric, so the we use cylindrical coordinates $(r,\phi,z)$ through the
paper with the $z$ axis perpendicular to the membrane through its center.

The basic theory and modeling for such a thin-film PZE-transducer-driven
acoustofluidic system, was developed in a perturbation scheme involving the
acoustic first-order fields and the steady time-averaged second-order fields
in our previous work [3], founded on the theory for bulk PZE-transducer-driven
systems [22] taking the acoustic boundary layers into account analytically
through effective boundary conditions [23]. Numerical simulations based on
this theory have been validated experimentally for several different
microscale acoustofluidic systems [22, 9, 10, 2]. In the following we briefly
summarize this basic theory and its numerical implementation and adapt the
previous cartesian-coordinate formulation into the cylindrical coordinates of
the present axisymmetric system. We do not study azimuthal variations, so the
full three-dimensional (3D) problem is independent of $\phi$ and is thus
reduced to a two-dimensional (2D) problem in the radial and axial coordinates
$r$ and $z$, respectively.

Figure 1:  The axisymmetric model in the $r$-$z$ plane consisting of a thin
silicon membrane (gray) clamped for $r>r_{\mathrm{mv}}$ and fully covered with
a thin-film (AlSc)N PZE transducer (beige) on its top surface. The grounded
transducer electrode (black) covers the entire transducer, but the excitation
electrode (red and blue) may be either fully covering or divided into several
sections each with the indicated excitation voltage. Below the membrane is air
(white, treated as vacuum) and above is a liquid (cyan). (a) The fundamental
resonance mode ($n=0$) of the membrane with a large amplitude as the fully-
covering excitation electrode (red) is compatible with the mode shape. (b) The
third-harmonic ($n=3$) resonance mode of the membrane with a small amplitude
as the fully covering excitation electrode (red) is incompatible with the mode
shape. (c) The third-harmonic ($n=3$) resonance mode of the membrane with a
large amplitude due to the specific patterning of the excitation electrode
(red and blue) into four pieces with specific excitation voltages.

###  II.1 First-order acoustic fields

Using first-order perturbation theory and complex-valued fields [3], the time-
harmonic electric potential $\tilde{\varphi}_{1}(\bm{r},t)$ applied to the PZE
thin-film excites a time-harmonic acoustic displacement field
$\tilde{\bm{u}}_{1}(\bm{r},t)$,

$\tilde{\varphi}_{1}(\bm{r},t)=\varphi_{1}(\bm{r})\>\mathrm{e}^{-{\mathrm{i}\omega
t}},\quad\tilde{\bm{u}}_{1}(\bm{r},t)=\bm{u}_{1}(\bm{r})\>\mathrm{e}^{-{\mathrm{i}\omega
t}},$ (1)

where $\omega=2\pi f$ is the angular frequency, $f$ the frequency, $\bm{r}$
the spatial coordinate, and $t$ the time.

We convert the cartesian coordinates $\bm{r}=(x,y,z)$ of Refs. [22, 3, 2] into
cylindrical coordinates $\bm{r}=(r\cos\phi,r\sin\phi,z)$, and assuming
$\phi$-independent excitation voltages, we have,

$\varphi_{1}(\bm{r})=\varphi_{1}(r,z),\quad\bm{u}_{1}(\bm{r})=u_{r}(r,z)\bm{e}_{r}+u_{z}(r,z)\bm{e}_{z}.$
(2)

The governing equations of $\bm{u}_{1}$ and $\varphi_{1}$ are the weakly
damped Cauchy elastodynamic equation and the quasi-static Gauss equation for a
charge-free linear dielectric,

$-\rho\omega^{2}(1+\mathrm{i}\Gamma_{\mathrm{sl}})\>\bm{u}_{1}=\bm{\nabla}\cdot\bm{\sigma}_{\mathrm{sl}},\quad\bm{\nabla}\cdot\bm{D}=0,$
(3)

where $\bm{\sigma}_{\mathrm{sl}}$ is the Cauchy stress tensor,
$\Gamma_{\mathrm{sl}}$ is the damping coefficient of the solid, and
$\bm{D}=-(1+\mathrm{i}\Gamma_{\varepsilon})\bm{\varepsilon}\cdot\bm{\nabla}\varphi_{1}$
is the electric displacement field with $\Gamma_{\varepsilon}$ and
$\bm{\varepsilon}$ being the dielectric damping coefficient and tensor,
respectively. For a purely mechanical solid, the stress tensor is defined in
Voigt notation as,

$\bm{\sigma}_{\mathrm{sl}}=\bm{C}\cdot\bm{S}_{\mathrm{sl}},$ (4)

where $\bm{C}$ is the elastic moduli tensor of rank two which relates the
strain tensor $\bm{S}_{\mathrm{sl}}$ to the stress tensor
$\bm{\sigma}_{\mathrm{sl}}$. In cylindrical coordinates, the Voigt notation-
form of $\bm{C}$, $\bm{S}_{\mathrm{sl}}$ and $\bm{\sigma}_{\mathrm{sl}}$ is,

$\displaystyle\bm{C}$
$\displaystyle=\left(\begin{array}[]{c@{\:}c@{\:}c@{\:}|c@{\:}c@{\:}c}C_{11}\hfil\>&C_{12}\hfil\>&C_{13}\hfil\>&0\hfil\>&0\hfil\>&0\\\
C_{12}\hfil\>&C_{11}\hfil\>&C_{13}\hfil\>&0\hfil\>&0\hfil\>&0\\\
C_{13}\hfil\>&C_{13}\hfil\>&C_{33}\hfil\>&0\hfil\>&0\hfil\>&0\\\ \hline\cr
0\hfil\>&0\hfil\>&0\hfil\>&C_{44}\hfil\>&0\hfil\>&0\\\
0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&C_{44}\hfil\>&0\\\
0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&0\hfil\>&C_{66}\\\ \end{array}\right),$
(5g) $\displaystyle\bm{S}_{\mathrm{sl}}$
$\displaystyle=\Big{(}\partial_{r}u_{1r},\frac{1}{r}u_{1r},\partial_{z}u_{1z},0,\partial_{r}u_{1z}+\partial_{z}u_{1r},0\Big{)}^{\mathrm{T}},$
(5h) $\displaystyle\bm{\sigma}_{\mathrm{sl}}$
$\displaystyle=\Big{(}\sigma_{rr},\sigma_{\phi\phi},\sigma_{zz},\sigma_{\phi
z},\sigma_{rz},\sigma_{r\phi}\Big{)}^{\mathrm{T}}.$ (5i)

Here, the column vectors $\bm{S}_{\mathrm{sl}}$ and
$\bm{\sigma}_{\mathrm{sl}}$ are written as the transpose “T” of the row
vectors. Moreover, $\bm{C}$ is written for the lowest symmetry in our problem,
the hexagonal crystal structure of (AlSc)N, and we note that it has the same
form in cartesian as in cylindrical coordinates when the polarization axis is
parallel to the $z$ axis. Whereas neither the displacement nor the electric
field has a component in the azimuthal $\phi$-direction, the strain and stress
do have a nonzero $\phi\phi$-component. Finally, we note that for silicon
(111) the above form for $\bm{C}$ still applies, but with the symmetry
constraints $C_{13}=C_{12}$, $C_{33}=C_{11}$, and $C_{66}=C_{44}$. Other
orientations of silicon crystal are not isotropic in the $r$-$\phi$ plane and
thus cannot be modeled as an axisymmetric system [24, 25].

In this work, our PZE transducer is taken to be (AlSc)N which has a hexagonal
piezoelectric crystal structure [26] with the material parameters listed in
Table 1. For PZE materials, the stress and electric displacement fields are
coupled to the electric and strain fields though the piezoelectric coupling
tensor,

$\displaystyle\bm{\sigma}_{\mathrm{sl}}$
$\displaystyle=\bm{C}\cdot\bm{S}_{\mathrm{sl}}-\bm{e}^{\mathrm{T}}\cdot\bm{E},$
(6a) $\displaystyle\bm{D}$
$\displaystyle=\bm{e}\cdot\bm{S}_{\mathrm{sl}}+\bm{\varepsilon}\cdot\bm{E},$
(6b)

where the electric field defined as $\bm{E}=-\bm{\nabla}\varphi_{1}$, and
where in the Voigt notation for the hexagonal crystal structure, $\bm{e}$ and
$\bm{\varepsilon}$ are the following tensors,

$\displaystyle\bm{e}$
$\displaystyle=\left(\begin{array}[]{ccc|ccc}0&0&0&0&0&e_{15}\\\
0&0&0&0&e_{15}&0\\\
e_{31}&e_{31}&e_{33}&0&0&0\end{array}\right)\\!,\;\bm{\varepsilon}=\left(\begin{array}[]{ccc}\varepsilon_{11}&0&0\\\
0&\varepsilon_{11}&0\\\ 0&0&\varepsilon_{33}\end{array}\right)\\!.$ (13)

The mechanical displacement in the solids couple into the adjacent liquid or
fluid as pressure waves. For a fluid with density $\rho_{\mathrm{fl}}$,
dynamic viscosity $\eta_{\mathrm{fl}}$, bulk viscosity
$\eta_{\mathrm{fl}}^{\mathrm{b}}$, compressibility $\kappa_{\mathrm{fl}}$, and
sound speed $c_{\mathrm{fl}}$, and damping coefficient $\Gamma_{\mathrm{fl}}$,
the acoustic pressure $p_{1}$ is governed by the weakly damped Helmholtz
equation, and the acoustic velocity $\bm{v}_{1}^{d}$ is derived from $p_{1}$,

$\displaystyle\tilde{p}_{1}$
$\displaystyle=p_{1}(\bm{r})\>\mathrm{e}^{-{\mathrm{i}\omega t}},\quad$
$\displaystyle\tilde{\bm{v}}_{1}^{d}$
$\displaystyle=\bm{v}_{1}^{d}(\bm{r})\>\mathrm{e}^{-{\mathrm{i}\omega t}},$
(14a) $\displaystyle\nabla^{2}p_{1}$
$\displaystyle=-\dfrac{\omega^{2}}{c_{\mathrm{fl}}^{2}}(1+\mathrm{i}\Gamma_{\mathrm{fl}})p_{1},\quad$
$\displaystyle\kappa_{\mathrm{fl}}$
$\displaystyle=(\rho_{\mathrm{fl}}c_{\mathrm{fl}}^{2})^{-1},$ (14b)
$\displaystyle\bm{v}_{1}^{d}$
$\displaystyle=-\mathrm{i}\dfrac{1-\mathrm{i}\Gamma_{\mathrm{fl}}}{\omega\rho_{\mathrm{fl}}}\bm{\nabla}p_{1},\quad$
$\displaystyle\Gamma_{\mathrm{fl}}$
$\displaystyle=\bigg{(}\frac{4}{3}\eta_{\mathrm{fl}}+\eta_{\mathrm{fl}}^{\mathrm{b}}\bigg{)}\omega\kappa_{\mathrm{fl}}.$
(14c)

In our $\phi$-independent axisymmetric case we have

$p_{1}(\bm{r})=p_{1}(r,z),\quad\bm{v}_{1}^{d}(\bm{r})=v_{1r}^{d}(r,z)\>\bm{e}_{r}+v_{1z}^{d}(r,z)\>\bm{e}_{z}.$
(15)

###  II.2 Second-order time-averaged fields

The nonlinearities in the fluid dynamics induce second-order terms in the
perturbation theory from products of first-order terms [3]. Here, we focus on
the time-averaged values of such fields,
$\bm{F}_{2}=\big{\langle}\tilde{\bm{F}}_{2}(\bm{r},t)\big{\rangle}=\dfrac{\omega}{2\pi}\int_{0}^{\frac{2\pi}{\omega}}\tilde{\bm{F}}_{2}(\bm{r},t)\>\mathrm{d}t$,
in particular the steady streaming velocity $\bm{v}_{2}$ and the acoustic
radiation force $\bm{F}^{\mathrm{rad}}$ on suspended particles. For the time-
averaged product of two first-order fields $A_{1}$ and $B_{1}$ we use the
well-known expression
$\big{\langle}\operatorname{Re}\big{[}\tilde{A}_{1}(\bm{r},t)\big{]}\operatorname{Re}\big{[}\tilde{B}_{1}(\bm{r},t)\big{]}\big{\rangle}=\frac{1}{2}\operatorname{Re}\big{[}A_{1}(\bm{r})B^{*}_{1}(\bm{r})\big{]}$.
The acoustic streaming $\bm{v}_{2}$ is an incompressible Stokes flow with a
non-negligible time-average Eckart-streaming bulk force

$\eta_{\mathrm{fl}}\nabla^{2}\bm{v}_{2}=\bm{\nabla}\mathcal{P}_{2}-\dfrac{\Gamma_{\mathrm{fl}}\omega}{2c_{\mathrm{fl}}^{2}}\operatorname{Re}\big{[}p_{1}^{*}\bm{v}_{1}^{d}\big{]},\quad\bm{\nabla}\cdot\bm{v}_{2}=0,$
(16)

where $\mathcal{P}_{2}$ is the pressure redefined to include the excessive
pressure [23]. The expression for the acoustic radiation force
$\bm{F}^{\mathrm{rad}}$ acting on a suspended particle of radius $a$, density
$\rho_{\mathrm{pt}}$, and compressibility $\kappa_{\mathrm{pt}}$ is [27],

$\displaystyle\bm{F}^{\mathrm{rad}}$ $\displaystyle=-\pi
a^{3}\bigg{\\{}\dfrac{2\kappa_{\mathrm{fl}}}{3}\operatorname{Re}\big{[}f_{0}^{*}p_{1}^{*}\bm{\nabla}p_{1}\big{]}-\rho_{\mathrm{fl}}\operatorname{Re}\big{[}f_{1}^{*}\bm{v}_{1}^{d*}\\!\cdot\\!\bm{\nabla}\bm{v}_{1}^{d}\big{]}\bigg{\\}},$
(17a) $\displaystyle f_{0}$
$\displaystyle=1-\dfrac{\kappa_{\mathrm{pt}}}{\kappa_{\mathrm{fl}}},\quad
f_{1}=\dfrac{2(1-\gamma)(\rho_{\mathrm{pt}}-\rho_{\mathrm{fl}})}{2\rho_{\mathrm{pt}}+(1-3\gamma)\rho_{\mathrm{fl}}},$
(17b)

where $f_{0}$ and $f_{1}$ are the monopole and dipole scattering coefficients,
$\delta$ is the width of the acoustic boundary layer, and $\gamma$ is a
correction factor,

$\delta=\sqrt{\frac{2\eta_{\mathrm{fl}}}{\omega\rho_{\mathrm{fl}}}},\quad\gamma=-\frac{3}{2}\bigg{[}1+\mathrm{i}\bigg{(}1+\frac{\delta}{a}\bigg{)}\bigg{]}\frac{\delta}{a}.$
(18)

The total force $\bm{F}^{\mathrm{pt}}$ acting on a suspended particle moving
with speed $\bm{v}_{\mathrm{pt}}$ at position $\bm{r}$ is the sum of the
radiation force $\bm{F}^{\mathrm{rad}}$, the Stokes drag force
$\bm{F}^{\mathrm{drag}}$, and the buoyancy-corrected gravitational force
$\bm{F}^{\mathrm{grav}}$ from the gravitational acceleration $\bm{g}$,

$\displaystyle\bm{F}^{\mathrm{pt}}$
$\displaystyle=\bm{F}^{\mathrm{rad}}+\bm{F}^{\mathrm{drag}}+\bm{F}^{\mathrm{grav}},$
(19a) $\displaystyle\bm{F}^{\mathrm{drag}}$
$\displaystyle=6\pi\eta_{\mathrm{fl}}a(\bm{v}_{2}-\bm{v}_{\mathrm{pt}}),\quad\bm{F}^{\mathrm{grav}}=\frac{4}{3}\pi
a^{3}(\rho_{\mathrm{pt}}-\rho_{\mathrm{fl}})\>\bm{g}.$ (19b)

###  II.3 Boundary conditions

In the following we state the boundary conditions that we apply to the first-
order acoustic and time-averaged second-order steady fields. On the axis
$r=0$, we have the usual axisymmetry conditions,

$\displaystyle\partial_{r}p_{1}$ $\displaystyle=0,\;\;$
$\displaystyle\partial_{r}u_{1z}$ $\displaystyle=0,\;\;$ $\displaystyle
u_{1r}$ $\displaystyle=0,\;\;$ $\displaystyle\partial_{r}\varphi_{1}$
$\displaystyle=0,$ (20a) $\displaystyle\partial_{r}p_{2}$
$\displaystyle=0,\;\;$ $\displaystyle\partial_{r}v_{2z}$
$\displaystyle=0,\;\;$ $\displaystyle v_{2r}$ $\displaystyle=0,$
$\displaystyle\text{at }\;r=0.$ (20b)

We control the electric potential $\varphi_{1}$ on the transducer electrodes.
At the ground electrode, $\varphi_{1}=0$ always, whereas the excited
electrodes have a peak-to-peak AC voltage amplitude of $\varphi_{0}$, so there
$\varphi=\frac{1}{2}\varphi_{0}\>\mathrm{e}^{\mathrm{i}\Delta}$ with $\Delta$
being a phase. Here, we work only with in-phase $\Delta=0$ (positive) and
anti-phase $\Delta=\pi$ (negative) excitation voltages, so at the electrode
surfaces we have

$\displaystyle\varphi_{1}$ $\displaystyle=0,\quad$ $\displaystyle\text{on
ground electrodes},$ (21a) $\displaystyle\varphi_{1}$
$\displaystyle=+\frac{1}{2}\varphi_{0},\quad$ $\displaystyle\text{on positive
electrodes},$ (21b) $\displaystyle\varphi_{1}$
$\displaystyle=-\frac{1}{2}\varphi_{0},,\quad$ $\displaystyle\text{on negative
electrodes}.$ (21c)

At interfaces away from the electrodes, we assume a zero-free-charge condition
on $\varphi_{1}$. On freely vibrating interfaces, we assume a zero-stress
condition on $\bm{u}_{1}$,

$\bm{D}\cdot\bm{n}=0,\quad\bm{\sigma}_{\mathrm{sl}}\cdot\bm{n}=\bm{0},\quad\text{at
the solid-air interface},$ (22)

At the solid-fluid interface, we use the effective boundary conditions for the
continuity of the first-order velocity and stress, developed in a coordinate-
free form by Bach and Bruus [23], where the acoustic boundary layer in the
fluid has been accounted for analytically. Here we state these boundary
conditions for our case with $\phi$-independent axisymmetric first-order
fields,

$\displaystyle\partial_{z}p_{1}$
$\displaystyle=\dfrac{\mathrm{i}\omega\rho_{0}}{1-\mathrm{i}\Gamma_{\mathrm{fl}}}\left[-\mathrm{i}\omega
u_{1z}-\dfrac{\mathrm{i}}{k_{\mathrm{s}}r}\partial_{r}\left(-\mathrm{i}\omega
ru_{1r}\right)\right]$
$\displaystyle\qquad-\dfrac{\mathrm{i}}{k_{\mathrm{s}}}\big{(}k_{\mathrm{c}}^{2}p_{1}+\partial_{z}^{2}p_{1}\big{)},$
(23a) $\displaystyle\bm{\sigma}_{\mathrm{sl}}\cdot\bm{e}_{z}$
$\displaystyle=-p_{1}\bm{e}_{z}+\mathrm{i}k_{\mathrm{s}}\eta_{0}\Big{(}-\mathrm{i}\omega\bm{u}_{1}+\dfrac{1}{\omega\rho_{0}}\bm{\nabla}p_{1}\Big{)},$
(23b) $\displaystyle k_{\mathrm{c}}^{2}$
$\displaystyle=(1+\mathrm{i}\Gamma_{\mathrm{fl}})\dfrac{\omega^{2}}{c_{\mathrm{fl}}^{2}},\qquad
k_{\mathrm{s}}=\dfrac{1+\mathrm{i}}{\delta}.$ (23c)

as well as the second-order streaming $\bm{v}_{2}$,

$\displaystyle v_{2r}^{d0}$
$\displaystyle=-\dfrac{1}{2\omega}\mathrm{Re}\bigg{\\{}\dfrac{1}{2}v_{1r}^{\delta
0\star}\partial_{r}v_{1r}^{\delta 0}+\dfrac{1}{2}v_{1z}^{\delta
0\star}\partial_{z}v_{1r}^{\delta 0}$
$\displaystyle\qquad\quad+\dfrac{2-\mathrm{i}}{2}\Big{(}\partial_{r}v_{1r}^{\delta
0\star}+\frac{1}{r}v_{1r}^{\delta 0\star}+\partial_{z}v_{1z}^{\delta
0\star}\Big{)}v_{1r}^{\delta 0}$
$\displaystyle\qquad\quad+\mathrm{i}\Big{(}\partial_{r}u^{0\star}_{1r}+\frac{1}{r}u^{0\star}_{1r}+\partial_{z}u^{0\star}_{1z}-\partial_{z}v_{1z}^{d\star}\Big{)}v_{1r}^{\delta
0}$ $\displaystyle\qquad\quad-\mathrm{i}\big{(}v_{1r}^{\delta
0\star}+v_{1z}^{\delta
0\star}\big{)}\partial_{z}u_{1r}-\mathrm{i}\bm{u}_{1}^{0\star}\cdot\bm{\nabla}v_{1r}^{d}\bigg{\\}},$
(24a) $\displaystyle v_{2z}^{d0}$
$\displaystyle=-\dfrac{1}{2\omega}\mathrm{Re}\Big{\\{}\mathrm{i}v_{1r}^{\delta
0\star}\partial_{r}u_{1z}^{0\star}+u_{1r}^{0\star}\partial_{r}\big{(}v_{1z}^{d}+v_{1z}^{\delta
0}\big{)}$
$\displaystyle\qquad\quad+u_{1z}^{0\star}\partial_{z}\big{(}v_{1z}^{d}+v_{1z}^{\delta
0}\big{)}\Big{\\}}$ (24b)

Here, $\bm{v}_{1}^{\delta 0}=\bm{u}_{1}^{0}-\bm{v}_{1}^{d0}$, which together
with all terms containing the factor $k_{\mathrm{s}}$ originates from the
boundary layer of width $\delta$.

Finally, when we consider systems without a lid, the acoustic waves
transmitted from the membrane into the fluid will propagate toward infinity
($r\rightarrow\infty$ and $z\rightarrow\infty$) until damped out. To keep a
finite-sized computational domain, we therefore follow Refs. [28, 22] and
limit the physical domain of the fluid to the region $0<r<r_{\mathrm{fl}}$ and
$0<z<z_{\mathrm{fl}}$, then add so-called perfectly matched layers (PML)
around this domain for $r_{\mathrm{fl}}<r<r_{\mathrm{PML}}$ or
$z_{\mathrm{fl}}<z<z_{\mathrm{PML}}$, in which all outgoing waves from the
membrane are damped out, and from which no incoming waves are sent toward the
membrane,

$\displaystyle\text{For }\;r_{\mathrm{fl}}<r<r_{\mathrm{PML}}\text{ or
}\;z_{\mathrm{fl}}<z<z_{\mathrm{PML}}:$ $\displaystyle\chi(r,z)=$ (25a)
$\displaystyle
K_{\mathrm{PML}}\bigg{(}\theta[r\\!-\\!r_{\mathrm{fl}}]\bigg{[}\dfrac{r\\!-\\!r_{\mathrm{fl}}}{r_{\mathrm{PML}}\\!-\\!r_{\mathrm{fl}}}\bigg{]}^{2}\\!+\theta[z\\!-\\!z_{\mathrm{fl}}]\bigg{[}\dfrac{z\\!-\\!z_{\mathrm{fl}}}{z_{\mathrm{PML}}\\!-\\!z_{\mathrm{fl}}}\bigg{]}^{2}\bigg{)},$
$\displaystyle\partial_{r}\rightarrow[1+\mathrm{i}\chi(r,z)]\partial_{r},\quad\partial_{z}\rightarrow[1+\mathrm{i}\chi(r,z)]\partial_{z}.$
(25b)
$\displaystyle\mathrm{d}r\rightarrow\frac{\mathrm{d}r}{1+\mathrm{i}\chi(r,z)},\qquad\mathrm{d}x\rightarrow\frac{\mathrm{d}z}{1+\mathrm{i}\chi(r,z)}.$
(25c)

Here, $\theta$ is the Heaviside step function, and $K_{\mathrm{PML}}$ is a
constant chosen, such that the outgoing waves are damped in the PML without
reflections at $r=r_{\mathrm{fl}}$ and $z=z_{\mathrm{fl}}$.

The choice of boundary conditions at the outer surface of the PML are not
crucial due to the strong damping provided in the PML. We use the Dirichlet
condition $p_{1}=0$ there, but changing to Neumann condition
$\bm{n}\cdot\bm{\nabla}p_{1}=0$ leads to negligible relative deviation in the
result, less than $10^{-5}$ measured in terms of the L2-norm.

###  II.4 Numerical implementation in COMSOL

As in Ref. [3], we implement the model in the commercial finite-element
software COMSOL 5.5 [21], closely following the implementation method
described in Ref. [22]. The fields $p_{1}$, $\bm{u}_{1}$, $\varphi_{1}$,
$p_{2}$, and $\bm{v}_{2}$ are all computed with quartic-order polynomial test
functions. The mesh is defined to ensure at least 20 mesh elements per
wavelength, and typically more, resulting in more than 80 nodal points per
wavelength in the physical domain, and 10 elements per wavelength in the PML
domain. In the bulk fluid we use a triangular mesh, and in both the membrane
and the thin-film transducer, a structured mesh is used with at least four
mesh elements in the thickness direction. At the symmetry axis $r=0$ is added
a boundary layer of 15 elements, growing up from a smallest element of 0.005
times the bulk mesh size. This mesh ensured a successful mesh-convergence
analysis similar to that presented in Ref. [29].

Table 1:  Material parameter values used in the model. The iodixanol solution is computed from the interpolation polynomials given in Ref. [30]. For the MCF-7 cells, the radius is chosen as the center value $10~{}\textrm{\textmu{}m}$ in the observed range of 8.1-12.2 µm given in Ref. [31], and the scattering coefficients $f_{0}$ and $f_{1}$ are computed from Eq. (17b) using data from Ref. [32]. Parameter | Value |  | Parameter | Value
---|---|---|---|---
_Thin-film aluminum scandium nitride_ ,
$\textrm{Al}_{0.6}\textrm{Sc}_{0.4}\textrm{N}$ [33, 34]
$\rho_{\mathrm{sl}}$ | 3300 $\textrm{kg}\>\textrm{m}^{-3}$ |  | $\Gamma_{\mathrm{sl}}$ | 0.0005
$C_{11}$ | 313.8 GPa |  | $C_{33}$ | 197.1 GPa
$C_{12}$ | 150.0 GPa |  | $C_{44}$ | 108.6 GPa
$C_{13}$ | 139.2 GPa |  | $C_{66}$ | 81.9 GPa
$e_{31,f}$ | $-2.65~{}\textrm{C}\>\textrm{m}^{-2}$ |  | $e_{15}$ | $-0.32~{}\textrm{C}\>\textrm{m}^{-2}$
$e_{33}$ | 2.73 $\textrm{C}\>\textrm{m}^{-2}$ |  | $\Gamma_{\varepsilon}$ | $0.0005$
$\varepsilon_{11}$ | 22 $\varepsilon_{0}$ |  | $\varepsilon_{33}$ | 22 $\varepsilon_{0}$
_Membrane silicon, Si (111)_ [25, 35, 36]
$\rho_{\mathrm{sl}}$ | 2329 $\textrm{kg}\>\textrm{m}^{-3}$ |  |  |
$E$ | 168.9 GPa |  | $s$ | 0.262
$C_{11}$ | 207.5 GPa |  | $C_{44}$ | 73.7 GPa
$C_{12}$ | 66.9 GPa |  | $\Gamma_{\mathrm{sl}}$ | 0.0001
$c_{\mathrm{lo}}$ | 5594 $\textrm{m}\,\textrm{s}^{-1}$ |  | $c_{\mathrm{tr}}$ | 3425 $\textrm{m}\,\textrm{s}^{-1}$
_Water_ [29]
$\rho_{\mathrm{fl}}$ | $997~{}\textrm{kg}\>\textrm{m${}^{-3}$}$ |  | $\eta_{\mathrm{fl}}$ |  $0.890~{}\textrm{mPa}\>\textrm{s}$
$c_{\mathrm{fl}}$ | $1497~{}\textrm{m}\,\textrm{s}^{-1}$ |  | $\eta^{{\mathrm{b}}}_{\mathrm{fl}}$ | $2.485~{}\textrm{mPa}\>\textrm{s}$
$\kappa_{\mathrm{fl}}$ | $448~{}\textrm{TPa}^{-1}$ |  | $\Gamma_{\mathrm{fl}}$ | $10.3~{}\mathrm{THz}^{-1}f$
_Iodixanol 60% solution_ [30]
$\rho_{\mathrm{fl}}^{\mathrm{Idx,}60\%}$ | $1320~{}\textrm{kg}\>\textrm{m${}^{-3}$}$ |  | $c_{\mathrm{fl}}^{\mathrm{Idx,}60\%}$ |  $1498~{}\textrm{m}\,\textrm{s}^{-1}$
$\eta_{\mathrm{fl}}^{\mathrm{Idx,}60\%}$ | $7.690~{}\textrm{mPa}\>\textrm{s}$ |  | $\kappa_{\mathrm{fl}}^{\mathrm{Idx,}60\%}$ | $338~{}\textrm{TPa}^{-1}$
_Cancer cell MCF-7_ [31, 32]
$\rho_{\mathrm{MCF\text{-}7}}$ | $1055~{}\textrm{kg}\>\textrm{m${}^{-3}$}$ |  | $\rho_{\mathrm{MCF\text{-}7}}$ |  $373~{}\textrm{TPa}^{-1}$
$a_{\mathrm{MCF\text{-}7}}$ | $10~{}\textrm{\textmu{}m}$ |  |  |
$f_{0}^{\mathrm{Wa}}$ | $0.167$ |  | $f_{1}^{\mathrm{Wa}}$ | $0.037+0.00002\mathrm{i}$
$f_{0}^{\mathrm{Idx,}60\%}$ | $-0.104$ |  | $f_{1}^{\mathrm{Idx,}60\%}$ | $-0.153+0.0010\mathrm{i}$
Figure 2:  Simulation of the resonance modes $n$ for a 10-µm-thick silicon Si
(111) membrane of radius $r_{\mathrm{mv}}=100~{}\textrm{\textmu{}m}$ placed
below a half-space filled with a fluid and actuated by an attached 1-µm-thick
$\textrm{Al}_{0.6}\textrm{Sc}_{0.4}\textrm{N}$ thin-film PZE transducer with a
peak-to-peak AC voltage $\varphi_{0}=1~{}\textrm{V}_{\textrm{pp}}$. The fluid
domain is surrounded by a 100-µm wide PML domain. (a) Color plot of the
pressure $\operatorname{Re}[\mathrm{i}p_{1}]$ from $-7.8$ (blue) to
$+7.8~{}\textrm{kPa}$ (red) and of the displacement amplitude $|\bm{u}_{1}|$
from 0 (blue) to $0.75~{}\textrm{nm}$ (yellow) for the $n=0$ fundamental mode
at 2.16 MHz, excited by the uniform excitation electrode shown in (d) the
$r$-$z$ plane and (g) the $x$-$y$ plane. (b) The same as (a) but for the $n=3$
harmonic mode at 52.0 MHz with pressure amplitude $\pm 14~{}\textrm{kPa}$ and
displacement amplitude 0.08 nm, excited by the uniform excitation electrode
shown in (e) the $r$-$z$ plane and (i) the $x$-$y$ plane. (c) The same as (b)
but with with pressure amplitude $\pm 513~{}\textrm{kPa}$ and displacement
amplitude 0.75 nm, and with the patterned excitation electrode shown in (f)
the $r$-$z$ plane and (k) the $x$-$y$ plane. (h), (j), and (l) Color plot of
the out-of-phase in-plane strain
$\operatorname{Re}[\mathrm{i}\partial_{r}u_{1r}]$ from min (cyan) to max
(magenta) for the membrane modes shown respectively in (a), (b), and (c) for
which the vertical displacement $\operatorname{Re}[\mathrm{i}u_{1z}]$ for
clarity has been enhanced by a factor 15000, 30000, and 6000, respectively.
Animations of panels (a) and (c) are given in the Supplemental Material [37].

The materials used in the model are as follows: the PZE thin-film transducer
is $\textrm{Al}_{0.6}\textrm{Sc}_{0.4}\textrm{N}$, the silicon membrane is
Si-(111), the liquid is either pure water or a 60% aqueous solution of
iodixanol, and the cell is a MCF-7 breast-cancer cell. All material parameter
values used in the model are listed in Table 1 including references to the
literature.

We carried out the numerical simulations of the model on a workstation with a
16-cores Intel i9-7960X processor at 3.70 GHz boost clock frequency and 128 GB
of random access memory (RAM). At any given frequency, a simulation of the
first-order fields $p_{1}$, $\bm{u}_{1}$, and $\varphi_{1}$ comprised of 1.4
million degrees of freedom, used 14 GB RAM, and took 50 s, whereas the second-
order fields $p_{2}$ and $\bm{v}_{2}$ comprised of 2.7 million degrees of
freedom, used 27 GB RAM, and took 180 s.

##  III Excitation of higher-harmonic membrane modes by electrode patterning

We now introduce the excitation of higher-harmonic membrane modes in
acoustofluidics. As already indicated in Fig. 1, we show that only by
patterning the excitation electrode to be compatible with the shape of the
target higher-harmonic membrane mode, this mode can be excited with a
sufficiently large amplitude.

###  III.1 Membrane modes and patterned electrodes

In Fig. 2 we show the main effect of patterning the excitation electrodes
appropriately for a 10-µm-thick Si-disk membrane covered uniformly by a 1-µm-
thick (AlSc)N thin-film transducer. The inner part
$r<r_{\mathrm{mv}}=100~{}\textrm{\textmu{}m}$ of the membrane is free to
vibrate, whereas the surrounding ring
$r_{\mathrm{mv}}<r<r_{\mathrm{mc}}=300~{}\textrm{\textmu{}m}$ is clamped.
Above the membrane is placed a cylindrical domain of water of radius
$r_{\mathrm{mc}}$ and height $z_{\mathrm{fl}}=200~{}\textrm{\textmu{}m}$, and
outside this domain is placed a PML domain of thickness
$100~{}\textrm{\textmu{}m}$.

In Fig. 2(a), the ground and excitation electrode are fully covering the
transducer, and a peak-to-peak voltage $\varphi_{0}=1~{}\textrm{V}$ is applied
to the latter. A frequency sweep reveals the $\phi$-independent vibration
resonance modes $n$ of the membrane disk, $n=0,1,2,3,\ldots$, with resonance
frequencies $f_{n}=2.16,10.5,19.7,52.0$ MHz and a monotonically decreasing
maximum displacement amplitude $u_{1z,n}^{\mathrm{max}}=0.75,0.24,0.08,0.02$
nm for $n=0,1,2,3$. Animations of modes $n=0$ and $n=3$ are given in the
Supplemental Material 111See Supplemental Material at http://bruus-
lab.dk/files/Steckel_membrane_Suppl.zip for animations of $p_{1}$ and
$\bm{u}_{1}$ in Fig. 2(a), Fig. 2(c), Fig. 3(a), Fig. 4, and Fig. 5..

This decreasing mode amplitude is partly explained by the uniformity across
the thin-film-transducer surface of the applied perpendicular voltage drop
$\varphi_{0}$, see Fig. 2(d,g), which by the electromechanical coupling matrix
$\bm{e}$ promotes a uniform stretching or compression of the thin film at any
given time. This works well for the fundamental mode $n=0$ at 2.16 MHz, where
the film is either stretching across the entire membrane at any given time, or
contracting, although not uniformly, resulting in a perpendicular displacement
$u_{1z}$ that is either positive across the entire membrane, or negative, as
shown in Fig. 2(a). For the third harmonic mode $n=3$ at 52.0 MHz, the
oscillating displacement $u_{1z}$ changes sign 3 times along the radial
direction as shown in Fig. 2(b), and this counteracts uniform
stretching/contraction promoted by the uniform potential shown in Fig. 2(e,i),
and as a result the excitation amplitude is decreased. Clearly, by arranging
for a non-uniform potential following the non-uniform mechanical displacement,
we could hope to restore a large excitation amplitude. But how?

The clue is the out-of-phase in-plane strain $\mathrm{i}\partial_{r}u_{1r}$
for mode $n=0$ and 3 shown in Fig. 2(h,j), respectively. This field indicates
how the membrane would move were it free to vibrate in a given resonance mode.
For $n=0$ at 2.16 MHz, the strain indicates a contraction (cyan) in a large
center part of the membrane, which the uniform applied potential supports,
only counteracting the expansion indicated by the strain in a minor peripheral
part (magenta). In contrast, for mode $n=3$ at 52 MHz, the strain indicates
four ring-shaped domains respectively exhibiting contraction, expansion,
contraction, expansion, which is not compatible with the applied uniform
potential, as each full contraction-expansion period is nearly canceled by the
potential.

To make the potential for $n=3$ compatible with the strain pattern, we
therefore split the excitation electrode into the same 3 ring-shaped domains
(plus the central disk) defined by the strain with alternating excitation
voltages $+\frac{1}{2}\varphi_{0}$, $-\frac{1}{2}\varphi_{0}$,
$+\frac{1}{2}\varphi_{0}$, and $-\frac{1}{2}\varphi_{0}$, shown in Fig.
2(f,k). If these rings should be contacted by wires in the same layer [20] or
by conducting vias to wires in a second layer, we leave to the experts of MEMS
technology to figure out, here we just assume it possible in practise. The
resulting acoustic response for the same applied voltage amplitude is shown in
Fig. 2(c). Here we see that the resonance frequency remains $f_{3}=52.0$ MHz,
and that the displacement amplitude for the patterned electrode is
$u_{1z,3}^{\mathrm{max,p}}=0.75$ nm, one and half order of magnitude higher
than that of the uniform electrode, $u_{1z,3}^{\mathrm{max,u}}=0.02$ nm, and
equal to the uniform-electrode mode-0 amplitude
$u_{1z,0}^{\mathrm{max,u}}=0.76$ nm.

###  III.2 Traveling and standing pressure waves induced by higher-harmonic
membrane modes

Let us now turn to the resulting acoustic pressure field $p_{1}$ in the liquid
above the membrane vibrating in its resonance modes $n=0$ and 3. In Fig. 2(a)
we see that a strong pressure wave is emitted from the antinode part of the
vibrating membrane at 2.16 MHz in mode $n=0$ and absorbed in the PML domain.
In Fig. 2(b), and more clearly in Fig. 2(c), we also see partial waves being
emitted from each of the four antinodes in mode $n=3$ at 52.0 MHz. There is an
important qualitative difference between the fundamental mode $n=0$ and any
higher-harmonic mode $n$. The former emits only one partial wave, and thus
$p_{1}$ is devoid of interference patterns, whereas the latter emits $n+1$
partial waves from $n+1$ ring-shaped parts of the membrane surface. These
partial waves interfere and give rise to a complex interference pattern
stretching out in the bulk of the liquid above the membrane, such as seen in
Fig. 2(c). This has profound implications for the acoustofluidic properties of
the system, which we investigate further in Section IV.

Figure 3:  (a) 3D contour plot of instantaneous isobars at $-250$, $-125$,
$125$, and $250~{}\textrm{kPa}$ equivalent to Fig. 2(c) of the acoustic
pressure $p_{1}$ in the liquid above the membrane excited to resonance mode
$n=3$ at 52.0 MHz by the patterned 3-ring electrode shown in Fig. 2(f,k). The
displacement amplitude is enhanced by a factor 20,000 to be visible. An
animation of $p_{1}$ and $\bm{u}_{1}$ is given in the Supplemental Material
[37]. (b) A color plot similar to Fig. 2(c), but for the absolute value
$|p_{1}|$.

In Fig. 3(a) is shown a 3D contour plot of instantaneous isobars equivalent to
Fig. 2(c) of the acoustic pressure $p_{1}$ in the liquid above the membrane
excited to resonance mode $n=3$ at 52.0 MHz by the patterned electrode shown
in Fig. 2(f,k). The pressure pattern consists of ring-shaped waves propagating
upward in the $z$ direction, while exhibiting standing waves in the radial
direction.

This particular mixture of standing radial and traveling axial waves is
clearly revealed by studying the color plot Fig. 3(b) of the absolute value
$|p_{1}|$ of the acoustic pressure. Following a outward radial path from $r=0$
at any given fixed height $z=z_{0}$, reveals regular oscillations in
$|p_{1}(r,z_{0})|$ between nearly zero minima and large maxima, a
characteristic of a standing radial wave. The most pronounced oscillation
happens at height $z_{0}=60~{}\textrm{\textmu{}m}$, and here we find the ratio
$\max|p_{1}|/\min|p_{1}|\approx 140$, which corresponds to a nearly ideal
standing wave.

When inspecting the pressure nodal surfaces (gray) in $|p_{1}|$ in Fig. 3(b),
we see that they form two nearly perfect and one deformed coaxial cylinder
centered on the $z$-axis, but with a tendency to widen out as the distance to
the membrane increases. Following any vertical path upward parallel to the
$z$-axis, $|p_{1}|$ is mainly monotonically decreasing, which is the
characteristic of a nearly ideal damped traveling wave. One exception to this
behavior is found along the $z$-axis, where a local maximum in $|p_{1}(0,z)|$
is found for $z\approx 60~{}\textrm{\textmu{}m}$ above the center of the
membrane. Here we estimate the ratio $\max|p_{1}|/\min|p_{1}|\approx 1.3$,
which corresponds to a 1.0:0.3 mixture of a traveling and a standing wave. We
shall see in Section IV that this small component of a standing wave along $z$
together with the large standing wave component along $r$ spawns a particle
trap in all three spatial directions at some distance above the center of the
vibrating silicon membrane.

###  III.3 The effect of introducing a rigid lid

So far, we have assumed ideally absorbing surroundings by introducing the PML.
As a result the pressure waves are traveling away from the membrane without
being reflected back. In actual acoustofluidic systems, such a behavior is
approximately realized by using soft PDMS rubber walls, as demonstrated
recently by Skov et al. in a combined experimental and numerical study [38].
They also showed how replacing the soft lid with the other extreme, a hard
glass lid, the traveling wave field was replaced by a standing wave field.
Here, we briefly examine this latter case, by replacing the top PML in our
model by a rigid lid placed at $z=z_{\mathrm{fl}}$ with a hard-wall boundary
condition, $\partial_{z}p_{1}=0$, for the acoustic pressure. We keep the PML
at the side to mimic an infinitely broad system.

The results of the rigid-lid simulation of $p_{1}$ is shown in Fig. 4, where a
vertical standing-wave behavior is clearly seen as horizontal nodal lines in
the color plot of $|p_{1}|$ in Fig. 4(b). However, since the PML at the side
still allows for the waves to propagate away from the membrane, now
predominantly in the radial direction, neither the membrane modes nor the
maximum pressure amplitudes are much affected by replacing the PML lid with
the rigid lid, as seen when comparing Figs. 3(b) and 4(b). The smooth behavior
of $|p(0,z)|$ is now overlaid with minor oscillations, but there is still a
global pressure maximum along the $z$-axis.

Close to the membrane, the admixture of the vertical standing wave component
is strong for the rigid-lid system, but still traveling waves are seen, as
revealed by the animations of Figs. 2(c) and 4(b) given in the Supplemental
Material [37]. As the characteristics of the pressure field are not strongly
affected by the lid boundary condition, we continue using the previous case of
an ideally absorbing lid in the following.

Figure 4:  The small membrane of Figs. 2 and 3 but with a hard-wall lid at
$z=z_{\mathrm{fl}}$ instead of PML. (a) Same as Fig. 2(c) for the pressure
$p_{1}$. An animation of $p_{1}$ and $\bm{u}_{1}$ is given in the Supplemental
Material [37]. (b) Same as Fig. 3(b) for the magnitude $|p_{1}|$ of the
pressure $p_{1}$.

##  IV Cell trapping by higher-harmonic membrane modes

We choose to illustrate the possibilities offered by the use of higher-
harmonic membrane modes in acoustofluidics application with the example of
trapping of a single suspended biological cell. Particle trapping is a central
problem in the field studied by many groups. Examples are microparticle
trapping in acoustic tweezers [39, 40, 41, 42, 43, 44] and in disposable
capillary tubes[45, 28], nanoparticle by using seed particles [46, 47, 48],
and short and long term trapping[49, 50]. Trapping and the associated focusing
has been created by various methods, including standing bulk [51, 52],
traveling bulk [53], and surface [54, 55, 38] acoustic waves. An important
focus area in the field is the trapping and focusing of biological cells in
general [56, 57, 58, 50] and of circulating tumor cells in particular [59, 60,
30], in some cases accomplished by tuning the acoustic properties of the
suspension medium tune relative to those of the given cells [30, 61].

Figure 5:  A membrane device as in Fig. 2(c), but with
$r_{\mathrm{mv}}=500~{}\textrm{\textmu{}m}$,
$r_{\mathrm{mc}}=1500~{}\textrm{\textmu{}m}$,
$z_{\mathrm{fl}}=1000~{}\textrm{\textmu{}m}$, and with the excitation
electrode is divided into 10 (instead of 3) ring-shaped segments each
$44~{}\textrm{\textmu{}m}$ wide, separated by 4-µm-wide gaps (not shown), and
excited with alternating voltages $\pm\frac{1}{2}\varphi_{0}$ as in Fig.
2(f,k). Color plots of the pressure $\operatorname{Re}[\mathrm{i}p_{1}]$ from
$-235$ (blue) to $+235~{}\textrm{kPa}$ (red) and of $|\bm{u}_{1}|$ from 0
(blue) to $0.49~{}\textrm{nm}$ (yellow) for the $n=10$ higher-harmonic mode at
$18.5~{}\textrm{MHz}$. The top inset is a color plot of $|p_{1}|$. The bottom
inset shows details of $p_{1}$ and $|\bm{u}_{1}|$. An animation of $p_{1}$ and
$\bm{u}_{1}$ is given in the Supplemental Material [37].

###  IV.1 Design of the membrane for cell trapping

In the following numerical simulation analysis of cell trapping, we choose as
our model cell the breast cancer cell MCF-7 with the known acoustic parameters
listed in Table 1 and here assumed to be spherical with a radius
$a_{\mathrm{MCF\text{-}7}}=10~{}\textrm{\textmu{}m}$. To facilitate acoustic
trapping of such a cell, we prefer to work with an acoustic wavelength in the
fluid larger than the cell, $\lambda_{0}\gtrsim
8a_{\mathrm{MCF\text{-}7}}\approx 80~{}\textrm{\textmu{}m}$, or $f_{0}\lesssim
20~{}\textrm{MHz}$, a limit in which the long-wavelength limit expression
(17a) of the acoustic radiation force on a suspended particle is valid.
Furthermore, for membrane mode $n$ with radial wavelength
$\lambda_{r}^{(n)}\approx 2r_{\mathrm{mv}}/(n+\frac{1}{2})$, the mode
frequency scales as $f^{(n)}\propto\big{[}\lambda_{r}^{(n)}\big{]}^{-2}$, and
since $f^{(n)}\propto\lambda_{0}^{-1}$ in the fluid, we can at membrane
resonance write $\big{[}\lambda_{r}^{(n)}\big{]}^{2}=\lambda_{0}L$, where $L$
is a parameter of dimension length characteristic for the given system. We
find numerically that $L\approx 110~{}\textrm{\textmu{}m}$ in our system with
a 10-µm-thick silicon-membrane pushing on water. Combined with the long-
wavelength criterion, we obtain
$\lambda^{(n)}_{r}=\sqrt{\lambda_{0}L}\gtrsim\sqrt{8a_{\mathrm{MCF\text{-}7}}L}\approx
94~{}\textrm{\textmu{}m}$. Finally, to ensure the formation of a trapping
point on the $z$-axis, we must demand that
$\lambda_{0}<\alpha\lambda_{r}^{(n)}$ (where $\alpha\approx 1$ is a constant
we have not been able to compute) because if not, the acoustic wave in the
fluid propagates at an angle away from the $z$-axis. We thus conclude that all
above criteria are satisfied if
$8a_{\mathrm{MCF\text{-}7}}\lesssim\lambda_{0}\lesssim\alpha^{2}L$ (equivalent
to $12~{}\textrm{MHz}\lesssim f_{0}\lesssim 20\alpha^{2}~{}\textrm{MHz}$) and
$\lambda_{r}^{(n)}\gtrsim 94~{}\textrm{\textmu{}m}$.

Consequently, as a proof of concept, we use the design shown in Fig. 5,
consisting of a 10-µm-thick silicon Si (111) disk membrane of radius
$r_{\mathrm{mv}}=500~{}\textrm{\textmu{}m}$ driven in its $n=10$ harmonic mode
with $\lambda_{r}^{(10)}\approx~{}95~{}\textrm{\textmu{}m}$ and an estimated
resonance frequency $f^{(10)}\approx
Lc_{0}/\big{[}\lambda_{r}^{(10)}\big{]}^{2}\approx 18.2~{}\textrm{MHz}$.
Similar to Fig. 2(f,k) of mode $n=3$, the mode $n=10$ is excited by patterning
the excitation electrode of the 1-µm-thick
$\textrm{Al}_{0.6}\textrm{Sc}_{0.4}\textrm{N}$ transducer into ten ring-shaped
segments. A $1~{}\textrm{V}_{\textrm{pp}}$ AC voltage is applied to the these
segments with alternating phases $+\frac{1}{2}\varphi_{0}$ and
$-\frac{1}{2}\varphi_{0}$, resulting in a strong excitation of resonance mode
$n=10$ at $f=18.5~{}\textrm{MHz}$, close the prediction $f^{(10)}$. All the
features seen in Fig. 5 of the disk-membrane displacement field $\bm{u}_{1}$
and of the pressure field $p_{1}$ in the liquid are the same as in Fig. 2(c),
and the amplitudes are relatively large being
$|\bm{u}_{1}|=0.49~{}\textrm{nm}$ and $|p_{1}|=235~{}\textrm{kPa}$. Of
particular interest for trapping, we notice a local maximum in $|p_{1}(0,z)|$
in the inset of Fig. 5 near $z=100~{}\textrm{\textmu{}m}$, similar to the one
at $z=60~{}\textrm{\textmu{}m}$ shown in Fig. 3(b).

###  IV.2 Characterizing the cell trap

Figure 6:  Contour plots of the simulated components $F^{\mathrm{pt}}_{r}$ and
$F^{\mathrm{pt}}_{z}$ of the force $\bm{F}^{\mathrm{pt}}$ acting on MCF-7
cells in the system shown in Fig. 5, but for different liquids and voltage
amplitudes: (a,b) Water at $1~{}\textrm{V}_{\textrm{pp}}$, (c,d) 60% iodixanol
at $1~{}\textrm{V}_{\textrm{pp}}$, and (e,f) 60% iodixanol at
$5~{}\textrm{V}_{\textrm{pp}}$, using the parameters in Table 1. The contour
lines for $F^{\mathrm{pt}}_{r}$ are from $-200$ to $100~{}\textrm{pN}$ in
steps of $50~{}\textrm{pN}$ in (a,c), but from $-5$ to $2.5~{}\textrm{nN}$ in
steps of $1.25~{}\textrm{nN}$ in (e). The contour lines for
$F^{\mathrm{pt}}_{z}$ are from $-20$ to $50~{}\textrm{pN}$ in steps of
$5~{}\textrm{pN}$ in (b,d), but $-0.5$ to $1.25~{}\textrm{nN}$ in steps of
$0.125~{}\textrm{nN}$ in (f). The black arrows indicate the components
direction. (g) Line plot of the normalized vertical force component
$\hat{F}_{z}^{\mathrm{pt}}$ along the $z$-axis for case (b), (d), and (f).

To characterize the cell trap, we study the total force $\bm{F}^{\mathrm{pt}}$
(19a) that acts on a MCF-7 cell suspended in water. In Fig. 6(a,b) is shown
the contour plots of the radial and axial components $F^{\mathrm{pt}}_{r}$ and
$F^{\mathrm{pt}}_{z}$ of $\bm{F}^{\mathrm{pt}}$. To set the scale, we note
that in this case, the buoyancy-corrected gravitational force is
$F^{\mathrm{grav}}_{\mathrm{wa}}=2.3~{}\textrm{pN}$ and the Stokes drag force
is
$F^{\mathrm{drag}}_{\mathrm{wa}}(10~{}\textrm{\textmu{}m}/\textrm{s})=1.7~{}\textrm{pN}$
assuming that the velocity of the cell relative to water is
$10~{}\textrm{\textmu{}m}/\textrm{s}$. Firstly, we note that
$\bm{F}^{\mathrm{pt}}$ with a magnitude above $70~{}\textrm{pN}$ is completely
dominated by the radiation force $\bm{F}^{\mathrm{rad}}$, secondly that
$F^{\mathrm{pt}}_{r}\gtrsim F^{\mathrm{pt}}_{z}$, and thirdly that using pure
water, $\bm{F}^{\mathrm{pt}}$ expels the cell from the primary nodal plane,
instead of trapping it, as indicated by the black arrows.

This anti-trapping is a well-known result in trapping theory for a particle
more heavy and rigid than the fluid in a single traveling wave [44]. To change
the system into a trap, we need to tune the acoustic properties of the fluid
to reverse the signs of the scattering coefficients $f_{0}$ and $f_{1}$. This
can be done by using the density modifier iodixanol as demonstrated in Refs.
[30, 62, 63, 64].

In Fig. 6(c,d) is shown that indeed by using a 60% aqueous iodixanol solution
that is heavier than pure water, the particle force $\bm{F}^{\mathrm{pt}}$ is
reversed: the orange regions of $F^{\mathrm{pt}}_{r}>0$ in Fig. 6(a) turns
into the cyan regions of $F^{\mathrm{pt}}_{r}<0$ in Fig. 6(c), and vice versa.
Also, the yellow regions of $F^{\mathrm{pt}}_{z}>0$ in Fig. 6(b) turns into
the magenta regions of $F^{\mathrm{pt}}_{z}<0$ in Fig. 6(d), and vice versa.
The characteristic forces for the solution are
$F^{\mathrm{grav}}_{I60\%}=-11~{}\textrm{pN}$ and
$F^{\mathrm{drag}}_{I60\%}(10~{}\textrm{\textmu{}m}/\textrm{s})=15~{}\textrm{pN}$.
In Fig. 6(c) we see that radial force $F^{\mathrm{pt}}_{r}$ points toward the
$z$ axis for $r_{\mathrm{mv}}\lesssim 20~{}\textrm{\textmu{}m}$, reaching a
maximum amplitude of 234 pN at the height $z_{\mathrm{trap}}\approx
70~{}\textrm{\textmu{}m}$. However, the axial component $F^{\mathrm{pt}}_{z}$
is an order of magnitude smaller, and when plotting $F^{\mathrm{pt}}_{z}(0,z)$
along the $z$-axis, we find that although its has the desired upward direction
$F^{\mathrm{pt}}_{z}(0,z)>0$ for $z<z_{\mathrm{trap}}$, it is still positive
at its minimum at $z=z_{\mathrm{trap}}$,
$F^{\mathrm{pt}}_{z}(0,z_{\mathrm{trap}})\approx+5~{}\textrm{pN}$. The setup
is thus not a trap at all, because as the numerical simulation reveals, at
$z_{\mathrm{trap}}$ we have $F_{z}^{\mathrm{rad}}=-20~{}\textrm{pN}$,
$F_{z}^{\mathrm{drag}}=+14~{}\textrm{pN}$, and
$F_{z}^{\mathrm{grav}}=+11~{}\textrm{pN}$, such that
$F^{\mathrm{pt}}_{z}=(-20+14)~{}\textrm{pN}+11~{}\textrm{pN}=5~{}\textrm{pN}>0~{}\textrm{pN}$.

A trap is achieved simply by enhancing the AC voltage amplitude from the
$1~{}\textrm{V}_{\textrm{pp}}$ to $5~{}\textrm{V}_{\textrm{pp}}$. This
enhances the magnitude of the voltage-dependent second-order forces
$\bm{F}^{\mathrm{rad}}$ and $\bm{F}^{\mathrm{drag}}$ by a factor $5^{2}=25$,
leaving the gravitation force unchanged. Changing from
$1~{}\textrm{V}_{\textrm{pp}}$ to $5~{}\textrm{V}_{\textrm{pp}}$, we thus
expect the weakest point in the trap to attain the value
$\bm{F}^{\mathrm{pt}}=\big{(}-234\times 25,(-20+14)\times
25+11\big{)}~{}\textrm{pN}=(-5850,-119)~{}\textrm{pN}$, and this is exactly
what is obtained, as shown in Fig. 6(e,f). So using an appropriate tuning of
density of the liquid and of the excitation voltage, our system is able to
trap the MCF-7 cell in a trapping point located on the $z$-axis approximately
$70~{}\textrm{\textmu{}m}$ above the membrane as revealed by the red thick
line in Fig. 6(g).

As shown in Fig. 6(e), radial trapping occurs once the cell enters the cyan
region where $F^{\mathrm{pt}}_{r}<0$. However, we notice that this region is
surrounded by barrier-like orange region with $F^{\mathrm{pt}}_{r}>0$, which
repels the cell. This barrier could make it difficult for the cell to enter
the trapping region in the first place, however once trapped, the barrier
prevents more particles to enter the trap. This aspect is discussed further in
Section V. To gain a better understanding of how the trap is loaded, we study
in the steady acoustic streaming in the following.

###  IV.3 Steady acoustic streaming in the trap

Figure 7:  Vector plot (cyan arrows) of the steady streaming $\bm{v}_{2}$ and
a color plot of $|\bm{v}_{2}|$ from 0 (black) to
$247~{}\textrm{\textmu{}m}/\textrm{s}$ (white) in a 60% iodixanol solution at
$f=18.5~{}\textrm{MHz}$ and $5~{}\textrm{V}_{\textrm{pp}}$ corresponding to
Fig. 6(e,f) of the membrane driven in mode $n=10$ . Bulk-driven Eckart
streaming dominates over boundary-driven Rayleigh streaming.

In the system with an acoustic wave propagating upward from a vibrating
membrane, the acoustic streaming is dominated by the bulk-driven Eckart
streaming due to the force density (16)
$\bm{f}=-\dfrac{\Gamma_{\mathrm{fl}}\omega}{2c_{\mathrm{fl}}^{2}}\operatorname{Re}\big{[}p_{1}^{*}\bm{v}_{1}^{d}\big{]}$
from the upward traveling acoustic wave, whereas the Rayleigh boundary-driven
streaming is negligible. In Fig. 7 is shown the simulated streaming for the
60% iodixanol solution at 18.5 MHz and $5~{}\textrm{V}_{\textrm{pp}}$
corresponding to Fig. 6(e,f). The dominant feature is the toroidal vortex
centered at $(r,z)=(200,140)~{}\textrm{\textmu{}m}$, which drags in particles
down toward the membrane near $r=300~{}\textrm{\textmu{}m}$ and then sends
them upward along the lower part of the $z$-axis, exactly where the trapping
point is located. The flow velocity there is
$247~{}\textrm{\textmu{}m}/\textrm{s}$ corresponding to a Stokes drag force on
a trapped cell of $F^{\mathrm{drag}}=358~{}\textrm{pN}$, which is a factor of
1.4 lower than the dominant radiation force,
$F_{z}^{\mathrm{rad}}=500~{}\textrm{pN}$. The trap thus work well for the
large cancer cells, but less so for smaller cells, and we can estimate the
critical cell radius $a_{\mathrm{cr}}$, below which the cell is not trapped as
$a_{\mathrm{cr}}=\sqrt{F_{z}^{\mathrm{drag}}(10~{}\textrm{\textmu{}m})/F_{z}^{\mathrm{rad}}(10~{}\textrm{\textmu{}m})}\>10~{}\textrm{\textmu{}m}=8.5~{}\textrm{\textmu{}m}$.
[65]

##  V Discussion of the results

In the following, we discuss various aspects of the higher-harmonic membrane
trap, its advantages and disadvantages. Some of these aspects are particular
to the specific type of trap, while other apply to single-cell traps in
general.

Arrays of traps. One distinct advantage of the membrane trap is its relatively
small size, which permits the fabrication of arrays of single-cell traps for
parallel analysis that can be turned on and off in a controlled manner. The
design with $(r_{\mathrm{mv}},n)=(500~{}\textrm{\textmu{}m},10)$ analyzed in
Fig. 5 allows for placing one trap per mm2, a density that can be increased by
identifying membrane modes with the desired properties for smaller values of
$r_{\mathrm{mv}}$.

Improved loading by turning the trap upside down. In Section IV.3, it was
mentioned how the toroidal streaming vortex of the membrane trap helped
loading the trap. Given that $\bm{F}^{\mathrm{grav}}$ points upward, a further
improvement of the loading process can be achieved by placing the membrane
trap above the liquid instead of below. In this case cells would sediment
upward toward the membrane, even when the acoustics is turned off. By
subsequent actuation of the acoustics, the cells would be in a good position
to be carried into the trap by the toroidal vortex. A further advantage of the
inverted trap geometry is that trapping can be obtained at lower excitation
voltages and for smaller particles due to a favorable interplay between
$\bm{F}^{\mathrm{grav}}$ and $\bm{F}^{\mathrm{drag}}$ [66]. Turning the trap
described in Section IV.2 upside down reduces the critical trapping radius
$a_{\mathrm{cr}}$ by 30%.

Reduction of the radial force barrier. If the repelling force barrier
surrounding the trap region poses a problem, as mentioned in Section IV.2, it
would be possible to reduce and perhaps even remove it by careful further
tuning of the liquid. As studied by Qiu et al. [64], a range of carrier fluids
with different acoustic properties can be created using solutions of Ficoll or
iodixanol. In Fig. 8 we show the result of a purely theoretical computation of
$\bm{F}^{\mathrm{pt}}$ for the 60% iodixanol solution of Fig. 6(c,d), were we
assume it possible to keep all parameters constant, while lowering the density
from its actual value,
$\rho_{\mathrm{fl}}^{\mathrm{Idx,}60\%}=1320~{}\textrm{kg}\>\textrm{m${}^{-3}$}$,
to match that of the MCF-7 cell,
$\rho_{\mathrm{MCF\text{-}7}}=1055~{}\textrm{kg}\>\textrm{m${}^{-3}$}$, and
even lower to $930~{}\textrm{kg}\>\textrm{m${}^{-3}$}$. Were it possible to
find a liquid with these parameters, we see in Fig. 8(a,c), how the force
barrier is strongly reduced in the neutral buoyant case, and nearly vanished
in the low-density case of in Fig. 8(b,d).

Figure 8:  Color plots of $F^{\mathrm{pt}}_{r}$ and $F^{\mathrm{pt}}_{z}$ in a
60% iodixanol solution at 18.5 MHz and $1~{}\textrm{V}_{\textrm{pp}}$ as in
Fig. 6(c,d), but assuming theoretically a reduction in the density from its
actual value
$\rho_{\mathrm{fl}}^{\mathrm{Idx,}60\%}=1320~{}\textrm{kg}\>\textrm{m${}^{-3}$}$,
while keeping all other parameters constant. (a) The neutral buoyant case
$\rho_{\mathrm{fl}}^{\mathrm{Idx,}60\%}=\rho_{\mathrm{MCF\text{-}7}}=1055~{}\textrm{kg}\>\textrm{m${}^{-3}$}$.
(b) A low-density case
$\rho_{\mathrm{fl}}^{\mathrm{Idx,}60\%}=930~{}\textrm{kg}\>\textrm{m${}^{-3}$}$.

Creating an axial pressure node by breaking the axisymmetry. The anti-trapping
of cells in water shown in Fig. 6(a,b) is removed in Bessel-beam traps by
breaking the axisymmetry and create a pressure node along the $z$-axis [44].
It would be interesting to study if a similar effect could be achieved in the
membrane trap. It may be done by creating segments in the excitation electrode
both in the radial and azimuthal direction, and then excite a rotating higher-
harmonic mode by running the excitation voltage with appropriate phase shifts
in the azimuthal direction as in Ref. [67].

Sign-change of the radiation force for shorter wave lengths. The anti-trapping
of the membrane trap may be circumvented in a different way by working in the
regime where the acoustic wavelength $\lambda_{0}$ is comparable to the cell
radius $a$. In this regime diffraction effects may change the sign of the
radiation force [68]. A recent numerical study on suspended white blood cell
has shown such a sign reversal to occur for $a\approx 0.4\lambda_{0}$ [69].
This points to a way of getting the membrane trap to work without tuning of
the suspension medium, and would perhaps therefore also work for a standard
isotonic saline solution.

Trapping of smaller-sized particles. Given the relatively strong bulk-driven
Eckart streaming in the membrane trap, it is not easy to reduce the critical
radius from the above-mentioned value of
$a_{\mathrm{cr}}=8.5~{}\textrm{\textmu{}m}$. Several of the methods proposed
in the literature are for systems dominated by boundary-driven Rayleigh
streaming [70, 53, 71]. Perhaps only the method of trapped large seed
particles would be a way to trap smaller particles in the membrane trap [46].

##  VI Conclusion

In this paper, we have shown in a numerical study, how the concept of
selective and efficient excitation of higher-harmonic disk membrane modes can
be applied to acoustofluidic systems. The physical mechanism is based on the
actuation of 1-µm-thick piezoelectric
$\textrm{Al}_{0.6}\textrm{Sc}_{0.4}\textrm{N}$ thin-film transducers by
patterning the excitation electrodes to match the out-of-phase strain pattern
in the transducer at specific mode. As a proof of concept, we have shown that
the $n=3$ higher-harmonic mode in an axisymmetric 200-µm-diameter, 10-µm-thick
silicon membrane actuated at a frequency 52 MHz, emits pressure waves form
each anti-node which results in interference patterns in the pressure field in
a liquid placed above the membrane. When the condition that the wavelength in
the liquid is smaller than that of the node spacing in the membrane mode, the
interference pattern in the liquid can create a global pressure maximum in
some distance above the membrane.

Our main example demonstrated how the global pressure maximum can be
associated with the appearance of a single-cell trap, if the suspending liquid
is tuned such that the cells have a negative acoustic contrast factor. In our
model, this tuning was achieved by using a 60% iodixanol solution. To lower
the operating frequency and obtain a wave length larger than the cell, the
$n=10$ higher-harmonic mode in 1000-µm-diameter, 10-µm-thick membrane was was
excited at 18.5 MHz using appropriately patterned excitation electrodes. We
showed numerically that MCF-7 cancer cells could be trapped approximately 70
µm above the center of the membrane. Several aspects of the trap was discussed
arrays of traps, an inverted setup, reduction of the radial force barrier in
the trap, the creation of an axial pressure mode, and the sign-change of the
radiation force at shorter wave lengths.

Our analysis demonstrate the higher-harmonic MHz acoustic modes can be excited
selectively and efficiently in acoustofluidic systems by appropriate
patterning of transducer electrodes in thin-film transducers. Such modes may
prove useful in practical applications, such as the single-cell trap analyzed
in our main example.

## Acknowledgements

This work was supported by the BioWings project funded by the European Union’s
Horizon 2020 Future and Emerging Technologies (FET) programme, grant No.
801267.

## References

  * Reichert  _et al._ [2018] P. Reichert, D. Deshmukh, L. Lebovitz, and J. Dual, Thin film piezoelectrics for bulk acoustic wave (BAW) acoustophoresis, Lab Chip 18, 3655 (2018).
  * Steckel  _et al._ [2021] A. G. Steckel, H. Bruus, P. Muralt, and R. Matloub, Fabrication, characterization, and simulation of glass devices with $\mathrm{Al}\mathrm{N}$ thin-film transducers for excitation of ultrasound resonances, Phys. Rev. Applied 16, 014014, 1 (2021).
  * Steckel and Bruus [2021] A. G. Steckel and H. Bruus, Numerical study of bulk acoustofluidic devices driven by thin-film transducers and whole-system resonance modes, J. Acoust. Soc. Am. 150, 634 (2021).
  * Brenner  _et al._ [2019] K. Brenner, A. S. Ergun, K. Firouzi, M. F. Rasmussen, Q. Stedman, and B. P. Khuri-Yakub, Advances in capacitive micromachined ultrasonic transducers, Micromachines 10, 152, (27 pp) (2019).
  * Cui  _et al._ [2016] W. Cui, H. Zhang, H. Zhang, Y. Yang, M. He, H. Qu, W. Pang, D. Zhang, and X. Duan, Localized ultrahigh frequency acoustic fields induced micro-vortices for submilliseconds microfluidic mixing, Appl. Phys. Lett. 109, 253503 (6 pp) (2016).
  * Cui  _et al._ [2019] W. Cui, L. Mu, X. Duan, W. Pang, and M. A. Reed, Trapping of sub-100 nm nanoparticles using gigahertz acoustofluidic tweezers for biosensing applications, Nanoscale 11, 14625 (2019).
  * Qian  _et al._ [2021] J. Qian, R. Yang, H. Begum, and J. E.-Y. Lee, Reconfigurable acoustofluidic manipulation of particles in ring-like rich patterns enabled on a bulk micromachined silicon chip, in  _2021 21st International Conference on Solid-State Sensors, Actuators and Microsystems (Transducers)_ (IEEE, 2021) pp. 365–368.
  * Delsing  _et al._ [2019] P. Delsing, A. N. Cleland, M. J. A. Schuetz, J. Knoerzer, G. Giedke, J. I. Cirac, K. Srinivasan, M. Wu, K. C. Balram, C. Bauerle, T. Meunier, C. J. B. Ford, P. V. Santos, E. Cerda-Mendez, H. Wang, H. J. Krenner, E. D. S. Nysten, M. Weiss, G. R. Nash, L. Thevenard, C. Gourdon, P. Rovillain, M. Marangolo, J.-Y. Duquesne, G. Fischerauer, W. Ruile, A. Reiner, B. Paschke, D. Denysenko, D. Volkmer, A. Wixforth, H. Bruus, M. Wiklund, J. Reboud, J. M. Cooper, Y. Fu, M. S. Brugger, F. Rehfeldt, and C. Westerhausen, The 2019 surface acoustic waves roadmap, J. Phys. D Appl. Phys. 52, 353001 (40pp) (2019).
  * Bodé  _et al._ [2020] W. N. Bodé, L. Jiang, T. Laurell, and H. Bruus, Microparticle acoustophoresis in aluminum-based acoustofluidic devices with PDMS covers, Micromachines 11, 292 (2020).
  * Lickert  _et al._ [2021] F. Lickert, M. Ohlin, H. Bruus, and P. Ohlsson, Acoustophoresis in polymer-based microfluidic devices: Modeling and experimental validation, J. Acoust. Soc. Am. 149, 4281 (2021).
  * Bora and Shusteff [2015] M. Bora and M. Shusteff, Efficient coupling of acoustic modes in microfluidic channel devices, Lab Chip 15, 3192 (2015).
  * Moiseyenko and Bruus [2019] R. P. Moiseyenko and H. Bruus, Whole-system ultrasound resonances as the basis for acoustophoresis in all-polymer microfluidic devices, Phys. Rev. Applied 11, 014014 (2019).
  * Tahmasebipour  _et al._ [2020] A. Tahmasebipour, L. Friedrich, M. Begley, H. Bruus, and C. Meinhart, Toward optimal acoustophoretic microparticle manipulation by exploiting asymmetry, J. Acoust. Soc. Am. 148, 359 (2020).
  * Hammarström  _et al._ [2021] B. Hammarström, N. R. Skov, K. Olofsson, H. Bruus, and M. Wiklund, Acoustic trapping based on surface displacement of resonance modes, J. Acoust. Soc. Am. 149, 1445 (2021).
  * Wein  _et al._ [2013] F. Wein, M. Kaltenbacher, and M. Stingl, Topology optimization of a cantilevered piezoelectric energy harvester using stress norm constraints, Struct. Multidiscip. O. 48, 173 (2013).
  * Du  _et al._ [2017] S. Du, Y. Jia, S.-T. Chen, C. Zhao, B. Sun, E. Arroyo, and A. A. Seshia, A new electrode design method in piezoelectric vibration energy harvesters to maximize output power, Sensor. Actuat. A-Phys. 263, 693 (2017).
  * Fu  _et al._ [2018] H. Fu, G. Chen, and N. Bai, Electrode coverage optimization for piezoelectric energy harvesting from tip excitation, SENSORS 18 (2018), 10.3390/s18030804.
  * Yang  _et al._ [2018] Z. Yang, S. Zhou, J. Zu, and D. Inman, High-performance piezoelectric energy harvesters and their applications, Joule 2, 642 (2018).
  * Luo  _et al._ [2021] A. Luo, Y. Zhang, X. Guo, Y. Lu, C. Lee, and F. Wang, Optimization of MEMS vibration energy harvester with perforated electrode, J. Microelectromech. S. 30, 299 (2021).
  * Pulskamp  _et al._ [2012] J. S. Pulskamp, S. S. Bedair, R. G. Polcawich, G. L. Smith, J. Martin, B. Power, and S. A. Bhave, Electrode-shaping for the excitation and detection of permitted arbitrary modes in arbitrary geometries in piezoelectric resonators, IEEE transactions on ultrasonics, ferroelectrics, and frequency control 59, 1043 (2012).
  * Com [2019] COMSOL Multiphysics 5.5 (2019), http://www.comsol.com.
  * Skov  _et al._ [2019a] N. R. Skov, J. S. Bach, B. G. Winckelmann, and H. Bruus, 3D modeling of acoustofluidics in a liquid-filled cavity including streaming, viscous boundary layers, surrounding solids, and a piezoelectric transducer, AIMS Mathematics 4, 99 (2019a).
  * Bach and Bruus [2018] J. S. Bach and H. Bruus, Theory of pressure acoustics with viscous boundary layers and streaming in curved elastic cavities, J. Acoust. Soc. Am. 144, 766 (2018).
  * Hopcroft  _et al._ [2010] M. A. Hopcroft, W. D. Nix, and T. W. Kenny, What is the Young’s modulus of silicon, J. Microelectromech. Syst 19, 229 (2010).
  * Thomsen  _et al._ [2014] E. V. Thomsen, K. Reck, G. Skands, C. Bertelsen, and O. Hansen, Silicon as an anisotropic mechanical material: Deflection of thin crystalline plates, Sensors and Actuators A: Physical 220, 347 (2014).
  * Trolier-McKinstry and Muralt [2004] S. Trolier-McKinstry and P. Muralt, Thin film piezoelectrics for MEMS, Journal of Electroceramics 12, 7 (2004).
  * Settnes and Bruus [2011] M. Settnes and H. Bruus, Theoretical analysis of viscous corrections to the acoustic radiation force on cells in microchannel acoustophoresis, in  _Proc. 15th MicroTAS, 2 - 6 October 2011, Seattle (WA), USA_ , edited by J. Landers, A. Herr, D. Juncker, N. Pamme, and J. Bienvenue (CBMS, 2011) pp. 160–162.
  * Ley and Bruus [2017] M. W. H. Ley and H. Bruus, Three-dimensional numerical modeling of acoustic trapping in glass capillaries, Phys. Rev. Applied 8, 024020 (2017).
  * Muller and Bruus [2014] P. B. Muller and H. Bruus, Numerical study of thermoviscous effects in ultrasound-induced acoustic streaming in microchannels, Phys. Rev. E 90, 043016 (2014).
  * Augustsson  _et al._ [2016] P. Augustsson, J. T. Karlsen, H.-W. Su, H. Bruus, and J. Voldman, Iso-acoustic focusing of cells for size-insensitive acousto-mechanical phenotyping, Nat. Commun. 7, 11556 (2016).
  * An  _et al._ [2009] J. An, J. Lee, S. H. Lee, J. Park, and B. Kim, Separation of malignant human breast cancer epithelial cells from healthy epithelial cells using an advanced dielectrophoresis-activated cell sorter (dacs), Analytical and bioanalytical chemistry 394, 801 (2009).
  * Cushing  _et al._ [2017] K. W. Cushing, F. Garofalo, C. Magnusson, L. Ekblad, H. Bruus, and T. Laurell, Ultrasound characterization of microbead and cell suspensions by speed of sound measurements of neutrally buoyant samples, Anal. Chem. 89, 8917 (2017).
  * Caro  _et al._ [2015] M. A. Caro, S. Zhang, T. Riekkinen, M. Ylilammi, M. A. Moram, O. Lopez-Acevedo, J. Molarius, and T. Laurila, Piezoelectric coefficients and spontaneous polarization of ScAlN, J. Phys.-Condens. Mat. 27, 245901 (2015).
  * Olsson  _et al._ [2020] R. H. Olsson, Z. Tang, and M. D’Agati, Doping of aluminum nitride and the impact on thin film piezoelectric and ferroelectric device performance, in  _2020 IEEE Custom Integrated Circuits Conference (CICC)_ (2020) pp. 1–6.
  * Kim  _et al._ [2001] J. Kim, D.-i. D. Cho, and R. S. Muller, Why is (111) silicon a better mechanical material for MEMS? in  _Transducers 01 Eurosensors XV_ (Springer, 2001) pp. 662–665.
  * Hahn and Dual [2015] P. Hahn and J. Dual, A numerically efficient damping model for acoustic resonances in microfluidic cavities, Physics of Fluids 27, 062005 (2015).
  * Note [1] See Supplemental Material at http://bruus-lab.dk/files/Steckel_membrane_Suppl.zip for animations of $p_{1}$ and $\bm{u}_{1}$ in Fig. 2(a), Fig. 2(c), Fig. 3(a), Fig. 4, and Fig. 5.
  * Skov  _et al._ [2019b] N. R. Skov, P. Sehgal, B. J. Kirby, and H. Bruus, Three-dimensional numerical modeling of surface-acoustic-wave devices: Acoustophoresis of micro- and nanoparticles including streaming, Phys. Rev. Applied 12, 044028 (2019b).
  * Baresch  _et al._ [2016] D. Baresch, J.-L. Thomas, and R. Marchiano, Observation of a single-beam gradient force acoustical trap for elastic particles: Acoustical tweezers, Phys. Rev. Lett. 116, 024301 (2016).
  * Karlsen and Bruus [2017] J. T. Karlsen and H. Bruus, Acoustic tweezing and patterning of concentration fields in microfluidics, Phys. Rev. Applied 7, 034017 (2017).
  * Riaud  _et al._ [2017] A. Riaud, M. Baudoin, O. Bou Matar, L. Becerra, and J.-L. Thomas, Selective manipulation of microscopic particles with precursor swirling Rayleigh waves, Phys. Rev. Applied 7, 024007 (2017).
  * Baudoin  _et al._ [2019] 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).
  * Gong and Baudoin [2019] Z. Gong and M. Baudoin, Particle assembly with synchronized acoustic tweezers, Phys. Rev. Applied 12, 024045 (2019).
  * Baudoin and Thomas [2020] M. Baudoin and J.-L. Thomas, Acoustic tweezers for particle and fluid micromanipulation, Annual Review of Fluid Mechanics 52, 205 (2020).
  * Hammarström  _et al._ [2010] B. Hammarström, M. Evander, H. Barbeau, M. Bruzelius, J. Larsson, T. Laurell, and J. Nillsson, Non-contact acoustic cell trapping in disposable glass capillaries, Lab Chip 10, 2251 (2010).
  * Hammarström  _et al._ [2014] B. Hammarström, B. Nilson, T. Laurell, J. Nilsson, and S. Ekström, Acoustic trapping for bacteria identification in positive blood cultures with maldi-tof ms, Anal. Chem. 86, 10560 (2014).
  * Evander  _et al._ [2015] M. Evander, O. Gidlof, B. Olde, D. Erlinge, and T. Laurell, Non-contact acoustic capture of microparticles from small plasma volumes, Lab Chip 15, 2588 (2015).
  * Hammarström  _et al._ [2012] B. Hammarström, T. Laurell, and J. Nilsson, Seed particle enabled acoustic trapping of bacteria and nanoparticles in continuous flow systems, Lab Chip 12, 4296 (2012).
  * Hammarström  _et al._ [2014] B. Hammarström, M. Evander, J. Wahlström, and J. Nilsson, Frequency tracking in acoustic trapping for improved performance stability and system surveillance, Lab Chip 14, 1005 (2014).
  * Olofsson  _et al._ [2021] K. Olofsson, V. Carannante, M. Takai, B. Onfelt, and M. Wiklund, Ultrasound-based scaffold-free core-shell multicellular tumor spheroid formation, Micromachines 12, 329 1 (2021).
  * Hagsäter  _et al._ [2008] S. M. Hagsäter, A. Lenshof, P. Skafte-Pedersen, J. P. Kutter, T. Laurell, and H. Bruus, Acoustic resonances in straight micro channels: Beyond the 1D-approximation, Lab Chip 8, 1178 (2008).
  * Augustsson  _et al._ [2011] P. Augustsson, R. Barnkob, S. T. Wereley, H. Bruus, and T. Laurell, Automated and temperature-controlled micro-PIV measurements enabling long-term-stable microchannel acoustophoresis characterization, Lab Chip 11, 4152 (2011).
  * Bach and Bruus [2020] J. S. Bach and H. Bruus, Suppression of acoustic streaming in shape-optimized channels, Phys. Rev. Lett. 124, 214501 (2020).
  * Shi  _et al._ [2009] J. Shi, D. Ahmed, X. Mao, S.-C. S. Lin, A. Lawit, and T. J. Huang, Acoustic tweezers: patterning cells and microparticles using standing surface acoustic waves (SSAW), Lab Chip 9, 2890 (2009).
  * Collins  _et al._ [2016] D. J. Collins, A. Neild, and Y. Ai, Highly focused high-frequency travelling surface acoustic waves (saw) for rapid single-particle sorting, Lab Chip 16, 471 (2016).
  * Gascoyne  _et al._ [2002] P. Gascoyne, C. Mahidol, M. Ruchirawat, J. Satayavivad, P. Watcharasit, and F. F. Becker, Microsample preparation by dielectrophoresis: isolation of malaria, Lab Chip 2, 70 (2002).
  * Nilsson  _et al._ [2009] J. Nilsson, M. Evander, B. Hammarström, and T. Laurell, Review of cell and particle trapping in microfluidic systems, Analytica Chimica Acta 649, 141 (2009).
  * Gustafson  _et al._ [2021] K. T. Gustafson, K. T. Huynh, D. Heineck, J. Bueno, A. Modestino, S. Kim, A. Gower, R. Armstrong, C. E. Schutt, and S. D. Ibsen, Automated fluorescence quantification of extracellular vesicles collected from blood plasma using dielectrophoresis dagger, Lab Chip 21 (2021), 10.1039/d0lc00940g.
  * Li  _et al._ [2015] P. Li, Z. Mao, Z. Peng, L. Zhou, Y. Chen, P.-H. Huang, C. I. Truica, J. J. Drabick, W. S. El-Deiry, M. Dao, S. Suresh, and T. J. Huang, Acoustic separation of circulating tumor cells, Proc. Natl. Acad. Sci. U.S.A. 112, 4970 (2015).
  * Low and Abas [2015] W. S. Low and W. A. B. W. Abas, Benchtop technologies for circulating tumor cells separation based on biophysical properties, Biomed. Res. Int. 2015, 239362 (2015).
  * Olofsson  _et al._ [2020] K. Olofsson, B. Hammarstrom, and M. Wiklund, Acoustic separation of living and dead cells using high density medium, Lab Chip 20, 1981 (2020).
  * Karlsen  _et al._ [2016] J. T. Karlsen, P. Augustsson, and H. Bruus, Acoustic force density acting on inhomogeneous fluids in acoustic fields, Phys. Rev. Lett. 117, 114504 (2016).
  * Karlsen  _et al._ [2018] J. T. Karlsen, W. Qiu, P. Augustsson, and H. Bruus, Acoustic streaming and its suppression in inhomogeneous fluids, Phys. Rev. Lett. 120, 054501 (2018).
  * Qiu  _et al._ [2019] 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).
  * Muller  _et al._ [2012] 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).
  * Li  _et al._ [2021] 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, 1 (2021).
  * Tran  _et al._ [2012] S. B. Q. Tran, P. Marmottant, and P. Thibault, Fast acoustic tweezers for the two-dimensional manipulation of individual particles in microfluidic channels, Appl. Phys. Lett. 101, 114103 (2012).
  * Hasegawa [1977] T. Hasegawa, Comparison of 2 solutions for acoustic radiation pressure on a sphere, J. Acoust. Soc. Am. 61, 1445 (1977).
  * Habibi  _et al._ [2017] R. Habibi, C. Devendran, and A. Neild, Trapping and patterning of large particles and cells in a 1d ultrasonic standing wave, Lab Chip 17, 33279 (2017).
  * Qiu  _et al._ [2020] W. Qiu, H. Bruus, and P. Augustsson, Particle-size-dependent acoustophoretic motion and depletion of micro-and nano-particles at long timescales, Phys. Rev. E 102, 013108 (2020).
  * Winckelmann and Bruus [2021] B. Winckelmann and H. Bruus, Theory and simulation of electroosmotic suppression of acoustic streaming, J. Acoust. Soc. Am. 149, 3917 (2021).