Controlled mean dynamics of red blood cells by two orthogonal ultrasonic standing waves

Abstract

A computational model is developed to study the time-averaged mean dynamics of red blood cells (RBCs) driven by the time-averaged mean stress generated by two phase-shifted orthogonal ultrasonic standing waves in a viscous fluid. The cell is modelled as an ellipsoidal viscoelastic membrane enclosing the viscous fluid cytoplasm, the motion of which is described by the inclination angle of the ellipsoidal cell shape and the phase angle of the potential membrane cycle. Based on the acoustic perturbation method, the acoustic field and acoustic streaming field are solved to obtain the time-averaged mean stress, and then the temporal evolution equations of the inclination and phase angles of the cell are determined considering the torque balance and energy conservation. At a small acoustic pressure amplitude, this model reproduces the experimentally observed features of cell motion in orthogonal standing waves: the transition from steady stationary orientation to unsteady tumbling with the increase of the phase difference between the two standing waves. By turning up the acoustic pressure amplitude above a critical value, it is further predicted that the previously observed motions can be accompanied by the membrane tank-treading rotation. Observations of these motions, combined with the present computational model, can help to evaluate the mechanical properties of RBC membranes in an automated and high-throughput manner by acoustic methods.

1. Introduction

The viscoelastic properties of red blood cells (RBCs) are fundamentally important for maintaining cell functions. These properties are critical for microvascular function and can be altered in various blood-related diseases, such as cardiovascular diseases (Fornal et al. 2008) and malaria (Suresh et al. 2005). The structure of RBCs is composed of an internal homogeneous fluid and its enclosing viscoelastic biological membrane, which plays a major role in the mechanical behaviour of cells. In this regard, many techniques have been developed to measure the viscoelastic modulus of the RBC membrane.

Previous studies have shown that the viscoelastic moduli of an RBC membrane can be assessed based on their mechanical behaviour in response to quantified forces. This can be achieved through traditional experimental techniques, such as micropipette aspiration (Henriksen & Ipsen 2004; Daza et al. 2019), optical tweezers (Dao, Lim & Suresh 2003; Mills et al. 2004) and atomic force microscopy (Haase & Pelling 2015; Efremov et al. 2017). These techniques provide precise measurements but have a slow test speed of approximately one cell per minute. In recent years, a variety of microfluidic techniques have been proposed to improve detection throughput. Typical microfluidic techniques use hydrodynamic forces to deform cells as they pass through microfluidic channels (Gossett et al. 2012; Otto et al. 2015; Prado et al. 2015). While these techniques offer high throughput, such as real-time deformability cytometry (RT-DC) (Otto et al. 2015) which can perform thousands of single-cell deformability analyses in minutes, the drawback is that deforming cells at high rates in a short period of time may alter the microstructure of cells and affect their mechanical response (Urbanska et al. 2018). Furthermore, unlike traditional techniques, which deform cells in situ, typical microfluidic techniques deform cells in continuous flow. This makes it difficult for typical microfluidic techniques to identify each individual cell before and after detecting cell deformation, which may limit their multifunctional integration. An alternative to achieve reasonable throughput and deform cells in situ is acoustic microfluidics (i.e. the fusion of acoustics and microfluidics). Specifically, acoustic microfluidics can trap cells in well-designed acoustic potential wells, and high throughput can be easily achieved through parallel cell manipulation (Silva et al. 2019; Xie, Bachman & Huang 2019).

In typical acoustic microfluidic devices, standing waves are generated within a microfluidic channel or cavity. The particles immersed in it experience not only time-harmonic acoustic pressure directly caused by acoustic excitation, but also time-averaged forces caused by the acoustic nonlinear effect (Xin & Lu 2016; Drinkwater 2020). This force is called the acoustic radiation force, which can trap biological cells at the acoustic pressure nodes or antinodes of standing waves depending on the acoustophoretic contrast between the cells and the host medium (Settnes & Bruus 2012; Baasch, Qiu & Laurell 2022). In general, the cells in water-based host media have a positive acoustophoretic contrast and are directed towards the acoustic pressure node (Hartono et al. 2011; Li et al. 2015). For two-dimensional (2-D) standing waves consisting of two orthogonal standing waves, the cells can be patterned at grid-like acoustic pressure nodes, so operations in 2-D standing waves facilitate parallel cell manipulation at the single-cell level (Collins et al. 2015). Bernard et al. (2017) experimentally investigated the rigid body tumbling and translation of particles and biological cells (including red and white blood cells) in two orthogonal ultrasonic standing waves with phase difference $\zeta $. In their experiments, they found that the cells were in a steady stationary state when the phase difference $\zeta = 0$, while they execute a tumbling motion when the phase difference $\zeta = {\rm \pi}/2$. Increasing the phase difference $\zeta $ from $0$ to ${\rm \pi} /2$, the cells were observed to undergo a transition from the steady stationary state to tumbling motion. These features are generally similar to those of RBCs in shear flow. Specifically, the shear flow can be divided into an extensional component and a rotational component, which lead to steady stationary orientation and unstable tumbling of RBCs immersed in it, respectively (Vlahovska, Podgorski & Misbah 2009; Sinha & Thaokar 2018). The deformable RBCs in shear flows have been studied extensively and are known to exhibit intriguing dynamical behaviour (Rezghi & Zhang 2022). In shear flows with high shear rate or external fluid viscosity, the RBC membrane circulates around the deformed cell contour, which drives a vortex-like circular flow in the internal cytoplasm. This motion is called a tank-treading motion, and it provides researchers with a convenient way to measure the mechanical properties of RBCs (Tsubota 2021; Rezghi & Zhang 2022). It is worth noting that the work of Bernard et al. (2017) is limited to small acoustic inputs, where RBCs behave like rigid bodies. Further studies treating RBCs as elastic bodies under relatively large acoustic input are required to understand the complete dynamics of cells in 2-D standing waves and to explore whether 2-D standing waves can excite the tank-treading motion of RBCs like shear flow.

Over the years, some theoretical work with simplified models has been developed to analyse RBC dynamics in linear shear flow. Keller & Skalak (1982) proposed the first simplified theoretical model in this field. Their model treated the RBC as a fluid ellipsoid with a fixed shape, with a membrane that allows circulation around the ellipsoid shape, and the axis of symmetry of the ellipsoid was constrained in the shear plane. The degree of freedom describing RBC motion was reduced to two: that describing the orientation of the ellipsoid in space and that describing the circulation of the membrane relative to the ellipsoidal shape. The evolution equations of these two degrees of freedom were established from the perspective of torque balance and energy conservation. Abkarian, Faivre & Viallat (2007) further developed this model to explain the shear elasticity of a cell membrane. The enriched model achieved a quantitative agreement with experimental observations. Efforts have been devoted to including the stress-free shape of the cell membrane (Dupire, Abkarian & Viallat 2015), cell shape deformation (Noguchi 2010), three-dimensional effect (i.e. the axisymmetric axis of the cell is not in the shear plane) (Mendez & Abkarian 2018; Mignon & Mendez 2021) and general linear flow beyond shear flow (Ye et al. 2016).

To analyse cell dynamics in an acoustic field, an important issue is the theoretical treatment of the acoustically induced time-averaged force. Mathematically, the time-averaged force arises from the nonlinearity of the compressible Navier–Stokes (N-S) equations describing fluid dynamics (Baudoin & Thomas 2020). The acoustic perturbation method has been widely used in the acoustic microfluidic community to solve the compressible N-S equations (Muller et al. 2012; Nama et al. 2015). Expanding the fluid variables, the first-order variables represent the time-harmonic acoustic response, and the second-order variables after time averaging represent the nonlinear acoustic phenomenon, including the acoustic streaming, which describes the steady swirling fluid motion generated by acoustic dissipation. Several works have analytically solved the acoustic streaming around spherical particles located at the pressure node of 2-D standing waves, and have shown that the acoustic streaming exerts a torque on the particles (Lee & Wang 1989; Rednikov, Riley & Sadhal 2003; Zhang & Marston 2014). The torque is found to have a phase difference dependent factor $\sin \zeta $, which means that for phase difference $\zeta = 0$, the torque is zero, while for phase difference $\zeta = {\rm \pi}/2$, the torque is maximum. Due to the geometric complexity, the analytical expression of torque on a non-spherical particle has not been reported. Instead, a three-dimensional finite element model has been developed to calculate the torque on spherical and non-spherical particles (Hahn, Lamprecht & Dual 2016). Previous work mainly focused on the manipulation of rigid particles in ultrasonic standing waves, so there is still a lack of research on the manipulation of deformable cells in two-dimensional ultrasonic standing waves.

Recently, we have developed a finite element model to analyse cell dynamics in acoustic fields (Liu & Xin 2023a) and investigated the shape dynamics of two-dimensional capsules (an elastic membrane enclosing a viscous fluid) in two-dimensional standing waves (Liu & Xin 2023b). The complex dynamics of two-dimensional capsules, including steady stationary state, tumbling and tank treading, are identified for different phase differences and acoustic pressure amplitudes of the two-dimensional standing waves. However, the findings of this study have not been generalized to three-dimensional dynamics due to the huge computational effort. In this work, a three-dimensional model for red blood cell (RBC) dynamics in two phase-shifted orthogonal ultrasonic standing waves is proposed. The cell is modelled as an ellipsoidal viscoelastic membrane enclosing the viscous fluid cytoplasm. Applying the acoustic perturbation method, the Navier–Stokes equation is divided into two sets of equations: the first-order equations for ultrasonic propagation and the time-averaged second-order equations for the mean dynamics. The two sets of equations are numerically solved to calculate the acoustic-induced mean stress that drives cell motion through the finite element method. Based on the torque balance and energy conservation, the governing equations of cell motion are established and numerically solved by applying the fourth-order Runge–Kutta method. The conditions that determine the transition of different cell dynamic states are identified and a phase diagram of different cell dynamic states as a function of phase difference and acoustic pressure amplitude is constructed. This computational model provides a comprehensive understanding of erythrocyte dynamics in ultrasonic standing waves and not only reproduces previous experimental observations, but also predicts possible tank-treading motion beyond previous experimental conditions.

2. Theoretical model

The time-averaged mean dynamics of RBCs driven by two orthogonal ultrasonic standing waves in a viscous fluid medium is investigated. The RBC is composed of a membrane enclosing an homogeneous fluid. It is assumed that the RBC can be represented as an oblate ellipsoid of a prescribed shape that never changes during the motion. Thus, the present work is valid when the RBC deformation remains small. Nevertheless, the membrane of the cell is allowed to circulate around the fixed ellipsoidal shape. The consideration of an oblate ellipsoid shape facilitates theoretical modelling and can well reproduce the dynamic behaviour of real biconcave erythrocytes in fluid experiments (Abkarian et al. 2007). Two-dimensional ultrasonic standing waves are generated by the oscillation of two pairs of orthogonal piezoelectric transducers (PZTs), as shown in figure 1(a). When the two-dimensional standing waves are established, the cells are subjected to an acoustic radiation force, trapping the cells in the grid-like acoustic pressure nodes. Due to the presence of cavity/channel walls, boundary-driven Rayleigh streaming is generated (Hamilton, Ilinskii & Zabolotskaya 2003). This streaming can influence trapping stability and cell dynamics in a cell size dependent manner. The influence is expected to be relatively small due to the large size of RBCs (Muller et al. 2012). Furthermore, the Rayleigh streaming can be suppressed by a shape-optimized cavity/channel (Bach & Bruus 2020). Therefore, the boundary-driven Rayleigh streaming is neglected in our analysis.

Figure 1. Time-averaged mean dynamics of RBCs driven by two orthogonal ultrasonic standing waves. (a) Experimental set-up for generating two-dimensional standing waves. The background contour shows the acoustic pressure in (2.1) with phase difference $\varphi = 0$, and the colour ranges from $- {p_{am}}$ (blue) to ${p_{am}}$ (red). (b) Coordinate systems and variables describing cell motion. Here, $(\hat{x},\hat{y},\hat{z})$ is the fixed coordinate frame, $(x,y,z)$ is the body coordinate frame, $\theta $ is the inclination angle describing cell rigid tumbling and $\phi $ is the phase angle describing membrane tank treading.

As shown in figure 1(b), let $(\hat{x},\hat{y},\hat{z})$ denote a fixed Cartesian coordinate system that never changes with time and is therefore referred to as a fixed frame. The incident wave is composed of two orthogonal ultrasonic standing waves with the same acoustic pressure amplitude ${p_{am}}$ and frequency f, and the phase difference between them is $\zeta $. In the fixed frame $(\hat{x},\hat{y},\hat{z})$, the incident wave is given by

(2.1)\begin{equation}{p^{in}} = \frac{{{p_{am}}}}{2}\textrm{[}\sin ({k^e}\hat{y}) + \sin ({k^e}\hat{z})\,{\textrm{e}^{\textrm{i}\zeta }}\textrm{]}\,{\textrm{e}^{\textrm{i}\omega t}} + \textrm{c}.\textrm{c}.,\end{equation}

where ${k^e} = \omega /c_0^e$ is the acoustic wavenumber and $c_0^e$ is the sound speed in the external viscous fluid, $\omega = 2{\rm \pi} f$ is the angular frequency, $\textrm{i} = \sqrt { - 1} $ is the imaginary unit, and $\textrm{c}.\textrm{c}.$ represents the complex conjugate of the first two terms.

As shown in figure 1(b), let $(x,y,z)$ denote another Cartesian coordinate system, whose axes are the principal axes of the ellipsoidal cell, hence the name body frame. In the body frame $(x,y,z)$, the surface of the ellipsoidal cell is defined as

(2.2)\begin{equation}{\left( {\frac{x}{{{a_x}}}} \right)^2} + {\left( {\frac{y}{{{a_y}}}} \right)^2} + {\left( {\frac{z}{{{a_z}}}} \right)^2} = 1,\end{equation}

where ${a_x}$, ${a_y}$ and ${a_z}$ are the three semi-axes. Here, the $x$-axis coincides with the $\hat{x}$-axis, while the $y$- and $z$-axes are rotated by an angle $\theta $ relative to the $\hat{y}$- and $\hat{z}$-axes, respectively, as shown in figure 1(b). The inclination angle $\theta $ describes the rotation of the ellipsoidal cell. Therefore, the inclination angle $\theta $ is a function of time, and $\dot{\theta } = \textrm{d}\theta /\textrm{d}t$ is the angular velocity of the ellipsoidal cell. In the following, the theoretical framework is conveniently formulated in the body frame $(x,y,z)$.

2.1. Mean dynamics of RBCs

The internal and external media of the cell are assumed to be compressible viscous fluids, whose fluid properties are characterized by shear viscosity ${\eta ^{i,e}}$, bulk viscosity $\eta _b^{i,e}$, mass density $\rho _0^{i,e}$ and speed of sound $c_0^{i,e}$. The superscripts ‘$i$’ and ‘$e$’ relate to the internal and external media of the cell, respectively. Under high frequency acoustic excitation, the fluid dynamics is a combination of the fast time scale of ultrasound propagation $({\sim} 1\ {\rm \mu}{\rm s})$ and the slow time scale of mean dynamics $({\sim} 0.1\ \textrm{s)}$ for typical ultrasonic frequencies ~1 MHz. Here, the mean dynamics is actually time-averaged dynamics. During ultrasonic propagation on fast time scales, its nonlinear second-order terms produce a time-averaged mean stress, which drives the motion of cells on slow time scales. This process is called the mean dynamics of cells driven by ultrasonic waves.

The acoustic perturbation method based on the generalized Lagrangian formulation is employed to study the mean dynamics of cells driven by ultrasound (Nama, Huang & Costanzo 2017). In contrast to the widely used fully Eulerian formulation (Bruus 2012), this formulation employs the perturbation expansion in the mean configuration, which is not disturbed by acoustic oscillation, thus having an exact definition of the boundary conditions at cell membranes. As shown in figure 2(a), the generalized Lagrangian formulation adopts separate definitions of the mean configuration and the deformed configuration. In the mean configuration, the particle motion follows the mean smooth trajectory, while in the deformed configuration, the particle motion follows an actual oscillating trajectory, which differs from the mean trajectory by acoustic oscillating displacement. The compressible Navier–Stokes equations in the deformed configuration are reformulated in the mean configuration, and then the perturbation expansion is applied to linearize the reformulated Navier–Stokes equations. Specifically, the fluid quantity g is decomposed as $g = {g_0} + {g_1} + {g_2} + \cdots $, where ${g_0}$ is the zeroth-order (background) field, ${g_1}$ is the first-order field, and ${g_2}$ is the second-order field. The magnitude of the perturbation is characterized by the dimensionless acoustic Mach number $Ma = ||{\boldsymbol{v}_1}||/{c_0} \approx \textrm{1}{\textrm{0}^{ - 4}}$ with $||{\boldsymbol{v}_1}||$ being the amplitude of acoustic oscillating velocity, such that ${g_n} \propto M{a^n}$. Introducing the expansions of the reformulated Navier–Stokes equations and assuming zeroth-order fluid rest, one obtains two sets of governing equations: the first-order equations describe the ultrasound propagation and the time-averaged second-order equations describe the time-averaged mean dynamics (see Appendix A for details).

Figure 2. (a) Schematic of the material configuration ${\mathrm{{\mathcal{B}}}_0}$ before acoustic excitation, the mean configuration $\mathrm{{\mathcal{B}}}$ of particle mean motion and the deformed configuration ${\mathrm{{\mathcal{B}}}_t}$ of particle actual motion. (b) Schematic of mean and actual trajectories with initial position $\boldsymbol{X}$. The difference between the mean position $\boldsymbol{x}$ and the actual position $\boldsymbol{y}$ of the particle is measured by the acoustic oscillation displacement $\boldsymbol{\xi }$.

In ultrasonic standing waves, there are viscous and thermal boundary layers around the cell, and the effects of viscosity and thermal conduction occur in the boundary layers. For most liquids, such as water, the thermal effect is relatively small because the ratio of thermal boundary layer thickness to viscous boundary layer thickness ${\delta _t}/\delta \approx 0.3$ (Muller et al. 2012). To simplify the theoretical treatment, the coupling between the thermal field and other fluid quantities is neglected. In such a case, the first-order equations for ultrasound propagation are expressed as

(2.3a,b)\begin{equation}{\partial _t}{p_\textrm{1}} + {\rho _0}c_0^2\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}_\textrm{1}} = 0\quad \textrm{and}\quad {\rho _0}{\partial _t}{\boldsymbol{v}_\textrm{1}} - \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{P}_\textrm{1}} = \textbf{0},\end{equation}

with

(2.4)\begin{equation}{\boldsymbol{P}_\textrm{1}} =- {p_\textrm{1}}\boldsymbol{I} + \eta (\boldsymbol{\nabla }{\boldsymbol{v}_\textrm{1}} + \boldsymbol{\nabla }\boldsymbol{v}_\textrm{1}^T) + {\eta _b}(\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}_\textrm{1}})\boldsymbol{I}.\end{equation}

Here, all variables with subscript 1 correspond to first-order variables, ${p_1}$ is the first-order acoustic pressure, ${\boldsymbol{v}_\textrm{1}}$ is the first-order acoustic oscillation velocity, ${\rho _0}$ is the fluid density, ${c_0}$ is the sound speed of the fluid, $\eta$ and ${\eta _b}$ are the shear the bulk viscosities of the fluid, respectively, and $\boldsymbol{I}$ is the identity matrix. In (2.3), the nonlinear convective acceleration term is neglected in the first-order acoustic wave equation because the second-order small variables are much smaller than the first-order variables.

The second-order equations for the time-averaged mean dynamics are expressed as

(2.5a,b)\begin{equation}\boldsymbol{\nabla }\boldsymbol{\cdot }\langle {\boldsymbol{v}_2}\rangle = 0\quad \textrm{and}\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\langle {\boldsymbol{P}_2}\rangle = \textbf{0},\end{equation}

with the second-order time-averaged Piola–Kirchhoff stress

(2.6)\begin{align} \langle {\boldsymbol{P}_2}\rangle & =- \langle {p_2}\rangle \boldsymbol{I} + \eta (\boldsymbol{\nabla }\langle {\boldsymbol{v}_2}\rangle + \boldsymbol{\nabla }{\langle {\boldsymbol{v}_2}\rangle ^T}) - \eta \langle (\boldsymbol{\nabla }{\boldsymbol{v}_1}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\xi } + \boldsymbol{\nabla }{\boldsymbol{\xi }^T}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_1^T)\rangle \notag\\ & \quad - {\eta _b}\langle \boldsymbol{\nabla }{\boldsymbol{\xi }^T}:\boldsymbol{\nabla }{\boldsymbol{v}_1}\rangle \boldsymbol{I} + \langle {\boldsymbol{P}_1}({\boldsymbol{v}_1}) \boldsymbol{\cdot} [(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\xi })\boldsymbol{I} - \boldsymbol{\nabla }{\boldsymbol{\xi }^T}]\rangle . \end{align}

Here, all variables with subscript 2 correspond to second-order variables, $\langle {p_2}\rangle $ is the time-averaged second-order pressure, $\langle {\boldsymbol{v}_2}\rangle $ is the time-averaged second-order velocity, $\boldsymbol{\xi } = {\partial _t}{\boldsymbol{v}_1}$ is the acoustic oscillation displacement and the operator $\langle \cdot \rangle$ represents the time average over an acoustic oscillating period. The first line in (2.6) represents the stress of incompressible fluid, while the second line of (2.6) consists of the products of two first-order acoustic quantities, representing the driving force for the time-averaged fluid flow and cell motion. It is worth noting that, in the context of this work, there is a time-averaged second-order flow that includes the acoustic streaming generated by the driving terms and the Stokes flows generated by the cell motion. For later use, (2.5) and (2.6) can be rewritten as

(2.7a,b)\begin{equation}\boldsymbol{\nabla }\boldsymbol{\cdot }\langle {\boldsymbol{v}_2}\rangle = 0\quad \textrm{and}\quad \boldsymbol{\nabla }\boldsymbol{\cdot }[ - \langle {p_2}\rangle \boldsymbol{I} + \eta (\boldsymbol{\nabla }\langle {\boldsymbol{v}_2}\rangle + \boldsymbol{\nabla }{\langle {\boldsymbol{v}_2}\rangle ^T})]\textrm{ + }\boldsymbol{\nabla }\boldsymbol{\cdot }\mathrm{{\mathcal{B}}}({\boldsymbol{v}_1},\boldsymbol{\xi }) = \textbf{0},\end{equation}

with the tensor operator $\mathrm{{\mathcal{B}}}$ defined as

(2.8)\begin{align}\mathrm{{\mathcal{B}}}({\boldsymbol{v}_1},\boldsymbol{\xi }) &=- \eta \langle (\boldsymbol{\nabla }{\boldsymbol{v}_1}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\xi } + \boldsymbol{\nabla }{\boldsymbol{\xi }^T}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v}_1^T)\rangle - {\eta _b}\langle \boldsymbol{\nabla }{\boldsymbol{\xi }^T}:\boldsymbol{\nabla }{\boldsymbol{v}_1}\rangle \boldsymbol{I}\notag\\ &\quad + \langle {\boldsymbol{P}_1}({\boldsymbol{v}_1})\boldsymbol{\cdot }[(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\xi })\boldsymbol{I} - \boldsymbol{\nabla }{\boldsymbol{\xi }^T}]\rangle .\end{align}

In (2.8), by writing ${\boldsymbol{P}_1}({\boldsymbol{v}_1})$, ${\boldsymbol{P}_1}$ is regarded as a function of ${\boldsymbol{v}_1}$ according to the first equation of (2.3).

To study the possible motion of the cell, the membrane is assumed to have a tank-treading motion. Let $\boldsymbol{x}_s^0 = x_s^0{\boldsymbol{e}_x} + y_s^0{\boldsymbol{e}_y} + z_s^0{\boldsymbol{e}_z}$ represent the position vector of the material point on the membrane at the initial moment with $({\boldsymbol{e}_x},{\boldsymbol{e}_y},{\boldsymbol{e}_z})$ being the basis vector triad of the body frame. According to Keller & Skalak (1982), the membrane elements are prescribed to circulate the cell surface relative to the body frame, and the relative velocity of a material point attached to the cell membrane relative to the body frame is derived as

(2.9)\begin{equation}{\boldsymbol{v}^m} = \dot{\phi }[ - ({a_y}/{a_z})z{\boldsymbol{e}_y} + ({a_z}/{a_y})y{\boldsymbol{e}_z}],\end{equation}

where $\phi $ is the phase angle of the cell membrane (see figure 1b) and is a function of time t. It can be shown that ${\boldsymbol{v}^m}\boldsymbol{\cdot }\boldsymbol{n} = 0$, where $\boldsymbol{n}$ is the unit outward normal to the cell membrane. Thus, ${\boldsymbol{v}^m}$ is tangent to the cell membrane everywhere and represents the tank-treading motion. Finally, relative to the body frame, the velocity boundary conditions imposed on the external and internal fluids are

(2.10)\begin{gather}\langle \boldsymbol{v}_2^i\rangle = {\boldsymbol{v}^m}\quad \textrm{for}\;\boldsymbol{x} = {\boldsymbol{x}_s},\end{gather}

(2.11)\begin{gather}\langle \boldsymbol{v}_2^e\rangle = {\boldsymbol{v}^m}\quad \textrm{for}\;\boldsymbol{x} = {\boldsymbol{x}_s},\end{gather}

(2.12)\begin{gather}\langle \boldsymbol{v}_2^e\rangle =- \dot{\theta }{\boldsymbol{e}_x} \times \boldsymbol{x}\quad \textrm{for}\;|\boldsymbol{x}|\to \infty .\end{gather}

Recall that the equations are formulated in the body frame, and the far field flow condition (2.12) is due to the cell rigid tumbling motion. Moreover, the velocity of the material point attached to the membrane relative to the fixed frame is introduced as

(2.13)\begin{equation}{\boldsymbol{v}_s} = \boldsymbol{v}_s^{tt} + \boldsymbol{v}_s^{tu}.\end{equation}

Here, ${\boldsymbol{v}_s}$ describes the synthetic motion of the membrane tank-treading motion $\boldsymbol{v}_s^{tt}$ and the cell rigid tumbling $\boldsymbol{v}_s^{tu} = \dot{\theta }{\boldsymbol{e}_x} \times {\boldsymbol{x}_s}$.

The time-averaged flow field can be obtained by successively solving the acoustic equation (2.3) with the incident wave in (2.1), and the time-averaged dynamic equation (2.7) together with the velocity boundary conditions (2.10)–(2.12). Since (2.7) has the structure of the Stokes equations and is linear with respect to the unknown variables $(\langle {\boldsymbol{v}_2}\rangle ,\langle {p_2}\rangle )$, the solution to the time-averaged second-order flow field can be obtained by superposition of three ‘Stokes flow’ components

(2.14)\begin{equation}\langle {\boldsymbol{v}_2}\rangle = \langle \boldsymbol{v}_2^{ac}\rangle + \langle \boldsymbol{v}_2^{tu}\rangle + \langle \boldsymbol{v}_2^{tt}\rangle ,\end{equation}

where $\langle \boldsymbol{v}_2^{ac}\rangle$ represents the ‘Stokes flow’ (i.e. acoustic streaming) generated by the acoustic dissipation in the viscous fluid, $\langle \boldsymbol{v}_2^{tu}\rangle$ represents the ‘Stokes flow’ generated by the rigid tumbling of the cell, and $\langle \boldsymbol{v}_2^{tt}\rangle$ represents the ‘Stokes flow’ generated by the tank treading of the cell membrane. Specifically, the acoustic streaming $\langle \boldsymbol{v}_2^{ac}\rangle$ is driven by the body force $\boldsymbol{\nabla }\boldsymbol{\cdot }\mathrm{{\mathcal{B}}}({\boldsymbol{v}_1},\boldsymbol{\xi })$ in (2.7) and is calculated by successively solving equations (2.3) and (2.7) independently of any second-order velocity excitation at the boundaries (i.e. $\dot{\theta } = \dot{\phi } = 0$, or equivalently, $\langle {\boldsymbol{v}^{ac}}\rangle = \textbf{0}$ in velocity boundary conditions (2.10)–(2.12)). Since the first-order acoustic oscillation velocity depends linearly on the input acoustic pressure amplitude, the second-order acoustic streaming velocity depends quadratically on the input acoustic pressure amplitude, i.e. $\langle \boldsymbol{v}_2^{ac}\rangle \propto p_{am}^2$. The ‘Stokes flows’ $\langle \boldsymbol{v}_2^{tu}\rangle$ and $\langle \boldsymbol{v}_2^{tt}\rangle$ are obtained by solving the Stokes equations (i.e. (2.7) without acoustic excitation) for different velocity boundary conditions: $\langle \boldsymbol{v}_2^{tu}\rangle \propto \dot{\theta }$ driven by velocity boundary condition (2.12) and $\langle \boldsymbol{v}_2^{tt}\rangle \propto \dot{\phi }$ driven by velocity boundary conditions (2.10) and (2.11).

Finally, the traction jump across the cell membrane can be obtained from the total time-averaged stress as

(2.15)\begin{equation}\boldsymbol{f} = (\langle \boldsymbol{P}_2^e\rangle - \langle \boldsymbol{P}_2^i\rangle )\boldsymbol{\cdot }\boldsymbol{n}.\end{equation}

Similar to the flow, the fluid traction jump can be extended to

(2.16)\begin{equation}\boldsymbol{f} = {\boldsymbol{f}^{ac}} + {\boldsymbol{f}^{tu}} + {\boldsymbol{f}^{tt}},\end{equation}

where ${\boldsymbol{f}^{ac}}$, ${\boldsymbol{f}^{tu}}$ and ${\boldsymbol{f}^{tt}}$ are membrane tractions contributed by the flows $\langle \boldsymbol{v}_2^{ac}\rangle$, $\langle \boldsymbol{v}_2^{tu}\rangle$ and $\langle \boldsymbol{v}_2^{tt}\rangle$, respectively.

2.2. Acoustic streaming flow

The acoustic streaming flow $\langle \boldsymbol{v}_2^{ac}\rangle$ around a fixed rigid ellipsoid membrane tilted at an angle $\theta $ can be obtained by directly solving the acoustic equation (2.3) with the incident wave in (2.1) and the time-averaged dynamic equation (2.5), similar to the work by Hahn et al. (2016). The phase difference $\zeta$ and inclination angle $\theta $ are implicitly included in the solution. Furthermore, in this work, an acoustic decomposition technique based on the superposition principle is used to obtain the acoustic streaming flow $\langle \boldsymbol{v}_2^{ac}\rangle$. As we will see, by using the acoustic decomposition technique, the phase difference $\zeta$ and inclination angle $\theta $ are separated from the acoustic streaming flow $\langle \boldsymbol{v}_2^{ac}\rangle$. This makes the acoustic decomposition technique more suitable for incorporation into the theoretical model of Keller & Skalak (1982) for the analysis of cell dynamics. In addition, the acoustic decomposition technique provides insight into the underlying physical effects compared with directly solving the acoustic streaming flow $\langle \boldsymbol{v}_2^{ac}\rangle$.

The solution of the acoustic streaming flow $\langle \boldsymbol{v}_2^{ac}\rangle$ is formulated in the body frame. Considering that the wavelength is much larger than the size of the cell, the input acoustic pressure in (2.1) can be approximately expressed by Taylor series expansion as

(2.17)\begin{equation}{p^{in}} \approx \frac{{{p_{am}}}}{2}\textrm{(}{k^e}\hat{y} + {k^e}\hat{z}\,{\textrm{e}^{\textrm{i}\zeta }}\textrm{)}\,{\textrm{e}^{\textrm{i}\omega t}} + \textrm{c}.\textrm{c}.\end{equation}

Rewriting equation (2.17) in the body frame, one obtains

(2.18)\begin{equation}{p^{in}} \approx \frac{{{p_{am}}}}{2}\textrm{[(}\cos \theta + {\textrm{e}^{\textrm{i}\zeta }}\sin \theta \textrm{)}{k^e}y + ({\textrm{e}^{\textrm{i}\zeta }}\cos \theta - \sin \theta ){k^e}z\textrm{]}\,{\textrm{e}^{\textrm{i}\omega t}} + \textrm{c}.\textrm{c}.\end{equation}

For the incident wave with vanishing viscous mode, the second equation of (2.3) becomes ${\rho _0}{\partial _t}{\boldsymbol{v}_1} =- \boldsymbol{\nabla }\boldsymbol{\cdot }{p_1}$. Thus, the corresponding acoustic oscillation velocity is

(2.19)\begin{equation}{\boldsymbol{v}^{in}} \approx \frac{{{v_{am}}}}{2}\textrm{[}(\cos \theta + {\textrm{e}^{\textrm{i}\zeta }}\sin \theta ){\boldsymbol{e}_y} + ({\textrm{e}^{\textrm{i}\zeta }}\cos \theta - \sin \theta ){\boldsymbol{e}_z}\textrm{]}\,{\textrm{e}^{\textrm{i}\omega t}} + \textrm{c}.\textrm{c}.,\end{equation}

where ${v_{am}} = {p_{am}}/(\textrm{i}\rho _0^ec_0^e)$ is the amplitude of the acoustic particle velocity. The solution to the acoustic equation (2.3) must match with the input acoustic excitation far from the cell (for $(x,y,z) \to \infty $).

Subjected to the acoustic excitation expressed by (2.18) and (2.19), the solution to the acoustic equations (2.3) can be assumed as

(2.20)\begin{equation}{\boldsymbol{v}_1} = {\tilde{p}_{am}}[\textrm{(}\cos \theta + {\textrm{e}^{\textrm{i}\zeta }}\sin \theta \textrm{)}\,{\textrm{e}^{\textrm{i}\omega t}}{\boldsymbol{v}_T} + ({\textrm{e}^{\textrm{i}\zeta }}\cos \theta - \sin \theta )\,{\textrm{e}^{\textrm{i}\omega t}}{\boldsymbol{v}_A}] + \textrm{c}.\textrm{c}.,\end{equation}

where ${\tilde{p}_{am}} = {p_{am}}/{p_0}$ is the normalized acoustic pressure amplitude with the reference pressure ${p_0} = 1\ \textrm{MPa}$. Here, ${\boldsymbol{v}_T}$ is the one-dimensional standing wave response to the transverse acoustic excitation along the $y$-direction, while ${\boldsymbol{v}_A}$ is the one-dimensional standing wave response to the axial acoustic excitation along the $z$-direction. The subscripts ‘T’ and ‘A’ correspond to the transverse and axial acoustic excitations with respect to the cell, respectively. The linear superposition of these components constitutes the general solution. The component ${\boldsymbol{v}_T}$ is the solution of the following linear system and associated boundary conditions:

(2.21a,b)\begin{equation}(\textrm{i}\omega ){p_T} + {\rho _0}c_0^2\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}_T} = 0\quad \textrm{and}\quad (\textrm{i}\omega ){\rho _0}{\boldsymbol{v}_T} - \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{P}_T} = \textbf{0},\end{equation}

with

(2.22a,b)\begin{equation}{p_T} = {p_0}\frac{{{k^e}y}}{2}\quad \textrm{and}\quad {\boldsymbol{v}_T} = \frac{{{p_0}}}{{2\textrm{i}\rho _0^ec_0^e}}{\boldsymbol{e}_y}\ \textrm{as}\ (x,y,z) \to \infty .\end{equation}

Similarly, the component ${\boldsymbol{v}_A}$ obeys

(2.23a,b)\begin{equation}(\textrm{i}\omega ){p_A} + {\rho _0}c_0^2\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}_A} = 0\quad \textrm{and}\quad (\textrm{i}\omega ){\rho _0}{\boldsymbol{v}_A} - \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{P}_A} = \textbf{0},\end{equation}

with

(2.24a,b)\begin{equation}{p_A} = {p_0}\frac{{{k^e}z}}{2}\quad \textrm{and}\quad {\boldsymbol{v}_A} = \frac{{{p_0}}}{{2\textrm{i}\rho _0^ec_0^e}}{\boldsymbol{e}_z}\ \textrm{as}\ (x,y,z) \to \infty .\end{equation}

The solution of the time-averaged dynamic problem is obtained by substituting the acoustic solution (2.20) into the time-averaged dynamic equation (2.5). Since the time-averaged dynamic equation (2.5) is linear, the solution of $\langle \boldsymbol{v}_2^{ac}\rangle$ can be assumed to be the sum of the solution components resulting from different terms in the expansion of $\mathrm{{\mathcal{B}}}({\boldsymbol{v}_1},\boldsymbol{\xi })$ as given in (B3) (see Appendix B). Formally, the solution of $\langle \boldsymbol{v}_2^{ac}\rangle$ can be expressed in the form

(2.25)\begin{align}\langle \boldsymbol{v}_2^{ac}\rangle & = (1 + \sin 2\theta \cos \zeta )\tilde{p}_{am}^2{\boldsymbol{v}_{TT}} + (1 - \sin 2\theta \cos \zeta )\tilde{p}_{am}^2{\boldsymbol{v}_{AA}}\notag\\ & \quad + \cos 2\theta \cos \zeta \tilde{p}_{am}^2\boldsymbol{v}_{AT}^0 + \sin \zeta \tilde{p}_{am}^2\boldsymbol{v}_{AT}^{{\rm \pi} /2}. \end{align}

The terms ${\boldsymbol{v}_{TT}}$, ${\boldsymbol{v}_{AA}}$, $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$ are the solution components from the different contributions to the acoustic streaming. These terms depend only on the geometry of the cell and on the visco-acoustic properties (density, sound speed and viscosity) of the internal and external fluids. Specifically, in the first row of (2.25), the term ${\boldsymbol{v}_{TT}}$ corresponds to the nonlinear interaction of the acoustic response ${\boldsymbol{v}_T}$ with itself, and similarly, the term ${\boldsymbol{v}_{AA}}$ corresponds to the nonlinear interaction of the acoustic response ${\boldsymbol{v}_A}$ with itself. They are the solutions of the following linear problems:

(2.26) \begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}_{TT}} = 0\quad \textrm{and}\quad \begin{array}{*{20}{@{}l@{}}}\boldsymbol{\nabla }\boldsymbol{\cdot }[ - {p_{TT}}\boldsymbol{I} + \eta (\boldsymbol{\nabla }{\boldsymbol{v}_{TT}} + \boldsymbol{\nabla }\boldsymbol{v}_{TT}^T)]\\ \quad +\, \boldsymbol{\nabla }\boldsymbol{\cdot }[\mathrm{{\mathcal{B}}}({\boldsymbol{v}_T},\boldsymbol{\xi }_T^ \ast ) + \mathrm{{\mathcal{B}}}(\boldsymbol{v}_T^ \ast ,{\boldsymbol{\xi }_T})] = \textbf{0} \end{array},\end{gather}

(2.27) \begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{v}_{AA}} = 0\quad \textrm{and}\quad \begin{array}{*{20}{@{}l@{}}}\boldsymbol{\nabla }\boldsymbol{\cdot }[ - {p_{AA}}\boldsymbol{I} + \eta (\boldsymbol{\nabla }{\boldsymbol{v}_{AA}} + \boldsymbol{\nabla }\boldsymbol{v}_{AA}^T)]\\ \quad +\, \boldsymbol{\nabla }\boldsymbol{\cdot }[\mathrm{{\mathcal{B}}}({\boldsymbol{v}_A},\boldsymbol{\xi }_A^ \ast ) + \mathrm{{\mathcal{B}}}(\boldsymbol{v}_A^ \ast ,{\boldsymbol{\xi }_A})] = \textbf{0} \end{array},\end{gather}

where the asterisk indicates complex conjugate. Physically, ${\boldsymbol{v}_{TT}}$ and ${\boldsymbol{v}_{AA}}$ can be understood as acoustic streaming flow driven by a one-dimensional standing wave propagating in the $y$-direction and $z$-direction, respectively. The terms $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$ in the second line of (2.25) correspond to the nonlinear interaction between the acoustic response ${\boldsymbol{v}_T}$ and the acoustic response ${\boldsymbol{v}_A}$. They are the solutions of the following linear problem:

(2.28) \begin{equation}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}_{AT}^Z = 0\quad \textrm{and}\quad \begin{array}{*{20}{@{}l@{}}} {\boldsymbol{\nabla }\boldsymbol{\cdot }[ - p_{AT}^Z\boldsymbol{I} + \eta (\boldsymbol{\nabla }\boldsymbol{v}_{AT}^Z + \boldsymbol{\nabla }{\bf (}\boldsymbol{v}_{AT}^Z{{\bf )}^T})] + \boldsymbol{\nabla }\boldsymbol{\cdot }[\mathrm{{\mathcal{B}}}({\boldsymbol{v}_T},{{({\textrm{e}^{\textrm{i}Z}}{\boldsymbol{\xi }_A})}^ \ast })}\\ {\quad +\, \mathrm{{\mathcal{B}}}({\textrm{e}^{\textrm{i}Z}}{\boldsymbol{v}_A},\boldsymbol{\xi }_T^ \ast ) + \mathrm{{\mathcal{B}}}(\boldsymbol{v}_T^ \ast ,{\textrm{e}^{\textrm{i}Z}}{\boldsymbol{\xi }_A}) + \mathrm{{\mathcal{B}}}({{({\textrm{e}^{\textrm{i}Z}}{\boldsymbol{v}_A})}^ \ast },{\boldsymbol{\xi }_T})] = \textbf{0}} \end{array},\end{equation}

with $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$ corresponding to $Z = \textrm{0}$ and $Z = {\rm \pi}/2$, respectively.

By successively solving the acoustic problem posed by (2.21)–(2.24) and the time-averaged dynamic problem posed by (2.26)–(2.28), one can obtain the acoustic streaming components ${\boldsymbol{v}_{TT}}$, ${\boldsymbol{v}_{AA}}$, $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$ (see Appendix C for detailed numerical calculations). According to (2.25), the acoustic-induced traction jump ${\boldsymbol{f}^{ac}}$ can be expanded as

(2.29) \begin{align} {\boldsymbol{f}^{ac}} & = (1 + \sin 2\theta \cos \zeta )\tilde{p}_{am}^2{\,\boldsymbol{f}_{TT}} + (1 - \sin 2\theta \cos \zeta )\tilde{p}_{am}^2{\,\boldsymbol{f}_{AA}}\notag\\ & \quad + \cos 2\theta \cos \zeta \tilde{p}_{am}^2\,\boldsymbol{f}_{AT}^0 + \sin \zeta \tilde{p}_{am}^2\,\boldsymbol{f}_{AT}^{{\rm \pi} /2}, \end{align}

where ${\boldsymbol{f}_{TT}}$, ${\boldsymbol{f}_{AA}}$, $\boldsymbol{f}_{AT}^0$ and $\boldsymbol{f}_{AT}^{{\rm \pi} /2}$ are tractions contributed by the acoustic streaming components ${\boldsymbol{v}_{TT}}$, ${\boldsymbol{v}_{AA}}$, $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$, respectively. By using the acoustic decomposition technique, (2.29) isolates the influence of the phase difference $\zeta$ and inclination angle $\theta $ on the acoustic-induced traction jump ${\boldsymbol{f}^{ac}}$.

2.3. Torque balance and energy conservation for mean dynamics

The motion of the cell is completely characterized by the inclination angle $\theta $ and the phase angle $\phi $, and their time evolution equations are established using the torque balance and energy conservation, respectively. The derivation of the time evolution equations follows the work by Keller & Skalak (1982), who considered red blood cells in shear flow. Specifically, the Stokes flows ($\langle \boldsymbol{v}_2^{tt}\rangle$ and $\langle \boldsymbol{v}_2^{tu}\rangle$) and their resulting membrane tractions ($\,{\boldsymbol{f}^{tt}}$ and ${\boldsymbol{f}^{tu}}$) can be found in the work of Keller & Skalak (1982) and references therein. The solutions of the acoustic streaming $\langle \boldsymbol{v}_2^{ac}\rangle$ and its resulting membrane traction ${\boldsymbol{f}^{ac}}$ have been formulated in § 2.2. With these solutions, the torque balance equation and the energy conservation equation for the red blood cell in a two-dimensional standing wave are given below.

The torque balance requires that the torque acting on the cell is zero, which can be expressed as

(2.30)\begin{equation}M = {\boldsymbol{e}_x}\boldsymbol{\cdot }\int_A {{\boldsymbol{x}_s} \times \boldsymbol{f}\,\textrm{d}A} = \int_A {(\kern 0.06em {y_s}\,{f_z} - {z_s}\,{f_y})\,\textrm{d}A} = 0,\end{equation}

where ${f_y}$ and ${f_z}$ denote the y and z components of membrane traction $\boldsymbol{f}$, respectively, and A denotes the surface area of the cell membrane. According to the traction decomposition in (2.16), the torque balance can be expressed as

(2.31)\begin{equation}M = {M^{ac}} + {M^{tu}} + {M^{tt}} = 0,\end{equation}

where ${M^{ac}}$, ${M^{tu}}$ and ${M^{tt}}$ are torques associated with the membrane tractions ${\boldsymbol{f}^{ac}}$, ${\boldsymbol{f}^{tu}}$ and ${\boldsymbol{f}^{tt}}$, respectively. The torques ${M^{tu}}$ and ${M^{tt}}$ are given by (see Appendix D for details)

(2.32a,b)\begin{equation}{M^{tu}} =- C\dot{\theta }(a_y^2 + a_z^2)\quad \textrm{and}\quad {M^{tt}} =- 2C{a_y}{a_z}\dot{\phi },\end{equation}

with

(2.33)\begin{equation}C = \frac{{4V}}{{a_y^2{g_y} + a_z^2{g_z}}}.\end{equation}

Here, V is the volume of the cell, and ${g_y}$ and ${g_z}$ are two geometric parameters of the ellipsoidal cell (see Appendix D for definitions). Due to the formulation of the acoustic streaming flow developed above, the corresponding torque ${M^{ac}}$ can be further expanded according to (2.29). Recalling that ${\boldsymbol{v}_{TT}}$ and ${\boldsymbol{v}_{AA}}$ are the acoustic streaming flow driven by the one-dimensional standing wave propagating in the $y$-direction and $z$-direction, respectively, the resulting membrane tractions ${\boldsymbol{f}_{TT}}$ and ${\boldsymbol{f}_{AA}}$ in (2.29) are symmetrical with respect to the three coordinate planes (i.e. the $xy$-, $yz$- and $xz$-planes), and therefore does not contribute to the torque acting on the cell. Consequently, the acoustic-induced torque is found to be

(2.34)\begin{equation}{M^{ac}} = \cos (2\theta )\cos \zeta \tilde{p}_{am}^2M_{AT}^0 + \sin \zeta \tilde{p}_{am}^2M_{AT}^{{\rm \pi} /2},\end{equation}

where $M_{AT}^0$ and $M_{AT}^{{\rm \pi} /2}$ are the torques contributed by the membrane tractions $\boldsymbol{f}_{AT}^0$ and $\boldsymbol{f}_{AT}^{{\rm \pi} /2}$, respectively. The first term in (2.34) depends on the inclination angle with a factor $\cos (2\theta )$ and hence alternates between positive and negative with respect to $\theta$, while the second term is independent of $\theta$. These two torques were also qualitatively identified by Bernard et al. (2017) based on the ideas of acoustic radiation torque and acoustic viscous torque, but their work did not give an accurate method to calculate the two torques. Substituting (2.32) and (2.34) into (2.31), the evolution equation of the inclination angle $\theta$ is

(2.35)\begin{equation}\dot{\theta } = \frac{{\cos \zeta \tilde{p}_{am}^2M_{AT}^0\cos (2\theta ) + \sin \zeta \tilde{p}_{am}^2M_{AT}^{{\rm \pi} /2} - 2C{a_y}{a_z}\dot{\phi }}}{{C(a_y^2 + a_z^2)}}.\end{equation}

In (2.35), the fluid viscosity appears explicitly in the parameter C originating from the torques ${M^{tu}}$ and ${M^{tt}}$ generated by the Stokes flows, while appears implicitly in the torques $M_{AT}^0$ and $M_{AT}^{{\rm \pi} /2}$ generated by the acoustic streaming.

The energy conservation provides a constraint on the possible motion of the cell, that is, the rate of work done by the surrounding fluid is equal to the elastic power stored or dissipated in the membrane. By employing (2.13), the rate of work done by the surrounding fluid can be written as

(2.36)\begin{equation}\dot{W} = \int_A {{\boldsymbol{v}_s}\boldsymbol{\cdot }\boldsymbol{f}\,\textrm{d}A} = \int_A {\boldsymbol{v}_s^{tt}\boldsymbol{\cdot }\boldsymbol{f}\,\textrm{d}A} + \dot{\theta }{\boldsymbol{e}_x}\boldsymbol{\cdot }\int_A {{\boldsymbol{x}_s} \times \boldsymbol{f}\,\textrm{d}A} .\end{equation}

The last integral in (2.36) disappears according to the torque balance, so

(2.37)\begin{equation}\dot{W} = \int_A {\boldsymbol{v}_s^{tt}\boldsymbol{\cdot }\boldsymbol{f}\,\textrm{d}A} .\end{equation}

Substituting the membrane tank-treading velocity $\boldsymbol{v}_s^{tt}$ given in (2.9) into (2.37), $\dot{W}$ can be expressed as

(2.38)\begin{equation}\dot{W} = \dot{\phi }W\quad \textrm{with}\;W = \int_A {\left( {\frac{{{a_z}}}{{{a_y}}}{y_s}{\,f_z} - \frac{{{a_y}}}{{{a_z}}}{z_s}{\,f_y}} \right)\,\textrm{d}A} .\end{equation}

Similar to the decomposition of the torque in (2.31), the work can be decomposed as

(2.39)\begin{equation}W = {W^{ac}} + {W^{tu}} + {W^{tt}},\end{equation}

where ${\dot{W}^{ac}}$, ${\dot{W}^{tu}}$ and ${\dot{W}^{tt}}$ are the rates of work associated with the membrane tractions ${\boldsymbol{f}^{ac}}$, ${\boldsymbol{f}^{tu}}$ and ${\boldsymbol{f}^{tt}}$, respectively. The works ${W^{tu}}$ and ${W^{tt}}$ are given by (see Appendix D for details)

(2.40a,b)\begin{equation}{\dot{W}^{tu}} =- 2{a_y}{a_z}C\dot{\theta }\quad \textrm{and}\quad {\dot{W}^{tt}} =- V{\eta ^i}{f_1}\dot{\phi } + V{\eta ^e}{f_2}\dot{\phi } - \frac{{4a_y^2a_z^2}}{{a_y^2 + a_z^2}}C\dot{\phi }.\end{equation}

Here, ${f_1}$ and ${f_2}$ are two geometric parameters of the ellipsoidal cell (see Appendix D for definitions). Similar to the torque ${M^{ac}}$, the work ${W^{ac}}$ in (2.39) is expanded according to (2.29). Also, based on the symmetry of membrane tractions ${\boldsymbol{f}_{TT}}$ and ${\boldsymbol{f}_{AA}}$, it can be clarified that they do not contribute to the work of the membrane tank-treading motion. Since the acoustic streamings ${\boldsymbol{v}_{TT}}$ and ${\boldsymbol{v}_{AA}}$ represent the acoustic streaming flow driven by one-dimensional standing wave propagation, this reflects the factor that one-dimensional standing waves will not lead to the membrane tank-treading motion. Therefore, the work ${W^{ac}}$ can be expanded as

(2.41)\begin{equation}{W^{ac}} = \cos (2\theta )\cos \zeta \tilde{p}_{am}^2W_{AT}^0 + \sin \zeta \tilde{p}_{am}^2W_{AT}^{{\rm \pi} /2},\end{equation}

where $W_{AT}^0$ and $W_{AT}^{{\rm \pi} /2}$ are understood as the works contributed by the terms $\boldsymbol{f}_{AT}^0$ and $\boldsymbol{f}_{AT}^{{\rm \pi} /2}$, respectively. However, the mechanical behaviour of the cell membrane is modelled by the Kelvin–Voigt viscoelastic model. The membrane stress is $\boldsymbol{\sigma } = 2{\eta ^m}\boldsymbol{D} + 2{\mu ^m}\boldsymbol{E}$, where ${\eta ^m}$ and ${\mu ^m}$ are the viscosity and shear modulus of the membrane, respectively. Here, the membrane viscosity ${\eta ^m}$ and shear modulus ${\mu ^m}$ are related to the two-dimensional membrane viscosity $\eta _{2D}^m$ and shear modulus $\mu _{2D}^m$ reported in other literature (e.g. Tsubota 2021; Rezghi & Zhang 2022) by $\eta _{2D}^m = {\eta ^m}e$ and $\mu _{2D}^m = {\mu ^m}e$, respectively, where e is the membrane thickness. Here, $\boldsymbol{D}$ and $\boldsymbol{E}$ are the Eulerian strain rate tensor and the Eulerian strain tensor of the membrane, respectively. The elastic energy stored and dissipated in the membrane elements during tank treading is (Abkarian et al. 2007)

(2.42)\begin{equation}\dot{E} = \int_\varOmega {\boldsymbol{\sigma }:\boldsymbol{D}\,\textrm{d}\varOmega } = \dot{\phi }\left[ {{f_1}{\eta^m}\dot{\phi } + \frac{1}{2}\dot{\phi }{f_1}{\mu^m}\sin (2\phi )} \right]\varOmega ,\end{equation}

where $\varOmega = Ae$ is the membrane volume. It is important to note that ${\mu ^m}$ in (2.42) should be understood as the effective shear modulus, which contains the effect of the initial deformation from the unstressed shape to the elliptical shape of the cell at rest (Dupire et al. 2015). Using energy conservation, which is translated as $\dot{W} = \dot{E}$, and employing (2.35) and (2.39)–(2.42), the evolution equation for the phase angle $\phi$ is

(2.43)\begin{align}\dot{\phi } = \frac{{\cos \zeta \tilde{p}_{am}^2(W_{AT}^0 - {f_3}M_{AT}^0)\cos (2\theta ) + \sin \zeta \tilde{p}_{am}^2(W_{AT}^{{\rm \pi} /2} - {f_3}M_{AT}^{{\rm \pi} /2}) - {\textstyle{1 \over 2}}\varOmega {\mu ^m}{f_1}\sin (2\phi )}}{{(V{\eta ^i}{f_1} - V{\eta ^e}{f_2} + \varOmega {\eta ^m}{f_1})}},\end{align}

where ${f_3} = (2{a_y}{a_z})/(a_y^2 + a_z^2)$ is a geometric parameter of the ellipsoidal cell. In (2.43), the fluid viscosity appears explicitly in the denominator originating from the works ${W^{tu}}$ and ${W^{tt}}$ generated by the Stokes flows, while appears implicitly in the torques $M_{AT}^0$, $M_{AT}^{{\rm \pi} /2}$, $W_{AT}^0$ and $W_{AT}^{{\rm \pi} /2}$ generated by the acoustic streaming.

Equations (2.35) and (2.43) are the extension of the theoretical model describing RBC dynamics in shear flow by Keller & Skalak (1982) to 2-D standing waves. Exactly, by replacing the torque and work due to 2-D standing waves with those due to shear flow in (2.35) and (2.43), one can recover the simplified theoretical model describing RBC dynamics in shear flow. After calculating $M_{AT}^0$, $M_{AT}^{{\rm \pi} /2}$, $W_{AT}^0$ and $W_{AT}^{{\rm \pi} /2}$ (see Appendix C), the solutions of the cell motion are obtained by numerically solving the two coupled first-order ordinary differential equations (2.35) and (2.43) using the fourth-order Runge–Kutta method. The results are independent of the initial positions including the initial inclination angle ${\theta _0}$ and initial phase angle ${\phi _0}$, that is, changing the initial position will only produce a time shift in the results. Therefore, they are fixed to be $({\theta _0},{\phi _0}) = (0,0)$.

3. Results and discussion

The geometric parameters of the ellipsoidal RBC are considered to be consistent with previous studies: ${a_1} = {a_2} = 4\ {\rm \mu}{\rm m}$, ${a_3} = 1.5\ {\rm \mu}{\rm m}$ and $e = 50\ \textrm{nm}$ (Abkarian et al. 2007; Dupire et al. 2015). The physical properties of the RBC cytoplasm are mass density $\rho _0^i = 1168\ \textrm{kg}\ {\textrm{m}^{ - 3}}$, sound speed $c_0^i = 1680\ \textrm{m}\ {\textrm{s}^{ - 1}}$ and shear viscosity ${\eta ^i} = 6 \times {10^{ - 3}}\ \textrm{Pa}\ \textrm{s}$ (Bagchi, Johnson & Popel 2005; Mishra, Hill & Glynne-Jones 2014), while the physical properties of the surrounding fluid (set to water) are $\rho _0^e = 998\ \textrm{kg}\ {\textrm{m}^{ - 3}}$, $c_0^e = 1498\ \textrm{m}\ {\textrm{s}^{ - 1}}$ and ${\eta ^e} = 0.89 \times {10^{ - 3}}\ \textrm{Pa}\ \textrm{s}$. For simple fluids, the bulk viscosity is ${\eta _b} =- 2\eta /3$. The operating frequency is fixed at $f = 375\ \textrm{kHz}$ (Hahn et al. 2016). In the typical acoustic microfluidic experiment, the range of phase difference is $\zeta = 0 \sim 0.5{\rm \pi} $, and the range of acoustic pressure amplitude is ${p_{am}} = 0 \sim 1\ \textrm{MPa}$ (Mishra et al. 2014). Using these parameters, the acoustic streaming and RBC dynamics are calculated.

3.1. Flow structure of acoustic streaming

In this model, the dynamics of the RBC is driven by the torques and works of the acoustic streaming. To understand the response of RBCs in two-dimensional standing waves, the acoustic flow is first investigated. In (2.25), the acoustic streaming $\langle \boldsymbol{v}_2^{ac}\rangle $ is expressed as a linear combination of four acoustic streaming components ${\boldsymbol{v}_{TT}}$, ${\boldsymbol{v}_{AA}}$, $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$. Figures 3(a–d) and 3(e–h) show the flow structures of the four acoustic streaming components and the corresponding time-averaged membrane traction forces generated by them, respectively. As mentioned earlier, ${\boldsymbol{v}_{TT}}$ and ${\boldsymbol{v}_{AA}}$ can be understood as the acoustic streamings driven by one-dimensional standing waves propagating along the $y$-axis and $z$-axis, respectively. Two classical streaming flow structures emerge (Chan et al. 2022): ${\boldsymbol{v}_{TT}}$ falls into a bilayer regime outside the cell, characterized by an internal small vortex structure and an external flow extending to infinity, as shown in figure 3(a); ${\boldsymbol{v}_{AA}}$ falls into a monolayer regime as the inner layer thickness increases and eventually diverges, as shown in figure 3(b). Due to the symmetry of the acoustic streaming components ${\boldsymbol{v}_{TT}}$ and ${\boldsymbol{v}_{AA}}$ with respect to the $y$-axis and $z$-axis, the membrane traction forces ${\boldsymbol{f}_{TT}}$ and ${\boldsymbol{f}_{AA}}$ are also symmetrical with respect to the $y$-axis and $z$-axis, as shown in figures 3(e) and 3(f). Therefore, ${\boldsymbol{f}_{TT}}$ and ${\boldsymbol{f}_{AA}}$ make no contribution to the torque and work. In contrast, $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$ do not have this symmetry, as shown in figures 3(c) and 3(d), and the same is true for $\boldsymbol{f}_{AT}^0$ and $\boldsymbol{f}_{AT}^{{\rm \pi} /2}$, as shown in figures 3(g) and 3(h). The contributions of $\boldsymbol{f}_{AT}^0$ and $\boldsymbol{f}_{AT}^{{\rm \pi} /2}$ to the torque and work are therefore presented in table 1. As shown in figure 3(e–h), the maximum traction acting on the RBC membrane is ${\sim} 1\ \textrm{Pa}$, This level of traction in the 2-D standing wave considered in this work is comparable with that in the shear flow in the work of Abkarian et al. (2007). Considering that the RBC is a fixed ellipsoid but allows membrane tank treading in the theoretical model, good agreement between the theoretical results and experimental observation is achieved in the work of Abkarian et al. (2007). Therefore, it is believed that the present theoretical model, also considering the RBC as a fixed ellipsoid but allowing membrane tank treading, can give reasonable predictions of RBC dynamics.

Figure 3. (a–d) Time-averaged second-order velocity field and corresponding streamlines in the $\hat{y}\hat{z}$ plane of the body frame. (e–f) Plots of the membrane traction: the arrow plots show the direction and relative magnitude of membrane traction; the line plots show the normal component ${f_n} = \boldsymbol{f}\boldsymbol{\cdot }\boldsymbol{n}$ and tangential component ${f_t} = \boldsymbol{f}\boldsymbol{\cdot }\boldsymbol{t}$ along the arclength s with $\boldsymbol{n}$ being the unit normal vector pointing outward the cell, and $\boldsymbol{t}$ being the unit tangential vector pointing clockwise. The first to fourth rows correspond to the four acoustic streaming components ${\boldsymbol{v}_{TT}}$, ${\boldsymbol{v}_{AA}}$, $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$, respectively. The red ellipse represents the RBC membrane.

Table 1. Values of torques and works.

To gain insight into the role of the acoustic streaming components $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$, two special cases of phase differences $\xi = 0$ and ${\rm \pi} /2$ are considered. In view of (2.25), the total acoustic streamings in the two cases are

(3.1)\begin{gather}\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = 0}} = (1 + \sin 2\theta ){\boldsymbol{v}_{TT}} + (1 - \sin 2\theta ){\boldsymbol{v}_{AA}} + \cos 2\theta \boldsymbol{v}_{AT}^0,\end{gather}

(3.2)\begin{gather}\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = {\rm \pi}/2}} = {\boldsymbol{v}_{TT}} + {\boldsymbol{v}_{AA}} + \boldsymbol{v}_{AT}^{{\rm \pi} /2},\end{gather}

where the acoustic pressure amplitude is set to ${p_{am}} = 1\ \textrm{MPa}$. According to (3.1) and (3.2), $\boldsymbol{v}_{AT}^0$ and $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$ can be understood as the asymmetric part of the acoustic streaming $\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = 0}}$ and $\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = {\rm \pi}/2}}$, respectively. At the phase difference $\xi = 0$, the two-dimensional standing wave approximates the propagating one-dimensional wave, as shown in figure 4(a). As the RBC is non-spherical, the acoustic streaming $\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = 0}}$ depends on the inclination angle $\theta $. This is reflected by the different structures of $\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = 0}}$ at the four different inclination angles $\theta = 0$, $- {\rm \pi}/4$, $- {\rm \pi}/2$ and $- 3{\rm \pi} /4$, as shown in figure 5(a–d). It can also be seen that $\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = 0}}$ at $\theta =- {\rm \pi}/4$ and $- 3{\rm \pi} /4$ has the same structure as the acoustic streaming components ${\boldsymbol{v}_{AA}}$ and ${\boldsymbol{v}_{TT}}$, while its magnitude is twice those of ${\boldsymbol{v}_{AA}}$ and ${\boldsymbol{v}_{TT}}$, respectively. The dependence of $\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = 0}}$ on the inclination angle $\theta $ is responsible for the factor $\cos 2\theta$ before $\boldsymbol{v}_{AT}^0$ in the acoustic streaming expression (2.25). At the phase difference $\xi = {\rm \pi}/2$, the two-dimensional standing wave approximates a vortex field, as shown in figure 4(b). This rotationally symmetrical acoustic field leads to the same acoustic streaming structure, regardless of the orientation of the RBC, as shown in figure 5(e–h). Thus, $\langle \boldsymbol{v}_2^{ac}\rangle {|_{\xi = {\rm \pi}/2}}$ is independent of the inclination angle $\theta $, which explains the fact that the factor before $\boldsymbol{v}_{AT}^{{\rm \pi} /2}$ in the acoustic streaming expression (2.25) is independent of the inclination angle $\theta $.

Figure 4. Contours of the root-mean-square acoustic pressure in the $\hat{y}\hat{z}$ plane of the fixed frame for incident waves with different phase differences: (a) $\xi = 0$ and (b) $\xi = {\rm \pi}/2$. The acoustic pressure amplitude is ${p_{am}} = 1\ \textrm{MPa}$. The red ellipse represents the RBC membrane.

Figure 5. Acoustic streaming velocity and corresponding streamlines around the RBC in 2-D standing waves in the $\hat{y}\hat{z}$ plane of the fixed frame for different phase differences: (a–d) $\xi = 0$ and (e–h) $\xi = {\rm \pi}/2$. The first to fourth rows correspond to the four inclination angles $\theta = 0$, $- {\rm \pi}/4$, $- {\rm \pi}/2$ and ${\rm \pi} /4$, respectively. The acoustic pressure amplitude is ${p_{am}} = 1\ \textrm{MPa}$. The red ellipse represents the RBC membrane.

3.2. Mean dynamics of red blood cells

The motion modes of the cell are controlled by two main conditions. First, (2.34) indicates that the cell is subjected to two types of acoustic-induced torques: the torque in the first term of (2.34) tends to align the major axis of the oblate cell at the specific angle $\theta =- {\rm \pi}/4$, the torque in the second term tends to rotate the cell continuously. When the phase difference $\zeta $ increases from $0$ to ${\rm \pi} /2$, the magnitude of the latter torque gradually increases, and finally can exceed the former torque. Thus, the cell is aligned in a specific direction at small phase differences and tumbles at large phase differences (see the upper branch in figure 6). Second, the time-averaged fluid shear stress acting on the cell membrane, which is proportional to the square of the input acoustic pressure amplitude ${p_{am}}$, tends to circulate the membrane elements around the cell contour, while the membrane elastic restoring force tends to keep them in their initial positions. The latter is generally understood as an energy barrier and corresponds to the second term in (2.42). Let us consider the ${\rm \pi} /2$-circulation of the membrane around its contour: the membrane elements initially in the short axis region enter the long axis region, and vice versa. In this process, the cell membrane strains and stores elastic energy, which should be provided by the external work done by the time-averaged fluid shear stress. When ${p_{am}}$ is small, the shear stress is too small to overcome the energy barrier, and the membrane cannot tank tread. When ${p_{am}}$ is large, the shear stress can overcome the energy barrier, enabling the membrane to perform the tank-treading motion (see the lower branch in figure 6).

Figure 6. Schematic of the motion transition of an RBC, with the blue and green particles representing two material points attached to the cell membrane. The upper branch shows that the cell tumbles when the phase difference exceeds the critical phase difference, while the lower branch shows that the cell tank treads when the acoustic pressure amplitude exceeds the critical acoustic pressure amplitude.

More specifically, for small phase differences and small acoustic pressure amplitudes, the steady stationary state of the cell corresponds to steady solutions of (2.35) and (2.43). Letting the right-hand side of the equations equal to zero, the steady solutions of the inclination angle ${\theta ^\# }$ and phase angle ${\phi ^\# }$ can be found as

(3.3)\begin{gather}{\theta ^\# } = \frac{1}{2}{\cos ^{ - 1}}\left( { - \frac{{M_{AT}^{{\rm \pi} /2}}}{{M_{AT}^0}}\tan \zeta } \right),\end{gather}

(3.4)\begin{gather}{\phi ^\# } = \frac{1}{2}{\sin ^{ - 1}}\left[ {\frac{{2\tilde{p}_{am}^2}}{{\varOmega {\mu^m}{f_1}}}\left( {W_{AT}^{{\rm \pi} /2} - W_{AT}^0\frac{{M_{AT}^{{\rm \pi} /2}}}{{M_{AT}^0}}} \right)\sin \zeta } \right].\end{gather}

The condition for the above solution is that the argument of the inverse trigonometric function should be in the range of [−1, 1], or equivalently,

(3.5)\begin{gather}\zeta < {\zeta ^c} \equiv {\tan ^{ - 1}}\left( {\frac{{M_{AT}^0}}{{M_{AT}^{{\rm \pi} /2}}}} \right),\end{gather}

(3.6)\begin{gather}{\tilde{p}_{am}} < \tilde{p}_{am}^c \equiv \sqrt {\frac{{\varOmega {\mu ^m}{f_1}}}{{2\sin \zeta }}{{\left( {W_{AT}^{{\rm \pi} /2} - W_{AT}^0\frac{{M_{AT}^{{\rm \pi} /2}}}{{M_{AT}^0}}} \right)}^{ - 1}}} .\end{gather}

Above the critical values of phase difference ${\zeta ^c}$ and acoustic pressure amplitude $p_{am}^c$, steady solutions no longer exist, and tumbling or tank treading occurs, as shown in figure 6.

Then, a phase diagram is constructed for the dimensionless acoustic force ${\tilde{P}^{ac}} = p_{am}^2\kappa _0^eR/(e{\mu ^m})$ and phase difference $\zeta $, as shown in figure 7(a). Here, the dimensionless acoustic force is introduced and defined as ${\tilde{P}^{ac}} = p_{am}^2\kappa _0^eR/(e{\mu ^m})$, where $\kappa _0^e = 1/[\rho _0^e{(c_0^e)^2}]$ is the adiabatic compressibility of the external fluid and $R = {({a_x}{a_y}{a_z})^{1/3}}$ is the effective radius of the cell. For the simplicity of the dimensionless formulation, the dimensionless acoustic force here is expressed as the ratio of the membrane traction force outside the cell generated by 2-D standing waves to the elastic force of the cell membrane, while the actual force that drives the deformation of the cell is the difference between the membrane traction force inside and outside the cell generated by 2-D standing waves. Due to the small difference in acoustic impedance of fluid inside and outside the cell, the difference in membrane traction force inside and outside the cell is also not significant, and the actual magnitude of this net membrane traction force is shown in figure 3. The effective shear modulus of the membrane is ${\mu ^m} = 0.5\ \textrm{Pa}$ and the viscosity of the membrane is ${\eta ^m} = 0\ \textrm{Pa}\ \textrm{s}$.

Figure 7. (a) Phase diagram of RBC motion in 2-D ultrasonic standing waves as a function of the dimensionless acoustic force ${\tilde{P}^{ac}} = p_{am}^2\kappa _0^eR/(e{\mu ^m})$ and phase difference $\zeta $, showing four different regimes. The solid lines are the numerical solutions of (2.35) and (2.43), while the circles are the analytical solutions of (3.5) and (3.6). (b) Phase diagram of RBC motion in general linear flow as a function of the capillary number $Ca = G{\eta ^e}R/(e{\mu ^m})$ and parameter $\varsigma $ measuring the ratio of the strain rate to vorticity, also showing four different regimes. The mechanical properties of the cell membrane are: effective shear modulus ${\mu ^m} = 0.5\ \textrm{Pa}$ and surface viscosity ${\eta ^m} = 0\ \textrm{Pa}\ \textrm{s}$. (a) 2-D standing waves and (b) general linear flow.

As shown in figure 7(a), the above two conditions of (3.5) and (3.6) divide the phase space into four regions: in the steady stationary state for small $\zeta $ and small ${\tilde{P}^{ac}}$, the inclination angle $\theta $ and phase angle $\phi $ asymptotically approach constant values ${\theta ^\# }$ and ${\phi ^\# }$, respectively (figure 8c); in the tank-treading motion state for small $\zeta $ and large ${\tilde{P}^{ac}}$, $\theta $ rotates and $\phi $ oscillates (figure 8d); in the tumbling motion state for large $\zeta $ and small ${\tilde{P}^{ac}}$, $\theta $ oscillates and $\phi $ rotates (figure 8a); in the intermittent rotation state for large $\zeta $ and large ${\tilde{P}^{ac}}$, $\theta $ and $\phi $ both rotates (figure 8b). The boundary of the steady stationary state is determined by the analytical expressions of (3.5) and (3.6), or the numerical solutions of (2.35) and (2.43). However, the boundaries between the other three motions need to be determined numerically because the motions of the inclination angle $\theta $ and phase angle $\phi $ are nonlinearly coupled in these states. At small dimensionless acoustic force ${\tilde{P}^{ac}}$, this phase diagram is consistent with the experimental observations of Bernard et al. (2017). They investigated the effect of phase difference $\zeta $ on the motion of RBCs in a surface acoustic wave device and found that cell tumbling could only be observed when $\zeta $ was sufficiently large at the input acoustic pressure amplitude ${p_{am}} \approx 0.1\ \textrm{MPa}$. However, due to limited experimental working conditions, Bernard et al. (2017) did not detect the tank-treading motion of the cells in their experiments. In contrast, at relatively large input acoustic pressure amplitude ${p_{am}} > 0.5\ \textrm{MPa}$, our theory and simulation have well predicted the tank-treading motion of the cells here, and have calculated the distributions of acoustic streaming and acoustic-induced stress in § 3.1, thus explaining the dynamic motion mechanism of cells at the level of the stress tensor. Furthermore, a similar phase diagram was observed in our previous 2-D finite element simulations, but the intermittent rotation regime has not yet been obtained (Liu & Xin 2023b). This is the new finding of the present 3-D finite element simulation.

Figure 8. Typical temporal evolution of inclination angle $\theta $ and phase angle $\phi $ in the four different regions in figure 7(a): (a) ${p_{am}} = 0.3\ \textrm{MPa}$ and $\zeta = 0.4{\rm \pi} $; (b) ${p_{am}} = 0.7\ \textrm{MPa}$ and $\zeta = 0.4{\rm \pi} $; (c) ${p_{am}} = 0.3\ \textrm{MPa}$ and $\zeta = 0.2{\rm \pi} $; and (d) ${p_{am}} = 0.9\ \textrm{MPa}$ and $\zeta = 0.29{\rm \pi} $. (a) Tumbling, (b) intermittent rotation, (c) steady stationary state and (d) tank treading.

For comparison, the RBC dynamics in a general linear flow beyond shear flow is presented. The velocity field of the flow is expressed as

(3.7)\begin{equation}{\boldsymbol{v}^g} = (s + w)\hat{z}{\hat{\boldsymbol{e}}_y} + (s - w)\hat{y}{\hat{\boldsymbol{e}}_z},\end{equation}

where s and w are the strain rate and vorticity, respectively. The theoretical model describing the dynamics of RBCs in general linear flow is given in Appendix D. Referring to the 2-D ultrasonic standing waves, s and w are rewritten as $s = G{\cos ^2}\varsigma $ and $w = G{\sin ^2}\varsigma $, where $G = s + w$ measures the intensity of the flow and $\varsigma $ measures the ratio of the strain rate to vorticity.

A phase diagram is constructed for the capillary number $Ca = G{\eta ^e}R/(e{\mu ^m})$ and parameter $\varsigma $ measuring the ratio of the strain rate to vorticity, as shown in figure 7(b). The phase diagrams in general linear flow and 2-D standing wave are globally similar, showing four dynamical regimes. The similarity can be analysed based on the two conditions governing the RBC dynamics mentioned above. For the first condition describing the competition between two types of torques, the two types of torque can be identified in both shear flow and 2-D standing wave. When $\zeta = \varsigma = 0$, (2.1) locally approximates a one-dimensional standing wave and (3.7) describes a purely extensional flow. The suspended cells are subjected to a pure extensional effect (Liu & Xin 2023b), and the extensional effect tends to align RBCs in a specific direction. When $\zeta = \varsigma = {\rm \pi}/2$, (2.1) locally approximates a vortex acoustic field and (3.7) describes a purely rotational flow. The suspended cells are subjected to a pure rotational effect (Liu & Xin 2023b), and the rotational effect tends to rotate the RBCs continuously. The $\varsigma $ and $\zeta $ control the ratio of the two types of torque in 2-D standing wave and general linear flow, respectively; in the same way, the larger the value of $\varsigma $ and $\zeta $, the larger the ratio of the torque tending to continuously rotate the RBCs to that tending to align the RBCs in a specific direction.

For the second condition describing the competition between the work done by the fluid shear stress and the energy barrier, the capillary number $Ca$ and dimensionless acoustic force ${\tilde{P}^{ac}}$ control the ratio of the input work to the energy barrier in a 2-D standing wave and general linear flow, respectively; in the same way, the larger the value of $Ca$ and ${\tilde{P}^{ac}}$, the larger the ratio of the work to the energy barrier. As for the difference between 2-D ultrasonic standing waves and general linear flow, two separate regimes of intermittent rotational states are observed in 2-D ultrasonic standing waves, while only one regime is observed in general linear flow. This difference comes from the competition between the work done by the fluid shear stress and the energy barrier, here the fluid shear stress can be generated by cell tumbling and input energy (i.e. general linear flow and 2-D standing wave). In general linear flow, the work associated with the rotational part and the cell tumbling cancel each other out and only a regime of intermittent rotational state is observed, where the energy barrier is overcome by the work associated with the extensional part. In 2-D standing waves, the work associated with the rotational part and the cell tumbling are different. In one regime of the intermittent rotational state (small acoustic pressure amplitudes), the energy barrier is overcome by the work associated with the 2-D standing waves, while in the other regime of the intermittent rotational state (large acoustic pressure amplitudes), the energy barrier is overcome by the work associated with cell tumbling. It is worth noting that for the special case of $\varsigma = {\rm \pi}/4$, (3.7) describes the shear flow. As indicated by the dashed line in figure 7(b), only the tumbling motion and intermittent rotation are observed for shear flow due to the fixed ratio of the extension to rotation (Deschamps et al. 2009).

As can be seen from (3.6), the transition between the steady stationary state and tank treading is related to the effective shear modulus of the membrane ${\mu ^m}$. It can therefore be anticipated that the experimental values of the critical acoustic pressure amplitude fitted with (3.6) can be used to determine the effective shear modulus of cell membrane ${\mu^m}$. However, this transition is independent of the membrane viscosity ${\eta ^m}$, as the energy barrier of the membrane tank-treading motion is independent of the membrane viscosity. To assess the membrane viscosity, the dynamical parameters of tank-treading motion are analysed. Figure 9 shows the tank-treading frequency ${f^{tt}}$ (i.e. the reciprocal of the time for $2{\rm \pi} $-rotation of the phase angle) and the average inclination ${\theta ^m}$ (i.e. the average inclination angle during oscillation) versus the membrane viscosity ${\eta ^m}$, where the values of ${\eta ^m}$ cover the literature value of RBC membranes ranging as $0.7 \sim 2\ \textrm{Pa}\ \textrm{s}$ (Abkarian et al. 2007). Figure 9 also shows the influence of cytoplasmic viscosity. The range of literature values for cytoplasmic viscosity is $5 \sim 55\ \textrm{mPa}\ \textrm{s}$. Four different values of cytoplasmic viscosity${\eta ^i}$ are considered: 6, 12, 24 and $48\ \textrm{mPa}\ \textrm{s}$. Figure 9(a) shows that the tank-treading frequency ${f^{tt}}$ decreases with increasing membrane viscosity and increasing cytoplasmic viscosity. This is because the greater the viscosity of the cell membrane and cytoplasm, the greater the viscous resistance impeding the cell membrane tank-treading motion, and the slower the angular velocity or frequency of the cell tank-treading motion. It can also be observed from figure 9(a) that the four lines representing the different cytoplasmic viscosities tend to merge. Thus, at high membrane viscosities, the tank-treading frequency is almost independent of the cytoplasmic viscosity. Figure 9(b) shows that increasing membrane viscosity and increasing cytoplasmic viscosity align the cell closer to the inclination angle ${\theta ^\# }$ of RBCs with a rigid membrane, indicated by the dashed lines. This is because increasing membrane viscosity and increasing cytoplasmic viscosity stiffen the cell membrane.

Figure 9. (a) Tank-treading frequency ${f^{tt}}$ and (b) average inclination ${\theta ^m}$ versus membrane viscosity for different cytoplasmic viscosity at the acoustic pressure amplitude ${p_{am}} = 0.8\ \textrm{MPa}$ and the phase difference $\zeta = 0.2{\rm \pi} $.

4. Concluding remarks

A computational model is developed to study the time-averaged dynamics of RBCs trapped in the acoustic pressure nodes of two orthogonal ultrasonic standing waves with phase difference. In the context of the generalized Lagrangian formulation, the time-averaged mean stress caused by acoustic nonlinear effect is obtained by the acoustic perturbation method. The cell is modelled as a viscoelastic membrane enclosing a homogeneous fluid, while simplifying its dynamics by constraining two degrees of freedom: the inclination angle of the ellipsoidal cell shape and the phase angle of the potential tank-treading motion of the cell membrane. Based on torque balance and energy conservation, the evolution equations of the inclination angle and the phase angle are derived, respectively.

In two-dimensional standing waves, RBCs have four different types of motion: the steady stationary state, tumbling motion, tank-treading motion and intermittent rotation, in which both inclination angle and phase angle rotate. The transition from the steady stationary state to tumbling motion is caused by the competition between two types of acoustic-induced torques: one tends to align the cell in a specific direction and the other tends to cause the cell to tumble. Consistent with experimental observations, the transition is found by increasing the phase difference until the latter torque overcomes the former. The tank-treading motion occurs when the acoustic-induced shear stresses overcome the energy barrier of the cell membrane. This requirement is met when the acoustic pressure amplitude is relatively large (larger than $0.5\ \textrm{MPa}$) compared to those used for only rigid-body translational or rotational manipulation of cells (approximately $0.1\ \textrm{MPa}$), but still within the cell-friendly range. It is hoped that the present results will help guide future experimental work to control the time-averaged dynamics of RBCs with acoustic techniques, and further help to extract the viscoelastic properties of cell membranes by fitting experimental observations with the present computational model.