Bulk-driven acoustic streaming at resonance in closed microcavities

TL;DR

This paper investigates how bulk-driven acoustic streaming, especially in double-mode resonances within closed microcavities, can generate rotating flow patterns, challenging the traditional boundary-driven streaming dominance in acoustofluidic devices.

Contribution

It derives analytical solutions for double modes in rectangular cavities, revealing how rotating bulk-driven streaming can be excited by non-rotating actuation and identifying optimal geometries for such effects.

Findings

Rotating bulk-driven streaming can be excited without rotating actuation.

High-frequency patterns are dominated by bulk-driven streaming.

Optimal cavity geometries maximize the excitation of rotating streaming patterns.

Abstract

Bulk-driven acoustic (Eckart) streaming is the steady flow resulting from the time-averaged acoustic energy flux density in the bulk of a viscous fluid. In simple cases, like the one-dimensional single standing-wave resonance, this energy flux is negligible, and therefore the bulk-driven streaming is often ignored relative to the boundary-driven (Rayleigh) streaming in the analysis of resonating acoustofluidic devices with length scales comparable to the acoustic wavelength. However, in closed acoustic microcavities with viscous dissipation, two overlapping resonances may be excited at the same frequency as a double mode. In contrast to single modes, the double modes can support a steady rotating acoustic energy flux density and thus a corresponding rotating bulk-driven acoustic streaming. We derive analytical solutions for the double modes in a rectangular box-shaped cavity including the viscous boundary layers, and use them to map out possible rotating patterns of bulk-driven acoustic streaming. Remarkably, the rotating bulk-driven streaming may be excited by a non-rotating actuation, and we determine the optimal geometry that maximizes this excitation. In the optimal geometry, we finally simulate a horizontal 2-by-2, 4-by-4, and 6-by-6 streaming-roll pattern in a shallow square cavity. We find that the high-frequency 6-by-6 streaming-roll pattern is dominated by the bulk-driven streaming as opposed to the low-frequency 2-by-2 streaming pattern, which is dominated by the boundary-driven streaming.

Tables

Table 1. Table 1: The damping coefficients Γ Γ \Gamma , the full-width-half-maximum (FWHM) line widths, and the mode separation Δ Δ \Delta .

Parameter	Symbol	Reference
Bulk damping	${\Gamma }_{\mathrm{fl}}$	Eq. (4)
Boundary-layer damping	${\Gamma }_{\mathrm{bl}}^{lmn}$	Eq. (10)
Total damping	${\Gamma }^{lmn}$	Eq. (11)
FWHM of $\left|F^{lmn}\right|$	$\sqrt{3}{k}_0^{lmn}{\Gamma }^{lmn}$	Fig. 1
FWHM of $E_{\mathrm{ac}}^{lmn}\propto {\left|F^{lmn}\right|}^2$	$k_0^{lmn}{\Gamma }^{lmn}$	Eqs. (6),(12)
Optimal mode separation ${\Delta }_{lm}^{\mathrm{opt}}$	$\frac{1}{\sqrt{3}}{k}_{lm}^{\star }{\Gamma }_{lm}^{\star }$	Eq. (29b)

Table 2. Table 2: Material parameters for water (Muller and Bruus, 2014 ) and pyridine (C 5 H 5 N) (Dukhin and Goetz, 2009 ) at 25 C ∘ used in the numerical simulations.

Parameter	Symbol	Water	Pyridine	Unit
Mass density	${\rho }_{\mathrm{fl}}$	997.05	982	kg m-3
Compressibility	${\kappa }_{\mathrm{fl}}$	448	507	TPa-1
Speed of sound	$c_{\mathrm{fl}}$	1496.7	1417	m s-1
Dynamic viscosity	${\eta }_{\mathrm{fl}}$	0.890	0.879	mPa s
Bulk viscosity	${\eta }_{\mathrm{fl}}^b$	2.485	62.4	mPa s

Table 3. Table 3: Parameters used in the numerical simulations of the double modes l ​ 00 + 0 ​ l ​ 0 𝑙 00 0 𝑙 0 l00+0l0 . The geometry parameters are optimized for maximum acoustic rotation in water, and the resulting dependent parameters (damping factors and mode amplitudes) are listed for both water and pyridine.

Parameter	See Eq.	Unit	$l=2$	$l=4$	$l=6$
Geometry (optimized for water)
$d_{l0}^{\mathrm{opt}}$	(35a)	$\mu \text{m}$	$-3.922$	$-2.785$	$-2.286$
$A_{l0}^{\mathrm{opt}}$	(35b)	1	0.9980	0.9986	0.9989
$L^{\star }$		$\mu \text{m}$	.2000	.2000	.2000
$L_z$		$\mu \text{m}$	2.200	2.200	2.200
Dependent parameters (water)
$u_{1z}^0$	(36)	nm	1.450	7.474	14.027
$f_{l0}^{\star }$	(34)	MHz	0.747	1.494	2.243
${\Gamma }_{\mathrm{fl}}^{l00}\approx {\Gamma }_{\mathrm{fl}}^{0l0}$	(33a)	${10}^{-3}$	0.008	0.015	0.023
${\Gamma }_{\mathrm{bl}}^{l00}\approx {\Gamma }_{\mathrm{bl}}^{0l0}$	(33b)	${10}^{-3}$	3.389	2.396	1.957
Dependent parameters (pyridine)
$u_{1z}^0$	(36)	nm	1.757	9.665	19.617
$f_{l0}^{\star }$	(34)	MHz	0.707	1.415	2.123
${\Gamma }_{\mathrm{fl}}^{l00}\approx {\Gamma }_{\mathrm{fl}}^{0l0}$	(33a)	${10}^{-3}$	0.144	0.287	0.431
${\Gamma }_{\mathrm{bl}}^{l00}\approx {\Gamma }_{\mathrm{bl}}^{0l0}$	(33b)	${10}^{-3}$	3.488	2.467	2.014

Equations

$$\displaystyle\bm{u}(\bm{r},t)=\mathrm{Re}\big{[}\bm{u}_{1}(\bm{r})\mathrm{e}^{-\mathrm{i}\omega t}\big{]},\qquad\bm{r}\in\partial\Omega,$$

$$p(\mathbf{r},t)$$

$$\mathbf{v}(\mathbf{r},t)$$

$$\rho (\mathbf{r},t)$$

$$\displaystyle\nabla^{2}p_{1}+k_{\mathrm{c}}^{2}p_{1}=0,\quad k_{\mathrm{c}}=\Big{(}1+\mathrm{i}\frac{\Gamma_{\mathrm{fl}}}{2}\Big{)}k_{0},$$

$$\displaystyle\partial_{\perp}p_{1}+\frac{\mathrm{i}}{k_{\mathrm{s}}}\big{(}k_{\mathrm{c}}^{2}p_{1}+\partial_{\perp}^{2}p_{1}\big{)}=\frac{k_{0}^{2}U_{1\perp}(\bm{r})}{\kappa_{\mathrm{fl}}},$$

$$\displaystyle U_{1\perp}(\bm{r})=\frac{1}{1-\mathrm{i}\Gamma_{\mathrm{fl}}}\Big{(}-\bm{n}-\frac{\mathrm{i}}{k_{\mathrm{s}}}\bm{\nabla}\Big{)}\cdot\bm{u}_{1},$$

$$\displaystyle\Gamma_{\mathrm{fl}}=\bigg{(}\frac{4}{3}+\frac{\eta^{\mathrm{b}}_{\mathrm{fl}}}{\eta_{\mathrm{fl}}}\bigg{)}\eta_{\mathrm{fl}}\kappa_{\mathrm{fl}}\omega\ll 1,\quad k_{\mathrm{s}}=\frac{1+\mathrm{i}}{\delta},\quad\delta=\sqrt{\frac{2\eta_{\mathrm{fl}}}{\rho_{\mathrm{fl}}\omega}}.$$

$${\mathbf{v}}_1$$

$${\rho }_1$$

$$E_{\mathrm{ac}}$$

$${\mathbf{S}}_{\mathrm{ac}}$$

$${\mathcal{L}}_{\mathrm{ac}}$$

$$0$$

$${\mathbf{v}}_2$$

$${\mathbf{f}}_{\mathrm{ac}}$$

$${\mathbf{v}}_2^{\mathrm{slip}}$$

$$\displaystyle\quad+\frac{3}{2\omega}\big{\langle}[\mathrm{i}\partial_{\perp}v_{1\perp}]\bm{v}_{1\parallel}\big{\rangle},$$

$$\nabla \times {\mathbf{S}}_{\mathrm{ac}}={\omega }^2{\mathcal{L}}_{\mathrm{ac}},$$

$$\oint {\mathbf{f}}_{\mathrm{ac}}\cdot \mathrm{d}\mathbf{l}=\frac{{\Gamma }_{\mathrm{fl}}\omega }{c_{\mathrm{fl}}^2}\oint {\mathbf{S}}_{\mathrm{ac}}\cdot \mathrm{d}\mathbf{l}=\frac{{\Gamma }_{\mathrm{fl}}{\omega }^3}{c_{\mathrm{fl}}^2}\int {\mathcal{L}}_{\mathrm{ac}}\cdot \mathbf{n}\mathrm{d}A.$$

$$\displaystyle p_{1}^{lmn}\propto R^{lmn}(\bm{r})=\cos\big{(}K_{x}^{l}x\big{)}\cos\big{(}K_{y}^{m}y\big{)}\cos\big{(}K_{z}^{n}z\big{)}.$$

$$K_x^l$$

$$K_0^{lmn}$$

$$k_0^{lmn}$$

$${\Gamma }_{\mathrm{bl}}^{lmn}$$

$$\displaystyle\hskip 4.2679pt+\bigg{(}\frac{K_{y}^{m}}{K_{0}^{lmn}}\bigg{)}^{2}\bigg{(}\frac{\delta}{L_{z}^{{n}}}+\frac{\delta}{L_{x}^{{l}}}\bigg{)}$$

$$\displaystyle\hskip 4.2679pt+\bigg{(}\frac{K_{z}^{n}}{K_{0}^{lmn}}\bigg{)}^{2}\bigg{(}\frac{\delta}{L_{x}^{{l}}}+\frac{\delta}{L_{y}^{{m}}}\bigg{)},$$

$${\Gamma }^{lmn}={\Gamma }_{\mathrm{bl}}^{lmn}+{\Gamma }_{\mathrm{fl}}^{lmn}.$$

$$p_1^{lmn}(k_0,\mathbf{r})$$

$$P_1^{lmn}$$

$$F^{lmn}(k_0)$$

$${\varphi }_{\mathrm{tot}}^{lmn}={\varphi }_{\mathrm{act}}^{lmn}+{\varphi }_F^{lmn}.$$

$$\bm{f}_{\mathrm{ac}}=\Gamma_{\mathrm{fl}}\kappa_{\mathrm{fl}}\big{\langle}[\mathrm{i}p_{1}]\bm{\nabla}p_{1}\big{\rangle}.$$

$$p_1=p_{l^{\prime }{m}^{\prime }{n}^{\prime }}^{lmn}=p_1^{lmn}+p_1^{l^{\prime }{m}^{\prime }{n}^{\prime }},$$

$${\mathbf{f}}_{l^{\prime }{m}^{\prime }{n}^{\prime }}^{lmn}$$

$$\displaystyle=\Gamma_{\mathrm{fl}}\kappa_{\mathrm{fl}}\big{[}\big{\langle}[\mathrm{i}p_{1}^{lmn}]\bm{\nabla}p_{1}^{{l^{\prime}m^{\prime}n^{\prime}}}\big{\rangle}-\big{\langle}p_{1}^{{l^{\prime}m^{\prime}n^{\prime}}}[\bm{\nabla}\mathrm{i}p_{1}^{lmn}]\big{\rangle}\big{]}.$$

$${\mathbf{f}}_{l^{\prime }{m}^{\prime }{n}^{\prime }}^{lmn}$$

$${\mathcal{F}}_{l^{\prime }{m}^{\prime }{n}^{\prime }}^{lmn}(k_0)$$

$${\mathcal{K}}_{l^{\prime }{m}^{\prime }{n}^{\prime }}^{lmn}(\mathbf{r})$$

$${\mathbf{S}}_{l^{\prime }{m}^{\prime }{n}^{\prime }}^{lmn}$$

$${\mathcal{L}}_{l^{\prime }{m}^{\prime }{n}^{\prime }}^{lmn}$$

Full text

I Introduction

Acoustophoresis, particle migration by ultrasound in a microfluidic setting, is an increasingly popular method for non-contact and label-free handling of microparticles suspended in a fluid. Examples include acoustic separation Petersson et al. (2007); Manneberg et al. (2008); Ding et al. (2014), trapping Hammarström et al. (2012); Evander et al. (2015), and tweezing Collins et al. (2016); Lim et al. (2016); Baresch et al. (2016), as well as enrichment of cancer cells Augustsson et al. (2012); Antfolk et al. (2015) and bacteria Antfolk et al. (2014); Li et al. (2017), and size-independent sorting of cells Augustsson et al. (2016).

The acoustic field induces the acoustic radiation force that acts on suspended particles and scales with the particle volume King (1934); Yosioka and Kawasima (1955); Doinikov (1997); Settnes and Bruus (2012); Karlsen and Bruus (2015), as well as the acoustic streaming in the fluid Lord Rayleigh (1884); Schlichting (1932); Westervelt (1953); Nyborg (1958) that gives rise to a viscous Stokes drag force on the particles which scale with the linear size of the particles Muller et al. (2013). Consequently, the streaming-induced drag force is the dominating acoustic force for particles smaller than a critical size. For example, for an aqueous suspension of spherical polystyrene particles in a 1-MHz ultrasound field, the critical radius has been determined to be around $2~{}\upmu\textrm{m}$ Muller et al. (2012); Barnkob et al. (2012). To extend the application of acoustofluidics into the regime of sub-micrometer particles, such as bacteria, viruses and exosomes, a deep understanding of acoustic streaming is therefore important.

The first theoretical analysis of acoustic streaming dates back to Lord Rayleigh Lord Rayleigh (1884), who investigated the streaming generated by the viscous boundary layers close to the walls in planar and cylindrical systems. This analysis was later generalized by Nyborg Nyborg (1958), who formulated a slip velocity condition for the acoustic streaming outside the viscous boundary layers at weakly curved boundaries oscillating in the normal direction. Further extensions to this slip-velocity theory involves curvilinear corrections for curved boundary Lee and Wang (1989), general motion of a flat boundary Vanneste and Bühler (2011), and general motion of a curved boundary Bach and Bruus (2018). This kind of streaming is caused by the velocity mismatch between the boundary and the acoustic wave in the bulk, and we refer to it as boundary-driven acoustic streaming.

Acoustic streaming can also be generated by attenuation of the acoustic wave in the bulk as shown by Eckart Eckart (1948), who established the basic theory and investigated the steady flow caused by a sound beam. This kind of streaming is generated by a force, which is proportional to the acoustic energy flux density (or intensity) such that, say, a sound beam generates acoustic streaming in its direction of propagation. We refer to this kind of streaming as bulk-driven acoustic streaming.

While bulk-driven streaming is known to have a clear effect in systems much larger than the acoustic wavelength Riaud et al. (2017), it is often neglected in studies of resonating devices with length scales comparable with the acoustic wavelength. However, the bulk-driven streaming was taken into account implicitly by Antfolk et al. in their study of a long straight microchannel with a nearly-square cross section Antfolk et al. (2014), as they made a direct numerical simulation of the acoustofluidics in the system. They showed numerically that the experimentally observed single-vortex streaming pattern, which remarkably appeared in their microchannel instead of the usual quadrupolar Rayleigh-vortex pattern, could be explained if two orthogonal acoustic resonances with nearly identical resonance frequencies were excited with a quarter-period phase difference at the same frequency.

Such orthogonal phase-shifted acoustic waves have been used in other experiments to generate an acoustic radiation torque acting on suspended particles Busse and Wang (1981); Tran et al. (2012); Dual et al. (2012); Lee and Wang (1989); Courtney et al. (2013); Mao et al. (2017). Also travelling Bessel beams Baresch et al. (2018); Marston (2006) as well as arrays of holographic acoustic elements Marzo et al. (2015) have been used for this purpose, but not for generating specific streaming patterns.

Inspired by the work of Antfolk et al. Antfolk et al. (2014), we study in this paper bulk-driven acoustic streaming at resonance in rectangular-box cavities with acoustically hard walls and side lengths comparable to the acoustic wavelength. The paper is organized as follows: In Section II we present the governing equations for the acoustic pressure and the acoustic streaming at resonance, where the viscous boundary layers are taken into account using our recent boundary-layer analysis Bach and Bruus (2018). We apply these equations in Section III to derive analytical solutions for the acoustic single-mode pressure resonances in the rectangular cavity. We also introduce the so-called double-mode resonances, where two overlapping single modes are excited simultaneously at the same frequency. In Section IV we establish the so-called overlap and phase conditions for double-mode resonances that lead to a strong acoustic body force, which causes the bulk-driven streaming. A main result is presented in Section V, where we show how a weak symmetry breaking of a perfect square cavity can lead to a rotating body force even for a non-rotating actuation. We validate our theory in Section VI by direct numerical simulation in the case of a shallow nearly-square cavity similar to the geometry studied experimentally by Hagsäter et al. Hagsäter et al. (2007). We predict that the bulk-driven streaming in closed microcavities is enhanced by increasing the frequency and the bulk viscosity, just as is the case for Eckart streaming in open systems Eckart (1948). We propose to test this prediction by replacing water with pyridine in the microcavity. Finally, we discuss our results in Section VII and present our conclusions in Section VIII.

II Governing equations

The time-harmonic acoustic fields in the fluid domain $\Omega$ are induced by the time-harmonic displacement field $\bm{u}(\bm{r},t)$ of the wall at the boundary $\partial\Omega$ of the domain,

[TABLE]

where $\bm{r}$ is position, $t$ is time, $\omega=2\pi f$ is the angular driving frequency, and “Re” is the real part of complex-valued fields. In the quiescent and homogeneous fluid of density $\rho_{\mathrm{fl}}$ and ambient pressure $p_{\mathrm{fl}}$, we apply standard perturbation theory to describe the fluid velocity $\bm{v}(\bm{r},t)$, pressure $p(\bm{r},t)$, and density $\rho(\bm{r},t)$,

[TABLE]

where the subscript “1” denotes the first-order, complex-valued, time-harmonic acoustic fields and “2” the second-order, real-valued, time-averaged steady fields. We do not compute the oscillating second-order fields containing $\mathrm{e}^{\pm\mathrm{i}2\omega t}$, as they are typically not observed.

II.1 Pressure acoustics with boundary layers

The complex-valued acoustic pressure $p_{1}(\bm{r})$ satisfies a Helmholtz equation with a complex-valued compressional wavenumber $k_{\mathrm{c}}$ having the real part $k_{0}=\frac{\omega}{c_{\mathrm{fl}}}$. We apply the effective boundary condition for $p_{1}(\bm{r})$ recently derived by Bach and Bruus Bach and Bruus (2018), which takes the viscous boundary layer into account at the domain boundary $\partial\Omega$. We introduce the inward normal derivative $\partial_{\perp}=-\bm{n}\cdot\bm{\nabla}$ and the effective inward normal displacement $U_{1\perp}(\bm{r})$ of the wall in terms of the actual displacement $\bm{u}_{1}(\bm{r})$,

[TABLE]

Here, $\Gamma_{\mathrm{fl}}$ is the weak bulk damping coefficient and $k_{\mathrm{s}}$ is the viscous boundary layer wavenumber, given in terms of the viscous boundary-layer width $\delta$, the dynamic viscosity $\eta_{\mathrm{fl}}$, the bulk viscosity $\eta^{\mathrm{b}}_{\mathrm{fl}}$, and the isentropic compressibility $\kappa_{\mathrm{fl}}=\frac{1}{\rho_{\mathrm{fl}}}(\frac{\partial\rho}{\partial p})_{S}=\frac{1}{\rho_{\mathrm{fl}}c_{\mathrm{fl}}^{2}}$ of the fluid,

[TABLE]

From the acoustic pressure $p_{1}$ we directly obtain the irrotational acoustic velocity $\bm{v}_{1}$ and the acoustic density $\rho_{1}$ in Eq. (2) through the relations Bach and Bruus (2018),

[TABLE]

Further, we define the following relevant time-averaged first-order products: the acoustic energy density $E_{\mathrm{ac}}$, kinetic energy density $E_{\mathrm{kin}}$, potential energy density $E_{\mathrm{pot}}$, energy flux density vector $\bm{S}_{\mathrm{ac}}$, as well as the time-averaged acoustic angular momentum density $\bm{\mathcal{L}}_{\mathrm{ac}}$ with respect to unperturbed fluid position $\bm{r}$,

[TABLE]

Here, $\bm{d}_{1}(\bm{r})$ is the fluid displacement at position $\bm{r}$, and \big{\langle}A_{1}(\bm{r},t)B_{1}(\bm{r},t)\big{\rangle}=\frac{1}{2}\mathrm{Re}[A_{1}(\bm{r})B_{1}^{*}(\bm{r})] is the usual time-average of products of time-harmonic fields $A_{1}$ and $B_{1}$, with the asterisk denoting complex conjugation.

II.2 Acoustic streaming at fluid resonance

Outside the narrow viscous boundary layers of width $\delta\lesssim 500\,\textrm{nm}$ (for water at MHz frequencies), the second-order acoustic streaming $\bm{v}_{2}$ is described as a Stokes flow driven by the acoustic body force $\bm{f}_{\mathrm{ac}}$ in the domain $\Omega$, see Eq. (8) in Ref. Karlsen et al. (2018) and Eq. (52b) in Ref. Bach and Bruus (2018), and by the slip velocity $\bm{v}_{2}^{\mathrm{slip}}$ at the domain boundary $\partial\Omega$, see Eq. (61) in Ref. Bach and Bruus (2018),

[TABLE]

where we have used the expression for the slip velocity $\bm{v}_{2}^{\mathrm{slip}}$ valid near a fluid resonance, in which case the magnitude $|\bm{v}_{1}|$ of the acoustic velocity $\bm{v}_{1}$ in the bulk is much larger than that of the wall velocity \big{|}\omega\bm{u}_{1}\big{|}\ll|\bm{v}_{1}|, see Ref. Bach and Bruus (2018). We call the streaming driven by the acoustic body force $\bm{f}_{\mathrm{ac}}$ in Eq. (7d) for bulk-driven acoustic streaming, and the streaming driven by the slip velocity $\bm{v}_{2}^{\mathrm{slip}}$ in Eq. (7) for boundary-driven acoustic streaming. Note that $\bm{f}_{\mathrm{ac}}$ increases with the frequency squared, as $\Gamma_{\mathrm{fl}}\propto\omega$, whereas $\bm{v}_{2}^{\mathrm{slip}}$ is independent of the frequency, as $\bm{\nabla}_{\parallel},\partial_{\perp}\sim k_{0}=\frac{1}{c_{\mathrm{fl}}}\>\omega$.

One central quantity in the governing equations (7) for the acoustic streaming, is the acoustic energy flux density $\bm{S}_{\mathrm{ac}}$, which enters both in the body force $\bm{f}_{\mathrm{ac}}$, Eq. (7d), and in the slip velocity $\bm{v}_{2}^{\mathrm{slip}}$, Eq. (7). This quantity is rotating if the acoustic fields are rotating, as seen by taking the curl of \bm{S}_{\mathrm{ac}}=\big{\langle}p_{1}\bm{v}_{1}\big{\rangle} and using Eq. (5a) with $\bm{\nabla}p_{1}\approx\mathrm{i}\omega\rho_{\mathrm{fl}}\bm{v}_{1}$ and $\bm{\nabla}\times\bm{v}_{1}=\bm{0}$,

[TABLE]

Hence, bulk-driven streaming is a consequence of rotating acoustics with angular momentum density $\bm{\mathcal{L}}_{\mathrm{ac}}$.

III Acoustic resonance modes in a box with viscous boundary layers

We consider the rectangular-box cavity of dimensions $L_{x}\times L_{y}\times L_{z},$ where $0\leq x\leq L_{x}$, $0\leq y\leq L_{y}$, and $0\leq z\leq L_{z}$. Below, we give the solutions to Eq. (3) in this geometry for the resonance modes of the acoustic pressure with viscous boundary layers. We refer the reader to Appendix A for the detailed derivations. We then evaluate the second-order quantities $\bm{S}_{\mathrm{ac}}$, $\bm{f}_{\mathrm{ac}}$, and $\bm{\mathcal{L}}_{\mathrm{ac}}$ which vanish for all single-mode resonances but not for double-mode resonances, where two single modes are excited simultaneously by the same frequency.

III.1 Single-mode resonances

First, ignoring viscosity, the resonance solution for $p_{1}$ is proportional to the eigenfunctions $R^{lmn}(\bm{r})$ of the Helmholtz equation (3a), which are the usual hard-wall standing resonance modes with integer number $l$, $m$, and $n$ half waves in the $x$, $y$, and $z$ direction, respectively,

[TABLE]

In the presence of viscosity and hence viscous boundary layers, the resonant wavenumbers $k_{0}^{lmn}$ and angular frequencies $\omega^{lmn}$ are down-shifted slightly relative to the inviscid values due to the boundary-layer damping coefficients $\Gamma_{\mathrm{bl}}^{lmn}$,

[TABLE]

where $L_{x}^{{l}}=\int_{0}^{L_{x}}\cos^{2}(K_{x}^{l}x)\,\mathrm{d}x=\frac{1}{2}(1+\delta_{0l})L_{x}$ and similarly for $L_{y}^{{m}}$ and $L_{z}^{{n}}$. As an example, a pure $x$-mode will have $K_{0}^{l00}=K_{x}^{l}$ and $K_{y}^{0}=K_{z}^{0}=0$ giving $\Gamma_{\mathrm{bl}}^{l00}=\frac{\delta}{L_{y}}+\frac{\delta}{L_{z}}$ corresponding to the boundary-layer damping from the four boundaries parallel to the $x$-axis. The total damping coefficient $\Gamma^{lmn}$ of the resonance mode $lmn$ includes both the boundary-layer-damping coefficients $\Gamma_{\mathrm{bl}}^{lmn}$ and the bulk damping coefficient $\Gamma_{\mathrm{fl}}^{lmn}$, given in Eq. (4) with $\omega=\omega^{lmn}$,

[TABLE]

For frequencies close to the resonance frequencies, we have $k_{0}\approx k_{0}^{lmn}$ or $\omega\approx\omega^{lmn}$, and the resonance pressure mode $p^{lmn}_{1}$ satisfying Eq. (3) is the product of a complex-valued amplitude $P_{1}^{lmn}$, an internal frequency dependency $F^{lmn}(k_{0})$, and the usual hard-wall spatial dependency $R^{lmn}(\bm{r})$ from Eq. (9a),

[TABLE]

Here, we have also given $P_{1}^{lmn}$ and $F^{lmn}$ in their polar form and introduced the external actuation phase $\phi_{\mathrm{act}}^{lmn}$ as well as the internal frequency-dependent phase $\phi_{F}^{lmn}$. The phase $\phi_{F}^{lmn}$ is plotted in Fig. 1 together with the modulus $|F^{lmn}|$. The total phase $\phi_{\mathrm{tot}}^{lmn}$ of the mode $lmn$ is the sum of these phases,

[TABLE]

III.2 Double-mode resonances

The bulk-driven acoustic streaming is driven by the body force $\bm{f}_{\mathrm{ac}}$, defined in Eq. (7d) and restated here solely in terms of the pressure $p_{1}$, by using $\bm{v}_{1}\approx\frac{-\mathrm{i}}{\omega\rho_{\mathrm{fl}}}\bm{\nabla}p_{1}$ as stated in Eq. (5a),

[TABLE]

For any single-mode resonance $lmn$ of the form of Eq. (12a), this expression will vanish since $\mathrm{i}p_{1}^{lmn}$ is exactly $\frac{\pi}{2}$ out of phase with $\bm{\nabla}p_{1}^{lmn}$. This motivates us to consider double-mode resonances, which occur if two single modes $lmn$ and ${l^{\prime}m^{\prime}n^{\prime}}$ overlap, such that they can be excited simultaneously at the same frequency $f=\frac{1}{2\pi}k_{0}c_{\mathrm{fl}}$, with $k_{0}^{lmn}\approx k_{0}\approx k_{0}^{{l^{\prime}m^{\prime}n^{\prime}}}$. In that case, the total pressure is the sum of two modes,

[TABLE]

and the body force $\bm{f}_{\mathrm{ac}}$ in Eq. (14) becomes $\bm{f}_{{l^{\prime}m^{\prime}n^{\prime}}}^{lmn}$,

[TABLE]

Inserting the resonance (12a) with $P_{1}^{lmn}$ in polar form yields,

[TABLE]

Here, $\mathcal{F}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}\,(k_{0})$ gives the essential frequency dependency (as $\Gamma_{\mathrm{fl}}$ is nearly constant over the width of the double mode), while $\bm{\mathcal{K}}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}\,(\bm{r})$ with the dimension of a wave vector gives the spatial dependency and direction of the double-mode body force. Furthermore, we obtain the double-mode energy flux density $\bm{S}_{lmn}^{l^{\prime}m^{\prime}n^{\prime}}$ and acoustic angular momentum density $\bm{\mathcal{L}}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}$ directly from Eq. (17a) by use of Eqs. (7d) and (8a),

[TABLE]

Note that $\mathcal{F}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}(k_{0})$ is normalized to be between $-1$ and $1$, whereas $\bm{\mathcal{K}}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}(\bm{r})$ and $\bm{\nabla}\times\bm{\mathcal{K}}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}(\bm{r})=2\bm{\nabla}R^{lmn}\times\bm{\nabla}R^{{l^{\prime}m^{\prime}n^{\prime}}}$ has the maximum amplitudes $k_{0}$ and $2k_{0}^{2}$, respectively.

In Fig. 2 we show results for the nine lowest-frequency double-modes $lm0+ml0$ in a square cavity. We plot the body force $\bm{f}^{lm0}_{ml0}\propto\bm{\mathcal{K}}^{lm0}_{ml0}$ and the $z$-component of the acoustic angular momentum $\bm{\mathcal{L}}^{lm0}_{ml0}\propto\bm{\nabla}\times\bm{\mathcal{K}}^{lm0}_{ml0}$. In the Supplemental Material 111See Supplemental Material at [URL] extending the results of the following three figures: (1) Fig. 2 turned into one pdf-file with 55 sub-figures of the acoustic body force $\bm{f}_{\mathrm{ac}}$ of the double-modes $lm0+l^{\prime}m^{\prime}0$ with $l,m,l^{\prime},m^{\prime}=0,1,2,3,4$ (the MATLAB source code is also provided); (2) Fig. 4 turned into five animated gif files of the pressure $p_{1}$ in 3D and as line plots of the double modes $200+020$, $400+040$, and $600+060$ both for the wide $\frac{1}{2}L^{\star}$-width and the narrow $\frac{1}{12}L^{\star}$-width $G_{1}$ actuation; and (3) Fig. 8 re-calculated using the narrow $\frac{1}{12}L^{\star}$-width $G_{1}$ actuation., we show $55$ examples of double-mode resonances $lm0+l^{\prime}m^{\prime}0$ with $l,m,l^{\prime},m^{\prime}=0,1,2,3,4$ in rectangular cavities with aspect ratios between $1$ and $4$, and we provide a Matlab code to compute these modes.

IV Conditions for obtaining a strong acoustic body force

The magnitude of the body force $\bm{f}_{{l^{\prime}m^{\prime}n^{\prime}}}^{lmn}$, the acoustic energy flux density $\bm{S}_{{l^{\prime}m^{\prime}n^{\prime}}}^{lmn}$, and the acoustic angular momentum density $\bm{\mathcal{L}}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}$ are all proportional to the frequency-dependent factor $\mathcal{F}^{lmn}_{l^{\prime}m^{\prime}n^{\prime}}(k_{0})$ in Eq. (17b). To study this factor, we rewrite it to emphasize the total phase difference $\phi_{\mathrm{tot}}^{{l^{\prime}m^{\prime}n^{\prime}}}-\phi_{\mathrm{tot}}^{lmn}$ between the single mode phases given in Eq. (13),

[TABLE]

Consequently, two conditions must be fulfilled to yield a strong body force: (i) the overlap condition, i.e. the two single-mode resonances must have a large overlap in frequency space to ensure a large value of the amplitude \big{|}F^{lmn}\big{|}\big{|}F^{{l^{\prime}m^{\prime}n^{\prime}}}\big{|}, and (ii) the phase condition, i.e. the difference $\phi_{\mathrm{tot}}^{{l^{\prime}m^{\prime}n^{\prime}}}-\phi_{\mathrm{tot}}^{lmn}$ in phases must be close to $\pm\frac{\pi}{2}$ to ensure a large value of the sine factor.

IV.1 The overlap condition: Aspect ratios

The condition that both $|F^{lmn}|$ and $|F^{{l^{\prime}m^{\prime}n^{\prime}}}|$ are large in Eq. (19) is satisfied if the modes overlap. This occurs if the resonance frequencies are nearly identical, $k_{0}^{lmn}\approx k_{0}^{{l^{\prime}m^{\prime}n^{\prime}}}$, or equivalently,

[TABLE]

where we have used that $\Gamma_{\mathrm{bl}}^{lmn}\ll 1$ in Eq. (10a) so $k_{0}^{lmn}\approx K_{0}^{lmn}$. In the general case, Eq. (20) gives a myriad of possible mode combinations. If we restrict ourselves to horizontal modes, where $n=n^{\prime}=0$, we can solve for the aspect ratio $A=\frac{L_{x}}{L_{y}}$ in Eq. (20),

[TABLE]

Listing the integers $l,m,l^{\prime},m^{\prime}$ for which $A$ is real, we find that the square cavity with $A\approx 1$ is the richest case having the highest number of allowed double-mode excitations. In Fig. 2, we show the first 9 horizontal double-mode excitations in a square cavity with $A=1$ and $n=n^{\prime}=0$. In this case, the condition in Eq. (20) corresponds to two points \big{(}\frac{l\pi}{L_{x}},\frac{m\pi}{L_{y}}\big{)} and \big{(}\frac{l^{\prime}\pi}{L_{x}},\frac{m^{\prime}\pi}{L_{y}}\big{)} both lying on the same quarter circle of radius $K_{0}$, as illustrated by the black dots in Fig. 2(a).

IV.2 The phase condition: Rotating versus non-rotating actuation

The straightforward way to obtain rotating acoustics is to actuate a given system with two transducers running with a phase difference of $\phi_{\mathrm{act}}^{{l^{\prime}m^{\prime}n^{\prime}}}-\phi_{\mathrm{act}}^{lmn}=\pm\frac{\pi}{2}$, a technique that has been used in many experiments to generate acoustic radiation torques Busse and Wang (1981); Tran et al. (2012); Dual et al. (2012); Lee and Wang (1989); Courtney et al. (2013); Mao et al. (2017). Remarkably, even in the standard case of microfluidic ultrasound experiments driven by a single piezoelectric transducer with a non-rotating actuation, $\phi_{\mathrm{act}}^{{l^{\prime}m^{\prime}n^{\prime}}}-\phi_{\mathrm{act}}^{lmn}=0$, a rotating acoustic field may nevertheless be generated because the total phase difference $\phi_{\mathrm{tot}}^{{l^{\prime}m^{\prime}n^{\prime}}}-\phi_{\mathrm{tot}}^{lmn}=\phi_{F}^{{l^{\prime}m^{\prime}n^{\prime}}}-\phi_{F}^{lmn}$ may be close to $\pm\frac{\pi}{2}$ due to the internal phases $\phi_{F}^{lmn}$ of the frequency-dependent factors $F^{lmn}(k_{0})$, see Eq. (12c) and Fig. 1.

V Non-rotating actuation of rotating double-mode resonances in a nearly-square cavity

We now show how a non-rotating actuation can lead to a rotating double-mode resonance in a nearly-square cavity characterized by the average side length $L^{\star}=\frac{1}{2}(L_{x}+L_{y})$ and the small symmetry breaking $d=L_{x}-L_{y}$,

[TABLE]

The magnitude of the resulting acoustic body force $\bm{f}_{\mathrm{ac}}$ depends strongly on the symmetry breaking $d$, and for a given double mode we determine how much the aspect ratio $A=L_{x}/L_{y}\approx 1+d/L^{\star}$ should deviate from unity to maximize $\bm{f}_{\mathrm{ac}}$. Specifically, we consider the symmetric double-mode excitations $lm0+ml0$ in the horizontal $xy$ plane of the type shown in Fig. 2. For such double-modes, the analysis simplifies significantly, not only because $z$ drops out, but also because the damping coefficients $\Gamma^{lm0}$ and $\Gamma^{ml0}$ from Eqs. (10) and (11) are almost identical. In short, we consider the situation,

[TABLE]

To analyse the strength of the acoustic body force, we examine the factor $\mathcal{F}_{lm0}^{ml0}(k_{0})$ in Eq. (17b) that describes the frequency-dependency of the acoustic body force. For the double mode $lm0+ml0$, we introduce the short two-subscript notation $\mathcal{F}_{lm}(k_{0})$ for this quantity,

[TABLE]

where the functions $F^{lm0}(k_{0})$ and $F^{ml0}(k_{0})$ behave as shown in Fig. 1. Note that if the modes $lm0$ and $ml0$ have equal resonance frequencies, they have full mode overlap, and the product $|F^{lm0}||F^{ml0}|$ is maximized, while \sin\big{(}\phi_{F}^{ml0}-\phi_{F}^{lm0}\big{)} vanishes, and there is no body force. To maximize $\mathcal{F}_{lm}(k_{0})$, the two modes must therefore be separated sufficiently to ensure that the phase difference $\phi_{F}^{ml0}-\phi_{F}^{lm0}$ deviates from zero, but not too much, as a decreased mode overlap reduces the modulus product $|F^{lm0}||F^{ml0}|$. Therefore, we cannot in this case perfectly fulfil both the overlap condition, Section IV.1, and the phase condition, Section IV.2, so the strongest body force is found at a certain mode separation $\Delta_{lm}=k_{0}^{lm0}-k_{0}^{ml0}$ between the resonance wavenumbers, which gives only partial mode overlap and a phase difference that deviates from the ideal $\pm\frac{\pi}{2}$.

The mode separation $\Delta_{lm}$ is induced by the symmetry breaking $d$, and similar to Eq. (22), we introduce the center wavenumber $k^{\star}_{lm}=\frac{1}{2}(k_{0}^{lm0}+k_{0}^{ml0})$ corresponding to the center frequency between the two modes and write,

[TABLE]

The relation between the mode separation $\Delta_{lm}$, the symmetry breaking $d$, and the aspect ratio $A=L_{x}/L_{y}$, is found by inserting Eq. (22) into $k_{0}^{lm0}=\sqrt{(\frac{l\pi}{L_{x}})^{2}+(\frac{m\pi}{L_{y}})^{2}}+\mathcal{O}({\Gamma^{lm0}})$, expanding in the small ratio $\frac{d}{L^{\star}}$, and comparing with Eq. (25a),

[TABLE]

Since the damping coefficients are nearly identical $\Gamma^{lm0}\approx\Gamma_{lm}^{\star}\approx\Gamma^{ml0}$, with $\Gamma_{lm}^{\star}=\frac{1}{2}(\Gamma^{lm0}+\Gamma^{ml0})$ being the average value, it is convenient to scale all wavenumber quantities by their common line width $\Gamma^{\star}_{lm}k_{lm}^{\star}$, see Table 1,

[TABLE]

In terms of these quantities, we write the frequency-dependencies $F^{lm0}(k_{0})$ and $F^{ml0}(k_{0})$ from Eq. (12c) as,

[TABLE]

In Fig. 3(a) and (b) we show as in Fig. 1 the modulus and phase, respectively, of $F^{lm0}$ and $F^{ml0}$ for four different mode separations $\tilde{\Delta}_{lm}=0.10,0.58,1.50,3.00$. In Fig. 3(c) we show for the same values of $\tilde{\Delta}_{lm}$ the frequency dependency $\mathcal{F}_{lm}(k_{0})$ from Eq. (V) of the acoustic body force. We note that there exist an optimal mode separation $\tilde{\Delta}_{lm}^{\mathrm{opt}}$ that yields the largest possible value $\mathcal{F}_{lm}^{\mathrm{opt}}$ of $\mathcal{F}_{lm}$ found at the optimal wave number $k_{0}^{\mathrm{opt}}$. In Eq. (84) of Appendix B we show that

[TABLE]

For all the double modes shown in Fig. 2, we have $l>m$, so for a non-rotating actuation, the rotation direction of the body force will therefore be as shown in the case of $d<0$ (corresponding to $L_{x}<L_{y}$ and $A<1$), and opposite for $d>0$. Note from Eq. (29e) that this asymmetry-induced rotating acoustics is extremely sensitive to $d$. In practise, it may be difficult to control this mechanism, as it requires a fabrication accuracy much smaller than $d/L^{\star}\sim\Gamma_{lm}^{\star}$, which in microfluidic devices typically is smaller than 1% Hahn and Dual (2015).

VI 3D simulation of acoustic streaming in a nearly-square cavity

In this section we validate by direct numerical simulation the theory for rotating acoustics induced by a non-rotating actuation in a weakly symmetry-broken geometry. We simulate the acoustic pressure $p_{1}$ and the acoustic streaming $\bm{v}_{2}$ in a shallow, nearly-square cavity in three dimensions (3D). We choose a cavity similar to that investigated experimentally by Hagsäter et al. Hagsäter et al. (2007), and simulated in 3D using the hard-wall approximation by Lei et al. Lei et al. (2014), and including the actuator and elastic walls by Skov et al. Skov et al. (2019). The cavity has the mean side length $L^{\star}=2000\,\upmu\textrm{m}$, $L_{x}=L^{\star}+\frac{1}{2}d$, $L_{y}=L^{\star}-\frac{1}{2}d$, and height $L_{z}=200\,\upmu\textrm{m}$. We use Eq. (29) to optimize the symmetry breaking $d$ and actuation frequency $f=\frac{1}{2\pi}k_{0}c_{\mathrm{fl}}$ to obtain the strongest possible rotating acoustic streaming. The coordinate system is the same as in Section III, with $0\leq x\leq L_{x}$, $0\leq y\leq L_{y}$, and $0\leq z\leq L_{z}$. The boundary is taken to be stationary everywhere, except at the bottom boundary $z=0$, where we apply the non-rotating Gaussian wall displacement in the normal direction $z$, with maximum amplitude $u_{1z}^{0}$ in the cavity center and minimum amplitude $\mathrm{e}^{-2}u_{1z}^{0}$ in the corners,

[TABLE]

In Section VI.3.1–VI.3.4 the fluid is water, and in Section VI.3.6 it is pyridine. The material parameters for these fluids are listed in Table 2.

In the following simulations, we excite the double-mode resonance $l00$+$0l0$ with $l=2,4,$ and $6$. Using our analytical expressions (17), we have in Fig. 2 shown the spatial dependency of the acoustic body force $\bm{f}_{\mathrm{ac}}=\bm{f}^{l00}_{0l0}$ for these particular double modes. As seen in Eq. (7), this body force is the source of the bulk-driven acoustic streaming, so we predict that the numerical simulation of the bulk-driven streaming will result in a 2-by-2 pattern as in Fig. 2(c) for $l=2$, a 4-by-4 pattern as in Fig. 2(h) for $l=4$, and similarly for $l=6$ (not shown in Fig. 2).

VI.1 Choosing the parameters

We derived in Section III analytical expressions for the acoustic resonance modes including the viscous boundary layers, and we derived in Section V the conditions on the geometry of a nearly-square cavity for which the body force $\bm{f}_{\mathrm{ac}}$ and the acoustic rotation are maximized for a non-rotating actuation. Remarkably, these analytic results allow an analytical determination of the parameter values that lead to the strongest rotating acoustic streaming in the $xy$-plane.

For the double mode $l00$+$0l0$, the small damping coefficients $\Gamma_{\mathrm{fl}}^{l00}\approx\Gamma_{\mathrm{fl}}^{0l0}$ in Eq. (4) and $\Gamma_{\mathrm{bl}}^{l00}\approx\Gamma_{\mathrm{bl}}^{0l0}$ in Eq. (10) are computed by approximating the angular frequency by the constant value $\omega=c_{\mathrm{fl}}k_{l0}^{\star}\approx c_{\mathrm{fl}}\frac{l\pi}{L^{\star}}$, see Eq. (25b),

[TABLE]

Combining Eqs. (10a), (29a), and (25b), we obtain an expression for the boundary-layer shifted, optimal double-mode actuation frequency, which is the center frequency of the two single modes,

[TABLE]

The optimal symmetry breaking $d_{l0}^{\mathrm{opt}}$ and the corresponding aspect ratio $A_{l0}^{\mathrm{opt}}$ then follow from Eqs. (29e) and (29f),

[TABLE]

Finally, we set the pressure amplitudes $P_{1}^{l00}=P_{1}^{0l0}=1~{}\textrm{MPa}$ in Eq. (12b), and compute the corresponding amplitude $u^{0}_{1z}$ of the actuation displacement in Eq. (32) as,

[TABLE]

For water and $l=2$, 4, and 6, the values of the expressions in Eq. (33)-(36) are listed in Table 3. For pyridine, we adjust the actuation amplitude and actuation frequency to maintain the resonance amplitudes of $P_{1}^{l00}=P_{1}^{0l0}=1~{}\textrm{MPa}$. Note that for pyridine, we use the same optimal symmetry breaking $d_{l0}^{\mathrm{opt}}$ as used for water. By combining Eqs. (V), (28) and (29e) with the material parameters in Table 2, this leads to the values $\mathcal{F}_{l0}=0.65$, 0.64, and 0.63 for $l=2$, 4, and 6, respectively, which are nearly the same as the optimal value $\mathcal{F}_{l0}^{\mathrm{opt}}=0.65$ given in Eq. (29c).

VI.2 Implementation in COMSOL

As in our previous work Muller et al. (2012); Ley and Bruus (2016); Skov et al. (2019), the simulation is performed in two steps using the Weak Form PDE module in the finite-element-method software COMSOL Multiphysics COM . In the first step, we solve for the acoustic pressure $p_{1}$ by implementing Eq. (3) in weak form, and in the second step we solve Eq. (7) for the steady streaming $\bm{v}_{2}$ and pressure $p_{2}$. The zero-level of $p_{2}$ is fixed by using the Global Constraint $\int_{\Omega}p_{2}\,\mathrm{d}V=0$. For numerical efficiency, we exploit the symmetry at the vertical planes $x=\frac{1}{2}L_{x}$ and $y=\frac{1}{2}L_{y}$, where we apply the symmetry conditions $\bm{n}\cdot\bm{\nabla}p_{1}=0$, $\bm{n}\cdot\bm{v}_{2}=0$, and $\bm{n}\cdot\bm{\nabla}\bm{v}_{2\parallel}=\bm{0}$. The computational domain is thus reduced to the quadrant $0\leq x\leq\frac{1}{2}L_{x}$ and $0\leq y\leq\frac{1}{2}L_{y}$. We use quartic-, cubic-, and quadratic-order Lagrangian shape functions for $p_{1}$, $\bm{v}_{2}$, and $p_{2}$, respectively, and with a tetrahedral finite-element mesh of mesh size $\frac{1}{6}L_{z}$, this leads to $1.3\times 10^{6}$ degrees of freedom and a relative accuracy better than 1 %. The simulations were performed on a workstation with a 3.5-GHz Intel Xeon CPU E5-1650 v2 dual-core processor, and with a memory of 128 GB RAM.

VI.3 Simulation results

VI.3.1 The acoustic pressure

In Fig. 4, we show the simulation results $p_{1}^{\mathrm{num}}$ of the acoustic pressure $p_{1}$ in water for $l=6$ in 3D and lineplots for $l=2,4$, and $6$. These three modes are examples of rotating double modes actuated by the wide, single, non-rotating Gaussian actuation $G_{1}$, Eqs. (32) and (36), at the bottom boundary $z=0$. In the Supplemental Material Note (1) we show the rotating dynamics of the pressure contours similar to Fig. 4(a) for $l=2,4,$ and $6$ as gif animations. In Fig. 4 we also show the the analytical result $p_{1}^{\mathrm{ana}}=p_{0l0}^{l00}=p_{1}^{l00}+p_{1}^{0l0}$, which is obtained from Eqs. (12) and (15) with the optimal mode separation $k_{0}^{l00}-k_{0}^{0l0}=\frac{1}{\sqrt{3}}k_{l0}^{\star}\Gamma_{l0}^{\star}$, Eq. (29b), and mode amplitudes $P_{1}^{l00}=P_{1}^{0l0}=1\,\textrm{MPa}$,

[TABLE]

In Fig. 4(b), we compare \operatorname{Im}\big{(}p_{1}^{\mathrm{num}}\big{)} and \operatorname{Im}\big{(}p_{1}^{\mathrm{ana}}\big{)} and find a quantitative agreement better than 1 %. The reason for choosing the imaginary part is to single out the resonant part of the field, which is phase shifted by the factor \exp\big{(}\mathrm{i}\frac{\pi}{2}\big{)}=\mathrm{i} relative to the actuation displacement $\bm{u}_{1}$ at the boundary, Eq. (32). All non-resonant contributions are in phase with $\bm{u}_{1}$. In the Supplemental Material Note (1) is shown an animated gif-file “p1_246_wide_G1.gif” of the full cycle in the time-domain of the line plots in Fig. 4(b) of $p_{1}$. These animations reveal that the deviation of the analytical solution from the numerical one is $\sim 1$ % for $l=2$ at 0.75 MHz, $\sim 15$ % for $l=4$ at 1.5 MHz, and $\sim 30$ % for $l=6$ at 2.2 MHz. The deviation increases for increasing mode number, because the wide vertical $G_{1}$-actuation couples strongly to the vertical 001-mode through the amplitude factor $P^{001}_{1}$ of Eq. (12b), while the frequency-dependency factor $F^{001}(k_{0})$ of Eq. (12c), which is very small away from resonance, increases significantly as the actuation frequency $f=c_{\mathrm{fl}}k_{0}$ approaches the vertical-mode resonance $f^{001}_{0}=3.7$ MHz. We can reduce the coupling $P^{001}_{1}$ to the vertical 001-mode by reducing the width of the $G_{1}$ actuation profile. This is demonstrated in the animated gif-file “p1_246_narrow_G1.gif” of the Supplemental Material Note (1), where the width of the $G_{1}$-actuation has been reduced by a factor of 6 from $\frac{1}{2}L^{\star}$ to $\frac{1}{12}L^{\star}$. This results in a decreased deviation, now less than $\sim 1$ %, of the analytical result from the numerical one for $p_{1}$ for all three modes $l=2,4,6$.

VI.3.2 The acoustic body force and angular momentum density

In Fig. 5, we validate our analysis in Section V of rotating double modes induced by a non-rotating actuation. For the double mode $200+020$, we plot the body force $\bm{f}_{\mathrm{ac}}=\bm{f}_{020}^{200}$ from Eq. (17), which according to Eq. (7) is the source of the bulk-driven acoustic streaming, and the acoustic angular momentum density $\bm{\mathcal{L}}_{\mathrm{ac}}=\bm{\mathcal{L}}_{020}^{200}$ from Eq. (18b) in the center plane $z=\frac{1}{2}L_{z}$. We show results for three different cavity aspect ratios $A=L_{x}/L_{y}$ for both the non-rotating wide Gaussian actuation $G_{1}(x,y)$, given in Eq. (32), and for the externally-controlled rotating actuation $G_{16}(x,y)$ consisting of four quadruples of narrow Gaussian actuations, see Eq. (89) of Appendix C. As expected, the non-rotating actuation $G_{1}$ used in Fig. 5(a)-(c) induces no acoustic rotation for a perfect square cavity $A=1$, Fig. 5(b), while it induces different rotation directions for $A$ slightly larger and smaller than unity, see respectively Fig. 5(a) and (c). In contrast, the rotating actuation $G_{16}$ used in Fig. 5(d)-(f) for the same three aspect ratios, induces an acoustic rotation, which is maximized in the perfect square cavity $A=1$, see Fig. 5(e). We show in Eq. (92) of Appendix C that this rotating actuation excites the modes $200$ and $020$ with a phase difference $\phi_{\mathrm{act}}^{020}-\phi_{\mathrm{act}}^{200}\approx 0.5\pi$ such that the factor $\mathcal{F}_{020}^{200}$ in Eq. (19) is positive, and consequently the rotation pattern in Fig. 5(d)-(f) is in the same direction as shown in Fig. 2(c).

In Fig. 6 we investigate the strength and rotation direction of the body force $\bm{f}_{020}^{200}$ as a function of the mode separation $\tilde{\Delta}_{20}$ or equivalently, the aspect ratio $A$. From Eq. (17a), this dependency is quantified through the factor $\mathcal{F}_{020}^{200}$ which we obtain numerically from the simulation of the double mode $200+020$ in Fig. 5 by calculating the acoustic angular momentum $\bm{\mathcal{L}}_{\mathrm{ac}}$ and using Eq. (18b),

[TABLE]

Here, for pressure amplitudes \big{|}{P_{1}^{200}}\big{|}=\big{|}{P_{1}^{020}}\big{|}=1 MPa, the denominator evaluates to $0.954\times 10^{-4}\,~{}\textrm{kg}\,~{}\textrm{m}^{-1}\textrm{s}^{-1}$ , which is indeed the maximum value found numerically in Fig. 5(e), such that $(\mathcal{F}_{020}^{200})_{\mathrm{num}}$ takes the maximum value of unity, as expected. In Fig. 6 we plot the largest possible value of $(\mathcal{F}_{020}^{200})_{\mathrm{num}}$ both for the non-rotating actuation used in Fig. 5(a)-(c) and for the rotating actuation used in Fig. 5(d)-(f). We also plot the analytical expressions for $(\mathcal{F}_{ml0}^{lm0})_{\mathrm{ana}}^{\mathrm{extr}}$ derived in Appendix B and find a quantitative agreement better than 1 %. For the non-rotating actuation (the dark blue curve in Fig. 6), the acoustic rotation direction is determined by the aspect ratio $A$, and it reverses as $A$ crosses over from $A<1$ to $A>1$ around the perfect square $A=1$.

For the rotating actuation (the light blue curves in Fig. 6), the rotation direction can be either co- and counter-rotating (at different frequencies) with respect to the rotation direction of the actuation. However, for small mode separations $|\tilde{\Delta}_{lm}|<1$, which for the double mode $200+020$ corresponds to $|A-1|<0.34\%$, the rotation direction of the acoustic field can only be co-rotating. For further details, see Appendix B.

VI.3.3 The acoustic streaming

In Fig. 7, we show the acoustic streaming $\bm{v}_{2}$ in 3D resulting from the double mode $600$+$060$ displayed in Fig. 4(a), which is excited by the non-rotating actuation $G_{1}$. As mentioned above, this situation is similar to that investigated experimentally by Hägsater et al. Hagsäter et al. (2007). Remarkably, in the optimized geometry of Fig. 7, the horizontal 6-by-6 streaming-roll pattern at any height $0\leq z\leq L_{z}$ is strongest in the center of the cavity, $z=\frac{1}{2}L_{z}$, thus indicating its bulk-driven origin. In the bottom plane $z=0$ (shown in Fig. 7) and in the top plane $z=L_{z}$ (not shown in Fig. 7), the slip-velocity $\bm{v}_{2}^{\mathrm{slip}}$ is dominated by the gradient term in the expression (7) for the slip velocity, and has only a weak rotating component.

In Fig. 8, we compare the results for the acoustic streaming $\bm{v}_{2}$, see Eq. (7), for the double modes $l00+0l0$ with $l=2,4,6$. The streaming is shown in the horizontal $x$-$y$ center plane at height $z=\frac{1}{2}L_{z}$ (dark blue edge in Fig. 7) and in the vertical $x$-$z$ center plane at $y=\frac{1}{2}L_{y}$ (light green edge in Fig. 7). In the first row, Fig. 8(a)-(c), we show the bulk-driven streaming alone ($\bm{f}_{\mathrm{ac}}$ on, $\bm{v}_{2}^{\mathrm{slip}}$ off); in the second row, Fig. 8(d)-(f), we show the boundary-driven streaming alone ($\bm{f}_{\mathrm{ac}}$ off, $\bm{v}_{2}^{\mathrm{slip}}$ on); in the third row, Fig. 8(g)-(i), we show the total streaming ($\bm{f}_{\mathrm{ac}}$ on, $\bm{v}_{2}^{\mathrm{slip}}$ on). Comparing the bottom-row and top-row panels in Fig. 8 strongly indicate that the respective 4-by-4 and 6-by-6 streaming-roll patterns in Figs. 8(h) and (i) are bulk-driven. Furthermore, Fig. 8(a)-(c) show that the bulk-driven acoustic streaming becomes stronger for higher mode number $l$ and thus for higher frequency $f$. In contrast, the boundary-driven acoustic streaming Fig. 8(d)-(f) stays almost constant in amplitude and shows no frequency dependency. Note that the bulk-driven streaming in Fig. 8(a)-(c) has its maximum velocity $v_{2}^{\mathrm{max}}$ in the bulk of the fluid whereas the boundary-driven streaming in Fig. 8(d)-(f) has its maximum velocity $v_{2}^{\mathrm{max}}$ at the boundaries.

VI.3.4 Comparison with the analytical expression for the bulk-driven streaming

When the no-slip condition on the vertical side walls is ignored, we can derive an approximate analytical expression for the bulk-driven streaming $(\bm{v}_{2}^{\mathrm{blk}})_{l0}$ for the double modes $l00+0l0$. The details are given in Appendix D, where we find $(\bm{v}_{2}^{\mathrm{blk}})_{l0}$ to be proportional to the acoustic energy flux density $\bm{S}_{l0}=\bm{S}_{0l0}^{l00}$ of Eq. (18a),

[TABLE]

Inserting here $k_{l0}^{\star}=\frac{l\pi}{L^{\star}}$ and expression (18a) for $\bm{S}_{l0}$ with $\mathcal{F}_{l0}=\frac{3\sqrt{3}}{8}$ from Eq. (29c), we obtain the maximum value of the bulk-driven streaming velocity in the center plane $z=\frac{1}{2}L_{z}$,

[TABLE]

For the three cases shown in Fig. 8(a)-(c), where we have used $P_{1}^{l00}=P_{1}^{0l0}=1\,\textrm{MPa}$ and $l=2,4,6$, respectively, we calculate from Eq. (40) the maximum values 0.018, 0.060, and $0.102~{}\textrm{mm}/\textrm{s}$, which deviate less than 7 % from the corresponding results 0.017, 0.060, and $0.110~{}\textrm{mm}/\textrm{s}$ from the numerical simulation, listed in Fig. 8(a)-(c). For $l=2$, the analytical prediction of $0.018~{}\textrm{mm}/\textrm{s}$ over estimates by 5 % the numerical value $0.017~{}\textrm{mm}/\textrm{s}$ because it neglects the significant viscous damping from the vertical side walls in this case. For $l=6$, the analytical prediction of $0.102~{}\textrm{mm}/\textrm{s}$ under estimates by 7 % the numerical value $0.110~{}\textrm{mm}/\textrm{s}$, because the latter includes the coupling to the vertical 001-mode through the wide $\frac{1}{2}L^{\star}$-width vertical $G_{1}$-actuation as discussed in Section VI.3.1. When changing to the narrow $\frac{1}{12}L^{\star}$-width vertical $G_{1}$-actuation, also discussed in Section VI.3.1, the coupling to the vertical 001-mode diminish so much that the numerical streaming amplitudes becomes 0.017, 0.060, and $0.103~{}\textrm{mm}/\textrm{s}$ for $l=2,4,6$, respectively, where the latter now deviates less than 1 % from the analytical prediction $0.102~{}\textrm{mm}/\textrm{s}$. These numerical results for the acoustic streaming is shown in the file “fig_08_narrow_G1.pdf” in the Supplemental Material Note (1).

VI.3.5 Bulk-driven versus boundary-driven streaming rolls

We address the question on whether the horizontal streaming rolls in the $xy$-plane are due to the bulk-driven streaming or the boundary-driven streaming. We restrict the discussion to horizontal double modes $l00+0l0$, for which we have found the analytical expressions for the bulk-driven streaming $(\bm{v}_{2}^{\mathrm{blk}})_{l0}$, Eq. (39), and the slip velocity $\bm{v}_{2}^{\mathrm{slip}}$, Eq. (7), at the top ($z=L_{z}$) and bottom ($z=0$) boundaries.

According to Eq. (7), the boundary-driven streaming velocity at the double-mode resonance $l00+0l0$ is driven by a slip-velocity, which contain the term $\frac{1}{2}\kappa_{\mathrm{fl}}\bm{S}_{\mathrm{ac}}$ that rotates in the direction of the acoustic energy flux density $\bm{S}_{\mathrm{ac}}$, just as the bulk-driven streaming, see Eq. (39). The other terms in the slip velocity have similar magnitude, but they are gradient terms which cannot drive a rotating streaming in the horizontal $xy$-plane. Therefore, as seen in Fig. 8(d)-(f), the boundary-driven streaming does not lead to horizontal flow rolls as clearly as the bulk-driven streaming does, see Fig. 8(a)-(c). There may, however, exist horizontal planes at a few specific heights $z$, where the horizontal flow-roll pattern appears more clearly as suggested by Lei et al. Lei et al. (2014) and analyzed further by Skov et al. Skov et al. (2019).

By using Eqs. (7) and (39), we compute the ratio of the horizontally rotating bulk-driven streaming \big{|}{(\bm{v}_{2}^{\mathrm{blk}}(z=\frac{L_{z}}{2}))_{l0}}\big{|} in the center of the cavity to the horizontally rotating boundary-driven streaming \big{|}{\frac{1}{2}\kappa_{\mathrm{fl}}\bm{S}_{\mathrm{ac}}(z=0)}\big{|} at the bottom boundary,

[TABLE]

This velocity ratio takes values from 0 (boundary-driven streaming dominates) for a flat cavity, $lL_{z}\ll L^{\star}$, to $\frac{4}{3}+\frac{\eta^{\mathrm{b}}_{\mathrm{fl}}}{\eta_{\mathrm{fl}}}$ (bulk-driven streaming dominates) for a high cavity, $lL_{z}\gg L^{\star}$. For water, the velocity ratio (41) lies between 0 and 4.13, and in the setup of Figs. 4-8 with $L_{z}=0.1L^{\star}$ and water as the fluid, we obtain the values of the velocity ratio to be 0.38, 1.22 and 2.09 for $l=2$, 4 and 6, respectively. This result indicates that, that the 6-by-6-streaming roll pattern with $l=6$ observed by Hagsäter et al. Hagsäter et al. (2007) is predominantly bulk-driven, as also seen by comparing Fig. 8(c) and (i). In contrast, for the lower mode $l=2$, the horizontal streaming rolls are predominantly boundary-driven, which is seen by comparing Fig. 8(d) and (g).

The cross-over from boundary- to bulk-driven horizontal streaming rolls may be characterized by the height-to-width ratio $B=lL_{z}/L^{\star}$: we compute the critical value $B^{\mathrm{crit}}$, where the velocity ratio (41) is unity, and thus,

[TABLE]

For water, Eq. (41) leads to the ratio $B^{\mathrm{crit}}_{\mathrm{water}}=0.35$, and since $L_{z}=0.1L^{\star}$ in the setup of Figs. 4-8, we find that the horizontal streaming rolls in the double modes $l00+0l0$ with $l\leq 3$ are predominantly boundary-driven, while for $l\geq 4$ they are bulk-driven.

VI.3.6 Enhancement of the bulk-driven streaming velocity by increasing the bulk viscosity

For a given resonance mode with a given acoustic pressure amplitude in a given geometry, we obtain from Eqs. (40) and (7) the following scaling laws for the bulk- and boundary-driven streaming, respectively:

[TABLE]

To study the consequences of these scaling laws, we have chosen to compare water with the heterocyclic organic compound pyridine having the chemical formula C5H5N. Table 2 reveals that pyridine has the same material parameters (within 5 %) as water, but a bulk viscosity $\eta^{\mathrm{b}}_{\mathrm{fl}}$ that is $25.1$ times larger than that of water. Using these parameter values in Eq. (43), we find a large enhancement of the bulk-driven streaming, while the boundary-driven streaming is nearly unchanged,

[TABLE]

We test this prediction in Fig. 9, where we show simulations of the total acoustic streaming in pyridine in the same geometry as in Fig. 8(g)-(i) but with adjusted actuation amplitude $u_{1z}^{0}$ and frequency $f$ given in Table 3 to maintain the same pressure amplitudes $P_{1}^{l00}=P_{1}^{0l0}=1\ \textrm{MPa}$. As expected, we find that the total acoustic streaming for the three double modes is enhanced compared to that in water shown in Fig. 8(g)-(i), and the enhancement factors are 2.8, 8.2, and 13.9, respectively. The increasing factors reflect the increasing weight of bulk-driven relative to boundary-driven streaming as the mode index $l$ increases. The enhancement factors for the bulk-driven streaming alone (not shown in Fig. 9) are 21.4, 21.1, and 21.2, as predicted by Eq. (44a), while they are 1.2, 1.2, and 1.3 for boundary-driven streaming alone (also not shown in Fig. 9) in good agreement with Eq. (44b).

For all three cases $l=2$, 4, and 6, the horizontally rotating bulk-driven streaming in the center $z=\frac{1}{2}L_{z}$ is larger than the horizontally rotating boundary-driven streaming at the top-bottom boundaries, $z=0$ and $z=L_{z}$. This is also predicted by Eq. (VI.3.5), since $l\frac{L_{z}}{L^{\star}}=0.2,0.4,0.6>B^{\mathrm{crit}}_{\mathrm{pyridine}}=0.08$ for all three modes in pyridine.

VII Discussion

In the following, we raise a few discussion points related to the results of the above analysis, which are based on the single- and double mode resonances occurring in a closed box-shaped cavity with narrow viscous boundary layers. We note that the idealized, hard-wall, numerical model for the acoustic streaming is strikingly sensitive to small differences of only $0.1\%$ between the cavity side lengths $L_{x}$ and $L_{y}$. Our results for the resonance frequency in Eq. (10) and the optimal aspect ratio in Eq. (29f) are therefore crucial in the investigation of the largest possible bulk-driven acoustic streaming. This optimization was not taken into account in the experimental study by Hagsäter et al. Hagsäter et al. (2007), and the corresponding numerical simulation by Skov et al. Skov et al. (2019) of the shallow, nearly-square cavity of aspect ratio of $A=0.990$. As this aspect ratio differs from the optimal value $A=0.999$ for the hard-wall device in Fig. 8(i), it may be possible to increase the bulk-driven acoustic streaming significantly by fine tuning the geometry.

A priori, one may argue that it is difficult to obtain these optimal conditions with a relative precision of the order of 0.1 % in experiments. However, comparing our Fig. 7 with the experimental results in Hagsäter et al. Hagsäter et al. (2007), we find good qualitative agreement even though none of the optimal conditions were considered when designing the nearly-square fluid cavity. Another example is provided by the straight microchannel with a square cross section and non-rotating actuation that was studied experimentally by Antfolk et al. Antfolk et al. (2014). There, a single large vertical streaming roll similar to Fig. 2(b) was observed experimentally and analyzed numerically in terms of a double-mode excitation.

There may be several reasons for why the bulk-driven acoustic streaming can be observed experimentally in spite of the low 0.1-%-tolerance level of the aspect ratio suggested by our simplified, hard-wall, micro-cavity analysis in Fig. 6. Firstly, there are other damping mechanisms in a real experiment Hahn and Dual (2015), besides the bulk and boundary-layer damping considered here, which may broaden the resonance peaks significantly and thereby increase the change for overlapping resonances. Secondly, the double-mode body force $\bm{f}_{{l^{\prime}m^{\prime}n^{\prime}}}^{lmn}$ in Eq. (17a), which is the source of the bulk-driven acoustic streaming, depends on the product \big{|}{P_{1}^{lmn}}\big{|}\big{|}{P_{1}^{{l^{\prime}m^{\prime}n^{\prime}}}}\big{|} of the mode amplitudes, and therefore one can obtain a rotating body force even if the pressure is dominated by one of the modes, say $lmn$, as long as the other mode ${l^{\prime}m^{\prime}n^{\prime}}$ is just weakly represented. Thirdly, in typical experiments, the driving frequency is scanned over a broad span of frequencies, say 100 kHz around 2 MHz, which increases the possibility for exciting double-mode resonances. Fourthly, unintentional asymmetries and imperfections in the experimental setup may cause modes to overlap locally in some regions at some frequencies and in other regions at other frequencies, rather than the idealized homogeneous pattern at a single frequency shown in Fig. 7. This point was also mentioned by Hagsäter et al.: “If the frequency is shifted slightly in the vicinity of 2.17 MHz, the same vortex pattern will still be visible, but the strength distribution between the vortices will be altered”. Moreover, in the 3D simulation by Skov et al. Skov et al. (2019) including the surrounding solid and the piezo-electric transducer, this inhomogeneous distribution of the streaming was indeed reproduced in their Fig. 2(b2). Revisiting the experiment by Antfolk et al. Antfolk et al. (2014) of a long straight capillary of square cross section, this effect may also be observed, because it is likely that the aspect ratio is slightly less than unity in some regions of the capillary and slightly above unity in other regions, thereby causing the rotation direction to alternate along the channel.

We have in this work mainly considered horizontal acoustics with $n=0$ half-waves in the vertical $z$-direction. A similar analysis can be done for perpendicular acoustics (with respect to the top and bottom boundaries). In Sections VI.3.1 and VI.3.4, we have already briefly discussed the role of the vertical single-mode 001 for the pressure $p_{1}$ and the streaming $\bm{v}_{2}$. It is of course also possible to consider vertical rotating double modes involving, say, a single half-wave in the vertical direction, $lm1+l^{\prime}m^{\prime}1$. For relatively shallow cavities, a vertical excitation implies a relatively high frequency and a correspondingly increased body force $\bm{f}_{\mathrm{ac}}$. Furthermore, for perpendicular acoustics where $\partial_{\perp}v_{1\perp}\approx\bm{\nabla}\!\cdot\bm{v}_{1}=\mathrm{i}\omega\kappa_{\mathrm{fl}}p_{1}$, the slip-velocity $\bm{v}_{2}^{\mathrm{slip}}$ on the top and bottom boundaries is dominated by the first and last term in Eq. (7), \bm{v}_{2}^{\mathrm{slip}}\approx\frac{1}{2}\kappa_{\mathrm{fl}}\bm{S}_{\mathrm{ac}}+\frac{3}{2\omega}\big{\langle}[\mathrm{i}\partial_{\perp}v_{1\perp}]\bm{v}_{1\parallel}\big{\rangle}\approx-\kappa_{\mathrm{fl}}\bm{S}_{\mathrm{ac}}. Remarkably, this slip velocity is in the opposite direction of the body force $\bm{f}_{\mathrm{ac}}\propto\bm{S}_{\mathrm{ac}}$, so we expect that the bulk- and boundary-driven streaming will be in opposite directions, and thus the cavity could be designed such these streaming flows tend to cancel. A study of this effect is left for future work.

VIII Conclusion

We have derived analytical expressions for the single-mode pressure resonances $p_{1}^{lmn}=P_{1}^{lmn}F^{lmn}R^{lmn}$ and the corresponding boundary-layer damping coefficients $\Gamma_{\mathrm{bl}}^{lmn}$ for a closed box-shaped cavity with narrow viscous boundary layers, see Section III.1 and Appendix A. Based on these expressions, we have shown that the body force $\bm{f}_{\mathrm{ac}}$ can drive a strong bulk-driven acoustic streaming, if two overlapping single modes are excited simultaneously at the same frequency, thereby forming a double-mode resonance $lmn+{l^{\prime}m^{\prime}n^{\prime}}$, see Sections III.2 and IV. In contrast to the conventional wisdom that the bulk-driven streaming can be ignored on the scale of a few acoustic wavelengths, these double modes constitute an important example where the bulk-driven streaming is significant.

We have shown in Eq. (8) that the appearance of bulk-driven streaming in closed microcavities is directly related to the amount of rotation (acoustic angular momentum density) of the underlying acoustic fields. Whereas there is no rotation in a single mode, we have demonstrated in Sections IV.2 and V that the double modes can be rotating not only due to an externally-controlled rotating actuation, but, remarkably, also due to a non-rotating actuation coupled with a weak asymmetry in the geometry of the cavity. Furthermore, in Section V we have derived analytical expressions for the optimal aspect ratios that maximizes the bulk-driven streaming in a nearly-square cavity with a non-rotating actuation. We found that the bulk-driven acoustic streaming is sensitive to small deviations between the cavity side lengths of only $0.1\%$.

We have validated the theory by direct numerical simulation in Section VI of the double modes $200+020$, $400+040$, and $600+060$ in a nearly-square cavity with side lengths $L_{x}\approx 2000\,\upmu\textrm{m}$ and $L_{y}\approx 2000\,\upmu\textrm{m}$ and height $L_{z}=200\,\upmu\textrm{m}$, similar to the device investigated experimentally by Hagsäter et al. Hagsäter et al. (2007). Using the analytically-known optimal conditions (29), for which the bulk-driven acoustic streaming is maximized, we have avoided the otherwise required time-consuming parametric sweeping in geometry and frequency to locate these optimal conditions. The numerical results shown in Figs. 4-6 agrees with the analytical expressions with a relative accuracy better than 1 %.

In Fig. 7, the 3D simulation of the acoustic streaming of the double-mode $600+060$ is seen to reproduce the horizontally rotating 6-by-6 streaming-roll pattern observed by Hagsäter et al. In Fig. 8 we showed that this streaming pattern at the higher frequency 2.24 MHz is dominated by the bulk-driven streaming, whereas the streaming of the double mode 200+020 at the lower frequency 0.75 MHz is dominated by the boundary-driven streaming. This frequency-dependent crossover occurs because the body force $\bm{f}_{\mathrm{ac}}$ increases as frequency to the power two, while the slip velocity $\bm{v}_{2}^{\mathrm{slip}}$ depends only weakly on frequency, see Eq. (7). An analytical expression for computing the crossover for the double mode $l00+0l0$ is given in Eq. (41).

Finally, in agreement with Eckart Eckart (1948), we have shown in Eq. (39) that the bulk-driven acoustic streaming increases with the ratio $\frac{\eta^{\mathrm{b}}_{\mathrm{fl}}}{\eta_{\mathrm{fl}}}$ between the bulk and the dynamic viscosity. This dependency on $\eta^{\mathrm{b}}_{\mathrm{fl}}$ was studied numerically in Section VI.3.6 by exchanging the water in the microcavity with pyridine. As shown in Fig. 9, we found that the bulk-driven streaming for pyridine was enhanced by a factor $\simeq 21$ relative to water. This prediction obviously calls for experimental validation.

In this paper we have pointed out the significance of the bulk-driven acoustic streaming in resonating acoustic cavities even when the geometry length scale is comparable to a few half-wave lengths. The fundamental requirement for obtaining this effect, is the existence of two overlapping resonance modes excited simultaneously at the same frequency. This requirement is easily fulfilled in integrable geometries in 3D, such as rectangular, cylindrical, and spherical cavities, or in 2D, such as rectangular and circular cross sections of long, straight capillary channels, all of which are frequently encountered in experiments reported in the literature. The insight provided by our analysis is therefore relevant and important for acoustofluidics in general, both for fundamental studies and for technological applications.

Appendix A The acoustic pressure in a box with boundary layers

In this section, we derive the solution for the acoustic pressure $p_{1}$ with boundary layers in a box-shaped cavity $\Omega$ for which $0\leq x\leq L_{x}$, $0\leq y\leq L_{y}$, and $0\leq z\leq L_{z}$. The governing equation and boundary conditions for $p_{1}$ is given in Eq. (3). We consider first an actuation $U_{1\perp}(x,y)$ which is only non-zero at the bottom boundary $z=0$. The solution for a general actuation on all six boundaries is then finally constructed by superposition. Using separation of variables to solve the Helmholtz equation (3a) gives $p_{1}$ in the following form,

[TABLE]

Inserting this form into boundary condition (3b) gives,

[TABLE]

where $\perp$ is the inward direction perpendicular to the boundary, and the quantity $\alpha_{\perp}$ is the boundary-layer correction to the inward-normal derivative of $p_{1}$. Note that $\partial_{\perp}$ has a sign change at two opposite boundaries, say $x=0$ and $x=L_{x}$, for which $\partial_{\perp}=\partial_{x}$ and $\partial_{\perp}=-\partial_{x}$, respectively, whereas $\alpha_{\perp}(k_{\perp})$ depends on $k_{\perp}^{2}=k_{x}^{2}$, without a sign change at opposite boundaries. We assume in the following that the length scale introduced by $\alpha_{\perp}$ is much longer than the cavity length scale $L_{\perp}$,

[TABLE]

According to Eq. (46), the boundary conditions at the stationary boundaries at $x=0$ and $x=L_{x}$ are

[TABLE]

where $\alpha_{\perp}=\alpha_{x}$ in both equations, whereas $\partial_{\perp}=\partial_{x}$ in Eq. (48a), while $\partial_{\perp}=-\partial_{x}$ in Eq. (48b). Non-trivial solutions with $C_{x}$ and $S_{x}$ different from zero must satisfy the usual criterion for the equation determinant $\mathcal{D}_{x}(k_{x})$,

[TABLE]

We write the solutions $k_{x}=k_{x}^{l}$ to Eq. (49a) as perturbations away from the inviscid solutions $K_{x}^{l}$,

[TABLE]

where $\epsilon_{x}^{l}L_{x}\ll 1$. By inserting this into Eq. (49b), expanding in $\epsilon_{x}^{l}$, and writing $\alpha^{l}_{x}=\alpha_{x}(k_{x}^{l})$, we find,

[TABLE]

From assumption (47), $({\alpha_{x}^{l}})^{2}\epsilon_{x}^{l}L_{x}=[{\alpha_{x}^{l}}L_{x}]\epsilon_{x}^{l}{\alpha_{x}^{l}}$ is much smaller than $2\epsilon_{x}^{l}{\alpha_{x}^{l}}$, so it can be ignored. We then factorize out $K_{x}^{l}+\epsilon_{x}^{l}=k_{x}^{l}$ and obtain,

[TABLE]

For $l=0$, we have $K_{x}^{0}=0$ and $k_{x}^{0}(k_{x}^{0}-K_{x}^{0})=(k_{x}^{0})^{2}$, while for $l>0$ we have $k_{x}^{l}(k_{x}^{l}-K_{x}^{l})\approx\frac{1}{2}(k_{x}^{l}+K_{x}^{l})(k_{x}^{l}-K_{x}^{l})$ to first order in $\epsilon_{x}^{l}$. In either case, by introducing $L_{x}^{{l}}$ as

[TABLE]

the resulting expression for $\mathcal{D}_{x}(k_{x}^{l})$ for all $l$ is written as,

[TABLE]

The wavenumbers $k_{x}^{l}$ thus fulfil the zero-determinant criterion (49a) $\mathcal{D}_{x}(k_{x}^{l})=0$ if

[TABLE]

and similarly for $(k_{y}^{m})^{2}$ for all $l$ and $m$. From Eqs. (A) and (48a) with $k_{x}=k_{x}^{l}$, we obtain the corresponding eigenfunctions $X^{l}(x)$,

[TABLE]

where the prefactor $C_{x}^{l}$ is absorbed into $Z(z)$ in Eq. (A). We examine the ratio $\frac{\alpha_{x}^{l}}{k_{x}^{l}}$ for all $l$ using Eq. (55),

[TABLE]

By assumption (47), we note that even for $l=0$ this ratio is small, but not as small as for $l>0$.

By inserting $k_{x}^{l}$ from Eq. (55) into Eq. (56) and Taylor expanding in ${\alpha_{x}^{l}}/k_{x}^{l}$, we recover the usual inviscid hard-wall eigenfunctions \cos\big{(}K_{x}^{l}x\big{)} plus a small correction of the order $\alpha_{x}^{l}L_{x}$, due to the boundary layers,

[TABLE]

Note that $(\partial_{x}+\alpha_{x}^{l})X^{l}|_{x=0}=(\partial_{-x}+\alpha_{x}^{l})X^{l}|_{x=L_{x}}=0$ as required by the boundary condition Eq. (46). Similar expressions are valid for the $y$-eigenfunctions $Y^{m}(y)$.

In general, the pressure is an infinite sum of the eigenfunctions,

[TABLE]

where $k_{z}^{lm}$ depends on the angular frequency $\omega$ of the actuation through $k_{\mathrm{c}}$. At $z=L_{z}$, the sets of coefficients $C_{z}^{lm}$ and $S_{z}^{lm}$ satisfy a condition similar to Eq. (48b), which leads to

[TABLE]

where $\alpha_{z}^{lm}=\alpha_{z}(k_{z}^{lm})$. Further, at $z=0$ we have a condition similar to Eq. (48a) with $k_{0}^{2}\kappa_{\mathrm{fl}}^{-1}U_{1z}(x,y)$ on the right-hand side. Combined with Eqs. (59) and (60), the boundary condition at $z=0$ becomes,

[TABLE]

where $\mathcal{D}_{z}$ is defined similar to Eq. (49b). To find the coefficients $C^{lm}_{z}$, we write the wall actuation function $U_{1z}(x,y)$ as a generalized Fourier series in the functions $X^{l}(x)$ and $Y^{m}(y)$ that form a complete basis set on the interval $0\leq x\leq L_{x}$ and $0\leq y\leq L_{y}$ to order $\alpha_{x}^{l}L_{x}$ and $\alpha_{y}^{m}L_{y}$, respectively,

[TABLE]

Inserting the expansion (62a) into Eq. (61), we obtain the amplitudes $C_{z}^{lm}$,

[TABLE]

Finally, using Eq. (60) for $S_{z}^{lm}$ and Eq. (63) for $C_{z}^{lm}$ in Eq. (59) for $p_{1}$, yields the expression for the pressure,

[TABLE]

The infinite sum in Eq. (64a) is general and applies to all frequencies and actuations at the bottom boundary. For a given actuation $U_{1z}$, this solution is largest for frequencies $\omega=c_{\mathrm{fl}}k_{0}$ where $\mathcal{D}_{z}(k_{z}^{lm})$ is smallest, which gives the hard-wall resonances $lmn$ studied below. In real systems, the surrounding solid will have its own resonance properties and therefore, the actuation $U_{1z}$ will depend strongly on the frequency. In this case, the largest value of the prefactor $\frac{k_{0}^{2}k_{z}^{lm}\hat{U}_{1z}^{lm}}{\kappa_{\mathrm{fl}}\mathcal{D}_{z}(k_{z}^{lm})}$ is not necessarily found where $\mathcal{D}_{z}(k_{z}^{lm})$ is smallest.

A.1 Single-mode resonances in a box

Resonance occurs when $k_{0}=\frac{1}{c_{\mathrm{fl}}}\omega$ equals one of the values $k_{0}^{lmn}$ that minimizes $\mathcal{D}_{z}(k_{z}^{lm})$ in Eq. (64a). We expect these values of $k_{0}^{lmn}$ to be near the inviscid values $K_{0}^{lmn}$, and write similar to Eq. (50),

[TABLE]

where $\epsilon_{0}^{lmn}$ should be chosen to minimize $\mathcal{D}_{z}(k_{z}^{lm})$. The modes $l$ and $m$ in the $x$ and $y$ direction together with the resonance condition $k_{0}=k_{0}^{lmn}$ fixes the wavenumber $k_{z}^{lmn}=k_{z}^{lm}(k_{0}^{lmn})$ in Eq. (64c) pertaining to the $z$ direction,

[TABLE]

Combining this expression with Eq. (55) for $(k_{x}^{l})^{2}$ and $(k_{y}^{m})^{2}$, we find,

[TABLE]

To evaluate $\mathcal{D}_{z}(k_{z}^{lmn})$ in Eq. (64a), we note that $k_{z}^{lmn}$ is close to $K_{z}^{n}$, so we can use Eq. (54) with $(k_{x}^{l})^{2}-(K_{x}^{l})^{2}\rightarrow(k_{z}^{lmn})^{2}-(K_{z}^{n})^{2}$, which is found from Eq. (A.1),

[TABLE]

The (real) value of $\epsilon_{0}^{lmn}$, which minimizes this expression and leads to the resonance wavenumber $k_{0}^{lmn}$, is then,

[TABLE]

where we have inserted the expressions for $\alpha_{x}^{l}$, $\alpha_{y}^{m}$, and $\alpha_{z}^{n}$ defined by Eq. (46b), and introduced the boundary layer damping coefficients $\Gamma_{\mathrm{bl}}^{lmn}$.

We are now in a position to determine the third and last wavenumber $k_{z}^{lmn}$, the one in the actuation direction. Using the value of $\epsilon_{0}^{lmn}$ from Eq. (69) in expression (A.1) for $k_{z}^{lmn}$, we find,

[TABLE]

with $k_{z}^{n}$ defined similar to $k_{x}^{l}$ in Eq. (55) and the total damping coefficient given by,

[TABLE]

Note that in comparison with $k_{x}^{l}$ and $k_{y}^{m}$, the wavenumber $k_{z}^{lmn}$ has an additional dependency on $\Gamma^{lmn}$. The corresponding $z$-dependent functions $Z^{lmn}(z)$ are computed to lowest order in the small parameters by inserting Eq. (A.1) into Eq. (64b),

[TABLE]

Here, $Z^{lmn}(z)$ is close to the eigenfunction $Z^{n}(z)$ for stationary boundaries, analogous to $X^{l}(x)$ and $Y^{m}(y)$, see Eq. (56), but it contains an extra term, which satisfies the boundary condition at the moving actuated boundary $z=0$.

Finally, we evaluate $\mathcal{D}_{z}(k_{z}^{lmn})$ in Eq. (A.1) at resonance $k_{0}=k_{0}^{lmn}$ where $\epsilon_{0}^{lmn}$ takes the value given in Eq. (69),

[TABLE]

and by inserting Eqs. (73) and (72) into Eq. (64a), we find an expression for each resonance mode $p_{1}^{lmn}$,

[TABLE]

Note that for any $l$, $m$, and $n$, the product term $X^{l}(x)Y^{m}(y)Z^{n}(z)$ only satisfies the boundary condition Eq. (46a) for a stationary boundary $U_{1\perp}=0$. The actual boundary condition is satisfied by the infinite sum (64a) over terms $X^{l}(x)Y^{m}(y)Z_{\mathrm{act}}^{n}(z)$, because $\partial_{z}Z_{\mathrm{act}}^{n}=\frac{L_{z}^{{n}}}{(L_{z})^{2}}$ at $z=0$ and $\partial_{z}Z_{\mathrm{act}}^{n}=0$ at $z=L_{z}$.

As Eq. (74) is only valid exactly at resonance $k_{0}=k_{0}^{lmn}$, it must be corrected to deal with frequencies $k_{0}\approx k_{0}^{lmn}$. This is done by adding $k_{0}-k_{0}^{lmn}$ to the right-hand side of Eq. (65). By Eq. (A.1), this procedure is seen to be equivalent to the substitution

[TABLE]

in Eq. (74).

A.2 Approximate solutions near resonance for actuation at all boundaries

At resonance, the expression (74) satisfies all boundary conditions to first order in $\alpha_{x}^{l}L_{x},\alpha_{y}^{m}L_{y},$ and $\alpha_{z}^{n}L_{z}$. Ignoring these first-order corrections as well as the small $Z^{n}_{\mathrm{act}}$-term in Eq. (74), the eigenfunctions are approximately equal to the usual hard-wall eigenfunctions $R^{lmn}$,

[TABLE]

To generalize our results from an actuation acting on only the bottom boundary $z=0$ to all six boundaries, we note first that $L_{z}^{{n}}\approx\int_{0}^{L_{z}}[Z^{n}(z))]^{2}\,\mathrm{d}z$, whereby Eq. (62b) can be rewritten as

[TABLE]

Using this expression together with the substitution (75) valid for $k_{0}\approx k_{0}^{lmn}$, we obtain the following expression for the resonance modes $lmn$ in a box $\Omega$ with viscous boundary layers and an inward normal displacement $U_{1\perp}(\bm{r})$ specified on all six boundaries $\partial\Omega$, Eq. (3c),

[TABLE]

The quality factor $Q^{lmn}$ for resonance mode $lmn$ is obtained from the linewidth of the acoustic energy density $E_{\mathrm{ac}}$ in Eq. (6), and since $E_{\mathrm{ac}}\propto|p_{1}^{lmn}|^{2}\propto|F^{lmn}|^{2}$, we find

[TABLE]

where $\Gamma^{lmn}$ is given in Eq. (71).

Appendix B Maximum acoustic rotation

In the following, we evaluate the frequency-dependent factor $\mathcal{F}_{ml0}^{lm0}$ in Eq. (17b) for rotating and non-rotating actuations of the nearly-square cavity considered in Section V. For convenience, we first introduce the double-mode quantities, $\tilde{\xi}_{lm}$ and $\phi_{lm}^{\mathrm{act}}$,

[TABLE]

where $\tilde{k}_{0}$ and $\tilde{k}_{lm}^{\star}$ are defined in Eq. (27), and $\phi_{lm}^{\mathrm{act}}$ is the difference in the phases by which the modes $lm0$ and $ml0$ are actuated, see Eq. (12b). By using these quantities and the expressions for $F^{lm0}$ and $F^{ml0}$ in Eq. (28), we calculate the frequency-dependency $\mathcal{F}_{ml0}^{lm0}$ from Eq. (17b),

[TABLE]

This equation expresses the magnitude of the acoustic rotation as a function of the actuation phase difference $\phi_{lm}^{\mathrm{act}}$, the actuation frequency $\tilde{\xi}_{lm}$, and the mode separation $\tilde{\Delta}_{lm}$, which by Eq. (26) is related to the aspect ratio $A$. In the following, we study the non-rotating actuation with $\phi_{lm}^{\mathrm{act}}=0$ and the rotating actuation with $\phi_{lm}^{\mathrm{act}}=\pm\frac{1}{2}\pi$, and we determine the optimal mode separation $\tilde{\Delta}_{lm}^{\mathrm{opt}}$ and frequency $\tilde{\xi}_{lm}^{\mathrm{opt}}$, for which $\mathcal{F}_{ml0}^{lm0}$ takes its largest positive or negative value $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}$.

B.1 Extremum values of $\mathcal{F}_{ml0}^{lm0}$ for a non-rotating actuation $\phi_{lm}^{\mathrm{act}}=0$

For a non-rotating actuation, we have $\phi_{lm}^{\mathrm{act}}=0$, and Eq. (81) becomes

[TABLE]

The extremum values of this expressions are

[TABLE]

In Fig. 6 (“non-rotating act.”) we plot the value $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}$ from this expression as a function of $\tilde{\Delta}_{lm}$. The optimal value $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{opt}}$ of $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}$ is found at the optimal mode separation $\tilde{\Delta}_{lm}^{\mathrm{opt}}$, and optimal rescaled frequency $\tilde{\xi}_{lm}^{\mathrm{opt}}$,

[TABLE]

B.2 Extremum values of $\mathcal{F}_{ml0}^{lm0}$ for a rotating actuation $\phi_{lm}^{\mathrm{act}}=\pm\frac{1}{2}\pi$

For an externally controlled rotating actuation, we have $\phi_{lm}^{\mathrm{act}}=\pm\frac{1}{2}\pi$ and Eq. (81) becomes,

[TABLE]

This expression is more complicated than Eq. (82), since it only has one extremum for a sufficiently small mode separation $\tilde{\Delta}_{lm}$, namely where the angular momentum and the actuation are co-rotating, $\mathrm{sign}(\mathcal{F}_{ml0}^{lm0})=\mathrm{sgn}_{\phi}^{lm}$. For larger values of $\tilde{\Delta}_{lm}$, the frequency can be tuned to two different extremum values, namely both the co- and counter-rotating cases, $\mathrm{sign}(\mathcal{F}_{ml0}^{lm0})=\pm\mathrm{sgn}_{\phi}^{lm}$. We therefore identify two branches of extremum values. The co-rotating branch $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}_{\mathrm{co}}$ with subscript “co” and the counter-rotating branch $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}_{\mathrm{cntr}}$ with subscript “cntr”.

For the co-rotating branch, $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}_{\mathrm{co}}$ is

[TABLE]

For the counter-rotating branch, $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}_{\mathrm{cntr}}$ is

[TABLE]

In Fig. 6 (Rotating act.), we plot for $\phi_{lm}^{\mathrm{act}}=+\frac{1}{2}\pi$ both the co- and counter rotating branches, $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}_{\mathrm{co}}$ and $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{extr}}_{\mathrm{cntr}}$.

From Eqs. (86a) and (87a), we find that the largest value $(\mathcal{F}_{ml0}^{lm0})^{\mathrm{opt}}$ of $(\mathcal{F}_{ml0}^{lm0})$ is the co-rotating double mode found at the optimal mode separation $\tilde{\Delta}_{lm}^{\mathrm{opt}}$, and the rescaled frequency $\tilde{\xi}_{lm}^{\mathrm{opt}}$,

[TABLE]

We see that in contrast to the non-rotating actuation (84), the rotating actuation optimizes the acoustic rotation for zero mode separation, $\tilde{\Delta}_{lm}^{\mathrm{opt}}=0$.

Appendix C The rotating actuation G16 of the double mode 200+020

In the simulation shown in Fig. 5(d-f) of Section VI.3.2 we excite the double mode $200+020$ by using the rotating actuation,

[TABLE]

To obtain the complex pressure amplitudes $P_{1}^{200}$ and $P_{1}^{020}$ resulting from $G_{16}$, we first evaluate the integral in the denominator of Eq. (36), where $G_{1}$ is substituted by $G_{16}$,

[TABLE]

We then use Eq. (36) to compute the desired amplitudes,

[TABLE]

We see that the actuation in Eq. (89a) excites the two modes $200$ and $020$ with a phase difference of

[TABLE]

Finally, we use Eq. (91) to choose the amplitude $u_{1z}^{0}$ of the actuation such that $|P_{1}^{200}|=|P_{1}^{020}|=1$ MPa,

[TABLE]

where we used the parameters for water given in Table 2.

Appendix D Approximate solution for the bulk-driven streaming from double modes $l00+0l0$

In this section, we calculate an approximate expression for the bulk-driven acoustic streaming velocity $(\bm{v}_{2}^{\mathrm{blk}})_{0l0}^{l00}=(\bm{v}_{2}^{\mathrm{blk}})_{l0}$ driven by the acoustic body force $\bm{f}_{\mathrm{ac}}=\bm{f}_{0l0}^{l00}=\bm{f}_{l0}$ resulting from the combination of the two perpendicular, horizontal modes $l00$ and $0l0$ in a square cavity. As in Eq. (V), the double subscript $l0$ refers to the double mode $l00+0l0$. The bulk-driven streaming velocity $(\bm{v}_{2}^{\mathrm{blk}})_{l0}$ satisfies an equation similar to Eq. (7) except for having no-slip at the boundaries,

[TABLE]

Here, $\bm{S}_{l0}$ is defined in Eq. (18a) combined with Eq. (V) for $\mathcal{F}_{l0}(k_{0})$, Eq. (17c) for $\bm{\mathcal{K}}_{l0}(\bm{r})$, and Eq. (9a) for $R_{l0}(\bm{r})$,

[TABLE]

The impact on the bulk-driven streaming from the no-slip condition at the side walls will decrease exponentially with the length scale $w$ away from the side walls, where $w$ is the minimum shear length scale, i.e.  either the vertical length scale $L_{z}$ or the horizontal length scale $(k_{l0}^{\star})^{-1}$,

[TABLE]

For the narrow channel considered in Section VI where $L_{z}/L^{\star}=0.1$, we have $w/L^{\star}=0.1$ for $l=2$, $w/L^{\star}=\frac{1}{4\pi}=0.08$ for $l=4$, and $w/L^{\star}=\frac{1}{6\pi}=0.05$ for $l=6$, and therefore the vertical side walls can be ignored in the majority of the cavity for these three double modes.

To compute the acoustic streaming velocity, we ignore the side walls and assume that $(\bm{v}_{2}^{\mathrm{blk}})_{l0}$ is proportional to the $xy$-dependent acoustic energy flux density $\bm{S}_{l0}(x,y)$ and some $z$-dependent function $\zeta(z)$ to be found,

[TABLE]

By using $\bm{S}_{l0}$ from Eq. (95a), we find the three identities,

[TABLE]

Consequently, by Eq. (98a), the ansatz (97) satisfies the continuity equation (94a). Moreover, inserting (97) into Eq. (94b) using \Gamma_{\mathrm{fl}}=\Big{(}\frac{4}{3}+\frac{\eta^{\mathrm{b}}_{\mathrm{fl}}}{\eta_{\mathrm{fl}}}\Big{)}\eta_{\mathrm{fl}}\kappa_{\mathrm{fl}}\omega from Eq. (4), and Eq. (98b) we find,

[TABLE]

Since $\bm{S}_{l0}(x,y)$ has no $z$-component, $(p_{2}^{\mathrm{blk}})_{l0}$ on the left-hand side is independent of $z$. Therefore, the square bracket on the right-hand side must be a constant. This constant must be zero, since otherwise taking the curl of Eq. (99) would lead to $\bm{0}=\bm{\nabla}\times\bm{S}_{l0}$ in conflict with Eq. (98c). The condition of a vanishing square bracket in Eq. (99) and the no-slip boundary conditions $\zeta(0)=\zeta(L_{z})=0$, lead to the following expression for $\zeta(z)$,

[TABLE]

The solution for $(\bm{v}_{2}^{\mathrm{blk}})_{l0}$ is then obtained by inserting Eqs. (95) and (100) in Eq. (97).