1. Introduction
The acoustic fields imposed on fluids exhibit several interesting nonlinear acoustic phenomena including an acoustic radiation force acting on particles or interfaces, and acoustic streaming (Friend & Yeo Reference Friend and Yeo2011). This subject has a long history, beginning with early investigations by Faraday (Reference Faraday1831), Rayleigh (Reference Rayleigh1884), King (Reference King1934) and Lighthill (Reference Lighthill1978). Over the last two decades, employing these acoustic forces in microscale flows has become a rapidly growing research field known as ‘acoustofluidics’. It has far-ranging applications in biological (Petersson et al. Reference Petersson, Åberg, Swärd-Nilsson and Laurell2007; Wiklund Reference Wiklund2012; Collins et al. Reference Collins, Morahan, Garcia-Bustos, Doerig, Plebanski and Neild2015; Lee et al. Reference Lee, Shao, Weissleder and Lee2015; Ahmed et al. Reference Ahmed, Ozcelik, Bojanala, Nama, Upadhyay, Chen, Hanna-Rose and Huang2016), medical (Augustsson et al. Reference Augustsson, Magnusson, Nordin, Lilja and Laurell2012; Li et al. Reference Li2015) and chemical sciences (Suslick et al. Reference Suslick, Didenko, Fang, Hyeon, Kolbeck, William B. Mcnamara, Mdleleni and Wong1999; Li & Huang Reference Li and Huang2019; Chen et al. Reference Chen, Pei, Zhao, Zhang, Wei and Hao2021). Recently, two interesting phenomena were observed in microchannel experiments: (i) inhibition of Rayleigh streaming vortices (bulk flow rolls outside the boundary layer due to interaction of the acoustic wave with a solid boundary) called acoustic streaming suppression and (ii) relocation or stabilization of inhomogeneous fluids to a stable configuration called acoustic relocation (Deshmukh et al. Reference Deshmukh, Brzozka, Laurell and Augustsson2014; Augustsson et al. Reference Augustsson, Karlsen, Su, Bruus and Voldman2016). This phenomenon opens the door to several applications including sorting of submicron particles such as bacteria and nanoparticles using acoustofluidic systems (Gautam et al. Reference Gautam, Gurung, Fencl and Piyasena2018; Van Assche et al. Reference Van Assche, Reithuber, Qiu, Laurell, Henriques-Normark, Mellroth, Ohlsson and Augustsson2020), tweezing and patterning of fluids (Karlsen & Bruus Reference Karlsen and Bruus2017; Baudoin & Thomas Reference Baudoin and Thomas2020) and on-demand stream to stream or stream to drop relocation of immiscible fluids (Hemachandran et al. Reference Hemachandran, Karthick, Laurell and Sen2019, Reference Hemachandran, Hoque, Laurell and Sen2021).
The main goal of this work is to develop a theory of nonlinear acoustics for inhomogeneous fluids that explains the phenomena of acoustic relocation and streaming suppression. Remarkably, our theory predicts that the acoustic relocation/stabilization of inhomogeneous fluids in a microchannel subjected to standing acoustic waves is possible only if there exists an impedance (
$Z=\rho c$) gradient, which agrees well with recent experiments. We demonstrate that the amplitude of the first-order fields is highly dependent on fluid configuration, thus acoustic energy density (
$E_{ac}$) varies significantly during the process of relocation and diffusion, whereas 
$E_{ac}$ is assumed to be constant in previous works (Karlsen, Augustsson & Bruus Reference Karlsen, Augustsson and Bruus2016; Karlsen Reference Karlsen2018). Also, we successfully separate the streaming term and acoustic relocation term from the generalized acoustic body force, which was previously claimed not to be possible by Karlsen et al. (Reference Karlsen, Qiu, Augustsson and Bruus2018). Furthermore, this theory also explains the recently observed relocation of immiscible fluids under acoustic fields (Hemachandran et al. Reference Hemachandran, Karthick, Laurell and Sen2019).
Although Karlsen et al. (Reference Karlsen, Augustsson and Bruus2016) explained these phenomena using the acoustic force density, 
$f_{ac}=-\boldsymbol {\nabla }\boldsymbol {\cdot }\langle p_2I+\rho _0{\boldsymbol {v}}_1{\boldsymbol {v}}_1\rangle$, it is not derived from first principles and involves the assumption on the time-averaged second-order mean Eulerian pressure as 
$\langle p_2\rangle =(1/2)[\kappa _0\langle |p_1|^2\rangle - \rho _0\langle |{\boldsymbol {v}}_1|^2\rangle ]$, which is claimed as the central hypothesis of their theory (Karlsen Reference Karlsen2018). Here, 
$\rho _0$ is the zeroth-order density, 
$\kappa _0$ is the zeroth-order compressibility, 
$p_1$ is the first-order pressure and 
${\boldsymbol {v}}_1$ is the first-order velocity. The following are our objections to this assumption: first, like any other pressure field, the mean Eulerian pressure has to be derived from the Navier–Stokes (N-S) equations. Second, if the assumed 
$p_2$ is the second-order Eulerian pressure then the pressure that results from the N-S equations after the substitution of 
$f_{ac}$ lacks clarity. In this work, we improve the clarity of the acoustic force density by deriving it from first principles without invoking any assumption on the time-averaged second-order mean Eulerian pressure 
$\langle p_2\rangle$.
2. Physics of the problem
The hydrodynamics of inhomogeneous fluids considered in this study is governed by the mass continuity, momentum and advection–diffusion equations (Landau & Lifshitz Reference Landau and Lifshitz1987)
\begin{gather} \partial_t \rho+\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho {\boldsymbol{v}})=0,\end{gather}
\begin{gather} \rho[\partial_t {\boldsymbol{v}} + ({\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla}){\boldsymbol{v}}]={-}\boldsymbol{\nabla} p +\eta\nabla^2 {\boldsymbol{v}}\nonumber\\\hspace{-3pc} +\beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}) + \rho{\boldsymbol{g}}, \end{gather}
\begin{gather} \partial_t s+{\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla} s = D \nabla^2{s}, \end{gather}
where 
$\rho$ is the density, 
${\boldsymbol {v}}$ is the velocity, 
$p$ is the pressure, 
$\eta$ is the dynamic viscosity of the fluid, 
$\xi$ is the bulk viscosity, 
$\beta =(\xi /\eta )+(1/3)$, g is acceleration due to gravity in negative y-direction 
$(\boldsymbol{g}=-g\boldsymbol{j})$, 
$s$ is the solute concentration and 
$D$ is the diffusivity. When the fluid is subjected to acoustic waves, the following thermodynamic pressure–density relation in terms of material derivative 
$({\rm d} / {\rm d} t)=\partial _t+({\boldsymbol {v}}\boldsymbol {\cdot }\boldsymbol {\nabla })$ is also required (Bergmann Reference Bergmann2005):
\begin{equation} \frac{{\rm d} \rho}{{\rm d} t}=\frac{1}{c^2}\frac{{\rm d} p}{{\rm d} t}, \end{equation}
where 
$c^2=(\partial p/\partial \rho )|_{S}$ and 
$c$ is the adiabatic local speed of sound.
According to perturbation theory (Eckart Reference Eckart1948), the dependent fields 
$f$ are decomposed as
\begin{equation} f\approx f_0({\boldsymbol{r}},\tau)+f_1({\boldsymbol{r}},\tau, t_f)+f_2({\boldsymbol{r}},\tau), \end{equation}
where 
$f_{0}$ are zeroth-order (background) fields, 
$f_1$ are first-order time-harmonic acoustic fields 
$f_1=f_a({\boldsymbol {r}},\tau )\, {\rm e}^{-{\rm i}\omega t_f}$, actuated at an angular frequency 
$\omega ({\sim }1\ {\rm MHz})$, and 
$f_{2}$ are second-order fields (in general 
$f_{2}\ll f_1$). In microscale flows, since the hydrostatic pressure 
$\rho g H ({\sim }1\ {\rm Pa}) \ll p_1 ({\sim }10^6\ {\rm Pa})$, the variation of pressure and velocity fields due to gravity is accounted in the second-order effects. Thus, in a quiescent fluid, we take the zeroth-order velocity 
${\boldsymbol {{v}}}_0=0$, and pressure 
$p_0 = {\rm const.} (\boldsymbol {\nabla } p_0 = 0)$.
The first-order acoustic fields vary on the fast time scale 
$t_f$ (
$t_f\sim {1}/{\omega } \sim 0.1\ \mathrm {\mu }{\rm s}$), whereas the second-order hydrodynamic fields vary on the slow time scale 
$\tau (\tau \gg t_f)$. Usually in perturbation theory, 
$f_0(\boldsymbol{r},\tau )$ is assumed to be constant for homogeneous fluids. Whereas for inhomogeneous fluids, the variation in the background fields (
$\rho _0$, 
$s_0$, 
$c_0$ and 
$\eta _0$) with space as well as the slow time scale has to be accounted for due to the gravity stratification and second-order acoustic effects. As the first-order acoustic fields (
$f_1$) are sensitive to the inhomogeneous configuration (
$\rho _0$, 
$s_0$, 
$c_0$ and 
$\eta _0$) which varies on the slow time scale, the amplitude of these acoustic fields are considered to be a function of the slow time scale (
$f_a({\boldsymbol {r}},\tau )$) unlike homogeneous fluids, where the amplitude of these fields is only a function of space (
$f_a({\boldsymbol {r}})$).
The variation of zeroth-order fields on the fast time scale (
$t_f$) can be neglected. Also, since 
$D\sim {O}(10^{-9})$, the diffusion term is negligible on the fast time scale, due to which the composition of any given fluid particle remains unchanged as it moves. Consequently, the governing equations up to first order (
$t\sim t_f$) reduce to (Bergmann Reference Bergmann2005)
$$\begin{gather} \partial_t \rho_1 +\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho_0 {\boldsymbol{{v}}}_1)=0, \end{gather}$$
$$\begin{gather}\rho_0 \partial_t {\boldsymbol{v}}_1 ={-}\boldsymbol{\nabla} p_1 + \eta\nabla^2 {\boldsymbol{v}}_1+\beta\eta \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1), \end{gather}$$
$$\begin{gather}\partial_ts_1+{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} {s}_0 =0, \end{gather}$$
$$\begin{gather}\partial_t\rho_1+({\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla})\rho_0=(1/c^2)[\partial_tp_1]. \end{gather}$$ The first-order fields have a harmonic time dependence, thus the time average of these fields is zero. Therefore, first-order terms in (2.3) cannot cause any bulk fluid motion. However, the N-S equation is nonlinear and the above linearized equations (2.3) are not exact. Hence, proceeding to solve (2.1) up to second order (
$t\sim \tau$),
\begin{gather} \langle\partial_t(\rho_0+\rho_2)\rangle+ \boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_1 {\boldsymbol{v}}_1\rangle + \boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0 {\boldsymbol{v}}_2\rangle = 0, \end{gather}
\begin{gather} \langle\rho_1 \partial_t {\boldsymbol{v}}_1\rangle + \rho_0\langle({\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla}){\boldsymbol{v}}_1\rangle-\langle(\rho_0 +\rho_2){\boldsymbol{g}}\rangle={-}\boldsymbol{\nabla}\langle p_2 \rangle \nonumber\\ \hspace{-3pc} + \eta\nabla^2 \langle{\boldsymbol{v}}_2\rangle+ \beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\langle{\boldsymbol{v}}_2 \rangle)-\langle\rho_0 \partial_t {\boldsymbol{v}}_2\rangle, \end{gather}
$$\begin{gather} \langle\partial_t (s_0+s_2)\rangle+\langle{\boldsymbol{v}}_2\boldsymbol{\cdot}\boldsymbol{\nabla} s_0\rangle+ \langle{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} s_1\rangle= D\nabla^2\langle(s_0+s_2)\rangle, \end{gather}$$
$$\begin{gather}\langle\partial_t(\rho_0+\rho_2)\rangle+\langle({\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla})\rho_1\rangle+\langle({\boldsymbol{v}}_2\boldsymbol{\cdot}\boldsymbol{\nabla}) \rho_0\rangle=(1/c^2)[\langle\partial_tp_2\rangle+ \langle({\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla})p_1\rangle], \end{gather}$$
where 
$\langle \cdots \rangle$ denotes time average over one oscillation period. Since the time average of the product of two first-order fields is non-zero (the time average of two first-order fields 
$\langle u_1v_1\rangle$ is defined as 
$\frac {1}{2}{\rm Real}(u_1^\star v_1)$, where 
$\star$ denotes complex conjugation), they act as a source term for second-order fields (slow hydrodynamic time scale). The first two terms on the left side of (2.4b) together comprise the divergence of the Reynolds stress tensor, 
$\boldsymbol {\nabla } \boldsymbol {\cdot} \langle \rho _0 {\boldsymbol {v}}_1\otimes {\boldsymbol {v}}_1\rangle =\boldsymbol {\nabla }\boldsymbol {\cdot} \langle \rho _0{\boldsymbol {v}}_1 {\boldsymbol {v}}_1\rangle$.
In microscale flows, the variation of second-order fields with respect to the slow time scale is negligible (Bruus Reference Bruus2011; Friend & Yeo Reference Friend and Yeo2011). From (2.3), we get 
$\rho _1\ll \rho _0$, 
$s_1\ll s_0$ (see Appendix A) and from general perturbation theory, 
$f_2\ll f_1$ (for any perturbed field 
$f$), thus 
$\rho _2\ll \rho _0$ and 
$s_2\ll s_0$. We can also neglect 
$\langle {\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla } s_1\rangle$, 
$\langle \rho _1\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_1\rangle$ since the first-order fields in both these terms are out of phase (see Appendix A). Using the above arguments, combining ((2.4a) and (2.4d)),
\begin{equation} \langle\rho_0\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_2\rangle={-}(1/c^2)\langle{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} p_1\rangle. \end{equation} Substituting the above relation in (2.4b) and analysing the order of magnitude, 
${O}(\beta \eta \boldsymbol {\nabla }({1}/{\rho _0 c_0^2})({\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla })p_1)\ll {O}(\boldsymbol {\nabla } \boldsymbol {\cdot }\langle \rho _0{\boldsymbol {v}}_1{\boldsymbol {v}}_1\rangle )$. As this term, 
$(1/c^2)\langle {\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla } p_1\rangle$ does not contribute to the momentum equation, it does not affect the second-order velocity and pressure fields and thus can be neglected. Hence, (2.5) becomes divergence free and second-order flow is incompressible i.e. 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\langle {\boldsymbol {v}}_2\rangle =0$ (Nyborg's approximation; Nyborg Reference Nyborg2005; Baasch, Doinikov & Dual Reference Baasch, Doinikov and Dual2020). Accounting for all the above arguments, the governing equations reduce to
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\langle{\boldsymbol{v}}_2\rangle=0, \end{gather}$$
$$\begin{gather}-\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle+ \langle\rho_{0}{\boldsymbol{g}}\rangle-\boldsymbol{\nabla}\langle p_2 \rangle + \eta\nabla^2 \langle{\boldsymbol{v}}_2\rangle=0, \end{gather}$$
$$\begin{gather}\langle\partial_t {s}_0\rangle+\langle{\boldsymbol{v}}_2\boldsymbol{\cdot}\boldsymbol{\nabla} {s}_0\rangle = D\nabla^2\langle{s}_0\rangle. \end{gather}$$The above (2.6) govern the dynamics of inhomogeneous fluids in microscale acoustofluidics and can also be derived by another approach, see Appendix B. From (2.6b), it is evident that second-order slow hydrodynamic flows are created due to the divergence of the Reynolds stress tensor consisting of the product of first-order fast acoustic fields. Thus, we introduce the body force due to acoustic fields as
\begin{equation} {\boldsymbol{f}}_{ac}={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle. \end{equation} It is well known that the above force is responsible for boundary-driven Rayleigh streaming and bulk-driven Eckart streaming in homogeneous fluids. In this work, we proceed to show that the same force is also responsible for recently observed streaming suppression and acoustic relocation of miscible as well as immiscible inhomogeneous fluids. A microchannel of width 
$w = 380\ \mathrm {\mu }{\rm m}$ and height 
$h = 160\ \mathrm {\mu }{\rm m}$ is chosen for study from Muller & Bruus (Reference Muller and Bruus2014) containing an inhomogeneous miscible solution, as shown in figure 1. Fluid properties are taken from Qiu et al. (Reference Qiu, Karlsen, Bruus and Augustsson2019). An acoustic standing half-wave along the width is imposed in the microchannel by actuating both the sidewalls in phase at a specific frequency ‘
$\nu$’ and wall displacement ‘
$d$’ in the 
$x$-direction. The first-order fields due to this standing wave are obtained by solving (2.3) in the frequency domain, see Appendix C. These first-order fields are then used to obtain the acoustic body force (2.7) to solve the time-averaged second-order fields (2.6). At each time step, the first-order fields are solved for the updated inhomogeneous fluid configuration (
$\rho _0, c_0, \eta _0$) from the previous time step. Hence, both the first-order and second-order equations are bidirectionally coupled and solved numerically at all slow time steps in COMSOL Multiphysics 5.6. The mesh convergence for the numerical simulation is discussed in Appendix G.
Figure 1. Sketch of acoustofluidic microchannel with an imposed half-wave acoustic pressure resonance, containing fluids whose density, speed of sound, dynamic viscosity and solute concentration are functions of space and time.
3. Results and discussion
From figure 2, the acoustic body force 
${\boldsymbol {f}}_{ac}$ tends to relocate high impedance fluid to nodes and low impedance to anti-nodes, called a stable configuration (whereas any other configuration with an impedance gradient is considered unstable). Once the stable configuration is reached, due to the existing impedance gradient, this 
${\boldsymbol {f}}_{ac}$ inhibits any fluid motion due to gravity and suppresses acoustic streaming, which tries to disturb the stable configuration, as seen in figures 2 and 3. As time progresses, due to diffusion, the fluid profile becomes homogeneous, where the same 
${\boldsymbol {f}}_{ac}$ induces boundary-driven Rayleigh streaming. From these results, it is evident that the acoustic body force 
${\boldsymbol {f}}_{ac}$ is responsible for acoustic relocation and streaming suppression in inhomogeneous fluids as well as acoustic streaming in homogeneous fluids. Remarkably, in the process of acoustic relocation and diffusion, as shown in figures 2(a) and 3(a), the amplitude of the first-order fields (
$p_a$ and 
$v_a$) varies significantly (see Appendix D) as the background 
$\rho _0$ and 
$c_0$ fields change on the slow time scale. Thus, the acoustic energy density 
$E_{ac}$, which is a function of 
$p_a$, 
$\rho _0$ and 
$c_0$, also changes significantly, in contrast to constant 
$E_{ac}$ assumed in previous studies (Karlsen et al. Reference Karlsen, Augustsson and Bruus2016; Karlsen Reference Karlsen2018; Karlsen et al. Reference Karlsen, Qiu, Augustsson and Bruus2018). Surprisingly, it is observed that, for the case of constant impedance (figure 4), relocation does not occur irrespective of the 
$\rho _0$ and 
$c_0$ configurations, thus any constant impedance configuration is called a neutral configuration. This demonstrates that the impedance gradient is the requisite and governing factor for acoustic relocation. In addition to the impedance gradient being a prerequisite, the sufficient conditions for acoustic relocation include: the fluid interface should not be at a node (See Appendix E) and the fluid configuration must be unstable.
Figure 2. (a) Relocation of unstable configuration to stable configuration with 
$d=0.261$ nm and 
$\nu=1.96$ MHz. Initially (
$t=0$ s), there is low impedance deionized (DI) water (blue) at the centre and high impedance 10 % Ficoll PM70 (red) at the sides. The green image indicates the homogeneous fluid profile after relocation and complete diffusion. (b) Second-order velocity 
$\langle {\boldsymbol {v}}_2\rangle$.
Figure 3. (a) Stable configuration with no relocation, with low impedance DI water (blue) at the sides and high impedance 10 % Ficoll PM70 (red) at the centre. The green image indicates the homogeneous fluid profile after complete diffusion. (b) Second-order velocity 
$\langle {\boldsymbol {v}}_2\rangle$. Here, 
$d=0.261$ nm and 
$\nu=1.96$ MHz. Relocation force (
${\boldsymbol {f}}_{rl}$) stabilizes this configuration against gravitational stratification.
Figure 4. For constant impedance fluids any configuration of 
$\rho _0$ and 
$c_0$ is a neutral configuration. (a) No relocation due to zero relocation force (
${\boldsymbol {f}}_{rl} = 0$). To clearly show the absence of a relocation force in the constant impedance case, gravity is neglected. The presence of gravity would stratify the fluids. The green image indicates the homogeneous fluid profile after complete diffusion. (b) Second-order velocity 
$\langle {\boldsymbol {v}}_2\rangle$. Here, 
$d=2$ nm and 
$\nu=1.96$ MHz.
In order to explain the above results mathematically, we analyse the acoustic body force, 
${\boldsymbol {f}}_{ac}=-\boldsymbol {\nabla }\boldsymbol {\cdot }\langle \rho _0 {\boldsymbol {v}}_1 {\boldsymbol {v}}_1\rangle$ in detail
\begin{equation} {\boldsymbol{f}}_{ac}=\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0 {\boldsymbol{v}}_1\otimes {\boldsymbol{v}}_1\rangle={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0 {\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle={-}\rho_0\langle{\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{v}}_1\rangle-\langle\rho_1 \partial_t {\boldsymbol{v}}_1\rangle. \end{equation} Using the following identity 
${\boldsymbol {A}}\boldsymbol {\cdot }\boldsymbol {\nabla }{\boldsymbol {A}}= \boldsymbol {\nabla }({\boldsymbol {A}}^2/2)-{\boldsymbol {A}}\times (\boldsymbol {\nabla }\times {\boldsymbol {A}})$ and substituting first-order fields from (2.3), (3.1) becomes
\begin{align} {\boldsymbol{f}}_{ac}&={-}\frac{1}{2\rho_0}\boldsymbol{\nabla}\langle|\rho_0 {\boldsymbol{v}}_1|^2\rangle+\langle{\boldsymbol{v}}_1\times\boldsymbol{\nabla}\times(\rho_0 {\boldsymbol{v}}_1)\rangle+\frac{1}{2}\kappa_0\boldsymbol{\nabla}\langle |p_1|^2\rangle\nonumber\\ &\quad -\langle(\kappa_0p_1)(\eta\nabla^2{\boldsymbol{v}}_1+ \beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle. \end{align}For a detailed derivation refer to Appendix F. From the scaling analysis, the last term in (3.2) can be neglected, thus reducing to
\begin{align} {\boldsymbol{f}}_{ac}&=\tfrac{1}{2}\boldsymbol{\nabla}(\kappa_0\langle |p_1|^2\rangle-\rho_0 \langle|{\boldsymbol{v}}_1|^2\rangle) +\langle{\boldsymbol{v}}_1\times\boldsymbol{\nabla} \times(\rho_0{\boldsymbol{v}}_1)\rangle\nonumber\\ &\quad -\tfrac{1}{2}[ \langle|p_1|^2\rangle\boldsymbol{\nabla}\kappa_0+ \langle|{\boldsymbol{v}}_1|^2\rangle\boldsymbol{\nabla}\rho_0]. \end{align} In the above (3.3), the first term 
$(1/2)\boldsymbol {\nabla }(\kappa _0\langle |p_1|^2\rangle -\rho _0 \langle |{\boldsymbol {v}}_1|^2\rangle )$ is a pure gradient term (conservative), and the curl of that term is zero, thus it does not cause relocation or streaming but only contributes to the second-order pressure. The third term 
$(-1/2)[\langle |p_1|^2\rangle \boldsymbol {\nabla }\kappa _0 +\langle |{\boldsymbol {v}}_1|^2\rangle \boldsymbol {\nabla }\rho _0]$ is responsible for relocation, since the curl of this term, in general, is non-zero. From (2.3b), the first-order fields in the inviscid region satisfy the following relation, 
$\boldsymbol {\nabla }\times (\rho _0 {\boldsymbol {v}}_1)=0$. Thus the second term 
$\langle {\boldsymbol {v}}_1\times \boldsymbol {\nabla }\times (\rho _0{\boldsymbol {v}}_1)\rangle$ is only significant inside the boundary layer (
$\delta \sim 1\ \mathrm {\mu }{\rm m}$) of first-order fields, which is responsible for boundary-driven streaming in inhomogeneous fluids. The competition between the relocation force (third term) and streaming force (second term) accounts for streaming suppression, as seen in figures 2 and 3. In (3.3), we have analytically separated relocation and streaming causing terms from the general body force term, 
${\boldsymbol {f}}_{ac}$, which was previously claimed not to be possible by Karlsen et al. (Reference Karlsen, Qiu, Augustsson and Bruus2018). For homogeneous fluids (3.3) reduces to (Bradley Reference Bradley1998; Friend & Yeo Reference Friend and Yeo2011)
\begin{equation} {\boldsymbol{f}}_{ac}={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle =\tfrac{1}{2}(\kappa_0\boldsymbol{\nabla}\langle |p_1|^2\rangle-\rho_0\boldsymbol{\nabla}\langle|{\boldsymbol{v}}_1|^2\rangle) +\{\rho_0\langle{\boldsymbol{v}}_1\times(\boldsymbol{\nabla}\times{\boldsymbol{v}}_1)\rangle\}, \end{equation}where the first term is the homogeneous second-order mean Eulerian pressure and the second term is responsible for acoustic streaming.
To understand the second-order pressure and relocation phenomenon in inhomogeneous fluids clearly, we study the acoustic body force outside the boundary layer of first-order fields (neglecting the streaming term). Thus, (3.3), after time averaging, reduces to
\begin{align} {\boldsymbol{f}}_{ac}&={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle =\tfrac{1}{4}\boldsymbol{\nabla}(\kappa_0 |p_1|^2-\rho_0 |{\boldsymbol{v}}_1|^2) [-\tfrac{1}{4}|p_1|^2\boldsymbol{\nabla}\kappa_0-\tfrac{1}{4}|{\boldsymbol{v}}_1|^2 \boldsymbol{\nabla}\rho_0]\nonumber\\ &= {\boldsymbol{f}}_{1}+[{\boldsymbol{f}}_{2}]. \end{align}
Considering the case of a stable configuration (figure 3), the second-order velocity will be zero as relocation does not take place due to which (2.6b) will reduce to 
$-\boldsymbol {\nabla }\boldsymbol {\cdot }\langle \rho _0{\boldsymbol {v}}_1 {\boldsymbol {v}}_1\rangle = \boldsymbol {\nabla }\langle p_2\rangle$. Now, it is evident that, in the inviscid case, both the terms 
$(\boldsymbol{f}_1, \boldsymbol{f}_2)$ in (3.5) contribute to the mean Eulerian pressure 
$\langle p_2\rangle$, whereas, in Karlsen et al. (Reference Karlsen, Augustsson and Bruus2016), Karlsen (Reference Karlsen2018) and Karlsen et al. (Reference Karlsen, Qiu, Augustsson and Bruus2018), it is hypothesized that only the first term in (3.5) contributes to the mean Eulerian pressure, 
$\langle p_2\rangle =\frac {1}{2}\kappa _0 \langle |p_1|^2\rangle -\frac {1}{2}\rho _0\langle |{\boldsymbol {v}}_1|^2\rangle$. Thus, our work clearly demonstrates that 
$\langle p_{2}\rangle$ can be derived from the governing equations and their assumption on 
$\langle p_{2}\rangle$ is incorrect and needless.
In the case of constant impedance (
$Z_0=\rho _0 c_0 ={\rm const.}$) inhomogeneous fluids, subjected to a one-dimensional acoustic standing half-wave (
$\lambda \approx 2w$), using the relations 
$p_1=p_a \sin (k{x})$ and 
$v_1=({p_a}/{i\rho _0 c_0})\cos (k{x})$ where 
$k= 2{\rm \pi} /\lambda$ is wavenumber, the relocation force term 
${\boldsymbol {f}}_{2}$ reduces to
\begin{equation} {\boldsymbol{f}}_{2|_{Z=C}}={-}\frac{1}{4} |p_1|^2\boldsymbol{\nabla}\kappa_0-\frac{1}{4}|{\boldsymbol{v}}_1|^2\boldsymbol{\nabla}\rho_0 ={-}\boldsymbol{\nabla}\left(\frac{p_a^2\rho_0}{4Z^2}\right). \end{equation} Since 
${\boldsymbol {f}}_{2}$ reduces to a purely gradient term in the case of constant impedance inhomogeneous fluids, it does not induce relocation or motion but only contributes to the second-order pressure. Thus, we prove that the impedance gradient is the necessary condition for acoustic relocation, which agrees with experimental results (Deshmukh et al. Reference Deshmukh, Brzozka, Laurell and Augustsson2014; Hemachandran et al. Reference Hemachandran, Karthick, Laurell and Sen2019).
For the case of variable impedance, the relocation term 
$\boldsymbol{f}_2$ in (3.5) is written as
$$\begin{gather} {\boldsymbol{f}}_{2}=\frac{p_a^2 \sin^2({kx}) }{4}\left(\frac{\boldsymbol{\nabla}\rho_0}{\rho_0^2 c_0^2}\right) + \frac{p_a^2 \sin^2({kx})}{2}\left(\frac{\boldsymbol{\nabla} c_0}{\rho_0 c_0^3}\right) - \frac{p_a^2 \cos^2({kx})}{4}\left(\frac{\boldsymbol{\nabla}\rho_0}{\rho_0^2 c_0^2}\right), \end{gather}$$
$$\begin{gather}{\boldsymbol{f}}_{2}={-}\frac{p_a^2\cos({2kx})}{4\rho_{0}^2 c_{0}^3}\boldsymbol{\nabla} Z_0 + \frac{p_a^2}{4\rho_{0} c_{0}^3}\boldsymbol{\nabla} c_0. \end{gather}$$ Analogous to the Boussinesq approximation, 
${\boldsymbol {\nabla } Z_0}/{\rho _0^2 c_0^3}\approx {\boldsymbol {\nabla } Z_0}/{\rho _{avg}^2 c_{avg}^3}, {\boldsymbol {\nabla } c_0}/{\rho _0 c_0^3}\approx {\boldsymbol {\nabla } c_0}/{\rho _{avg} c_{avg}^3}$ and substituting (3.8) in (3.5), separating gradient and non-gradient terms,
\begin{align} {\boldsymbol{f}}_{ac}&={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle ={-}E_{ac}\cos({2kx})\boldsymbol{\nabla}\hat Z_0\nonumber\\ &\quad -\boldsymbol{\nabla}(\tfrac{1}{4}\rho_0 |{\boldsymbol{v}}_1|^2-\tfrac{1}{4}\kappa_0 |p_1|^2-E_{ac}\hat c_0), \end{align}
where 
$E_{ac}=p_a^2/(4\rho _{avg} c_{avg}^2)$, 
$\hat c_0 = c_0/ c_{avg}$, 
$\hat \rho _0 = \rho _0/ \rho _{avg}$ and 
$\hat Z_0 = \hat \rho _0\hat c_0$.
It is clear that only the first term in (3.9) is responsible for acoustic relocation in inhomogeneous fluids, whereas the second term resembles a conservative or purely gradient term, thus inducing only pressure. Thus, in addition to the impedance gradient being the necessary factor for relocation, now we write the relocation force in terms of the impedance gradient as follows:
\begin{equation} {\boldsymbol{f}}_{rl}={-}E_{ac}\cos({2kx})\boldsymbol{\nabla}\hat Z_0. \end{equation} The above force term, which is a part of the generalized force (2.7), is responsible for the relocation of the unstable configuration in figure 2 and maintaining the stable configuration by inhibiting acoustic streaming as well as gravity stratification in figure 3. Whereas, for constant impedance fluids (figure 4), irrespective of any fluid configuration, the relocation force is always absent (
${\boldsymbol {f}}_{rl} = 0$).
In several previous studies (Karlsen & Bruus Reference Karlsen and Bruus2017; Nath & Sen Reference Nath and Sen2019; Qiu et al. Reference Qiu, Karlsen, Bruus and Augustsson2019), including ours (Pothuri, Azharudeen & Subramani Reference Pothuri, Azharudeen and Subramani2019; Kumar et al. Reference Kumar, Azharudeen, Pothuri and Subramani2021), (3.8) has been used to study relocation without realizing that it contains implicit gradient or conservative terms which do not contribute to relocation. This leads to an incorrect scaling analysis and thus it is necessary to ignore this gradient term in such circumstances. Furthermore, (3.10) captures all the aspects of relocation and does not cause scaling issues like (3.8). It must be noted that (3.10) is valid only outside the boundary layer of first-order fields, whereas inside the boundary layer, it cannot be reduced from 
${\boldsymbol {f}}_{2}$ of (3.5). This is because 
$p_1=p_a \sin (kx)$ remains the same inside and outside the boundary layer (
$\delta \sim 1\ \mathrm {\mu }{\rm m}$), but 
$v_1\neq ({p_a}/{i\rho _0 c_0})\cos (kx)$ and 
$\partial _y v_1$ is significant inside the boundary layer due to the no-slip condition. Thus, in the case of constant impedance fluids, the relocation force is absent in the bulk but non-zero within the boundary layer, which is responsible for the disturbance of homogeneous streaming, as in figure 4.
Immiscible fluids. The phenomenon of acoustic relocation of immiscible fluids due to a standing acoustic wave is also governed by (2.7). This theory predicts that, for the relocation of immiscible fluids to occur, the applied acoustic energy density 
$E_{ac}$ must be greater than the threshold 
$E_{cr}$ to overcome the interfacial tension force. The above prediction is in agreement with the experimental studies by Hemachandran et al. (Reference Hemachandran, Karthick, Laurell and Sen2019). However, in their study, relocation is achieved irrespective of the interface location with respect to the node. Also, the frequency employed is very different from the resonant frequency (
$\nu=c_{avg}/2w$) that corresponds to the standing half-wave actuated only along the width of the microchannel, as shown in figure 1. Contrary to their study, when we actuated only sidewalls, relocation seemed to be highly dependent on the position of the interface and was absent when the interface and node coincide. Whereas, when all the walls were actuated in the direction of their surface normal, we achieved relocation, independent of the interface location, (figure 5), which is in agreement with their experiments (Hemachandran et al. Reference Hemachandran, Karthick, Laurell and Sen2019). It is evident that relocation can be achieved through two modes; one is due to a one-dimensional standing wave resulting from the actuation of the sidewalls at the one-dimensional resonant frequency (
$\nu=c_{avg}/2w$) while the other is due to a two-dimensional standing wave resulting from the actuation of all walls at the two-dimensional resonant frequency (between 
$\nu=c_{avg}/2w$ and 
$\nu=c_{avg}/2h$). Since the relocation in Hemachandran et al. (Reference Hemachandran, Karthick, Laurell and Sen2019) is due to a two-dimensional standing wave, the frequency required (
${\sim }2.1\text {--}2.4$ MHz) for relocation (for silicone–mineral oil) is different from the one-dimensional resonant frequency (
$1.66$ MHz). The velocity and pressure fields of the relocation of immiscible fluids due to a two-dimensional standing wave are shown in Appendix E.
Figure 5. Relocation of immiscible fluids by actuating all walls of a microchannel (two-dimensional), consisting of high impedance mineral oil (red) and low impedance silicone oil (blue) with a surface tension of 
$1$ mN m
$^{-1}$ actuated at a frequency of 
$2.1$ MHz for unstable configurations. (a) No relocation for 
$E_{ac} < E_{cr}$, slight fluid displacement is attributed to slanted node alignment. (b) Relocation for 
$E_{ac} > E_{cr}$. (c) Stable configuration, no relocation for any 
$E_{ac}$, slight fluid displacement is attributed to slanted node alignment.
4. Conclusion
We have put forward a theory of nonlinear acoustics that governs the acoustic phenomena of relocation and streaming suppression in inhomogeneous miscible and immiscible fluids (including streaming in homogeneous fluids). This theory also confirms the fact that the divergence of the time-averaged Reynolds tensor 
$-\boldsymbol {\nabla }\boldsymbol {\cdot }\langle \rho _0{\boldsymbol {v}}_1{\boldsymbol {v}}_1\rangle$ is alone responsible for all the above processes. We showed that first-order fields 
$p_1, {\boldsymbol {v}}_1$ and energy density 
$E_{ac}$ vary significantly during the process of acoustic relocation and diffusion. Importantly, we have proved that an impedance gradient is the necessary condition for relocation. The other conditions for acoustic relocation in one- and two-dimensional modes will be addressed in detail using stability analysis for both miscible and immiscible fluids in an upcoming paper. The fundamental understanding from this study can give new insights into particle (cells/drops/beads) and inhomogeneous fluid handling in microchannels under acoustic fields.
Funding
This work is supported by IIITDM Kancheepuram via Grant No: IRG/2019-20/KS/006.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Order of magnitude analysis
The order of magnitude analysis for different terms considered at various steps in the derivations is given below. The proof of the stated assumption 
$\rho _1\ll \rho _0$ and 
$s_1\ll s_0$ is as follows. We shall first obtain the order for 
${\boldsymbol {v}}_1$ from (2.3b). The viscous effects are dominant only inside the thin viscous boundary layer (
$\sim 1\ \mathrm {\mu }{\rm m}$), due to which we can take only the first term, i.e. 
$\rho _0\partial _t{\boldsymbol {v}}_1$, which has an effect on the entire medium. Thus comparing its magnitude with 
$\boldsymbol {\nabla } p_1$, we get
\begin{equation} {O}(\rho_0\partial_t{\boldsymbol{v}}_1)\sim{O}(\boldsymbol{\nabla} p_1). \end{equation}
It is known that the first-order pressure (
$p_1$) in the domain in our study is of the order of 
$\sim 10^6$, 
$\rho _0 \sim {O}(10^3)$, wavenumber 
$k\sim {O}(10^4)$ and 
$\omega \sim {O}(10^6)$ (the time derivative can be written in terms of the angular frequency as 
$-{\rm i} \omega$, where 
${\rm i}$ represents a complex root)
\begin{equation} {\boldsymbol{v}}_1\sim\frac{k p_1}{\rho_0 \,{\rm i}\omega}\sim\frac{1}{0.76}\sim 1. \end{equation}
In the above simplification, wavelength 
$\lambda$ is considered as 
$760\ \mathrm {\mu }{\rm m}$ (it is chosen as twice the width of the channel in order to get a one-dimensional half-standing wave along this dimension of the channel). Using this order of 
${\boldsymbol {v}}_1$ in ((2.3c) and (2.3d)), we have
$$\begin{gather} s_1\sim\frac{{\boldsymbol{v}}_1k s_0}{{\rm i}\omega}\sim10^{{-}2}(s_0), \end{gather}$$
$$\begin{gather}\rho_1\sim\frac{{\boldsymbol{v}}_1 k\rho_0}{{\rm i}\omega}\sim10^{{-}4}(\rho_0). \end{gather}$$
Thus, it is conclusive that 
$s_1\ll s_0$ and 
$\rho _1\ll \rho _0$.
On simplifying to the second-order equations, we have stated that 
${O}(\beta \eta \boldsymbol {\nabla }({1}/{\rho _0 c_0^2})\langle ({\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla })p_1\rangle )\ll {O}(\boldsymbol {\nabla }\boldsymbol {\cdot }\langle \rho _0{\boldsymbol {v}}_1{\boldsymbol {v}}_1\rangle )$ and the proof for this follows. Using the order of magnitudes as 
$\beta \sim {O}(10^{-1})$, 
$\eta \sim {O}(10^{-3})$, 
$\delta \sim {O}(10^{-6})$, 
$c\sim {O}(10^3)$ and the magnitudes for other variables as considered above, we have
$$\begin{gather} \beta\eta\boldsymbol{\nabla}\left(\frac{1}{\rho_0 c_0^2}\langle({\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla})p_1\rangle\right)\sim\frac{\beta\eta}{\delta} \left(\frac{1}{\rho_0 c_0^2}({\boldsymbol{v}}_1 k)p_1\right) \sim\frac{10^{{-}1}10^{{-}3}}{10^{{-}6}}\left(\frac{10^4 10^6}{10^9}\right)\sim 10^3, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1{\boldsymbol{v}}_1\rangle \sim k \rho_0{\boldsymbol{v}}_1{\boldsymbol{v}}_1\sim 10^4 10^3 \sim 10^7. \end{gather}$$Thus, 
$\beta \eta \boldsymbol {\nabla }({1}/{\rho _0 c_0^2})\langle ({\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla })p_1\rangle$ is 
$10^4$ orders less than 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\langle \rho _0{\boldsymbol {v}}_1{\boldsymbol {v}}_1\rangle$.
We have also stated that the first-order fields in the terms 
$\langle {\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla } s_1\rangle$, 
$\langle \rho _1\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_1\rangle$ are out of phase. The below is the proof for this. From (C1a) and 
${\boldsymbol {v}}_1={\boldsymbol {v}}_a(r,\tau ))\, {\rm e}^{-{\rm i}\omega t_f}$ (first-order time-harmonic acoustic field), we have
$$\begin{gather} -{\rm i}\omega\rho_1=\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0{\boldsymbol{v}}_1)= \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0{\boldsymbol{v}}_a(r,\tau))\, {\rm e}^{-{\rm i}\omega t_f}, \end{gather}$$
$$\begin{gather}\rho_1=Re\left[\frac{{\rm i} \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0{\boldsymbol{v}}_a(r,\tau))\, {\rm e}^{-{\rm i}\omega t_f}}{\omega}\right], \end{gather}$$
where 
${\boldsymbol {v}}_a$ is the velocity amplitude. The complex ‘
${\rm i}$’ in the 
$\rho _1$ expression indicates a phase shift between the velocity and density oscillations. Thus, 
$\rho _1$ and 
${\boldsymbol {v}}_1$ are out of phase (Muller & Bruus Reference Muller and Bruus2014). Similarly, from (C1c), we have
$$\begin{gather} -{\rm i}\omega s_1={-}{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} {s}_0={-}{\boldsymbol{v}}_a(r,\tau)\, {\rm e}^{-{\rm i}\omega t_f}\boldsymbol{\cdot}\boldsymbol{\nabla} s_0, \end{gather}$$
$$\begin{gather}s_1=Re\left[\frac{-{\rm i}{\boldsymbol{v}}_a(r,\tau)\, {\rm e}^{-{\rm i}\omega t_f}\boldsymbol{\cdot}\boldsymbol{\nabla} s_0}{\omega }\right]. \end{gather}$$
Thus 
${\boldsymbol {v}}_1$ and 
$s_1$ are out of phase.
Appendix B. Derivation of governing equations for inhomogeneous fluids – alternate approach
The hydrodynamics of inhomogeneous fluids considered in this study is governed by the mass continuity, momentum and advection–diffusion equations, which are given below
$$\begin{gather} \partial_t \rho+\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho {\boldsymbol{v}})=0, \end{gather}$$
$$\begin{gather}\rho[\partial_t {\boldsymbol{v}} + ({\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla}){\boldsymbol{v}}]={-}\boldsymbol{\nabla} p +\eta\nabla^2 {\boldsymbol{v}} +\beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}) + \rho{\boldsymbol{g}}, \end{gather}$$
$$\begin{gather}\partial_t {s}+{\boldsymbol{v}}\boldsymbol{\cdot}\boldsymbol{\nabla} {s} = D\nabla^2{s}, \end{gather}$$
where 
$\rho$ is the density, 
${\boldsymbol {v}}$ is the velocity, 
$p$ is the pressure, 
$\eta$ is the dynamic viscosity of the fluid, 
$\xi$ is the volume fluid viscosity, 
$\beta =(\xi /\eta )+(1/3)$, 
$s$ is the solute concentration and 
$D$ is the diffusivity. When the fluid is subjected to acoustic waves, the following thermodynamic pressure density relation is required,
\begin{equation} \frac{{\rm d} \rho}{{\rm d} t}=\frac{1}{c^2}\frac{{\rm d} p}{{\rm d} t}. \end{equation}Say the fields are decomposed based on their variation on different time scales
\begin{equation} f=f_0({\boldsymbol{r}},\tau)+f_1({\boldsymbol{r}},\tau)\, {\rm e}^{-{\rm i}\omega t_f}, \end{equation}
where 
$f_0$ varies due to gravity and nonlinear effects on the slow time scale and 
$f_1$ varies on the fast time scale. The time scale of acoustics is of the order of microseconds (
$t_f\sim {1}/{\omega } \sim 0.1\ \mathrm {\mu }{\rm s}$), whereas the hydrodynamic time scale of the order of milliseconds (
$t_f \ll \tau$).
Substituting (B2) in the governing equations, we have
\begin{gather} \partial_t (\rho_0+\rho_1) +\boldsymbol{\nabla} \boldsymbol{\cdot}[(\rho_0+\rho_1) ({\boldsymbol{{v}}}_0+{\boldsymbol{{v}}}_1)]=0, \end{gather}
\begin{gather} \hspace{-8pc}(\rho_0+\rho_1) \partial_t ({\boldsymbol{{v}}}_0+{\boldsymbol{v}}_1)+(\rho_0+\rho_1)(({\boldsymbol{v}}_0 +{\boldsymbol{v}}_1)\boldsymbol{\cdot}\boldsymbol{\nabla})({\boldsymbol{v}}_0+{\boldsymbol{v}}_1)\nonumber\\ \hspace{1pc}={-}\boldsymbol{\nabla} (p_0+p_1) + \eta\nabla^2({\boldsymbol{{v}}}_0+ {\boldsymbol{v}}_1)+\beta\eta \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}({\boldsymbol{{v}}}_0+{\boldsymbol{v}}_1))+(\rho_0+\rho_1)g, \end{gather}
$$\begin{gather} \partial_t(s_0+s_1)+({\boldsymbol{v}}_0+{\boldsymbol{v}}_1)\boldsymbol{\cdot}\boldsymbol{\nabla} (s_0+s_1) =D\nabla^2 (s_0+s_1), \end{gather}$$
$$\begin{gather}\partial_t(\rho_0+\rho_1)+(({\boldsymbol{v}}_0+{\boldsymbol{v}}_1) \boldsymbol{\cdot}\boldsymbol{\nabla})(\rho_0+\rho_1)=(1/c^2)[\partial_t(p_0+p_1)]. \end{gather}$$It is important to mention here that the pressure field on the slow time scale is decomposed as, 
$p_0=p_c+p_g+p_s$, where 
$p_c$ is the reference pressure, 
$p_g$ is the pressure due to gravity and 
$p_s$ is the second-order mean Eulerian pressure. Neglecting the influence of gravity and accounting for the fact that the reference pressure remains constant, the variation in 
$p_0$ is only due to the variation in 
$p_s$.
Before we proceed to the fast time scale equations, we propose the validity of this theory, i.e. it holds good only when 
${\boldsymbol {{v}}}_0\ll {\boldsymbol {{v}}}_1$ and 
$\boldsymbol {\nabla } p_0\ll \boldsymbol {\nabla } p_1$, which is true in our case. On the fast time scale (
$t_f$), the variation of the 
$f_0$ fields with respect to time can be neglected and since 
$D\sim {O}(10^{-9})$, the diffusion term is negligible, due to which the composition of any given fluid particle will remain unchanged as it moves. Also, from (B3c) and (B3d), it is evident that 
$s_1\ll s_0$ and 
$\rho _1\ll \rho _0$. Accounting for all the above arguments, the fast time scale equations reduce to
$$\begin{gather} \partial_t\rho_1+\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0{\boldsymbol{{v}}}_1)=0, \end{gather}$$
$$\begin{gather}\rho_0\partial_t{\boldsymbol{v}}_1={-}\boldsymbol{\nabla} p_1+ \eta\nabla^2{\boldsymbol{v}}_1+\beta\eta \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1), \end{gather}$$
$$\begin{gather}\partial_ts_1+{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} s_0=0, \end{gather}$$
$$\begin{gather}\partial_t\rho_1+{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_0=(1/c^2)[\partial_tp_1]. \end{gather}$$ However, the time average of these 
$f_0$ fields over one complete oscillation is zero, thus these fields cannot cause any bulk fluid motion. But the N-S equation is nonlinear and the above linearized equation is not exact. Thus, we must proceed to solve the fields on the slow time scale, which is responsible for the nonlinear effects. Since 
$\rho _1 \ll \rho _0$ and 
${\boldsymbol {{v}}}_0\ll {\boldsymbol {{v}}}_1$, the fast time-averaged, slow time scale equations reduce to
$$\begin{gather} \langle\partial_t\rho_0\rangle +\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{{v}}}_0+ \rho_1{\boldsymbol{{v}}}_1\rangle=0, \end{gather}$$
$$\begin{gather}\rho_0\langle\partial_t {\boldsymbol{{v}}}_0\rangle+\rho_1\langle\partial_t {\boldsymbol{{v}}}_1\rangle+\rho_0\langle({\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla}){\boldsymbol{v}}_1\rangle={-}\boldsymbol{\nabla} \langle p_0\rangle + \eta\nabla^2\langle{\boldsymbol{{v}}}_0\rangle+\beta\eta \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\langle{\boldsymbol{{v}}}_0\rangle)+\rho_0 g, \end{gather}$$
$$\begin{gather}\langle\partial_ts_0\rangle+\langle{\boldsymbol{v}}_0\boldsymbol{\cdot}\boldsymbol{\nabla} s_0\rangle+\langle{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} s_1\rangle=D\nabla^2 \langle s_0\rangle, \end{gather}$$
$$\begin{gather}\langle\partial_t\rho_0\rangle+\langle{\boldsymbol{v}}_0\boldsymbol{\cdot} \boldsymbol{\nabla}\rho_0\rangle+\langle{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} \rho_1\rangle=(1/c^2)\langle{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} p_1\rangle, \end{gather}$$
where 
$\langle \cdots \rangle$ denotes the time average for one oscillation period. From (B5b), it can be seen that the product of the fast time scale varying fields along with gravity act as source terms responsible for second-order fields, which cause a change in the fluid profile, due to which the amplitudes of fast scale fields also change accordingly at all time instants on the slow scale. Also, in inhomogeneous fluids, changing the background velocity field 
$v_0$, results in variation of 
$\rho _0$ and 
$s_0$ with respect to the slow time scale (
$t\sim \tau$), thus 
$\langle \partial _t\rho _0\rangle \neq 0$ and 
$\langle \partial _t s_0\rangle \neq 0$. Along with 
$\rho _{1} {g} \ll \rho _{0} {g}$, 
$D\boldsymbol {\nabla } s_{1} \ll D\boldsymbol {\nabla } s_{0}$, we also neglect 
$\langle {\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla } s_1\rangle$, 
$\langle \rho _1\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_1\rangle$ since both these terms are out of phase. Using the above arguments, combining equations (B5a) and (B5d) we get
\begin{equation} \langle\rho_0\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1\rangle={-}(1/c^2) \langle{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} p_1\rangle. \end{equation}
Substituting the above equation in (B6), since 
$\beta \eta \boldsymbol {\nabla }({1}/{\rho _0c^2})({\boldsymbol {v}}_1\boldsymbol {\cdot } \boldsymbol {\nabla })p_1\sim {O}(10^3)\ll \boldsymbol {\nabla }\boldsymbol {\cdot }\langle \rho _0{\boldsymbol {v}}_1 {\boldsymbol {v}}_1\rangle \sim {O}(10^7)$ we can neglect the term 
$\beta \eta \boldsymbol {\nabla }(\boldsymbol {\nabla }\boldsymbol {\cdot }\langle {\boldsymbol {v}}_1 \rangle )$. The term 
$(1/c^2)\langle {\boldsymbol {v}}_1\boldsymbol {\cdot }\boldsymbol {\nabla } p_1\rangle$ does not contribute to any acoustic phenomena and can be neglected, thus we have 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\langle {\boldsymbol {v}}_1\rangle =0$, which indicates the flow is incompressible. This is stated as Nyborg's divergence-free approximation, which states that compressibility effects can be safely ignored. In microscale flows, the term 
$\rho _0\langle \partial _t {\boldsymbol {{v}}}_0\rangle$ in (B5b) can be neglected. Consequently, our governing equations reduce to
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\langle{\boldsymbol{v}}_0\rangle=0, \end{gather}$$
$$\begin{gather}-\boldsymbol{\nabla}\langle p_0 \rangle+\eta\nabla^2 \langle{\boldsymbol{v}}_0 \rangle-\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1 {\boldsymbol{v}}_1 \rangle+\langle\rho_{0}{\boldsymbol{g}}\rangle = 0, \end{gather}$$
$$\begin{gather}\langle\partial_t {s}_0\rangle+\langle{\boldsymbol{v}}_0\boldsymbol{\cdot} \boldsymbol{\nabla} {s}_0\rangle = D\nabla^2\langle{s}_0\rangle. \end{gather}$$Hence, the body force due to acoustic fields is construed as
\begin{equation} \boldsymbol{f}_{ac}={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0{\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle. \end{equation}Appendix C. First-order equations in frequency domain
The first-order equations, which give rise to the first-order fields, ((2.3) in the main text) are specified to be solved in the frequency domain on the fast time scale. As we are solving the governing equations at a particular frequency for all acoustic phenomena, they can be rewritten as
$$\begin{gather} -{\rm i}\omega\rho_1=-\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0{\boldsymbol{v}}_1), \end{gather}$$
$$\begin{gather}-{\rm i}\omega\rho_0{\boldsymbol{v}}_1={-}\boldsymbol{\nabla} p_1 + \eta\nabla^2 {\boldsymbol{v}}_1+\beta\eta \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1), \end{gather}$$
$$\begin{gather}-{\rm i}\omega s_1+{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla} {s}_0=0, \end{gather}$$
$$\begin{gather}-{\rm i}\omega\rho_0\kappa_0p_1={-}{\rm i}\omega\rho_1+{\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_0. \end{gather}$$Also, combining (C1a) and (C1d) we get
\begin{equation} -{\rm i}\omega\kappa_0 p_1={-}\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1,\end{equation}
where 
$p_1$ is the first-order pressure field, 
$\rho _1$ refers to the first-order density field, 
${\boldsymbol {v}}_1$ is the first-order velocity field, 
$\omega$ is the angular frequency, 
$\eta$ is the dynamic viscosity of the fluid, 
$\xi$ is the volume fluid viscosity, 
$\beta =(\xi /\eta )+(1/3)$, 
$s$ is the solute concentration and 
$D$ is the diffusivity. In the time-dependent domain (C1e) can be written as 
$\kappa _0\partial _t p_1=-\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_1$.
Appendix D. Variation of first-order and second-order fields for different configurations
In support of the analytical calculations and numerical simulations for the theory proposed, the profiles of the first-order velocity and pressure fields are shown below.
From figure 6 it is observed that the first-order fields (
${\boldsymbol {v}}_1$ and 
$p_1$) are highly configuration dependent, thus they change with respect to time, which supports the fact that the assumption of 
$E_{ac}$ being constant is incorrect. In figure 7, it is seen that the magnitude of the second-order velocity 
$\langle {\boldsymbol {v}}_2\rangle$ is very high as the fluids evolve by an advection process to settle to a stable configuration. Also, it is noticed that the second-order pressure 
$\langle p_2\rangle$ is also significantly high. This is because of the large magnitude of 
$p_1$ at this instant that increases 
$E_{ac}$ drastically, which in turn affects 
$\langle p_2\rangle$ and 
$\langle {\boldsymbol {v}}_2\rangle$. The second-order streaming velocity in this study refers to the time-averaged Eulerian streaming velocity. It can also be better represented by the Lagrangian streaming velocity. The Eulerian and Lagrangian velocities are related using Stokes drift (Hamilton, Ilinskii & Zabolotskaya Reference Hamilton, Ilinskii and Zabolotskaya2003; Pavlic & Dual Reference Pavlic and Dual2021).
Figure 6. (a) Impedance profile for relocation of an unstable configuration to a stable configuration with 
$d=0.261$ nm and 
$\nu=1.96$ MHz. Initially (t=0s), there is low impedance DI water (blue) at the centre and high impedance 10 % Ficoll PM70 (red) at the sides. (b) Corresponding first-order velocity 
$v_1$. (c) Corresponding first-order pressure 
$p_1$. The microchannel width 
$w=380\ \mathrm {\mu }{\rm m}$ and height 
$h=160\ \mathrm {\mu }{\rm m}$.
Figure 7. Corresponding (a) second-order velocity (
$\langle {\boldsymbol {v}}_2\rangle$) and (b) second-order pressure (
$\langle p_2\rangle$) profile at 
$t=0.64$ s for the fluid configuration considered in figure (1).
From figure 8, it is evident that the first-order velocity has a large variation within the thin viscous boundary layer due to the no-slip condition whereas the first-order pressure remains the same inside and outside the boundary layer. Thus the acoustic standing wave relation 
$p_1=p_a \sin (kx)$ and 
${\boldsymbol {v}}_1=({p_a}/{i\rho _0 c_0})\cos (kx)$ (where 
$k= 2{\rm \pi} /\lambda$ is wavenumber and 
$\lambda \approx 2w$) holds good only outside the boundary layer.
Figure 8. (a) Variation of first-order velocity at t=0 s near the boundary layer. Plots for variation of (b) first-order pressure 
$p_1$ and (c) first-order velocity 
$v_1$ in the direction of height.
Figure 9 shows the relocation of immiscible fluids under acoustic fields. Here also, the first-order fields (
${\boldsymbol {v}}_1$ and 
$p_1$) are highly configuration dependent. Due to the application of a two-dimensional standing wave, the pressure nodes are not formed at the centre of the channel as a vertical line (figure 9c) unlike for miscible fluids (figure 6c).
Figure 9. Results of relocation of immiscible fluids by actuating all walls of a microchannel (two-dimensional), consisting of high impedance mineral oil (red) and low impedance silicone oil (blue) with a surface tension of 
$1$ mN m
$^{-1}$. The wall is actuated with a displacement 
$d= 15$ nm at a frequency of 2.1 MHz. (a) Impedance profile. (b) First-order velocity profile 
${\boldsymbol {v}}_1$. (c) First-order pressure profile 
$p_1$. The microchannel is of width 
$w=360\ \mathrm {\mu }{\rm m}$ and height 
$h=160\ \mathrm {\mu }{\rm m}$.
Appendix E. Fluid interface at the node for inhomogeneous fluids
In order to study the effect of the acoustic body force term alone, the influence of gravity has been excluded to avoid stratification. From figure 10, we see that acoustic relocation does not occur in the case of inhomogeneous fluids when the fluid interface is at the node. Acoustic streaming can disturb the interface location, as shown at 
$t=1.4$ s, after which the acoustic body force tries to stabilize the configuration.
Figure 10. No relocation of inhomogeneous fluid when the interface is at node.
Appendix F. Analysing the body force term 
${f}_{ac}$ in detail
The body force 
${\boldsymbol {f}}_{ac}=-\boldsymbol {\nabla }\boldsymbol {\cdot }\langle \rho _0{\boldsymbol {v}}_1 {\boldsymbol {v}}_1\rangle$ derived in the main text is analysed in detail by splitting it as
\begin{equation} -\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0 {\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle={-}\rho_0\langle{\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{v}}_1\rangle-\langle\rho_1 \partial_t {\boldsymbol{v}}_1\rangle. \end{equation} Consider the first term 
$-\rho _0\langle {\boldsymbol {v}}_1 \boldsymbol {\cdot }\boldsymbol {\nabla }{\boldsymbol {v}}_1\rangle$. Using the identity 
${\boldsymbol {A}}$, 
${\boldsymbol {A}}\boldsymbol {\cdot }\boldsymbol {\nabla }{\boldsymbol {A}}=\boldsymbol {\nabla }({\boldsymbol {A}}^2/2)-{\boldsymbol {A}}\times (\boldsymbol {\nabla }\times {\boldsymbol {A}})$, it can be written as
$$\begin{gather} -\rho_0\langle{\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{v}}_1\rangle={-}\frac{1}{2}\rho_0\boldsymbol{\nabla} \langle|{\boldsymbol{v}}_1^2|\rangle+\rho_0\langle{\boldsymbol{v}}_1 \times(\boldsymbol{\nabla}\times{\boldsymbol{v}}_1)\rangle, \end{gather}$$
$$\begin{gather}-\rho_0\langle{\boldsymbol{v}}_1 \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{v}}_1\rangle=\frac{-1}{2\rho_0}\boldsymbol{\nabla}\langle|\rho_0 {\boldsymbol{v}}_1|^2\rangle+\langle|{\boldsymbol{v}}_1|^2\rangle\boldsymbol{\nabla}\rho_0+ \rho_0\langle{\boldsymbol{v}}_1\times(\boldsymbol{\nabla}\times{\boldsymbol{v}}_1)\rangle. \end{gather}$$ Now consider the second term on the right-hand side of (F1), 
$-\langle \rho _1 \partial _t {\boldsymbol {v}}_1\rangle$. Using 
$\partial _t(\rho _1 {\boldsymbol {v}}_1)=0$, we have 
$-\langle \rho _1 \partial _t {\boldsymbol {v}}_1\rangle =\langle {\boldsymbol {v}}_1\partial _t\rho _1\rangle$ and from the first-order continuity equation, we have 
$\partial _t\rho _1=-\boldsymbol {\nabla }\boldsymbol {\cdot }(\rho _0 {\boldsymbol {v}}_1)$. Substituting these and using the identity 
$\boldsymbol {\nabla }\boldsymbol {\cdot }(AB)=B\boldsymbol {\cdot }\boldsymbol {\nabla }(A)+A\boldsymbol {\nabla }\boldsymbol {\cdot }(B)$, (where 
$A$ and 
$B$ are two vectors), we get
\begin{equation} -\langle{\boldsymbol{v}}_1\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0 {\boldsymbol{v}}_1)\rangle={-}\langle{\boldsymbol{v}}_1( {\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_0)\rangle-\langle \rho_0{\boldsymbol{v}}_1\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1\rangle. \end{equation} From 
$\kappa _0\partial _t p_1=-\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_1$, the second term on the right-hand side of (F3) can be simplified as 
$-\langle \rho _0{\boldsymbol {v}}_1\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_1\rangle =\langle \rho _0{\boldsymbol {v}}_1(\kappa _0\partial _t p_1)\rangle =-\langle \kappa _0 p_1(\rho _0\partial _t{\boldsymbol {v}}_1)\rangle$ (since 
$\langle f_1(\partial _t g_1)\rangle =-\langle g_1(\partial _t f_1)\rangle$ is valid for time-harmonic fields 
$f_1$ and 
$g_1$). Now, using the first-order momentum equation (
$\rho _0\partial _t{\boldsymbol {v}}_1=-\boldsymbol {\nabla } p_1 +\eta \nabla ^2{\boldsymbol {v}}_1+\beta \eta \boldsymbol {\nabla }(\boldsymbol {\nabla }\boldsymbol {\cdot }{\boldsymbol {v}}_1)$), (F3) reduces to
$$\begin{gather} -\langle{\boldsymbol{v}}_1\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0 {\boldsymbol{v}}_1)\rangle =\langle-{\boldsymbol{v}}_1({\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_0) \rangle-\langle\kappa_0 p_1(-\boldsymbol{\nabla} p_1 +\eta\nabla^2{\boldsymbol{v}}_1 +\beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle, \end{gather}$$
$$\begin{gather}-\langle{\boldsymbol{v}}_1\boldsymbol{\nabla}\boldsymbol{\cdot}(\rho_0 {\boldsymbol{v}}_1)\rangle =\tfrac{1}{2}\kappa_0\boldsymbol{\nabla}\langle|p_1^2|\rangle-\langle{\boldsymbol{v}}_1( {\boldsymbol{v}}_1\boldsymbol{\cdot}\boldsymbol{\nabla}\rho_0)\rangle-\langle\kappa_0 p_1(\eta\nabla^2{\boldsymbol{v}}_1+\beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle. \end{gather}$$
\begin{gather} {\boldsymbol{f}}_{ac}=\frac{1}{2}\kappa_0\boldsymbol{\nabla}\langle|p_1^2|\rangle -\frac{1}{2\rho_0}\boldsymbol{\nabla}\langle|\rho_0 {\boldsymbol{v}}_1|^2\rangle+\langle|{\boldsymbol{v}}_1|^2 \rangle\boldsymbol{\nabla}\rho_0+\langle{\boldsymbol{v}}_1\times \boldsymbol{\nabla}\times(\rho_0{\boldsymbol{v}}_1)\rangle\nonumber\\ \hspace{-3pc} -\langle{\boldsymbol{v}}_1({\boldsymbol{v}}_1\boldsymbol{\cdot} \boldsymbol{\nabla}\rho_0)\rangle-\langle\kappa_0p_1(\eta\nabla^2{\boldsymbol{v}}_1+ \beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle, \end{gather}
\begin{gather} =\frac{1}{2}\kappa_0\boldsymbol{\nabla}\langle|p_1^2|\rangle-\frac{1}{2\rho_0}\boldsymbol{\nabla}\langle|\rho_0 {\boldsymbol{v}}_1|^2\rangle+\langle{\boldsymbol{v}}_1\times \boldsymbol{\nabla}\times(\rho_0{\boldsymbol{v}}_1)\rangle-\langle \kappa_0p_1(\eta\nabla^2{\boldsymbol{v}}_1+\beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle,\end{gather}
\begin{gather} \hspace{-7pc}=\tfrac{1}{2}\boldsymbol{\nabla}(\kappa_0\langle |p_1|^2\rangle-\rho_0 \langle|{\boldsymbol{v}}_1|^2\rangle) +\{\langle{\boldsymbol{v}}_1\times\boldsymbol{\nabla}\times(\rho_0{\boldsymbol{v}}_1)\rangle\}\nonumber\\ -\tfrac{1}{2}[ \langle|p_1|^2\rangle\boldsymbol{\nabla}\kappa_0+\langle|{\boldsymbol{v}}_1|^2 \rangle\boldsymbol{\nabla}\rho_0]-\langle(\kappa_0p_1)(\eta\nabla^2{\boldsymbol{v}}_1+ \beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle. \end{gather} For an acoustic wave propagating in the medium with frequency in the range of MHz, the orders of magnitude of individual parameters are wavenumber 
$k\sim {O}(10^4), \rho _0\sim {O}(10^3)$, 
$\delta \sim {O}(10^{-6})$, 
$c_0\sim {O}(10^3)$, 
$\eta \sim {O}(10^{-3})$ and 
$p_1\sim \rho _0 c_0 {\boldsymbol {v}}_1$ (where 
$u_1$ and 
$v_1$ denote 
$x$ and 
$y$ components of first-order velocity). The order of magnitude analysis results in, outside the boundary layer,
$$\begin{gather} \frac{1}{2}\boldsymbol{\nabla}(\kappa_0\langle |p_1|^2\rangle-\rho_0 \langle|{\boldsymbol{v}}_1|^2\rangle\sim \frac{\rho_0 k u_1^2}{2}\sim{O}(10^6)(u_1^2), \end{gather}$$
$$\begin{gather}-\langle(\kappa_0p_1)(\eta\nabla^2{\boldsymbol{v}}_1+\beta\eta \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle\sim \frac{\eta k^2 u_1^2}{c_0}\sim{O}(10^2)(u_1^2). \end{gather}$$Inside the boundary layer,
$$\begin{gather} \{\langle{\boldsymbol{v}}_1\times\boldsymbol{\nabla}\times(\rho_0{\boldsymbol{v}}_1) \rangle\}\sim\frac{\rho_0u_1^2}{\delta}\sim{O}(10^9)(u_1^2), \end{gather}$$
$$\begin{gather}-\langle(\kappa_0p_1)(\eta\nabla^2{\boldsymbol{v}}_1+\beta\eta \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{v}}_1))\rangle\sim \frac{\eta u_1^2} {c_0 \delta^2}\sim{O}(10^6)(u_1^2). \end{gather}$$From the above arguments, we can conclude that the last term in (F5c) is insignificant inside as well as outside the boundary layer on comparing the magnitudes of acoustic parameters with those of other terms, thus we can certainly neglect this term to proceed further. Hence, the equation reduces to
\begin{align} {\boldsymbol{f}}_{ac}&={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\langle\rho_0 {\boldsymbol{v}}_1 {\boldsymbol{v}}_1\rangle=\tfrac{1}{2}\boldsymbol{\nabla}(\kappa_0\langle |p_1|^2\rangle-\rho_0 \langle|{\boldsymbol{v}}_1|^2\rangle)+ \{\langle{\boldsymbol{v}}_1\times\boldsymbol{\nabla}\times(\rho_0{\boldsymbol{v}}_1)\rangle\}\nonumber\\ &\quad -\tfrac{1}{2}[ \langle|p_1|^2\rangle\boldsymbol{\nabla}\kappa_0 +\langle|{\boldsymbol{v}}_1|^2\rangle\boldsymbol{\nabla}\rho_0]. \end{align}Appendix G. Mesh convergence and numerical uncertainty
The mesh convergence analysis is performed as shown in figures 11 and 12 to determine the number of grids large enough for all dependent variables to converge, ensuring a mesh-independent result. The dark blue line, which indicates a mesh constituting of 24 794 triangular grid elements (maximum element size of about 2.655 
$\mathrm {\mu }$m), is chosen as a trade-off between computational time and accuracy.
Figure 11. Mesh convergence analysis for first-order pressure amplitude 
$p_a$.
Figure 12. Mesh convergence analysis for maximum second-order velocity 
$\langle {\boldsymbol {v}}_2\rangle$.
The discretization error and numerical uncertainty for this mesh is estimated using the Richardson extrapolation method. The discretization error is given as (Freitas Reference Freitas2002; Phillips & Roy Reference Phillips and Roy2014)
\begin{equation} \varepsilon\approx[(f_f-f_c)/f_f]/[r^p-1], \end{equation}
where 
$p=\ln [(f_m-f_c)/(f_f-f_m)]/ \ln [r]$ is the order of the convergence rate, subscript 
$f$ refers to the fine grid solution (24 794 elements), 
$c$ refers to the course grid solution (6070 elements), 
$m$ refers to the medium grid solution (12 308 elements) and 
$r$ is the ratio of grid spacings. For numerical solutions on the fine grid, grid convergence indices for the pressure amplitude and second-order velocity are calculated to be 1.028 % and 2.793 % respectively. Thus, the numerical uncertainty for pressure amplitude is 
$p_a\pm 0.01p_a$, and the second-order velocity is 
${\boldsymbol {v}}_2\pm 0.027{\boldsymbol {v}}_2$.
