2.1. Fluid dynamics equations

Consider a cylindrical microcavity resonator of radius $R$ and height $H$ filled with a quiescent viscous liquid. The rigid boundary walls that enclose the cavity are represented by $\mathcal{S}=S_{{T}}\cup S_{{B}}\cup S_{{L}}$ , with $S_{{T}}$ , $S_{{B}}$ and $S_{{L}}$ denoting the top, bottom and lateral boundary, respectively. The Eulerian reference frame is placed in the centre of $S_{{B}}$ at rest. The fluid flow is described by a mass density $\rho$ , pressure $p$ and fluid velocity $\tilde {\boldsymbol{v}}$ as a function of position $\boldsymbol{r}$ and time $t$ . To take advantage of the cavity’s shape, cylindrical coordinates are adopted $(\varrho , \varphi ,z)$ , with position vector $\boldsymbol{r} = \varrho \, \boldsymbol{e}_\varrho + \varphi \, \boldsymbol{e}_\varphi + z\, \boldsymbol{e}_z$ , in which $\boldsymbol{e}_i$ $(i=\varrho ,\varphi ,z)$ represents the transverse radial, azimuthal and axial unit vectors. The origin of the coordinate system is set at the centre of the bottom of the cavity.

The unperturbed state of the liquid is described by an ambient density $\rho _0$ , pressure $p_0$ , adiabatic speed of sound $c_0$ , isentropic compressibility $\kappa _0=1/\rho _0 c_0^2$ , and dynamic $\eta _0$ and bulk $\zeta _0$ viscosity. These parameters remain constant for the analysis under consideration. The physical parameters of water are listed in table 1.

Table 1.

Physical properties of water at ${25}\,{^{\circ }\mathrm{C}}$ (Holmes, Parker & Povey Reference Holmes, Parker and Povey2011).

In the investigation of acoustic wave propagation within the enclosed cavity, the process is assumed to be isentropic, with both internal heat generation and conduction effects being considered negligible. This assumption allows for a simplification in complex interactions between thermodynamic and acoustic phenomena. A detailed analysis of adopting the isentropic hypothesis is presented in Appendix A.

The fluid dynamics equations are described by the mass and linear momentum conservation equations along with the density-pressure constitutive relation,

(2.1a) \begin{align} \partial _t \rho + \boldsymbol{\nabla } \cdot (\rho {\boldsymbol{v}}) &= 0, \end{align}

(2.1b) \begin{align} \partial _t (\rho {\boldsymbol{v}}) -\boldsymbol{\nabla } \cdot \left (\boldsymbol{\sigma } - \rho {\boldsymbol{v}}{\boldsymbol{v}} \right ) &= \boldsymbol{0}, \end{align}

(2.1c) \begin{align} p &= c_0^2 \rho . \end{align}

The stress tensor is given by

(2.2) \begin{equation} \boldsymbol{\sigma } = \left [- p + \left ( \zeta _0 -\dfrac {2}{3}\eta _0 \right ) (\boldsymbol{\nabla } \cdot {\boldsymbol{v}})\right ] \unicode{x1D644} + \eta _0\! \left [\boldsymbol{\nabla } {\boldsymbol{v}} + \left (\boldsymbol{\nabla } {\boldsymbol{v}}\right )^{\text{T}} \right ], \end{equation}

where $\unicode{x1D644}$ is the unit tensor and the superscript $^{\text{T}}$ denotes the transpose of a tensor. The energy conservation equation is not needed to describe acoustic phenomenon within the isentropic hypothesis.

2.2. Perturbation expansion

An acoustic wave can be injected into the microcavity through a time-harmonic oscillation with angular frequency $\omega$ , at the bottom surface $S_{{B}}$ at $z=0$ . By considering that $S_{{B}}$ moves back and forth along the $z$ axis, we express the surface displacement as

(2.3) \begin{equation} S(\varrho , \varphi ; t) = d(\varrho , \varphi )\, \textrm {e}^{-\textrm {i} \omega t + \textrm {i} \frac {\pi }{2}}, \quad \varrho \lt R, \end{equation}

where $d$ is the displacement function. The phase $\pi /2$ is added here for convenience. Note that no lateral motion of $S_{{B}}$ will be taken into account in our analysis. The corresponding excitation velocity is

(2.4) \begin{equation} v_{{ex}}(\varrho ,\varphi ; t) =\partial _t{S} = \omega \, d(\varrho ,\varphi )\, \textrm {e}^{-\textrm {i} \omega t}. \end{equation}

In the present analysis, the peak velocity of $S_{{B}}$ is assumed to be much smaller than the speed of sound in the medium. In terms of the Mach number, this implies $\textit{Ma}= \omega d_0/c_0\ll 1$ , in which $d_0$ is the given peak displacement at $S_{{B}}$ . To solve the fluid dynamics equations given in (2.1) in this approximation, the fields are conveniently expanded as

(2.5a) \begin{align} \rho (\boldsymbol{r},t) &= \rho _0[1 + \textit{Ma} \,\tilde {\rho }_1(\boldsymbol{r}) \textrm {e}^{-\textrm {i} \omega t}], \end{align}

(2.5b) \begin{align} p(\boldsymbol{r},t) &= p_0 + \textit{Ma} \kappa _0^{-1}\, \tilde {p}_1(\boldsymbol{r}) \textrm {e}^{-\textrm {i} \omega t}, \end{align}

(2.5c) \begin{align} {\boldsymbol{v}}(\boldsymbol{r},t) &= \boldsymbol{0} + \textit{Ma} c_0 \,\tilde {\boldsymbol{v}}_1(\boldsymbol{r}) \textrm {e}^{-\textrm {i} \omega t}. \end{align}

The dimensionless amplitude functions $\tilde {p}_1$ , $\tilde {\rho }_1$ and $\tilde {\boldsymbol{v}}_1$ represent respectively the pressure, density and fluid velocity in the linear approximation. Hereafter, a variable with a tilde accent is a dimensionless quantity. It is convenient to express the excitation in terms of a dimensionless function $\tilde {v}_{{b}}$ ,

(2.6) \begin{equation} v_{{ex}} = \textit{Ma} c_0 \tilde {v}_{{b}}\,\textrm {e}^{-\textrm {i} \omega t},\quad \text{with} \ \tilde {v}_{{b}} = \frac {d(\varrho ,\varphi )}{d_0} . \end{equation}

In acoustofluidic cavities based on piezoelectric ceramics, the peak displacement at the bottom is generally $d_0\approx {1}\, {\mathrm{nm}}$ in the low-megahertz range (Tran-Huu-Hue et al. Reference Tran-Huu-Hue, Levassort, Felix, Damjanovic, Wolny and Lethiecq2000). Hence, the expected Mach number for a water-filled cavity at ${1}\, {\mathrm{MHz}}$ is $\textit{Ma}\sim 10^{-6}$ .

2.3. Linear acoustics equations

A viscous fluid supports both compressional and shear flows. So we introduce a decomposition of the fluid velocity as a compressional part ( $\tilde {\boldsymbol{v}}^{\textit{c}}$ ) plus a shear part ( $\tilde {\boldsymbol{v}}^{\textit{s}}$ ),

(2.7) \begin{equation} \tilde {\boldsymbol{v}}_1 = \tilde {\boldsymbol{v}}_{{c}} + \tilde {\boldsymbol{v}}_{{s}}. \end{equation}

The compressional component is irrotational, while the shear part is incompressible,

(2.8) \begin{equation} \boldsymbol{\nabla } \times \tilde {\boldsymbol{v}}_{{c}} = \boldsymbol{0} \quad \text{and}\quad \boldsymbol{\nabla } \cdot \tilde {\boldsymbol{v}}_{{s}} =0. \end{equation}

Replacing (2.5a ), (2.5c ) and (2.7) into (2.1a ), while expressing the dimensionless mass density with (2.1c ) as $\tilde {\rho }_1=\tilde {p}_1$ , we obtain the first-order dimensionless pressure as

(2.9) \begin{equation} \tilde {p}_1 = - \frac {\textrm {i} c_0}{\omega } \boldsymbol{\nabla } \cdot \tilde {\boldsymbol{v}}_{{c}}. \end{equation}

Inserting (2.5c ) with (2.7) into (2.1b ) and using (2.9) followed by the vector operator identity $\boldsymbol{\nabla }(\boldsymbol{\nabla }\cdot \square ) = \nabla ^2\square + \boldsymbol{\nabla }\times \boldsymbol{\nabla }\times \square$ , we arrive at the compressional and shear velocities, respectively

(2.10a) \begin{align} \tilde {\boldsymbol{v}}_{{c}} &= - \frac {\textrm {i} c_0 }{\omega } \left (1- \textrm {i} \Gamma _{{c}}\right ) \boldsymbol{\nabla } \tilde {p}_1, \end{align}

(2.10b) \begin{align} \tilde {\boldsymbol{v}}_{{s}} &= \frac {\textrm {i} \eta _0 }{\omega \rho _0}\nabla ^2 \tilde {\boldsymbol{v}}_{{s}}. \end{align}

The compressional damping factor (in the fluid bulk) is given by

(2.11) \begin{equation} \Gamma _{{c}} =\kappa _0 \omega \left (\zeta _0 + \frac {4}{3}\eta _0 \right ). \end{equation}

Taking the divergence of (2.10a ) and using (2.9) yields the Helmhotz equation for pressure:

(2.12) \begin{align} \big (\nabla ^2 + k_{{c}}^2\big) \tilde {p}_1 = 0, \quad \text{with} \ k_{{c}} = \left ( 1 + \frac {\textrm {i}}{2}\Gamma _{{c}}\right )k. \end{align}

The real part of the compressional wavenumber is $k=\omega /c_0$ , with the associated acoustic wavelength being $\lambda = 2\pi /k$ . Equation (2.10b ) can be recast in the form

(2.13) \begin{align} \big(\nabla ^2 + k_{{s}}^2\big) \tilde {\boldsymbol{v}}_{{s}} = \boldsymbol{0},\quad k_{{s}} =\frac {1+\textrm {i}}{\delta }, \end{align}

where

(2.14) \begin{equation} \delta =\sqrt {\frac {2 \eta _0}{\rho _0 \omega }} \end{equation}

is the thickness of the viscous boundary layer. According to table 1, this thickness is $\delta ={0.53}\, {{\mu} \textrm{m}}$ for water at room temperature ( ${25}\,{^{\circ }\mathrm{C}}$ ) and ${1}\, {\mathrm{MHz}}$ .

3.1. Eigenfunctions of the Helmholtz equation

The regular eigenfunctions, i.e. independent solutions to (2.12) in cylindrical coordinates, for a rigid cavity resonator can be expressed by (Morse Reference Morse1948, p. 398)

(3.1) \begin{equation} \tilde {p}_1 = J_n (k_\varrho \varrho ) \big[a_1\cos \!\left (k_z z\right ) + a_2 \sin (k_z z)\big] \big(b_1\textrm {e}^{\textrm {i} n \varphi } + b_2 \textrm {e}^{-\textrm {i} n \varphi }\big), \end{equation}

where $J_n$ represents the $n$ th-order Bessel function with $n\in \mathbb{Z}$ . The coefficients $a_i$ and $b_i$ ( $i=1,2$ ) will be determined from the boundary conditions (4.2). Furthermore, from (2.12), we find that the axial ( $k_z$ ) and circular ( $k_\varrho$ ) wavenumbers are related by

(3.2) \begin{equation} k_z^2 + k_\varrho ^2 = (1+\textrm {i} \Gamma _{{c}})\, k^2. \end{equation}

Table 2.

Roots of equation $\textrm {d} J_n(x)/\textrm {d} x =0$ for $n,m\leqslant 4$ .

3.2. Natural frequency modes

Prior to commencing the analysis of a viscous fluid-filled cavity resonator, it is beneficial to present the natural frequency of the system. For a rigid cavity with boundary surface denoted by $S$ , the pressure satisfies the Neumann boundary condition at the walls,

(3.3) \begin{equation} \boldsymbol{n}\cdot \boldsymbol{\nabla } \tilde {p}_1{\kern-3pt}\left |_{\boldsymbol{r}\in S } = 0,\right . \end{equation}

where $\boldsymbol{n}$ is the normal unit vector of the surface $S$ pointing outwardly. In this case, the axial and circular wavenumbers become

(3.4) \begin{equation} k_{z}^l = \frac {l \pi }{H}, \quad k_{\varrho }^{nm} = \frac {j_{nm}'}{R}, \end{equation}

with $n,m,l=0, 1, 2, \ldots$ The function $j_{nm}'$ represents a root of equation $\textrm {d} J_n(x)/\textrm {d} x =0$ . Note that $k_\varrho ^{nm}$ is the eigenvalue of both natural modes with a positive or negative $n$ . Table 2 shows the roots considering mode numbers in the interval $0\leqslant n,m\leqslant 5$ . According to (3.2), the compression wavenumber followed by the angular frequency of the natural mode $(n\,m\,l)$ are given by

(3.5a) \begin{align} k^{nml}_0 &= \sqrt {\left (k_{z}^l\right )^2 + \big(k_{\varrho }^{nm}\big)^2}, \end{align}

(3.5b) \begin{align} \omega ^{nml}_0 &= c_0 k^{nml}_0. \end{align}

Purely radial waves of mode $(0\,m\,0)$ appear along the radial vector $\boldsymbol{e}_\varrho$ . Likewise, waves aligned with the $z$ axis are present in mode $(0\,0\,l)$ . Waves corresponding to displacement along the azimuthal direction, represented by $\boldsymbol{e}_\varphi$ , manifest in mode $(n\,0\,0)$ . Finally, modes characterised by $(n\,m\,0)$ are designated as circular modes.

3.3. Scaling parameter

Viscous effects near solid boundaries are significant due to the steep velocity gradients imposed by the no-slip condition. For acoustofluidic experimental settings in the megahertz range and above, the following condition involving the viscous boundary layer thickness is assumed to hold:

(3.6) \begin{equation} \delta \ll R, H, \lambda . \end{equation}

We define the scaling parameter for the viscous-boundary thickness as

(3.7) \begin{equation} \epsilon \equiv \frac { 2\pi \delta }{\lambda } = k \delta \ll 1. \end{equation}

Note that this parameter scales with frequency as $\epsilon \sim \sqrt {\omega }$ . At a frequency of ${1}\,{\mathrm{MHz}}$ , the scaling parameter for water is $\epsilon =0.0022$ , whereas at ${1}\,{\mathrm{GHz}}$ , it increases to $\epsilon =0.069$ . It is further assume that the peak displacement at the bottom is much smaller than the thickness of the viscous boundary layer, $d_0\ll \delta$ .

The compressional damping factor $\Gamma _{{c}}$ can also be expressed in terms of $\epsilon$ by using (2.11),

(3.8) \begin{equation} \Gamma _{{c}} = \left (\frac {2}{3} + \frac {\zeta _0}{2\eta _0}\right )\epsilon ^2. \end{equation}

Again, for water at the same conditions, we have $\Gamma _{{c}} = {5.06 \times 10^{-6}}$ . Equation (3.8) reveals that damping in the cavity bulk is generally very small.

6.1. Vortex excitation

In accordance with (2.4), the velocity excitation required to generate a mode with positive vorticity in the cavity is described as

(6.1) \begin{equation} \tilde {v}_{{b}} = J_{n}\big(k_{\varrho }^{nm}\varrho \big) \textrm {e}^{ \textrm {i} n \varphi }. \end{equation}

The induced pressure amplitude rotates counterclockwise when observed from above the $x{-}y$ plane – see figure 1. In this case, the azimuthal coefficients in (3.1) become $b_1=1$ and $b_2=0$ . Conversely, when the angular factor is $\textrm {e}^{-\textrm {i} n\varphi }$ , the induced pressure rotates clockwise from the same viewpoint, and $b_1=0$ and $b_2=1$ . For the case where $n=0$ , the excitation produces a non-vortex mode within the resonator, with $b_1=b_2=1/2$ .

6.2. Pressure coefficients

Substituting (3.1) and (6.1) into (4.2a ) and (4.2b ), we arrive at a system of linear equations for the pressure coefficients $a_1$ and $a_2$ :

(6.2a) \begin{align} \boldsymbol{A} \left ( \begin{matrix} a_1 \\ a_2 \end{matrix} \right ) &= \left ( \begin{matrix} \textrm {i} (k H)^2 \\ 0 \end{matrix} \right ), \end{align}

(6.2b) \begin{align} \boldsymbol{A} &= \left ( \begin{matrix} \epsilon \chi & k_z k H^2 \\[4pt] k_z k H^2 \sin (k_z H) + \epsilon \chi \cos (k_z H)& - k k_z H^2 \cos (k_z H)+ \epsilon \chi \sin ( k_z H) \end{matrix} \right ), \end{align}

with $\chi = (1+\textrm {i}) (k^2-k_z^2 ) H^2/2$ . Terms of order $\epsilon \Gamma _{{c}}$ were discarded. Also, we have taken the approximation $J_n(k_\varrho ^{nml} \varrho )\approx J_n(k_{\varrho }^{nm} \varrho )$ . According to Cramer’s rule, the solution of (6.2) is given by

(6.3a) \begin{align} a_1 &= \frac {\textrm {i} k}{k_z} \frac { k k_z H^2 \cos (k_z H) - \epsilon \chi \sin (k_z H)} { k k_z H^2 \sin (k_z H) + 2 \epsilon \chi \cos (k_z H)},\\[-11pt] \nonumber \end{align}

(6.3b) \begin{align} a_2 &= \frac {\textrm {i} k}{k_z} \frac {k k_z H^2 \sin (k_z H) + \epsilon \chi \cos (k_z H)} {k k_z H^2 \sin (k_z H) + 2 \epsilon \chi \cos (k_z H)}. \end{align}

It is important to compare the magnitude of the pressure coefficients. Referring to (5.2b ), we observe that $\cos (k_z z) \approx \cos (k_z^l z)$ and $\sin (k_z z) \approx \sin (k_z^l z)$ . This implies that $a_2/a_1 \sim \epsilon$ , allowing us to safely neglect the term $a_2 \sin (k_z z)$ in (3.1). Therefore, the pressure near the resonance is assumed to be

(6.4) \begin{equation} \tilde {p}_1^{nml}(\boldsymbol{r}) = a_1 J_n\big(k_\varrho ^{nm} \varrho \big) \cos \big (k_z^l z\big) \textrm {e}^{\textrm {i} n \varphi }. \end{equation}

The amplitude coefficient $a_1$ and the boundary-layer factor $\Gamma _\delta ^{nml}$ are determined in Appendix D. Gathering the result, we have

(6.5a) \begin{align} a_1 &= \frac {\textrm {i} (2-\delta _{l,0})}{k_0^{nml}H} \frac {\dfrac {1}{2}\omega ^{nml}}{\omega - \omega ^{n m l} + \dfrac {1}{2} \textrm {i}\,\omega ^{nml} \Gamma ^{n m l} }, \end{align}

(6.5b) \begin{align} \Gamma ^{nml} &= \Gamma _\delta ^{nml} + \Gamma _{{c}}^{nml}, \end{align}

(6.5c) \begin{align} \Gamma _\delta ^{nml} &= \epsilon ^{nml} \left \{ \left [\frac {2-\delta _{l,0}}{k_0^{nml} H} + \frac {1}{k_0^{nml} R} \frac {\big(k_{z}^l R\big)^2 + n^2 }{\big(k_{\varrho }^{nm} R\big)^2-n^2} \right ] \left (\frac {k_{\varrho }^{nm}}{k_0^{nml}}\right )^2 + \frac { \delta _{n,0} \delta _{m,0}}{ k_0^{nml} R} \left (\frac {k_{z}^{l}}{k_0^{nml}}\right )^2 \right \}. \end{align}

The total damping factor $\Gamma ^{nml}$ is determined by using the wavenumbers related to the natural frequency $\omega _0^{nml}$ .

6.3. Pressure field

We now replace (6.5a ) into (6.4) to obtain the pressure amplitude of the single-vortex mode $(n\,m\,l)$ ,

(6.6a) \begin{align} {p}_1^{nml} &= {A}^{nml} R^{nml} \Psi ^{nml},\quad (n,m,l=0, 1, 2, \ldots ), \end{align}

(6.6b) \begin{align} {A}^{nml} &= \frac {2-\delta _{l,0}}{\Gamma ^{nml} k_0^{nml}H}, \end{align}

(6.6c) \begin{align} {R}^{nml} &=\frac {1}{2}\frac { \omega ^{nml}\Gamma ^{nml} } { \omega - \omega ^{nml}\left (1 + \frac {1}{2}\textrm {i} \, \Gamma ^{nml}\right )}, \end{align}

(6.6d) \begin{align} \Psi ^{nml} &= J_n\big (k_\varrho ^{nm} \varrho \big) \cos \big(k_z^l z\big) \textrm {e}^{\textrm {i} \left (n \varphi +\frac {\pi }{2}\right )}. \end{align}

where ${A}^{nml}$ is the pressure amplitude, ${R}^{nml}$ is the resonance function and $\Psi ^{nml}$ is the cavity’s wavefunction. The pressure amplitude scales inversely with the total damping factor $(\Gamma ^{nml})^{-1}$ and the height-to-wavelength ratio $(k_0^{nml} H)$ . At the resonance frequency, we have $ R^{nml}|_{\omega =\omega ^{nml}} = -\textrm {i}$ . Hence, the pressure becomes

(6.7) \begin{equation} p_1^{nml}|_{\omega =\omega ^{nml}}= -\textrm {i} {A}^{nml} \Psi ^{nml}. \end{equation}

We note also that the pressure phase is composed of three components, namely, the resonance $(\varphi _{{r}}^{nml})$ , spatial ( $\varphi _{{s}} ^{nml}$ ) and vortex $(n\varphi )$ phase:

(6.8) \begin{equation} \varphi ^{nml} = \varphi _{{r}}^{nml} + \varphi _{{s}}^{nml}+ n\varphi + \frac {\pi }{2}. \end{equation}

6.4. Fluid velocity fields

We can determine the fluid velocity in the fluid bulk by using (2.10a ),

(6.9) \begin{equation} \tilde {\boldsymbol{v}}_{{c}}^{nml} = -\textrm {i} {A}^{nml} {R}^{nml}\, \tilde {\boldsymbol{\nabla }}\Psi ^{nml}, \end{equation}

where $\tilde {\boldsymbol{\nabla }}=(k_0^{nml})^{-1}\boldsymbol{\nabla }$ is the scaled gradient to the natural wavenumber. Simultaneously, the shear velocity near the boundaries can be calculated by applying (6.9) in combination with (B1) and (B7),

(6.10) \begin{equation} \tilde {\boldsymbol{v}}_{{s}}^{nml} = -\textrm {i} {A}^{nml} {R}^{nml} \Big[ \Big ( \textrm {e}^{\textrm {i} k_{{s}} z} + (-1)^l \textrm {e}^{\textrm {i} k_{{s}} (H-z)} \Big) \tilde {\boldsymbol{\nabla }}_{\!\boldsymbol{e}_z} + \textrm {e}^{\textrm {i} k_{{s}} (R-\varrho )} \tilde {\boldsymbol{\nabla }}_{\!\boldsymbol{e}_\varrho } \Big]\Psi ^{nml}, \end{equation}

with $\tilde {\boldsymbol{\nabla }}_{\!\boldsymbol{n}}= (\unicode{x1D644} - \boldsymbol{n}\boldsymbol{n} )\cdot \tilde {\boldsymbol{\nabla }}$ being the scaled surface gradient. The shear velocity is confined in the viscous boundary as it decays exponentially away from the cavity walls.

8.1. Acoustic energy density

The time-averaged acoustic energy density, often referred to as the energy landscape, offers a spatial representation of energy distribution within the cavity. This provides valuable insights into the system’s behaviour and performance. The energy landscape ( $E$ ) arises from the contributions of both kinetic ( $E_{{k}}$ ) and potential ( $E_{{p}}$ ) acoustic energy densities,

(8.2a) \begin{align} {E} &= {E}_{{k}} +{E}_{{p}}, \end{align}

(8.2b) \begin{align} {E}_{{k}}&= E_0\left |\tilde {\boldsymbol{v}}_{{c}}\right |^2, \end{align}

(8.2c) \begin{align} {E}_{{p}}&=E_0 \left |\tilde {p}_1\right |^2, \end{align}

where $E_0 = \textit{Ma}^2 \rho _0 c_0^2/4$ is the characteristic energy density of the cavity. Using (2.10a ) allows us to express the kinetic energy density as

(8.3) \begin{equation} {E}_{{k}}= \frac {E_0}{k^2} \boldsymbol{\nabla }\! \tilde {p}_1 \cdot \boldsymbol{\nabla }\! \tilde {p}_1^*= \frac {E_0}{k^2}\left (\frac {1}{2}\nabla ^2|\tilde {p}_1|^2-{\mathrm{Re}}\big[\tilde {p}_1^*\nabla ^2 \tilde {p}_1\big]\right )= E_0\left (\frac {1}{2k^2}\nabla ^2 + 1\right ) |\tilde {p}_1|^2. \end{equation}

Here, we have used (2.12). Hence, the energy density reads

(8.4) \begin{align} {E} = {E}_{{k}} + {E}_{{p}} = 2 E_0 \left (\frac {1}{4k^2}\nabla ^2 + 1\right ) |\tilde {p}_1|^2. \end{align}

We define the dimensionless energy density as $\tilde {E}=E/E_0$ . Using (6.6a ) and (7.5), we obtain the single-mode energy density as

(8.5) \begin{equation} \tilde {E}^{nml} = 2\big|A^{nml}R^{nml} \big |^2 \left (\frac {1}{4}\tilde {\nabla }^2 + 1\right ) \big|\Psi ^{nml}\big|^2, \end{equation}

where $\tilde {\nabla }^2=k^{-2}\nabla ^2$ is the scaled Laplacian. Inserting (7.5) into (8.4), we obtain the energy density of the dual mode as

(8.6) \begin{equation} \tilde {E}^{n_1 m_1l_1}_{nml} = \tilde {E}^{nml} + \tilde {E}^{ n_1 m_1l_1} + \tilde {E}_{{cross}}. \end{equation}

The cross term of the energy density is given by

(8.7) \begin{equation} \tilde {E}_{{cross}}= 4 A^{nml} A^{n_1m_1 l_1} {\mathrm{Re}} \left[ \big (R^{nml}\big)^* R^{n_1m_1 l_1} \left (\frac {1}{4}\tilde {\nabla }^2 + 1\right ) \big({\Psi ^{nml}}\big)^* \Psi ^{n_1 m_1l_1}\right ]. \end{equation}

While the single-mode energy density does not depend on the azimuthal angle $\varphi$ , the crossed energy does.

8.2. Dissipated power

The time-averaged dissipated power $P^{nml}$ for a specific mode $(n\,m\,l)$ is linked with the energy stored in the cavity $(E_{{stored}})$ through the concept of the cavity’s quality factor, or $Q$ factor. As discussed by Hahn et al. (Reference Hahn, Leibacher, Baasch and Dual2015), this relationship is expressed as

(8.8) \begin{equation} Q = \frac {\omega ^{nml}}{P^{nml}} \int _{V_0} E^{nml} \, \textrm {d} V = \frac {\omega ^{n m l}}{\Delta \omega ^{n m l}}. \end{equation}

Here, the stored energy is derived from the integration of the energy density over the cavity volume $V_0$ . The term $\Delta \omega ^{n ml}$ represents the ${-3}\,{\rm dB}$ bandwidth of the resonance line, which is determined by solving the equation $|R^{n m l}(\omega )|^2 = ({1}/{2})$ , subsequently leading to $\Delta \omega ^{n m l} = \omega ^{n m l} \Gamma ^{n m l}.$ Replacing the bandwidth into (8.8) gives the $Q$ factor as

(8.9) \begin{equation} Q =\frac {1}{\Gamma ^{nml}}. \end{equation}

Consequently, the dissipated power by the cavity is

(8.10) \begin{equation} P^{nml}= \Gamma ^{nml}\omega ^{nml}\int _{V_0} E^{nml} \, \textrm {d} V. \end{equation}

8.3. Acoustic radiation force on a small particle

An important application of acoustofluidic microcavities is the contactless manipulations of small particles. We consider a spherical particle of volume $V_{{p}}$ , isentropic compressibility $\kappa _{{p}}$ and density $\rho _{{p}}$ . Moreover, the particle in consideration is assumed to be much smaller than the wavelength. In this case, the acoustic radiation force exerted on a particle at position $\boldsymbol{r}$ is given by (Gor’kov Reference Gor’kov1962)

(8.11a) \begin{align} \boldsymbol{F}&= -V_{{p}}\boldsymbol{\nabla } U, \end{align}

(8.11b) \begin{align} U &= f_\kappa {E}_{{p}} -\frac {3}{2}f_\rho {E}_{{k}}, \end{align}

(8.11c) \begin{align} f_\kappa &= 1 - \frac {\kappa _{{p}}}{\kappa _0}, \quad f_\rho = 2 \frac {\rho _{{p}} - \rho _0}{2\rho _{{p}}+ \rho _0}, \end{align}

where $U$ is the Gor’kov potential density, and $f_\kappa$ and $f_\rho$ are respectively the compressional and density contrast factors. After inserting the kinetic (8.2b ) and potential (8.2c ) energy densities into (8.11b ), we find the dimensionless Gor’kov potential density $(\tilde {U}=U/E_0)$ as

(8.12) \begin{equation} \tilde {U} = \left [ f_\kappa - \frac {3}{2} f_\rho \left ( \frac {1 }{2 } \tilde {\nabla }^2 +1 \right )\right ] \left |\tilde {p}_1\right |^2. \end{equation}

Substituting (6.6a ) into this equation yields the single-mode Gor’kov potential

(8.13) \begin{equation} \tilde {U}^{nml} = \big |A^{nml} R^{nml} \big|^2 \left [ f_\kappa - \frac {3}{2} f_\rho \left (\frac {1 }{2 } \tilde {\nabla }^2+1 \right )\right ] \big|\Psi ^{nml}\big|^2. \end{equation}

The Gor’kov potential of a dual mode is obtained by replacing (7.4) into (8.12),

(8.14) \begin{equation} \tilde {U}^{n_1 m_1l_1}_{nml} = \tilde {U}^{nml} + \tilde {U}^{ n_1 m_1 l_1} + \tilde {U}_{{cross}}. \end{equation}

The cross term of the dual-mode potential is

(8.15) \begin{align} \tilde {U}_{{cross}} = 2 A^{nml} A^{n_1m_1 l_1} {\mathrm{Re}} \! \left \{\big (R^{nml}\big )^* R^{n_1m_1 l_1} \left [ f_\kappa - \frac {3}{2} f_\rho \left ( \frac {1 }{2 } \tilde {\nabla }^2 +1 \right )\right ] {\Psi ^{nml}}^* \Psi ^{n_1 m_1l_1}\right \}. \end{align}

8.4. Spin angular momentum

The spin angular momentum density of an acoustic wave can be expressed by (Bliokh & Nori Reference Bliokh and Nori2019; Lopes et al. Reference Lopes, Lima, Leão Neto and Silva2020)

(8.16) \begin{equation} \boldsymbol{S} = \frac {\rho _0}{2 \omega }{\mathrm{Im}}\!\left [{\boldsymbol{v}}_{{c}}^* \times {\boldsymbol{v}}_{{c}}\right ], \end{equation}

which is given in the SI units of $\mathrm{J\,s\,m^{-3}}$ . Substituting (2.10a ) into (8.16) allows us to express the acoustic spin as

(8.17) \begin{align} \boldsymbol{{S}} &=\frac {1}{2\rho _0\omega ^3}{\mathrm{Im}} \!\left [ \boldsymbol{\nabla }\! {p}_1^* \times \boldsymbol{\nabla }\! {p}_1\right ]= \frac {1}{2\rho _0\omega ^3}{\mathrm{Im}} \!\left [\boldsymbol{\nabla } \times \left ({p}_1^*\boldsymbol{\nabla } \!{p}_1\right ) \right ]= \frac {1}{\omega ^2}\boldsymbol{\nabla } \times \dfrac {1}{2}\mathrm{Re}\! \left [ {p}_1^* {\boldsymbol{v}}_c\right ]\nonumber\\&= \frac {1}{\omega ^2}\boldsymbol{\nabla } \times \boldsymbol{I}. \end{align}

It is important to note that the spin angular momentum is a measure of the circulation of acoustic intensity, defined as $\boldsymbol{I} = ({1}/{2}) \text{Re} [ {p}_1^* \boldsymbol{v}_c ]$ . The acoustic spin density forms a solenoidal field, meaning $\nabla \cdot \boldsymbol{S} = 0$ .

We define the dimensionless acoustic spin density as $\tilde {\boldsymbol{S}}= \omega \boldsymbol{S}/E_0$ . To obtain the spin of a single-vortex mode at the resonance, $\omega =\omega ^{nml}$ , we substitute (6.6a ) into (8.17),

(8.18) \begin{equation} \tilde {\boldsymbol{S}}^{nml} = 2 \big|{A}^{nml}R^{nml} \big |^2{\mathrm{Im}}\! \left [ \tilde {\boldsymbol{\nabla }} \times \big( {\Psi ^{nml} }^* \tilde {\boldsymbol{\nabla }}\Psi ^{nml} \big) \right ]. \end{equation}

For the dual-mode configuration, one has

(8.19) \begin{align} \tilde {\boldsymbol{S}}_{nml}^{n_1m_1l_1} = \tilde {\boldsymbol{S}}^{nml} + \tilde {\boldsymbol{S}}^{ n_1 m_1l_1} + \tilde {\boldsymbol{S}}_{{cross}}, \end{align}

with the spin cross term being

(8.20) \begin{align} \nonumber \tilde {\boldsymbol{S}}_{{cross}} &= 2 A^{nml} A^{n_1m_1l_1} {\mathrm{Im}}\big \{ \tilde {\boldsymbol{\nabla }} \times \big [ \big (R^{nml}\big )^* R^{n_1m_1 l_1} \big (\Psi ^{nml} \big )^* \tilde {\boldsymbol{\nabla }}\Psi ^{n_1m_1l_1}\big . \big .\\ &\quad \big . \big .+ \big (R^{n_1m_1l_1}\big )^* R^{nm l}\big (\Psi ^{n_1 m_1l_1}\big )^* \tilde {\boldsymbol{\nabla }}{\Psi ^{nml} } \big ]\big \}. \end{align}

8.5. Acoustic spin and radiation torque connection

When a particle is immersed in an acoustic field, it may experience a radiation torque resulting from the transfer of angular momentum from the field to the particle. In the long-wavelength scattering regime, the acoustic radiation torque $\boldsymbol{\tau }$ exerted on an absorptive spherical particle (by a vortex beam or vortex standing wave) is connected to the spin density $\boldsymbol{S}$ through the relation (Silva Reference Silva2014, (29))

(8.21a) \begin{align} \boldsymbol{\tau } &= \gamma \boldsymbol{S}, \end{align}

(8.21b) \begin{align} \gamma &= -\frac {12 \pi c_0}{k^2} \big( {\mathrm{Re}}\!\left [s_1\right ]+ \left |s_1\right |^2{\kern-2pt} \big ), \end{align}

where $\gamma$ is the gyroacoustic factor and $s_1$ is the dipole scattering coefficient, which depends on energy dissipation by the particle. The derivation of (8.21a ) is based on the angular momentum flux formalism (Maidanik Reference Maidanik1958).

Using (8.21a ), we recover the expression for the axial radiation torque exerted by a Bessel vortex beam on a spherical particle of radius $ a$ , located at position $ \boldsymbol{r}$ , as previously derived by Zhang (Reference Zhang2018) and on-axis by Zhang & Marston (Reference Zhang and Marston2011 b). In the notation of Zhang (Reference Zhang2018), the dipole scattering coefficient is given by $ s_1 = 1 + \textrm {i} (ka)^3 (f_{2\textit{r}} + \textrm {i} f_{2\textit{i}})/6,$ where $ f_{2\textit{r}} = f_\kappa$ – as defined in (8.11c ) – and $ f_{2\textit{i}} = (3 \delta _{\textit{p}}/2a) f_\kappa ^2$ , with $ \delta _{\mathrm{p}}$ denoting the thickness of the viscous boundary layer surrounding the particle. This leads to the expression $ \gamma = 2\pi a^3 \omega f_{2\textit{i}}$ . The axial spin density $ S_z$ of an $ n$ th-order Bessel vortex beam is given by Lopes et al. (Reference Lopes, Lima, Leão Neto and Silva2020) as $S_z = 2 n (I_0/c_0) \sin \beta \, J_n(k \varrho \sin \beta ) J_n'(k \varrho \sin \beta ) / \omega k \varrho$ , where $ I_0 = p_0^2/2\rho _0 c_0$ is the time-averaged beam intensity, $ p_0$ is the pressure amplitude, $ \beta$ is the conical angle of the beam’s wave vector, and the prime symbol denotes differentiation with respect to the function’s argument. Substituting $ \gamma$ and $ S_z$ into (8.21a ), we obtain the axial radiation torque,

(8.22) \begin{equation} \tau _z = 6\pi n f_\kappa ^2 \delta _{{p}} a^2 \frac {I_0}{c_0} \sin ^2\beta \, J_n'(k \varrho \sin \beta ) \frac {J_n(k \varrho \sin \beta )}{k \varrho \sin \beta }. \end{equation}

This expression exactly matches the axial radiation torque derived by Zhang (Reference Zhang2018, (14a)) and Zhang & Marston (Reference Zhang and Marston2011b , (10)). We thus confirm the general validity of the spin-torque relation (8.21a ), which applies to spherical particles undergoing power absorption due to either internal dissipation or viscous losses in the surrounding fluid.

It is also important to mention that the radiation torque exerted by an acoustic wave (acoustic beam or standing wave) on a spherical particle of arbitrary size has also been derived using the angular momentum flux formalism (Zhang & Marston Reference Zhang and Marston2011 b; Silva et al. Reference Silva, Lobo and Mitri2012; Zhang & Marston Reference Zhang and Marston2014). These works demonstrate that torque arises from energy dissipation during scattering, which can occur either within the particle itself or in the surrounding viscous boundary layer. In the on-axis case, the proportionality of torque on power absorption for spheres in standing waves or beams (Zhang & Marston Reference Zhang and Marston2011 b, (10)) is not limited to the case of absorption associated with the dipole response emphasised in (8.22).

Recently, Bliokh & Nori (Reference Bliokh and Nori2019) formulated the concept of acoustic spin density within the framework of canonical field theory, drawing parallels with its electromagnetic analogue. Building on this formalism, Toftul et al. (Reference Toftul, Bliokh, Petrov and Nori2019) derived the spin-torque relationship given in (8.21a ) for isotropic particles. Subsequently, Lima et al. (Reference Lima, Leão-Neto, Marques, Silva, Lopes and Silva2020)and Lopes et al. (Reference Lopes, Lima, Leão Neto and Silva2020) extended the spin-torque relation to encompass small axisymmetric particles, such as spheroids, thereby broadening its applicability to more complex geometries.

9.1. Acoustic single modes

Two modes of particular interest are the axial and the transverse modes. The energy landscape of these modes near the resonance frequency $(\omega \approx \omega ^{nml})$ are derived by substituting (6.6) into (8.5). The result yields respectively

(9.1a) \begin{align} \tilde {E}^{l} &= \big |A^l R^{l} \big|^2, \end{align}

(9.1b) \begin{align} \tilde {E}^{nm} &= \left |A^{nm} R^{nm} \right |^2 \left [J_n\!\left (\tilde {\varrho }^{nm} \right )\left [\frac {J_n^{\prime}\!\left (\tilde {\varrho }^{nm}\right )}{\tilde {\varrho }^{nm}} + J_n^{\prime\prime}\!\left (\tilde {\varrho }^{nm}\right )\right ] +\left [J_n'\!\left (\tilde {\varrho }^{nm}\right )\right ]^2 + 2 J_n^2\!\left (\tilde {\varrho }^{nm}\right ) \right ], \end{align}

where $A^l=A^{00l}$ , $A^{nm}=A^{nm0}$ , $R^{l}=R^{00l}$ and $R^{nm}=R^{nm0}$ , and $\tilde {\varrho }^{nm}=k_\varrho ^{nm}\varrho$ . Note that the energy density of the vortex mode exhibits circular symmetry, while it remains spatially constant for the axial mode.

Now, we turn our analysis to the Gor’kov potential near the resonance frequency. For the axial mode $(0\,0\,l)$ , this potential is obtained from (8.13). Accordingly, we have

(9.2) \begin{equation} \tilde {U}^{l} = \tfrac {1}{2} \big |A^l R^l\big |^2 \big [\big (2f_\kappa + 3f_\rho \big ) \cos ^2 \! \tilde {z}^l - 3f_\rho \big ], \end{equation}

where $\tilde {z}^l=k_z^l z$ . The particles will be trapped at the Gor’kov potential minimum by the action of the acoustic radiation force. When $2f_\kappa + 3f_\rho \gt 0$ , the potential minima is in the axial pressure nodes at $z_{{min}}=(2\ell +1)H/2l$ , with $\ell =0,1,\ldots ,l-1$ . If $2f_\kappa + 3f_\rho \lt 0$ , the minima is at the pressure antinodes in $z_{{min}} = \ell H/l$ , where $\ell =0,1,\ldots ,l$ . The axial radiation force per unit length is obtained by taking the negative gradient of (9.2):

(9.3) \begin{equation} {F}_z^l =\frac {E_0 V_{{p}}}{2} \big |A^l R^l\big |^2 \big (2f_\kappa + 3f_\rho \big ) \sin \!\big ( 2\tilde {z}^l \big ). \end{equation}

Considering the transverse mode $(n\,m\,0)$ , we find from (8.13) that the Gor’kov potential is

(9.4) \begin{align} \nonumber \tilde {U}^{nm} &= \frac {1}{2} \left |A^{nm} R^{nm} \right |^2 \left \{ \left (2f_\kappa - 3f_\rho \right ) \left [J_n\!\left (\tilde {\varrho }^{nm}\right )\right ]^2 \vphantom{\frac {J_n'\!\left (\tilde {\varrho }^{nm}\right )}{\tilde {\varrho }^{nm}}} \right .\\ &\quad\left . - 3 f_\rho \left [J_n'\!\left (\tilde {\varrho }^{nm}\right )\right ]^2 - 3 f_\rho J_n\!\left (\tilde {\varrho }^{nm} \right ) \left [ J_n^{\prime\prime}\!\left (\tilde {\varrho }^{nm} \right ) + \frac {J_n'\!\left (\tilde {\varrho }^{nm}\right )}{\tilde {\varrho }^{nm}} \right ]\right \}. \end{align}

The Gor’kov potential here has no azimuthal dependence on $\varphi$ . For the sake of brevity, we will analyse the radial trapping around the symmetry axis of the cavity only. To do so, we use the asymptotic form of the Bessel function for small arguments, $J_n(x) \approx 2^{-n} x^n ( 1/n! - 1/4[(n+1)!])$ , with $x \rightarrow 0$ . The radial radiation force of the non-vortex and first vortex mode ( $n=0,1$ ) are given by

(9.5a) \begin{align} {F}_\varrho ^{0m} &= f_\kappa E_0 V_{{p}}\big |A^{0m} R^{0m} \big |^2 k_\varrho ^{0m} \varrho + \textrm {O}\big (\varrho ^3\big ), \end{align}

(9.5b) \begin{align} {F}^{1m}_\varrho &= -\big ( 2f_\kappa + 3f_\rho \big ) E_0 V_{{p}} \big |A^{1m} R^{1m} \big |^2 \varrho + \textrm {O}\big (\varrho ^3\big ). \end{align}

Some conclusions can be drawn from the analysis here. Particles in the non-vortex mode $(0\,m)$ can either be repelled from $(f_\kappa \gt 0)$ or attracted to $(f_\kappa \lt 0)$ the cavity axis. In the first-vortex mode $(1\,m)$ , particles will be trapped at the cavity axis if $2f_\kappa + 3f_\rho \gt 0$ . Conversely, they will be repelled if $2f_\kappa + 3f_\rho \lt 0$ .

We proceed with the derivation of the spin angular momentum for modes with $|n|\gt 0$ . Substituting (6.6d ) into (8.18), we find the only non-zero spin component is the axial one, which is given by

(9.6) \begin{align} \tilde {S}_z^{nm} = 2n \big|A^{nm}R^{nm} \big|^2 \frac { \big[J_n^2\!\big (k_\varrho ^{nm} \varrho \big )\big]^{\prime}}{k_\varrho ^{nm}\varrho }. \end{align}

Modes with no vorticity ( $n=0$ ) lack spin. Note that the acoustic radiation torque exerted on a small spherical particle by a Bessel vortex wave beam or standing wave) derived by Zhang (Reference Zhang2018, (12)) is proportional to the axial spin given in (9.6). Moreover, a detailed examination of (9.6), using the expression $\lim _{x \rightarrow 0} [{J_n^2(x)} ]'/x = {\delta _{n,\pm 1}}/{2}$ , indicates that only the $(\pm 1\,m)$ mode exhibits spin along the $z$ axis. Hence,

(9.7) \begin{equation} \lim _{\varrho \rightarrow 0}\tilde {S}_z^{\pm 1m} = \pm \big |A^{1m}R^{1m} \big|^2. \end{equation}

The amplitudes of the energy density and Gor’kov potential as well as the spin density can be obtained using (5.4b ), (6.5b ) and (6.6b ). It can be shown that

(9.8) \begin{equation} \big|A^l\big|^2, \big|A^{nm}\big|^2 \sim \textit{Sh}^2= \left (R/\delta \right )^2. \end{equation}

The parameter $\textit{Sh}=R/\delta$ is known as the shear wave number of the cavity (Tijdeman Reference Tijdeman1975). This suggests that for decoupled axial and vortex resonant modes, increasing the shear wave number ( $\textit{Sh}$ ) reduces the relative influence of boundary-layer-induced viscous dissipation on the bulk resonance response.

9.2. Acoustic dual modes

One potential application of dual modes in cylindrical cavities is to levitate cells or microparticles, providing an additional degree of freedom: rotation around an axis of symmetry. The dual mode, composed of an axial mode $(0\,0\,l)$ , ensures the formation of $l$ levitation planes at the mid-height of the cavity, while circular modes $(n\,m\,0)$ are responsible for clustering particles and inducing a spin-torque on the trapped particles. From (7.2), we note that the cavity’s aspect ratio to host the dual mode $ (0\,0\,l) + ( n\, m\, 0)$ is

(9.9) \begin{equation} \frac {R}{H} = \frac {j^{\prime}_{nm}}{l \pi }. \end{equation}

Moreover, according to (8.14) and (8.19), the Gor’kov potential and spin density of this mode near the resonance $(\omega \approx \omega ^{l}\approx \omega ^{nm0})$ can be expressed by

(9.10a) \begin{align}& \tilde {U}_{l}^{nm} = \tilde {U}^{l} + \tilde {U}^{nm}+ \tilde {U}_{{cross}}, \end{align}

(9.10b) \begin{align}&\quad \tilde {\boldsymbol{S}}_{l}^{nm} = \tilde {{S}}^{nm}_z \boldsymbol{e}_z + \tilde {\boldsymbol{S}}_{{cross}}. \end{align}

The cross terms are calculated with (8.15) and (8.20). Accordingly, we have

(9.11a) \begin{align} \nonumber \tilde {U}_{{cross}}&= A^lA^{nm}\big|R^{l}R^{nm}\big | J_{n_1}\big(\tilde {\varrho }^{nm}\big) \bigg [\big(2f_\kappa - 3 f_\rho \big) \cos \!\left (n\varphi - \Delta \varphi _l^{nm}\right )\vphantom{\frac {k_z^l}{k}} \\ &\quad + \;3f_\rho \bigg(\frac {k_z^l}{k}\bigg)^2\cos \!\left (n\varphi + \Delta \varphi _l^{nm}\right ) \bigg ]\cos \tilde {z}^l, \end{align}

(9.11b) \begin{align} \nonumber \tilde {\boldsymbol{S}}_{{cross}} &= 4 A^l A^{nm} \big|R^{l}R^{nm} \big| \bigg[ n \cos \big(n \varphi +\Delta \varphi _l^{nm}\big) \frac {J_{n}\big(\tilde {\varrho }^{nm}\big)}{\tilde {\varrho }^{nm}} \boldsymbol{e}_\varrho \\ & \quad -\sin\!\big(n \varphi +\Delta \varphi _l^{nm}\big) \, J_{n}'\big(\tilde {\varrho }^{nm}\big) \,\boldsymbol{e}_\varphi \bigg]\sin ^2 \tilde {z}^l, \end{align}

where the additional phase is $\Delta \varphi _l^{nm}=\angle (R^l)^*R^{n m}$ .

The spin component in (9.11b ) arises from the nonlinear interaction between two orthogonal modes, $(0\,0\,l)$ and $(n\,m\,0)$ . This component can be transferred to a suspended particle much smaller than the wavelength, resulting in an acoustic radiation torque on the particle (Zhang & Marston Reference Zhang and Marston2014). Consequently, the particle may rotate around one of its principal axes. Notably, even when $n=0$ , a spin component is still present, which is given by

(9.12) \begin{align} \tilde {\boldsymbol{S}}_{{cross}} &= 4 A^l A^{0m}\big|R^{l}R^{0m}\big| \sin \!\big(\Delta \varphi _l^{0m}\big) \, J_{1}\big(\tilde {\varrho }^{0m}\big) \sin ^2 \tilde {z}^l \,\boldsymbol{e}_\varphi . \end{align}

This result demonstrates that, even in the absence of vortex excitation, the spin angular momentum can still be generated in acoustofluidic resonators. It is also worth noting that metallic nanorods have been observed to rotate in acoustofluidic cylindrical cavities with a non-vortex excitation at the base (Wang et al. Reference Wang, Castro, Hoyos and Mallouk2012; Balk et al. Reference Balk, Mair, Mathai, Patrone, Wang, Ahmed, Mallouk, Liddle and Stavis2014). This rotational behaviour can be attributed to the spin-torque generated by the nonlinear interaction of orthogonal modes within the cavity.

Finally, we note that the cross terms also scale as $A^l A^{nm}\sim \textit{Sh}^2$ . We thus conclude that the amplitude of the Gor’kov potential and spin density of a dual mode also scales with shear wave number squared.

10.1. Single modes

In the results for single modes, a water-filled microcavity with dimensions $H=R={250}\, {{\mu} \textrm{m}}$ and a volume of $V_0={49}\, {\mathrm{nL}}$ is considered, representing a typical cavity in acoustofluidics devices (Rocha et al. Reference Rocha, Silva, Sales, D’Amato, Leite and Silva2023).

Figure 3.

Here and in the subsequent figures, the $(n\, m\, l)$ notation is defined in (6.6). The $(0\,0\,1)$ and $(0\,1\,0)$ mode induced in a microcavity with $H=R={250}\, {{\mu} \textrm{m}}$ . (a,c) The pressure-frequency plot (magnitude and phase in blue and red line) using (6.6a ), and the corresponding resonance frequency and $Q$ factor. Blue and red dots indicate numerical results. (b) Gor’kov potential in the $x{-}z$ plane (background image) obtained from (9.2), with the related radiation force field (white arrows) acting on a ${10}\, {{\mu} \textrm{m}}$ -polystyrene probe particle. (d) The energy density computed with (9.1b) and (e) the Gor’kov potential for the probe particle calculate with (9.4). All fields are obtained at the corresponding resonance frequency. The vertical white lines delimit the theoretical (left) and numerical (right) results.

Figure 3 presents the $(0\,0\,1)$ and $(0\,1\,0)$ mode induced in the microcavity. The physical parameters of water used in this analysis are given in table 1, while panels (a) and (c) plot pressure versus frequency as given in (6.6a ), with $Q=796.10$ and $622.01$ , respectively. The resonance frequencies are ${2.9917}\, {\mathrm{MHz}}$ and ${3.6494}\, {\mathrm{MHz}}$ . Solid lines denote the theoretical result, while dots represent numerical data. Panel (b) depicts the Gor’kov potential density (background image) of the $(0\,0\,1)$ mode at the $xz$ plane, with white arrows illustrating the radiation force field. Probe particles will collect in the pressure node at the cavity’s midheight. Panels (d) and (e) show the energy density and Gor’kov potential of the $(0\,1\,0)$ mode. The related radiation force will trap particles in an annular region (in dark brown) in panel (e). The vertical white lines separate the theoretical (left) and numerical (right) results. The numerical data and the theoretical predictions show excellent agreement with $\mathcal{R}^2\sim 1$ . It is important to mention that in the numerical simulations, the degrees of freedom (DOFs) for the two-dimensional axisymmetric model used for the $(0\,0\,1)$ and $(0\,1\,0)$ modes is $288\,045$ . Larger numbers of DOF typically require more computational resources, both in terms of memory and processing power.

In figure 3(e), we note that the radiation force is zero in a ring with radius, denoted by $\varrho _{{r}}$ , which is a solution of the transcendent equation derived by inserting (9.4) into (8.11a ):

(10.1) \begin{equation} 3 f_\rho J_2\big(k_\varrho ^{01} \varrho _{{r}}\big )- (4 f_\kappa + 3 f_\rho ) J_0\big(k_\varrho ^{01} \varrho _{{r}}\big ) = 0. \end{equation}

The numerical solution to this equation yields $\varrho _{{r}} = 0.6117 R$ or $\varrho _{{r}}={152}\,{{\mu \textrm{m}}}$ , which is compatible with the particle trap region seen in panel (e).

Figure 4.

Pure vortex mode $(1\,0\,0)$ and $(2\,0\,0)$ cases are considered here for the cavity described in figure 3. (a,e) Acoustic pressure plotted against frequency (magnitude and phase in blue and red line) using (6.6a ), along with the associated $Q$ factor and resonance frequency. The dots represent the numerical data. (b,f) Energy density calculated with (9.1b ). (c,g) Illustration of the Gor’kov potential for a ${10}\,{{\mu \textrm{m}}}$ polystyrene probe particle using (9.4). (d,h) Spin angular momentum according to (9.6). These fields are computed at the corresponding resonance frequency of each mode. The vertical white lines separate the theoretical (left) and numerical (right) results.

In figure 4, we illustrate the $(1\,0\,0)$ and $(2\,0\,0)$ mode within the same cavity introduced in figure 3. Panels (a) and (c) present the pressure plotted against frequency using (6.6a ). The associated $Q$ factor and resonance frequency are $436.97$ and ${1.7523}\,{\mathrm{MHz}}$ , and $453.37$ and ${2.9070}\,{\mathrm{MHz}}$ , respectively. The energy density given in (9.1b ) is shown in panels (b) and (f), while the Gor’kov potential from (9.4) for a ${10}\,{{\mu \textrm{m}}}$ polystyrene particle are seen in panels (c) and (g). The spin angular momentum, as given by (9.6), is shown in panels (d) and (h). These modes exhibit no levitation plane, indicating their purely vortex nature. Notably, the inward radiation force field is responsible for concentrating particles in the central portion of the resonator without levitation. The spatial pattern of the axial spin differs significantly between these modes. In the $(1\,0\,0)$ mode, collected particles experience a positive spin along the $z$ axis, peaking at $\varrho =0$ . Conversely, for the $(2\,0\,0)$ mode, the maximum axial spin, according to (9.6), occurs in a ring of radius $\varrho _{{s}}$ , which is determined by the solution of

(10.2) \begin{equation} \big [J_2^2\big(k_\varrho ^{20} \varrho _{{s}}\big )\big ]^{\prime\prime} =0 . \end{equation}

We find $\varrho _{{s}}={0.6759}R$ , yielding $\varrho _{{s}} = {169}\,{{\mu} \textrm{m}}$ . The axial spin potentially may induce a counterclockwise rotation on particles mostly depending on their geometric asymmetry. Note that the spin vanishes at $\varrho =0$ . In the numerical simulations, it was observed that the DOF of the three-dimensional model for the $(1\,0\,0)$ mode was $2\,705\,168$ . However, the DOF for the $(2\,0\,0)$ mode is $2\,475\,000$ . According to the coefficient of determination, $\mathcal{R}^2\sim 1$ , and visual inspection, the theoretical and numerical data show excellent agreement.

10.2. Dual modes

Here we consider the $(0\,1\,0)+(0\,0\,1)$ , $(1\,0\,0)+(0\,0\,1)$ and $(2\,0\,0)+(0\,0\,1)$ modes, as they correspond to the simplest modes for trapping and inducing torque on suspended particles. To develop these modes, the cavity’s aspect ratio should be respectively $R/H={1.2197}, {0.5860}$ , $0.9721$ . The height of all cavities is fixed to $H={250}\, {{\mu} \textrm{m}}$ . Also, the cavities are filled with water.

Figure 5.

The $(0\,1\,0)+(0\,0\,1)$ mode for a cavity with aspect ratio $R/H={1.2197}$ . Panel (a) shows the magnitude (in blue) and phase (in red) of the acoustic pressure around the resonance frequency ${2.992}\, {\mathrm{MHz}}$ , plotted using (6.6a ). The numerical data are represented by dots. Panel (b) presents the Gor’kov potential density for the probe particle, with white arrows depicting the corresponding acoustic radiation force. The theoretical part of the potential is calculated using (9.10a ). Panel (c) exhibits the Gor’kov potential at the cavity mid-height, while panel (d) illustrates the magnitude of the transverse spin density, with white arrows denoting the related vector field. The theoretical part of the spin is computed with (9.10b ). The acoustic fields in panels (b), (c) and (d) are computed at the resonance frequency. The vertical white lines separate the theoretical (left) and numerical (right) results.

Figure 6.

The $(1\,0\,0)+(0\,0\,1)$ mode for a cavity with $R/H={0.5860}$ and $R={250}\,{{\mu} \textrm{m}}$ . Panel (a) show the acoustic pressure (magnitude in blue and phase in red) around the resonance frequency ${2.991}\, {\mathrm{MHz}}$ , and plotted against frequency using (6.6a ). The dots represent the numerical data. Panel (b) presents the background map of the Gor’kov potential density for the probe particle. The theoretical part is computed using (9.10a ). Panel (c) exhibit the Gor’kov potential at the mid-height, while panels (d) and (e) illustrate the magnitude of the transverse and axial components of the spin density at the mid-height, respectively. The theoretical part of these components are calculated with (9.10b ). The white arrows in panel (d) correspond to the spin vector field. The acoustic fields are computed at the resonance frequency. The vertical white lines separate theoretical (left) and numerical (right) results.

Figure 7.

The $(2\,0\,0)+(0\,0\,1)$ mode for a cavity with an aspect ratio of $R/H = {0.9721}$ and radius $ R = {250}\,{{\mu} \textrm{m}}$ . Panels (a) and (e) display the acoustic pressure around the resonance at ${2.992}\,{\mathrm{MHz}}$ , plotted as a function of frequency using (6.6a ), with numerical data shown as dots. Panel (b) depicts the Gor’kov potential density for a probe particle, calculated using (9.10a ). Panel (c) shows this potential at mid-height with the corresponding radiation force (white arrows), while panels (d) and (e) present the transverse and axial spin density magnitudes, calculated via (9.10b ). White arrows in panel (d) represent the spin vector field. The presented acoustic fields are computed at resonance, with vertical white lines dividing the theoretical (left) and numerical (right) results.

Figure 5(a) displays the acoustic pressure for the $(0\,1\,0)+(0\,0\,1)$ mode near the resonance frequency of ${2.992}\, {\mathrm{MHz}}$ , based on numerical simulations with 331 164 DOF. Panel (b) depicts the Gor’kov potential density in the $xz$ plane, along with the radiation force field shown by the white arrows. This potential is axisymmetric. Panel (c) shows the Gor’kov potential in the $x{-}y$ plane at the mid-height of the cavity, where the radiation force field (white arrows) traps particles in an annular region around the cavity central axis. Panel (d) presents the transverse spin component, with white arrows representing the related vector field. This spin arises due to the interaction between the transverse and axial velocity fields in the fluid. The spin exhibits characteristics of a solenoidal field (circulation). The theoretical and numerical results are in excellent agreement, with $\mathcal{R}^2 \sim 1$ , as seen by visual inspection.

Figure 6(a) shows the acoustic pressure for the $(1\,0\,0)+(0\,0\,1)$ mode around the resonance frequency of ${2.991}\,{\mathrm{MHz}}$ . A significant increase in the DOF (1 900 422) is noted, since this simulation requires a three-dimensional (3-D) mesh. Panel (b) presents the Gor’kov potential density in the $x{-}z$ plane, along with the radiation force field (white arrows). The potential in this plane is asymmetric due to the cross term in (9.11a ), which depends on the azimuthal angle $\varphi$ . In contrast, panel (c) shows a symmetric Gor’kov potential (with $U_{{cross}}=0)$ in the mid-height plane, where the radiation force field directs particles towards the cavity’s central axis. Panels (d) and (e) illustrate the transverse and axial spin components, respectively. The white arrows in panel (d) indicate the transverse spin field at mid-height, which results from the interaction between the transverse and axial fluid velocity fields. The axial spin component is most pronounced in the central region of the cavity, but diminishes near the lateral walls, while the transverse spin is predominantly oriented left-to-right and is approximately four times stronger than the axial spin. A circulating pattern is observed in the transverse spin field, consistent with the behaviour expected of a solenoidal field. We see that theoretical and numerical results show excellent agreement, with $\mathcal{R}^2 \sim 1$ , confirmed by visual inspection.

Figure 7(a) displays the acoustic pressure for the $(2\,0\,0)+(0\,0\,1)$ mode near the resonance frequency of ${2.992}\, {\mathrm{MHz}}$ , with 1 900 422 DOF in the numerical simulations. Panel (b) shows the Gor’kov potential density and radiation force field (white arrows) in the $x{-}z$ plane. Conversely, panel (c) illustrates the Gor’kov potential in the $x{-}y$ plane at mid-height, with the radiation force (white arrows) pointing towards the cavity’s axis of symmetry. The potential well is more shallow compared with the one in figure 6. Moreover, the simulation results exhibit some numerical fluctuations in the computed radiation force field. These fluctuations can, in principle, be mitigated by increasing the mesh resolution, though this would require greater computational memory. Panels (d) and (e) depict the transverse and axial spin components, respectively, with white arrows in panel (d) representing the transverse spin field at mid-height. The axial spin is more dominant in the annular region around the central axis. The transverse spin shows a circulating pattern around the dark spots at the edge of the cavity. Finally, theoretical and numerical results align closely, with $\mathcal{R}^2 \sim 1$ .