2. Problem definition and governing equations

Figure 1 illustrates the typical geometry of a microbubble used to drive streaming flows in a microfluidic device: the semi-cylindrical bubble of radius $a$ protrudes from a blind side channel of width $2a$ branching off a main channel, whose opposite wall (in the $y$ -direction) is at a distance of $H=ha$ . The entire structure has a uniform depth $D=da$ in the $z$ -direction for ease of manufacture, so that the bubble axis spans the depth of the microchannel, see figure 1(a). The bubble contact lines at ( $x=\pm a,y=0$ ) are pinned to the side channel corners. The channel is filled with a fluid of kinematic viscosity ${\it\nu}$ and density  ${\it\rho}$ . The application of an acoustic driving pressure at an angular frequency ${\it\omega}=2{\rm\pi}f$ (via a piezoelectric transducer) drives harmonic oscillations of the bubble surface at a characteristic amplitude ${\it\epsilon}a$ , where ${\it\epsilon}\ll 1$ . Note that for the typical sizes ( $a\sim 50~{\rm\mu}\text{m}$ ) and frequencies ( $f\sim 20~\text{kHz}$ ) the acoustic wavelength is much greater than the bubble size. Using $a$ and ${\it\omega}^{-1}$ to normalize length and time respectively, the interface shape for small-amplitude oscillations can in general be written as

(2.1) $$\begin{eqnarray}R({\it\theta},z,t)=1-\text{i}{\it\epsilon}\text{e}^{\text{i}t}{\it\zeta}({\it\theta},z),\end{eqnarray}$$

where ${\it\zeta}({\it\theta},z)$ is a complex function and it is understood that only the real part of any complex expression has physical significance. Here, ${\it\theta}=\tan ^{-1}(y/x)$ is the polar angle. Oscillations of the interface drive a flow field which, in the limit of small oscillation amplitude, is conveniently solved for by an expansion in powers of ${\it\epsilon}$ (Riley Reference Riley1967; Davidson & Riley Reference Davidson and Riley1971; Longuet-Higgins Reference Longuet-Higgins1998; Riley Reference Riley2001; Sadhal Reference Sadhal2012; Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014).

The primary oscillatory flow $\boldsymbol{u}_{0}$ (characteristic scale ${\it\epsilon}a{\it\omega}$ ) is irrotational except in viscous boundary layers of dimensionless thickness ${\it\delta}\equiv a^{-1}\sqrt{2{\it\nu}/{\it\omega}}\ll 1$ , near the boundaries of the domain (Longuet-Higgins Reference Longuet-Higgins1953; Nyborg Reference Nyborg1958; Longuet-Higgins Reference Longuet-Higgins1998). Steady Reynolds stresses due to the oscillatory flow drive a secondary steady Eulerian flow $\langle \boldsymbol{u}_{1}\rangle$ ( $\langle \cdot \rangle$ denotes the time average over an oscillation cycle), with a characteristic velocity scale $u_{s}={\it\epsilon}^{2}a{\it\omega}$ and a streaming Reynolds number $\mathit{Re}_{s}=u_{s}a/{\it\nu}$ , which is $\ll 1$ in many practical situations (Nyborg Reference Nyborg1958; Wu & Du Reference Wu and Du1997; Longuet-Higgins Reference Longuet-Higgins1998; Riley Reference Riley2001), including microfluidics applications (Wang et al. Reference Wang, Jalikop and Hilgenfeldt2011, Reference Wang, Rallabandi and Hilgenfeldt2013; Sadhal Reference Sadhal2012). The time-averaged motion of fluid elements under the combination of the primary and secondary flows can be described in terms of a steady Lagrangian velocity field $\boldsymbol{U}=\langle \boldsymbol{u}_{1}\rangle +\boldsymbol{u}_{d}$ , where $\boldsymbol{u}_{d}$ is the Stokes drift, defined by

(2.2) $$\begin{eqnarray}\boldsymbol{u}_{d}=\left\langle \left(\int \boldsymbol{u}_{0}\,\text{d}t\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)\boldsymbol{u}_{0}\right\rangle .\end{eqnarray}$$

It is this time-averaged Lagrangian motion of the fluid that is practically relevant for fluid transport and is normally referred to as the ‘steady streaming’. Figure 1(b) shows the typical appearance of the steady streaming flow from a sessile bubble as viewed along the axis of the bubble. It was observed experimentally (Wang Reference Wang2013; Wang et al. Reference Wang, Rallabandi and Hilgenfeldt2013) and shown theoretically (Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014) that this pattern of a vortex pair in ‘fountain’ orientation (flow upwards at $x=0$ ) persists virtually unchanged over a range of frequencies. It was also observed (Wang Reference Wang2013) that the pattern does not vary much with the depth $z$ (compare out-of-focus tracks with the in-focus tracks closer to $z=0$ , in figure 1 b). However, the latter observations of the two-dimensional character of the flow were limited by the microscopic depth of field ( ${\gtrsim}10~{\rm\mu}\text{m}$ ) and the general inability to follow individual particles over long times. We will show here that new flow features appear beyond two-dimensional flow, and that they can be verified by more sophisticated experimental techniques.

For a theoretical description of the streaming $\boldsymbol{U}$ , it is convenient to first write the Reynolds stress term ${\bf\sigma}\equiv \langle \boldsymbol{u}_{0}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{0}\rangle$ as ${\bf\sigma}=\langle \boldsymbol{{\rm\nabla}}|\boldsymbol{u}_{0}|^{2}/2\rangle +\langle (\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}_{0})\times \boldsymbol{u}_{0}\rangle$ . Using the viscous stress scale ${\it\rho}{\it\nu}u_{s}/a$ to non-dimensionalize the time-averaged pressure, the dimensionless governing equations for the Lagrangian steady flow in the limit $\mathit{Re}_{s}\ll 1$ are

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{U}=0,\\ \displaystyle {\rm\nabla}^{2}\boldsymbol{U}={\rm\nabla}^{2}\boldsymbol{u}_{d}+\boldsymbol{{\rm\nabla}}p+\boldsymbol{f},\end{array}\right\}\end{eqnarray}$$

where

(2.4a,b ) $$\begin{eqnarray}p=\langle p_{1}\rangle +\frac{2}{{\it\delta}^{2}}\left\langle \frac{|\boldsymbol{u}_{0}|^{2}}{2}\right\rangle ,\quad \text{and}\quad \boldsymbol{f}=\frac{2}{{\it\delta}^{2}}\langle (\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}_{0})\times \boldsymbol{u}_{0}\rangle .\end{eqnarray}$$

Here, $p$ is a dimensionless effective pressure that incorporates both the time-averaged fluid pressure $\langle p_{1}\rangle$ and the kinetic energy of the oscillatory flow, and $\boldsymbol{f}$ is a body force related to the advection of oscillatory vorticity, which is confined to viscous boundary layers (Riley Reference Riley2001). The Lagrangian mean flow is no-penetration with a vanishing tangential stress at the mean bubble surface $\partial {\it\Omega}_{b}$ , and satisfies no-slip boundary conditions at the rigid wall surfaces $\partial {\it\Omega}_{w}$ ,

(2.5) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{U}=0\quad \text{on }\partial {\it\Omega}_{w}\\ \displaystyle U_{r}=0\quad \text{on }\partial {\it\Omega}_{b}\\ \displaystyle \unicode[STIX]{x1D64E}_{r{\it\theta}}=\unicode[STIX]{x1D64E}_{rz}=0,\quad \text{on }\partial {\it\Omega}_{b},\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D64E}$ is the symmetric part of $\boldsymbol{{\rm\nabla}}\boldsymbol{U}$ . Outside viscous boundary layers, where $\boldsymbol{f}$ is exponentially small, the Eulerian steady flow $\langle \boldsymbol{u}_{1}\rangle =\boldsymbol{U}-\boldsymbol{u}_{d}$ satisfies the Stokes equations, supported by the effective pressure  $p$ .

3. Axial streaming

To motivate the streaming solutions that will be developed in this section, we first provide the general form of the oscillatory flow $\boldsymbol{u}_{0}$ , which ultimately drives (2.3) through the Reynolds stress and the Stokes drift. For typical practical applications, ${\it\delta}\ll 1$ , allowing $\boldsymbol{u}_{0}$ to be evaluated using matched asymptotic expansions, with an ‘outer’ potential flow outside boundary layers and ‘inner’ (viscous boundary layer) solutions described by the unsteady Stokes equations (Longuet-Higgins Reference Longuet-Higgins1953; Nyborg Reference Nyborg1958; Riley Reference Riley2001). The oscillatory outer flow may in general be expressed as a modal decomposition

(3.1) $$\begin{eqnarray}\boldsymbol{u}_{0}=\mathop{\sum }_{n=0}^{\infty }\mathop{\sum }_{m=0}^{\infty }\boldsymbol{u}_{0}^{mn},\end{eqnarray}$$

where

(3.2) $$\begin{eqnarray}\boldsymbol{u}_{0}^{mn}=\text{e}^{\text{i}t}A_{mn}(2n{\it\alpha}\text{K}_{2m}^{\prime }(2n{\it\alpha}))^{-1}\boldsymbol{{\rm\nabla}}{\it\varphi}_{mn},\end{eqnarray}$$

and

(3.3) $$\begin{eqnarray}{\it\varphi}_{mn}=\text{K}_{2m}(2n{\it\alpha}r)\cos 2m{\it\theta}\cos 2n{\it\alpha}z\end{eqnarray}$$

is a cylindrical harmonic function ( ${\rm\nabla}^{2}{\it\varphi}_{mn}=0$ ). Here, $\text{K}_{{\it\nu}}$ represents a modified Bessel function of the second kind, with $\text{K}_{{\it\nu}}^{\prime }(z)=\text{d}/\text{d}z(\text{K}_{{\it\nu}}(z))$ and ${\it\alpha}\equiv {\rm\pi}/d$ . $A_{mn}$ are the (generally complex) coefficients of the Fourier decomposition of the interface deformation ${\it\zeta}({\it\theta},z)$ and are determined by bubble interface dynamics and the pinning of contact lines (Wang et al. Reference Wang, Rallabandi and Hilgenfeldt2013).

Figure 2. (a) Schematic of the flow superposition of modes in (3.18), indicating the planar fountain $\boldsymbol{u}$ , the axisymmetric fountain $c_{01}\boldsymbol{v}_{01}$ , and the mixed axial–azimuthal mode $c_{11}\boldsymbol{v}_{11}$ . (b) Streamline portrait of the planar flow $\boldsymbol{u}$ for a channel height $h=6.25$ , indicating limiting streamlines ${\it\psi}$ : ${\it\psi}={\it\psi}_{m}$ corresponds to the vortex centre of the two-dimensional flow (smallest streamline) and ${\it\psi}=0$ , the largest streamline. (c) Flow lines of the axisymmetric flow in an $r{-}z$ plane, showing a fountain orientation (radially outward velocities at $z=0$ ). The shaded region $r<1$ indicates the bubble.

Of this general combination of modes, planar modes ( $n=0$ ) have the largest amplitudes due to their relatively small damping, driving flows that decay algebraically away from the interface. On the other hand, axial modes ( $n>0$ ), which drive flows that decay exponentially ( ${\sim}\text{e}^{-2n{\it\alpha}r}$ ) away from the bubble surface are only weakly excited due to greater damping (see e.g. Sanz & Diez Reference Sanz and Diez1989; Wang et al. Reference Wang, Rallabandi and Hilgenfeldt2013). Therefore, to leading order, only planar modes of $\boldsymbol{u}_{0}$ contribute to (2.3), driving two-dimensional streaming in the absence of the confining top and bottom walls (Wang et al. Reference Wang, Rallabandi and Hilgenfeldt2013; Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014), see figure 2(b). Such a two-dimensional description of the flow has been used successfully to explain the experimentally observed streaming patterns and velocity fields in the $x{-}y$ plane in the case of oscillating bubbles (Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014), cf. figure 1(b), and solid cylinders (Lutz et al. Reference Lutz, Chen and Schwartz2006). In general, however, steady axial flow components (in the $z$ -direction) should be present.

Such axial velocity components of the steady streaming may be excited by two separate effects that are neglected in the two-dimensional theory: (i) the top and bottom walls necessitate that the steady fluid velocity vanish identically at the top and bottom walls ( $z=\pm d/2$ ), and (ii) axial oscillation modes of the bubble introduce Reynolds stress components that drive axial motion of the fluid.

We focus our attention on solutions of (2.3) that are valid outside steady boundary layers of thickness ${\it\delta}\ll 1$ near the boundaries of the domain. Formally, we are interested in the outer steady solution of a matched asymptotic expansion; this will be understood henceforth unless otherwise stated. Linearity of the governing equations makes it convenient to express the (three-dimensional) Lagrangian steady velocity field $\boldsymbol{U}$ as a superposition of a two-dimensional steady solution $\boldsymbol{u}$ (cf. figure 2 a,b) and a secondary solution $\boldsymbol{v}$ with non-zero axial velocity components, i.e.

(3.4) $$\begin{eqnarray}\boldsymbol{U}(r,{\it\theta},z)=\boldsymbol{u}(r,{\it\theta})+\boldsymbol{v}(r,{\it\theta},z),\end{eqnarray}$$

so that $\boldsymbol{v}$ satisfies (2.3), with $\boldsymbol{u}_{d}$ , $p$ and $\boldsymbol{f}$ now taken expressly to be axially non-uniform, while planar components are addressed by  $\boldsymbol{u}$ .

The outer flow must also satisfy appropriate matching conditions for proper asymptotic matching with the boundary layer solutions at the mean bubble surface ( $r=1$ ) and at the top and bottom walls ( $z=\pm d/2$ ). Since normal velocity variations are small over the boundary layer thickness, kinematic boundary conditions are directly applicable to the outer solution at all boundaries $\partial {\it\Omega}=\partial {\it\Omega}_{w}+\partial {\it\Omega}_{b}$ ,

(3.5) $$\begin{eqnarray}\displaystyle \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}_{\partial {\it\Omega}}=0\quad \text{on }\partial {\it\Omega},\end{eqnarray}$$

where $\boldsymbol{n}_{\partial {\it\Omega}}$ is the unit normal to $\partial {\it\Omega}$ .

No-slip conditions at the top and bottom walls are enforced with the help of steady boundary layers, which in general support large axial gradients of planar velocity components near $z=\pm d/2$ . This establishes effective slip conditions on the Lagrangian outer flow at the top and bottom walls (generally a slip velocity of $O(1)$ ) that depend on the particulars of the oscillatory flow $\boldsymbol{u}_{0}$ (see e.g Longuet-Higgins Reference Longuet-Higgins1953; Nyborg Reference Nyborg1958). Consider now the mixed boundary layer region near the contact line at the bottom wall, defined by $r=1+O({\it\delta})$ , $z=-d/2+O({\it\delta})$ , where the structure of the no-slip wall boundary layer must conform to that of the no-stress bubble boundary layer. Within this region, the pinning of contact lines, reflected in (3.5), ensures that $U_{r}\sim U_{z}\sim O({\it\delta})$ , which along with no-stress conditions at the bubble ( $\partial _{r}(r^{-1}U_{{\it\theta}})+r^{-2}\partial _{{\it\theta}}U_{r}=0$ ), requires that $U_{{\it\theta}}\sim O({\it\delta}^{2})$ . By symmetry, the same arguments are applicable to the flow in the vicinity of the contact line between the bubble and the top wall $r=1,z=d/2$ . All boundary layer velocity components are therefore $O({\it\delta})$ or smaller in the vicinity of the contact lines at the top and bottom walls. For proper asymptotic matching, the Lagrangian outer flow must therefore, to leading order in ${\it\delta}$ , satisfy no-slip conditions at $r=1$ , $z=\pm d/2$ , which may be written as

(3.6) $$\begin{eqnarray}\boldsymbol{v}=-\boldsymbol{u}\quad \text{on }r=1,z=\pm d/2.\end{eqnarray}$$

The bubble boundary conditions are more involved, since Reynolds stresses at the bubble may in general depend on a large number of mode combinations that are sensitive to the oscillatory boundary layer structure (see e.g. Davidson & Riley Reference Davidson and Riley1971; Longuet-Higgins Reference Longuet-Higgins1998; Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014). Rather than compute these conditions exactly, we analyse the azimuthal and axial symmetries of these terms, which are ultimately manifested in the streaming through the bubble boundary conditions (2.5) enforced at $r=1$ . Let us consider the Reynolds stress term engendered by the product of two arbitrary oscillatory flow modes $\boldsymbol{u}_{0}^{kl}$ and $\boldsymbol{u}_{0}^{mn}$ of (3.2), defined by ${\bf\sigma}^{kl,mn}=\langle \boldsymbol{u}_{0}^{mn}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{0}^{kl}+\boldsymbol{u}_{0}^{kl}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{0}^{mn}\rangle$ . Outside steady boundary layers, we have

(3.7) $$\begin{eqnarray}\displaystyle {\bf\sigma}^{kl,mn} & = & \displaystyle \frac{A_{mn}^{\ast }A_{kl}}{8nl{\it\alpha}^{2}\,\text{K}_{2m}^{\prime }(2n{\it\alpha})\text{K}_{2k}^{\prime }(2l{\it\alpha})}\boldsymbol{{\rm\nabla}}(\boldsymbol{{\rm\nabla}}{\it\varphi}_{mn}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\varphi}_{kl})\nonumber\\ \displaystyle & = & \displaystyle \frac{A_{mn}^{\ast }A_{kl}}{8nl{\it\alpha}^{2}\,\text{K}_{2m}^{\prime }(2n{\it\alpha})\text{K}_{2k}^{\prime }(2l{\it\alpha})}\boldsymbol{{\rm\nabla}}\mathop{\sum }_{i=1}^{2}\mathop{\sum }_{j=1}^{2}{\it\lambda}_{ij}^{kl,mn}(r)\cos [2(m+(-1)^{i}k){\it\theta}]\nonumber\\ \displaystyle & & \displaystyle \times \,\cos [2(n+(-1)^{j}l){\it\alpha}z],\end{eqnarray}$$

where the asterisk denotes the complex conjugate and

(3.8) $$\begin{eqnarray}{\it\lambda}_{ij}^{kl,mn}={\textstyle \frac{1}{4}}g_{mn}^{\prime }(r)g_{kl}^{\prime }(r)-[(-1)^{i}mk+(-1)^{j}{\it\alpha}^{2}nl]g_{mn}(r)g_{kl}(r),\end{eqnarray}$$

with $g_{mn}(r)\equiv \text{K}_{2m}(2{\it\alpha}nr)$ . The same azimuthal and axial symmetries also appear in analogous modes of the Stokes drift terms $\boldsymbol{u}_{d}^{kl,mn}$ defined by

(3.9) $$\begin{eqnarray}\boldsymbol{u}_{d}^{kl,mn}\equiv \left\langle \left(\int \boldsymbol{u}_{0}^{mn}\,\text{d}t\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)\boldsymbol{u}_{0}^{kl}+\left(\int \boldsymbol{u}_{0}^{kl}\,\text{d}t\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)\boldsymbol{u}_{0}^{mn}\right\rangle .\end{eqnarray}$$

Through (2.3), these axial and azimuthal dependences of ${\bf\sigma}^{kl,mn}$ and $\boldsymbol{u}_{d}^{kl,mn}$ are reflected directly in the steady outer solution (see e.g. Wang Reference Wang1968; Davidson & Riley Reference Davidson and Riley1971; Longuet-Higgins Reference Longuet-Higgins1998).

To first approximation, due to their large amplitudes, products of planar modes ( $n=0$ ) contribute most strongly to (3.7), driving two-dimensional flow. Non-trivial axial dependences of ${\bf\sigma}$ arise predominantly from the product of the most strongly excited axial mode $m=0,n=1$ (oscillating axisymmetric quadrupole) with the dominant two-dimensional mode $k=l=0$ (oscillating line source). The resulting forcing ${\bf\sigma}^{00,01}$ , by (3.7), is axisymmetric, with an axial wavelength equal to the channel depth $d$ , forming the most strongly excited mode of the forcing term ${\bf\sigma}$ . Decomposed into its components in cylindrical coordinates, and neglecting higher-order mode combinations, ${\bf\sigma}$ therefore takes the approximate form

(3.10) $$\begin{eqnarray}{\bf\sigma}\approx {\bf\sigma}^{00,01}=A_{00}A_{01}^{\ast }\{{\it\sigma}_{r}(r)\cos 2{\it\alpha}z,0,{\it\sigma}_{z}(r)\sin 2{\it\alpha}z\},\end{eqnarray}$$

where ${\it\sigma}_{r}(r)$ and ${\it\sigma}_{z}(r)$ are determined by the oscillatory flow functions. Note also that retaining the same axial modes $(0,1)$ , but including higher-order planar modes $(k>0,l=0)$ results in Reynolds stresses with both axial and azimuthal variation. Reynolds stress terms quadratic in axial mode amplitudes, as well as contributions involving higher-order axial modes, are subdominant due to their weak excitation and will be ignored here. As we shall show presently, flows excited due to these modes are further suppressed by their strong radial decay away from the bubble surface.

3.1. Separable Stokes solutions

Motivated by the axial and azimuthal symmetries of the Reynolds stresses at the bubble (3.7), we seek separable solutions of (2.3). In general, the velocity field $\boldsymbol{v}$ may be decomposed into its irrotational (conservative) part $\boldsymbol{v}^{C}$ and a rotational (non-conservative) part $\boldsymbol{v}^{N}$ (Doinikov Reference Doinikov1997; Doinikov & Bouakaz Reference Doinikov and Bouakaz2010). The irrotational part of the flow has a form similar to (3.2) and is expressible in terms of the cylindrical harmonics ${\it\varphi}_{mn}$ given by (3.3), allowing $\boldsymbol{v}_{mn}^{C}\equiv \boldsymbol{{\rm\nabla}}{\it\varphi}_{mn}$ to satisfy the kinematic condition at all solid walls, but not at the bubble. For completeness, we give the expression for the components of $\boldsymbol{v}_{mn}^{C}$ in cylindrical coordinates $(r,{\it\theta},z)$ :

(3.11) $$\begin{eqnarray}\boldsymbol{v}_{mn}^{C}=\{u_{mn}^{C}(r)\cos 2m{\it\theta}\cos 2n{\it\alpha}z,v_{mn}^{C}(r)\sin 2m{\it\theta}\cos 2n{\it\alpha}z,w_{mn}^{C}(r)\cos 2m{\it\theta}\sin 2n{\it\alpha}z\},\end{eqnarray}$$

where

(3.12) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle u_{mn}^{C}(r)=-\frac{2n{\it\alpha}}{2}(\text{K}_{2m+1}(2n{\it\alpha}r)+\text{K}_{2m-1}(2n{\it\alpha}r))\\ \displaystyle v_{mn}^{C}(r)=-\frac{2m}{r}\text{K}_{2m}(2n{\it\alpha}r)\\ \displaystyle w_{mn}^{C}(r)=-2n{\it\alpha}\text{K}_{2m}(2n{\it\alpha}r).\end{array}\right\}\end{eqnarray}$$

The no-penetration condition at the bubble is satisfied by making use of the rotational, pressure-dependent solutions of the Stokes equations ( $\boldsymbol{v}^{C}$ is, by definition, pressure independent). The effective pressure $p$ itself is harmonic outside boundary layers ( ${\rm\nabla}^{2}p=0$ ) due to (2.3), and may be decomposed into cylindrical harmonics (3.3). For the dimensionless pressure mode $p_{mn}={\it\varphi}_{mn}$ , rotational solutions $\boldsymbol{v}_{mn}^{N}$ of the Stokes equations ( ${\rm\nabla}^{2}\boldsymbol{v}_{mn}^{N}=\boldsymbol{{\rm\nabla}}p_{mn}$ ) take the form

(3.13) $$\begin{eqnarray}\boldsymbol{v}_{mn}^{N}=\{u_{mn}^{N}(r)\cos 2m{\it\theta}\cos 2n{\it\alpha}z,v_{mn}^{N}(r)\sin 2m{\it\theta}\cos 2n{\it\alpha}z,w_{mn}^{N}(r)\cos 2m{\it\theta}\sin 2n{\it\alpha}z\},\end{eqnarray}$$

where

(3.14) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle u_{mn}^{N}(r)=-\frac{m}{2n{\it\alpha}}\text{K}_{2m+1}(2n{\it\alpha}r)+\frac{r}{2}\text{K}_{2m+2}(2n{\it\alpha}r)\\ \displaystyle v_{mn}^{N}(r)=\frac{m+1}{2n{\it\alpha}}\text{K}_{2m+1}(2n{\it\alpha}r)\\ \displaystyle w_{mn}^{N}(r)=\frac{r}{2}\text{K}_{2m+1}(2n{\it\alpha}r).\end{array}\right\}\end{eqnarray}$$

Neglecting any axial Stokes drift contributions to the Lagrangian flow, the corresponding Lagrangian streaming mode $\boldsymbol{v}_{mn}$ that satisfies the no-penetration condition at the interface $r=1$ and at the walls $\partial {\it\Omega}_{w}$ , normalized to have unit peak axial velocity at $\partial {\it\Omega}_{b}$ , may be written as

(3.15) $$\begin{eqnarray}\boldsymbol{v}_{mn}=\left(\frac{w_{mn}^{C}(1)}{u_{mn}^{C}(1)}-\frac{w_{mn}^{N}(1)}{u_{mn}^{N}(1)}\right)^{-1}\left(\frac{\boldsymbol{v}_{mn}^{C}}{u_{mn}^{C}(1)}-\frac{\boldsymbol{v}_{mn}^{N}}{u_{mn}^{N}(1)}\right).\end{eqnarray}$$

A general class of no-penetration Stokes solutions may be constructed using linear combinations of $\boldsymbol{v}_{mn}$ (over  $m$ and  $n$ ), determined simultaneously by boundary conditions at the bubble and at the walls. Note, however, that $\boldsymbol{v}_{mn}$ are characterized by exponential radial decays ( $\propto \text{e}^{-2n{\it\alpha}r}$ ). The dominance of the lowest axial mode in the Reynolds stress, together with the exponential radial decays of (3.15) into the fluid away from the bubble surface, ensures that the lowest (non-trivial) axial modes are both (i) the most strongly forced, and (ii) the predominant axial streaming solutions in the bulk of the flow. It is therefore appropriate to consider only the lowest axial modes ( $n=1$ ) of (3.15), so that

(3.16) $$\begin{eqnarray}\boldsymbol{v}\approx \mathop{\sum }_{m=0}^{\infty }c_{m1}\boldsymbol{v}_{m1},\end{eqnarray}$$

where the $c_{m1}$ are determined by the slip condition (3.6) at the top and bottom wall. Since (3.16) is only a truncated representation of the flow, rather than satisfying the slip condition rigorously for all $r>1$ , we will focus on satisfying it only at $r=1$ . This is an appropriate simplification since the axial solutions presented here modify the flow significantly only near $\partial {\it\Omega}_{b}$ , where the slip at walls due to the planar solution $\boldsymbol{u}$ is the greatest (Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014). At $r=1$ , the slip condition (3.6), along with (3.16) when decomposed into a Fourier series in ${\it\theta}$ , yields the coefficients

(3.17) $$\begin{eqnarray}c_{m1}=-\frac{\displaystyle \int _{0}^{{\rm\pi}}\boldsymbol{u}(1,{\it\theta})\boldsymbol{\cdot }\hat{{\bf\theta}}\sin 2m{\it\theta}\,\text{d}{\it\theta}}{\displaystyle \int _{0}^{{\rm\pi}}\boldsymbol{v}(1,{\it\theta},d/2)\boldsymbol{\cdot }\hat{{\bf\theta}}\sin 2m{\it\theta}\,\text{d}{\it\theta}},\end{eqnarray}$$

where $\hat{{\bf\theta}}$ is the unit vector in the azimuthal direction.

Note that (3.17) does not determine the strength of the axisymmetric solution $c_{01}$ , since (3.17) is specifically a statement of vanishing azimuthal velocity at the interface. It should be expected, however, that such a flow component is in general excited, in particular due to bubble oscillations that set up axisymmetric Reynolds stresses (3.10). The structure of the axisymmetric flows $(r,z)$ is entirely analogous to planar flows $(r,{\it\theta})$ in the $x{-}y$ plane with respect to the governing equations and boundary conditions, albeit weaker in overall magnitude due to weaker forcing (smaller axial oscillation amplitudes). This suggests that the axisymmetric solution also inherits the ‘fountain’ structure of the two-dimensional solution due to the pinning of contact lines at the top and bottom walls (see e.g. Wang et al. Reference Wang, Rallabandi and Hilgenfeldt2013; Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014), with radially outward velocities at the mid-plane ( $z=0$ ) and radially inward velocities near the top and bottom walls (cf. figure 2 b,c).

In principle $c_{01}$ may be evaluated from a detailed calculation of the Reynolds stresses. This however requires both (i) a resolution of the axial oscillation mode amplitudes $A_{mn}$ as functions of frequency (see e.g. Gelderblom et al. Reference Gelderblom, Zijlstra, van Wijngaarden and Prosperetti2012; Wang et al. Reference Wang, Rallabandi and Hilgenfeldt2013), and (ii) a detailed description of both oscillatory and steady no-stress boundary layers due to these flow modes, accurate up to $O({\it\delta}^{2})$ (see e.g. Davidson & Riley Reference Davidson and Riley1971; Longuet-Higgins Reference Longuet-Higgins1998; Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014). We do not attempt to carry out this analysis rigorously here, but rather estimate $c_{01}$ by considering the ratio between the damping coefficient ${\it\gamma}_{00}$ of the monopole ( $m=0,n=0$ ) and that of the lowest axial mode ( $m=0,n=1$ ), denoted by ${\it\gamma}_{01}$ ; cf. equation (21) of Wang et al. (Reference Wang, Rallabandi and Hilgenfeldt2013). In the limit of small damping, the relative mode amplitude scales inversely with the relative damping coefficient; for the channel aspect ratio used in the present experiments ( $d/a=2.5$ ), this results in $|A_{01}|/|A_{00}|\approx {\it\gamma}_{00}/{\it\gamma}_{01}\approx 0.32$ . Under the present normalization ( $|A_{00}|=1$ ), the velocity scale due to the product of the monopole and the axisymmetric quadrupole is therefore ${\sim}0.32u_{s}$ , resulting in $|c_{01}|\sim 0.32$ . The value of $c_{01}$ inferred from experimental measurements, described in more detail in § 5, agrees closely with this estimate.

For $m\geqslant 1$ , the coefficients $c_{m1}$ may be computed directly from the two-dimensional theory, using (3.17). The ‘fountain’ structure of the planar flow $\boldsymbol{u}$ (one vortex per quadrant, cf. figure 2 b), naturally leads by (3.17) to $|c_{11}|\gg |c_{m1}|$ for $m>1$ . Retaining therefore only terms of (3.16) proportional to $c_{01}$ and $c_{11}$ , the leading contributions to the three-dimensional streaming flow may be written as

(3.18) $$\begin{eqnarray}\boldsymbol{U}=\boldsymbol{u}+c_{01}\boldsymbol{v}_{01}+c_{11}\boldsymbol{v}_{11}\quad (\boldsymbol{v}=c_{01}\boldsymbol{v}_{01}+c_{11}\boldsymbol{v}_{11}).\end{eqnarray}$$

The two axial contributions $c_{01}\boldsymbol{v}_{01}$ and $c_{11}\boldsymbol{v}_{11}$ represent, respectively, the lowest axisymmetric and axial–azimuthal modes, and are depicted in figure 2(a). These contributions address, respectively, the excitation of steady flow due to axial bubble oscillations (axial oscillatory flow gradients) and the observance of no-slip conditions at axially confining walls.

4. Results and discussion

Before we discuss the resultant flow $\boldsymbol{U}$ in detail, it will prove useful to write the planar flow $\boldsymbol{u}=\{u_{r},u_{{\it\theta}},0\}$ in terms of a stream function ${\it\psi}(r,{\it\theta})$ so that $u_{r}=r^{-1}\partial _{{\it\theta}}{\it\psi}$ and $u_{{\it\theta}}=-\partial _{r}{\it\psi}$ . The steady two-dimensional flow is in general composed of a dominant fountain vortex and a sub-dominant anti-fountain vortex near the wall at $y=0$ (Wang et al. Reference Wang, Rallabandi and Hilgenfeldt2013; Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014). Here, we neglect the effect of the anti-fountain near the wall, which is typically confined to a region comparable in extent to the boundary layer thickness  ${\it\delta}$ , and is therefore greatly suppressed by the (no-slip) wall boundary layer structure. The stream function ${\it\psi}$ describing the two-dimensional Lagrangian streaming flow in a channel of finite height $h$ is in general obtained by a modification of the solution obtained in the half-space ( $h\rightarrow \infty$ ), by the procedure outlined in appendix A. Focusing on the quadrant $(x>0,y>0)$ , the two-dimensional steady flow is characterized by a single vortex system (closed streamlines) governed by ${\it\psi}\in ({\it\psi}_{m},0)$ , where ${\it\psi}=0$ delineates the boundaries of the domain (largest enclosed area), and ${\it\psi}={\it\psi}_{m}<0$ identifies the two-dimensional vortex centre (vanishing enclosed area, see figure 2 b). The fountain vortex flow of figure 2(b), with its outward flow near the symmetry plane ${\it\theta}={\rm\pi}/2$ (Wang et al. Reference Wang, Jalikop and Hilgenfeldt2012, Reference Wang, Rallabandi and Hilgenfeldt2013; Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014) has a simple far-field approximation to the stream function for large channel height ( $h\rightarrow \infty$ ), namely,

(4.1) $$\begin{eqnarray}{\it\psi}\sim \frac{4e_{1}}{r}\sin ^{2}{\it\theta}\cos {\it\theta},\quad \text{as }r\rightarrow \infty ,\end{eqnarray}$$

where $e_{1}$ is an $O(1)$ constant related to planar oscillation mode amplitudes and phases and is negative for fountain-mode streaming (Rallabandi et al. Reference Rallabandi, Wang and Hilgenfeldt2014).

The superposition of the two-dimensional Lagrangian streaming with the flows computed in the previous section yields a three-dimensional steady flow that satisfies (i) the Stokes equations, (ii) no-penetration conditions at all boundaries and (iii) no-slip conditions at the confining walls close to the interface. We will focus on material trajectories in this flow, given by solutions $\boldsymbol{x}(t)$ of $\dot{\boldsymbol{x}}(t)=\boldsymbol{U}(\boldsymbol{x}(t))$ , which are integrated numerically using a fourth-order Runge–Kutta scheme. Due to the symmetry of the problem, we focus our results on the region $x>0,0<y<h,0<z<d/2$ . The smallest time scale relevant to the steady flow is the time scale of fluid advection along the bubble interface, which we define as $T_{s}\equiv a/u_{s}=1/({\it\epsilon}^{2}{\it\omega})$ , and use to non-dimensionalize time for all the results following below.

The two-dimensional vortex motion occurs over the nominal orbit time along closed two-dimensional streamlines (under the action of $\boldsymbol{u}$ alone), which is much greater than $T_{s}$ due to the algebraic decay of the flow velocity away from the bubble. In units of $T_{s}$ , the dimensionless planar orbit time on a streamline $T({\it\psi})$ is in general defined by

(4.2) $$\begin{eqnarray}T({\it\psi})=\oint _{{\it\psi}}\frac{1}{|\boldsymbol{u}|}\,\text{d}s=\oint _{{\it\psi}}\frac{1}{|\boldsymbol{{\rm\nabla}}{\it\psi}|}\,\text{d}s,\end{eqnarray}$$

where $\text{d}s=|\boldsymbol{u}|\,\text{d}t$ is a differential arclength element along the streamline and the path of integration along ${\it\psi}$ is counterclockwise. Using the far-field solution (4.1), one may estimate $T({\it\psi})$ as

(4.3) $$\begin{eqnarray}T_{h\rightarrow \infty }({\it\psi})\approx \int _{{\rm\pi}/2}^{0}\frac{r({\it\psi},{\it\theta})}{u_{{\it\theta}}({\it\psi},{\it\theta})}\,\text{d}{\it\theta}=-\frac{{\rm\pi}e_{1}^{2}}{2{\it\psi}^{3}}.\end{eqnarray}$$

This estimate, valid for $h\rightarrow \infty$ , becomes quantitatively modified for flow fields with finite $h$ (see appendix A), but the order of magnitude of $T({\it\psi})\sim 10T_{s}{-}100T_{s}$ obtained from (4.3) remains valid.

Figure 3. (a) Axial motion $z(t)$ of a typical fluid particle trajectory, showing the strong bursts of axial displacement (when the trajectory is close to the bubble surface) and the net axial displacement ${\rm\Delta}z$ over one cycle of radial motion. (b) Fluid particle trajectories for two different initial conditions, showing the nested torus structure: light red is a trajectory of large axial extent, while black indicates a trajectory close to the ‘core’ of the torus. The solid black line at $(x=1,y=0)$ indicates the bubble contact line with the bounding wall, dashed lines indicate the locations of symmetry planes intersecting the bubble surface at $(z=0,r=1)$ and ( $x=0,y=1$ ).

4.1. Three-dimensional effects on short and long time scales

The axial flows modify the two-dimensional motion in at least two significant ways. Perhaps most strikingly, a fluid particle trajectory close to the bubble ( $r\lesssim 1+d/(2{\rm\pi})$ ) shows strong bursts of axial displacement on the short time scale ${\sim}T_{s}$ , much shorter than the two-dimensional orbit time. Figure 3(a) shows an example of such a trajectory: as axial flow components of $\boldsymbol{v}$ very close to the bubble are in general comparable to planar velocity components due to the coupling (3.6), the fluid element approaches the bubble and is advected along the interface (in the direction ${\it\theta}\rightarrow {\rm\pi}/2$ ), suffering $O(1)$ axial displacements over the time scale  $T_{s}$ . The axial confinement necessitates a spatial reversal of the axial velocity components, so that a fluid element experiences axial displacements both towards and away from the mid-plane during the course of its motion near the bubble. As a consequence, the net axial displacement ( ${\rm\Delta}z$ ) of a fluid element over the intermediate time scales $T({\it\psi})$ , averaged over its nominal two-dimensional orbit, is typically small compared to the bubble radius (figure 3 a).

Figure 3(b) illustrates the long-time effect of these net displacements: similar to the axial ${\rm\Delta}z$ , planar components of $\boldsymbol{v}$ result in small deviations of the planar projections of fluid trajectories from nominally closed two-dimensional orbits. If the depth of the channel is on the order of the bubble radius (or greater), the small axial and planar drifts result in a motion over time scales longer still than $T({\it\psi})$ , over which a fluid element explores a significant fraction of the half-channel by systematically sampling a succession of planar orbits of varying planar extents (varying  ${\it\psi}$ ) at different depths $z$ (cf. figure 8). The superimposed axial and radial motions result in quasi-periodic trajectories that fill toroidal structures in three dimensions (figure 3 b). It can be seen that fluid elements at different initial positions follow toroidal trajectories of considerably varying spatial extents, none of which intersect. At the core of this nested system of tori is a singly periodic trajectory of vanishing toroidal volume.

The long-time toroidal motion of a fluid element is further illustrated in figure 4(a): its quasi-period ${\it\tau}$ depends on the initial position of the fluid particle, with net motion towards the mid-plane $z=0$ occurring over the widest planar orbits ( ${\it\psi}\sim 0$ ), followed by motion towards the wall ( $z=d/2$ ) occurring over ‘tight’ planar orbits that are closer to the two-dimensional vortex centre ( ${\it\psi}\sim {\it\psi}_{m}$ ). Reversals in the direction of the slow axial motion occur near the wall and the mid-plane, involving relatively large changes in the planar extent of the quasi-planar orbits, i.e. a widening of tight orbits near the wall, and a tightening of wide orbits close to the mid-plane (see figure 4). As discussed in § 5, all elements of this motion are observed in experiment as well (figure 4 b).

Figure 4. Typical toroidal trajectory, as (a) predicted by the theory for a fluid element, and (b) measured experimentally for an approximately passive tracer by APTV, cf. Marin et al. (Reference Marin, Rossi, Rallabandi, Wang, Hilgenfeldt and Kähler2015). The motion is towards the channel walls on the tight orbits (close to the nominal two-dimensional vortex centre), indicated in red, and towards the channel mid-plane during wider planar orbits, coloured grey. Arrows indicate these axial motions as well as the two-dimensional fountain orientation (blue).

4.2. Generalized Hamiltonian theory of three-dimensional effects

The wide spatial range explored by different material points underscores the importance of understanding the kinematics of these trajectories, particularly in the context of applications such as fluid mixing and manipulation of microparticles in these flows.

By definition, the displacement in the $z$ coordinate of a fluid element is simply the time integral of the axial component of the velocity field as a function of its instantaneous position $\boldsymbol{x}(t)$ . However, the drift in the $z$ coordinate of a fluid element over a time $T({\it\psi})$ , denoted by ${\rm\Delta}z$ , is typically much smaller than the channel depth $d$ (see figure 3 a), with significant axial motion only occurring over much longer time scales. In the limit ${\rm\Delta}z/d\ll 1$ (or otherwise $|\boldsymbol{v}|\ll |\boldsymbol{u}|$ ), the motion of the particle over a time $T({\it\psi})$ is therefore approximately planar to leading order in ${\rm\Delta}z$ , i.e.  $\boldsymbol{x}(t)\approx \tilde{\boldsymbol{x}}(t)=\{r(t),{\it\theta}(t),z\}$ such that $z$ and ${\it\psi}={\it\psi}(r(t),{\it\theta}(t))$ are approximately constant over the motion. For ${\rm\Delta}z/d\ll 1$ , one may therefore compute ${\rm\Delta}z$ as the integral axial displacement of a fluid element as it traverses a nominally closed two-dimensional orbit (constant  $z$ and  ${\it\psi}$ ), over the time $T({\it\psi})$ . To first approximation, we may write

(4.4) $$\begin{eqnarray}{\rm\Delta}z({\it\psi},z)=\int _{0}^{T({\it\psi})}v_{z}(\boldsymbol{x}(t))\,\text{d}t\approx \int _{0}^{T({\it\psi})}v_{z}(\tilde{\boldsymbol{x}}(t))\,\text{d}t.\end{eqnarray}$$

Note that the absolute values of the limits in the time integral are unimportant as long as their difference is $T({\it\psi})$ – this ensures an integral over the entire closed orbit  ${\it\psi}$ , valid for all material points on  ${\it\psi}$ . The ${\rm\Delta}z$ defined by (4.4) therefore represents the deviation of a material orbit from planarity over the time $T({\it\psi})$ .

The deviation from closedness of orbits over a time $T({\it\psi})$ , when projected on the $x{-}y$ plane, may be similarly quantified by a drift in ${\it\psi}$ (denoted by ${\rm\Delta}{\it\psi}$ ) of a fluid element. The time rate of change of the stream function value ${\it\psi}$ sampled by a fluid element as it is advected by the three-dimensional flow field $\boldsymbol{U}$ is simply the material derivative of ${\it\psi}$ , which for steady flow is given by $\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi}=\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi}$ . If we again assume small deviations from closedness over a planar orbit time scale, ${\rm\Delta}{\it\psi}$ over the time $T({\it\psi})$ is computed to first approximation as

(4.5) $$\begin{eqnarray}{\rm\Delta}{\it\psi}({\it\psi},z)=\int _{0}^{T({\it\psi})}\boldsymbol{v}(\boldsymbol{x}(t))\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi}(\boldsymbol{x}(t))\,\text{d}t\approx \int _{0}^{T({\it\psi})}\boldsymbol{v}(\tilde{\boldsymbol{x}}(t))\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\psi}(\tilde{\boldsymbol{x}}(t))\,\text{d}t.\end{eqnarray}$$

It may be shown (see appendix B) that for an arbitrary incompressible velocity field $\boldsymbol{v}$ , the approximate drifts ${\rm\Delta}z$ and ${\rm\Delta}{\it\psi}$ (using $\tilde{\boldsymbol{x}}(t)$ ) defined respectively by (4.4) and (4.5), identically satisfy the differential continuity equation

(4.6) $$\begin{eqnarray}\displaystyle \frac{\partial {\rm\Delta}{\it\psi}}{\partial {\it\psi}}+\frac{\partial {\rm\Delta}z}{\partial z}=0.\end{eqnarray}$$

This automatically allows them to be expressed in terms of a Hamiltonian $\mathscr{H}({\it\psi},z)$ as

(4.7a,b ) $$\begin{eqnarray}{\rm\Delta}{\it\psi}=\frac{\partial \mathscr{H}}{\partial z},\quad {\rm\Delta}z=-\frac{\partial \mathscr{H}}{\partial {\it\psi}},\end{eqnarray}$$

where $\mathscr{H}({\it\psi},z)$ may be written as

(4.8) $$\begin{eqnarray}\mathscr{H}({\it\psi},z)\equiv -\int \oint _{{\it\psi}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}s\,\text{d}z,\end{eqnarray}$$

and $\boldsymbol{n}\equiv -\boldsymbol{{\rm\nabla}}{\it\psi}/|\boldsymbol{{\rm\nabla}}{\it\psi}|$ is the local unit normal to streamlines in the $x{-}y$ plane; see appendix B for a more detailed discussion. Since these drifts occur over a planar orbit time $T({\it\psi})$ , one may define mean drift velocities as

(4.9a,b ) $$\begin{eqnarray}\overline{v}_{{\it\psi}}({\it\psi},z)\equiv \frac{1}{T}\frac{\partial \mathscr{H}}{\partial z},\quad \overline{v}_{z}({\it\psi},z)\equiv -\frac{1}{T}\frac{\partial \mathscr{H}}{\partial {\it\psi}}.\end{eqnarray}$$

The system (4.9) constitutes precisely an incompressible two-dimensional flow field in the abstract ${\it\psi}{-}z$ space, governed by the Hamiltonian $\mathscr{H}({\it\psi},z)$ , which plays the role of a stream function and remains invariant over the long-time three-dimensional steady motion of a fluid element. Since both $\boldsymbol{u}$ and $\boldsymbol{v}$ involve vanishing normal velocities on $\partial {\it\Omega}$ , it follows from (4.8) and (4.9) that the flow in the abstract space does not penetrate its boundaries, i.e.

(4.10) $$\begin{eqnarray}\overline{v}_{{\it\psi}}({\it\psi}_{m},z)=\overline{v}_{{\it\psi}}(0,z)=\overline{v}_{z}({\it\psi},0)=\overline{v}_{z}({\it\psi},d/2)=0.\end{eqnarray}$$

This automatically ensures the existence of at least one elliptic point ( ${\it\psi}_{c},z_{c}$ ) at which $\partial _{{\it\psi}}\mathscr{H}=\partial _{z}\mathscr{H}=0$ , which represents a two-dimensional curve in real space given by ${\it\psi}(r,{\it\theta})={\it\psi}_{c},z=z_{c}$ . Any material point on this curve, by virtue of lying on a (stable) stagnation point of the ${\it\psi}{-}z$ flow, remains confined to it under the action of the three-dimensional flow field, i.e. the curve in real space defined by $({\it\psi}_{c},z_{c})$ approximates the core of the system of nested toroidal fluid trajectories (cf. figure 3 b). Curves of constant $\mathscr{H}$ in the vicinity of the core are concentric about it, spanning finite ranges of axial coordinates $z$ and streamlines ${\it\psi}$ , and therefore delineate surfaces of the tori on which fluid elements move in real space. Figure 5(a) shows the ${\it\psi}{-}z$ phase portrait that provides this simplified description of the three-dimensional motion. It is noteworthy that this description assumes no particular forms for $\boldsymbol{u}$ and $\boldsymbol{v}$ , only requiring that the two-dimensional flow is steady and characterized by closed streamlines, and is therefore in principle applicable to any perturbed two-dimensional flow under confinement.

Figure 5. (a) Phase portrait of the Hamiltonian system in ${\it\psi}{-}z$ space, computed for $h\rightarrow \infty$ . Closed dashed lines are curves of constant $\mathscr{H}$ representing three-dimensional tori in real space, with the light red and black solid-line trajectories corresponding to the respective tori in figure 3(b). (b) Functions ${\it\eta}({\it\psi})$ at $f=20~\text{kHz}$ and $d=2.5$ for different channel heights  $h$ .

For $\boldsymbol{v}$ given by (3.18) (single axial wavenumber), the Hamiltonian is separable and takes the form

(4.11) $$\begin{eqnarray}\mathscr{H}={\it\eta}({\it\psi})\sin 2{\it\alpha}z,\end{eqnarray}$$

so that

(4.12a,b ) $$\begin{eqnarray}\overline{v}_{{\it\psi}}=\frac{2{\it\alpha}}{T({\it\psi})}{\it\eta}({\it\psi})\cos 2{\it\alpha}z,\quad \overline{v}_{z}=-\frac{1}{T({\it\psi})}{\it\eta}^{\prime }({\it\psi})\sin 2{\it\alpha}z.\end{eqnarray}$$

The core trajectory is therefore located at $z_{c}=d/4$ and ${\it\eta}^{\prime }({\it\psi}_{c})=0$ ; note that ${\it\eta}^{\prime }({\it\psi})$ is guaranteed to have at least one zero within the relevant range of stream function values $({\it\psi}_{m},0)$ . The temporal evolution of the stream function is therefore given by

(4.13) $$\begin{eqnarray}\frac{\text{d}{\it\psi}}{\text{d}t}=\pm \frac{2{\it\alpha}}{T({\it\psi})}\sqrt{{\it\eta}^{2}-\mathscr{ H}^{2}}.\end{eqnarray}$$

The characteristic function ${\it\eta}({\it\psi})$ itself is obtained from evaluating (4.8) and (4.11), and displayed in figure 5(b). Its shape is relatively insensitive to the channel width for $h$ not too small.

Equation (4.13) may be solved numerically for given initial conditions $(z_{0},{\it\psi}_{0})$ . For trajectories that approach the core, $\mathscr{H}\sim \mathscr{H}({\it\psi}_{c},z_{c})$ , orbits in the ${\it\psi}{-}z$ plane are elliptical so that ${\it\psi}(t)$ and $z(t)$ both vary sinusoidally with a single frequency. For arbitrary $\mathscr{H}$ , the temporal motion is asymmetric over a three-dimensional period ${\it\tau}$ due to the dependence of the drift velocities on $1/T({\it\psi})$ , so that the axial dynamics are faster for smaller orbits ${\it\psi}/{\it\psi}_{m}\sim 1$ , and slower for larger planar orbits ${\it\psi}/{\it\psi}_{m}\sim 0$ . In real space, this corresponds to the faster axial motion towards the top wall ( $z=d/2$ ) over tight planar orbits, and the slower axial drift towards the mid-plane ( $z=0$ ) over wider orbits, as depicted in figure 6(a), where the red line indicating the Hamiltonian theory captures the asymmetry of the long-time motion. The ${\it\psi}$ -motion depicted in figure 6(b) is much more symmetric, but is also subject to abrupt jumps on the short time scale.

Figure 6. (a) Axial coordinate $z$ and (b) normalized stream function ${\it\psi}/{\it\psi}_{m}$ as functions of time, over approximately one three-dimensional quasi-period  ${\it\tau}$ . Blue lines represent the numerical integration of particle trajectories, while red lines show the predictions of the Hamiltonian theory, which captures the long-time kinematics while averaging over the shorter time scales.

5. Experimental data for three-dimensional streaming flows

A microbubble streaming set-up of the type depicted in figure 1(a) was used to obtain three-dimensional trajectory data for approximately passive tracer particles (polystyrene beads of diameter $2~{\rm\mu}\text{m}$ ) utilizing APTV (Cierpka et al. Reference Cierpka, Rossi, Segura and Kähler2011; Cierpka & Kähler Reference Cierpka and Kähler2012; Rossi & Kähler Reference Rossi and Kähler2014). In this optical imaging technique, the image of a spherical tracer particle is distorted by a cylindrical lens into elliptic shapes whose eccentricity depends on the depth coordinate of the particle along the optical axis. This allows the quantification of all three velocity components of moving tracer particles in a three-dimensional domain with a single camera. The experiments were conducted using $a=40~{\rm\mu}\text{m}$ , $f=20~\text{kHz}$ , $d=2.5$ , and $h=25$ , with particle images taken by a high-speed camera connected to an inverted microscope, at a recording speed of 500 f.p.s. The measurement uncertainty of this technique arises primarily from errors in particle detection and tracking, and the degree to which the elliptical shapes of the particle images can be accurately determined (see Rossi & Kähler Reference Rossi and Kähler2014 for details). For the present experiments, the measurement uncertainty in the position of a particle is ${\approx}\pm 0.1~{\rm\mu}\text{m}$ in the $x$ and $y$ directions and ${\approx}\pm 1~{\rm\mu}\text{m}$ in the $z$ direction, much smaller than typical length scales of interest ( ${\sim}10~{\rm\mu}\text{m}$ ). More details about this use of the technique in a microfluidic context and more experimental data on bubble microstreaming are presented in Marin et al. (Reference Marin, Rossi, Rallabandi, Wang, Hilgenfeldt and Kähler2015). Figure 4(b) displays an example trajectory that shows all the qualitative features predicted by the three-dimensional streaming theory: one trajectory typically occupies half the channel (between one wall and the mid-plane) and consists of approximately two-dimensional trajectory loops whose radial extent is small as the particle is displaced towards the wall, then widens to a greater radial extent as the $z$ -motion reverses and the particle moves towards the mid-plane.

In order to model the axial motion quantitatively, we obtain the value of $c_{01}$ from the experimentally measured azimuthal profile of the axial velocity component $U_{z}({\it\theta})$ near the bubble ( $r\approx 1$ ), which determines the ratio $c_{01}/c_{11}$ . This yields $c_{01}$ , since $c_{11}\approx -1.41$ is given automatically by (3.17) and the theoretical $\boldsymbol{u}$ , thereby fixing the shape of the flow. The value of $c_{01}\approx -0.34$ obtained by this procedure confirms that the axisymmetric component flow in the experiment (i) has the ‘fountain’ orientation ( $c_{01}<0$ ) and (ii) is weaker than the planar component, with a magnitude that agrees well with the theoretical estimate (using damping coefficients) of $|c_{01}|\sim 0.32$ (cf. § 3). In order to non-dimensionalize the experimental times, we first extrapolate the experimentally measured velocity onto the bubble surface, which yields a maximum slip velocity $u_{\mathit{max}}$ at the interface, which is then translated into an experimental oscillation amplitude ${\it\epsilon}$ using the prediction (from two-dimensional theory) of $u_{\mathit{max}}/({\it\epsilon}^{2}a{\it\omega})\approx 1.54$ . The computed value of ${\it\epsilon}$ determines the time scale $T_{s}=1/({\it\epsilon}^{2}{\it\omega})$ , which is used to scale time in the experiments for a comparison with theory. It should be stressed that, therefore, all theoretical figures and computations in the present work (unless explicitly stated otherwise) directly use parameters employed in practical bubble microstreaming experiments: $d=2.5,h=6.25$ for the channel geometry, $a=40~{\rm\mu}\text{m}$ , $f=20~\text{kHz}$ , and the appropriate $c_{01}$ inferred above.

Figure 7. (a) Axial motion of a particle over short time scales in the vicinity of the bubble, predicted by the theory (solid line) and measured experimentally (markers). Note the large $O(1)$ axial excursion of the particle on the shortest time scales ${\sim}T_{s}$ as it passes close to the bubble, which ultimately results in a smaller net axial displacement ${\rm\Delta}z$ , informing the long-time three-dimensional motion. (b) Projection of a part of a three-dimensional particle trajectory in the $x{-}y$ plane computed from $\boldsymbol{U}$ (red), showing the radial drift as a systematic spiralling deviating from the two-dimensional periodic orbits (blue), with deviations concentrated at locations close to the bubble surface.

We show here that the experiment confirms the theoretical results on all three distinct time scales: (i) the fast time scale ${\sim}T_{s}$ of strong axial displacement; (ii) the intermediate time scale ${\sim}10T_{s}{-}100T_{s}$ of the approximately two-dimensional vortex motion; and (iii) the very long time scale ${\it\tau}\sim 1000T_{s}$ involving the toroidal motion with systematic variation of the axial position of the particle, accompanied by a slow variation of the extent of the planar motion. Figure 7(a) demonstrates that the strong, fast displacements along the $z$ -axis are temporally resolved in experiments and in very good agreement with the theoretical prediction, including the characteristic asymmetry between the positive and negative portions of the $z$ -displacement motion, which ultimately leads to the drifts ${\rm\Delta}z$ . Figure 7(b) shows that the motion projected into the $x{-}y$ plane still reproduces the ‘typical’ two-dimensional fountain dynamics. Limitations in the experimental time resolution prevent us from verifying directly the predicted deviations of the projected three-dimensional trajectories (red) near the bubble from the streamlines resulting from the two-dimensional theory (blue). It is interesting to note, though, that the $xy$ -projection of a tracer particle on such a trajectory shows a spiralling motion towards the two-dimensional vortex centre, an observation commonly made in cylinder streaming (Lutz, Chen & Schwartz Reference Lutz, Chen and Schwartz2005; Lutz et al. Reference Lutz, Chen and Schwartz2006) and ascribed to inertial forces on the particles (as supported by numerical computations, cf. Chong et al. Reference Chong, Kelly, Smith and Eldredge2013). Finally, figure 8 confirms that the pattern of motion over long time scales ${\it\tau}\sim 1000T_{s}$ is also modelled successfully: both the pattern of slow displacement along the $z$ -axis and the pattern of widening and tightening radial extent of trajectories are confirmed in experiments, as is the strong correlation between the two motion patterns. Figure 8 thus allows a quantitative verification of the toroidal motion pattern of figure 4(b).

Figure 8. Axial and radial positions of a fluid element determined by numerical integration (a) and of an approximately passive tracer in experiment (b). Both show that the long-time three-dimensional motion is organized into simpler two-dimensional orbits (constant  $z$ ). The slow axial motion is towards the wall $z=d/2$ on tight orbits (small maximum $r$ in regions of constant  $z$ ) and in the opposite direction for wide orbits (large maximum $r$ in regions of constant  $z$ ), quantifying the toroidal motion seen in figure 4(b).

Prior to these new APTV experiments, there has only been one report of three-dimensional boundary streaming flow, for a cylinder oscillated along a diameter (Lutz et al. Reference Lutz, Chen and Schwartz2005), cf. figure 9(a). The patterns of streaming along the cylinder axis show significant differences to the bubble streaming, but it must be realized that even a strictly two-dimensional theory of cylinder streaming (Holtsmark et al. Reference Holtsmark, Johnsen, Sikkeland and Skavlem1954; Wang Reference Wang1968; Bertelsen, Svardal & Tjøtta Reference Bertelsen, Svardal and Tjøtta1973) reveals very different patterns from two-dimensional bubble streaming. In particular, a ‘DC boundary layer’ is typically present, i.e. every quadrant of the streaming flow in the $x{-}y$ plane has two radially stacked vortices. Thanks to the general nature of our approach to describing three-dimensional flows as a perturbation of the two-dimensional patterns, however, it is possible to obtain a three-dimensional solution for cylinder streaming as well. This is helped by the absence of axial oscillations of the (solid) cylinder, so that $c_{01}$ is identically zero. The only flow mode magnitude $c_{11}$ is determined again by (3.17), with $\boldsymbol{u}$ now representing the steady streaming flow outside the no-slip cylinder boundary layer, as given by a leading-order expansion in powers of ${\it\delta}$ of the solution derived by Bertelsen et al. (Reference Bertelsen, Svardal and Tjøtta1973).

The qualitative experimental image of figure 9(a) from Lutz et al. (Reference Lutz, Chen and Schwartz2005) does not contain three-dimensional information, but is a projection with finite depth of field showing multiple flow lines. Likewise, figure 9(b) shows characteristic flow lines projected onto the $x{-}z$ plane ( $x$  being the axis of oscillation of the cylinder), the same plane as in figure 9(a). Streamlines point away from the bubble and towards the walls outside the nominal ‘DC boundary layer’, whose radial extent is marked by the ‘box’ outlines in the experimental and theoretical figures. Close to the mid-plane, the flow appears two-dimensional, with axial motion being restricted to regions close to the wall; these flow features are in qualitative agreement with the results of Lutz et al. (Reference Lutz, Chen and Schwartz2005) and demonstrate the wider applicability of our formalism to three-dimensional streaming flows.

Figure 9. (a) Fluid trajectories under three-dimensional streaming due to an oscillating solid cylinder under axial confinement, projected in a plane containing the cylinder axis $z$ and the axis of oscillation  $x$ . The right-hand part of the figure sketches the direction of motion of some tracer particles (reprinted with permission from Lutz et al. (Reference Lutz, Chen and Schwartz2005), AIP Publishing LLC). (b) Analogous figure utilizing the present theory, showing qualitative agreement with the flow directions and the flow organization, including the ‘DC boundary layer’ region indicated by boxed outlines in experiment and theory. The length of the cylinder axis relative to its radius in both the experimental image of Lutz et al. (Reference Lutz, Chen and Schwartz2005) and the present theory is $d=3$ .