# Eckart streaming with nonlinear high-order harmonics: an example at
gigahertz

Shiyu Li  State Key Laboratory of Ocean Engineering, School of Ocean and Civil
Engineering, Shanghai Jiao Tong University, Shanghai, 200240, China  Weiwei
Cui  State Key Laboratory of Precision Measuring Technology and Instruments,
Tianjin University, Tianjin 300072, China  Thierry Baasch  Lund University
Bin Wang  State Key Laboratory of Ocean Engineering, School of Ocean and Civil
Engineering, Shanghai Jiao Tong University, Shanghai, 200240, China  Key
Laboratory of Marine Intelligent Equipment and System, Ministry of Education,
China  Zhixiong Gong  Corresponding author<EMAIL_ADDRESS>URL:
https://en.naoce.sjtu.edu.cn/teachers/GongZhixiong.html State Key Laboratory
of Ocean Engineering, School of Ocean and Civil Engineering, Shanghai Jiao
Tong University, Shanghai, 200240, China  Key Laboratory of Marine Intelligent
Equipment and System, Ministry of Education, China

###### Abstract

Acoustic streaming shows great potential in applications such as bubble
dynamics, cell aggregation, and nano-sized particle isolation in the
biomedical and drug industries. As the acoustic shock distance decreases with
the increase of incident frequency, the nonlinear propagation effect will play
a role in acoustic streaming, e.g., Eckart (bulk) streaming at a few gigahertz
(GHz). However, the theory of source terms of bulk streaming is still missing
at this stage when high-order acoustic harmonics play a role. In this paper,
we derive the source term including the contribution of higher-order
harmonics. The streaming-induced hydrodynamic flow is assumed to be
incompressible and no shock wave occurs during the nonlinear acoustic
propagation as restricted by the traditional Goldberg number $\Gamma<1$ or
$\Gamma\approx 1$ which indicates the importance of nonlinearity relative to
dissipation. The derived force terms allow evaluating bulk streaming with
high-order harmonics at GHz and provide an exact expression compared to the
existing empirical formulas. Numerical results show that the contribution of
higher-order harmonics increases the streaming flow velocity by more than
$20\%$. We show that the expression introduced by Nyborg should be avoided in
numerical computations as it includes part of the acoustic radiation force
that does not lead to acoustic streaming.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

Gigahertz (GHz) acoustics have recently been used in experiments for
nanoparticle trapping, enrichment, and separation based on the acoustic
streaming effect [1, 2, 3]. Acoustic manipulation in the gigahertz range has
potential applications including nano-sized biosensors [2], nanoliter
microreactors [4], and microfluid jet producer [5]. However, the study of GHz
streaming is just at the beginning due to the challenges of fabrication
techniques of ultrahigh frequency resonators [6] and the huge computational
costs with direct numerical simulations at such small wavelengths especially
in three dimensions [7]. For a typical GHz tweezer, the wavelength is 1.5
$\mu$m at the frequency of 1 GHz in water, which is much smaller than typical
microchannel sizes, e.g., a few tens or hundreds of micrometers. In addition,
it is easy to induce nonlinear propagation at GHz since the acoustic shock
distance depends on the working frequency and the large vibration velocity on
the transducer surface. These make the theoretical and experimental studies of
GHz acoustical tweezers more challenging the same technique in the frequency
regime of megahertz (MHz).

The sound waves are produced based on mechanical vibration and are coupled
into the fluid medium for microfluidic applications. The highest efficiency of
energy conversion from electrical to mechanical energies typically occurs at
the mechanical resonance, which depends on the desired acoustic wavelength (or
driving frequency). The wavelengths of typical GHz transducers in the coupling
medium range from around 100 nanometers to a few micrometers, leading to the
difficulties of device fabrication. Typical piezoelectric transducers (PZT)
use the vibration of planar sources to produce acoustics with a frequency
regime of 1 - 10 MHz in the field of acoustofluidics. However, it is a
challenge to fabricate PZT at the thickness of nanometers for the frequency at
GHz since the resonance depends on the selected vibration mode depending on
the piezo thickness. Considering interdigitated transducers [8] at GHz, the
distances between electrode fingers are too small to fabricate with the
commonly-used fabrication process and cannot withstand high power [6]. This is
partly solved with the successful fabrication of the high-tone bulk acoustic
resonators on four substrates with very high Q factor (up to 48000) at 1GHz
[9]. Then, Cui et al. combined the fabrication of FBAR (Film Bulk Acoustic
Wave Resonator) technique with micro/nanofluidics and developed the field of
GHz acoustofluidics [1].

Compared with the recent development of fabrication techniques and
experimental demonstrations of GHz acoustical tweezers, the theories of
acoustic bulk streaming at such high frequencies are not well studied. For
most of the published experimental works at GHz, empirical expressions of the
source term of acoustic streaming are used in the numerical simulations, and
no one seems to verify the streaming flow velocities between the simulation
and experiment results [2, 3]. Indeed, it is easy to understand that the
streaming-induced flow patterns are similar in the confined microchannels
because of the mass conservation of the steady flow. Since the seminal works
of Eckart [10], Nyborg [11, 12] and Lighthill [13], the systematic theories of
the bulk (Eckart) streaming are built and developed. A good historical
perspective of bulk acoustic streaming can be found in Ref. [14]. It should be
noteworthy that Nyborg derived an expression of the source term of bulk
streaming to solve the Stokes equation for the hydrodynamic flow velocity [11,
12] which is widely used in the research community of this field. However, as
pointed out by Lighthill, the source term by Nyborg contains a gradient term
that makes a contribution to the acoustic radiation pressure instead of
streaming [13, 14]. The source term of acoustic streaming was recast recently
by Riaud et al. specifying the sole source of bulk streaming without the
gradient of acoustic Lagrangian [15]. In their work, they consider the bulk
streaming inside sessile droplets of size 1 mm under the activation of surface
acoustic waves at a frequency of around 20 MHz neglecting the nonlinear
propagation since Gol’dberg number [defined in Eq. (1) below] is much smaller
than 1 and the droplet size is much smaller than the shock wave distance.
Hence, the source term is limited to linear propagation when there are no
high-order acoustic harmonics. However, the nonlinear effect of acoustic
propagation can not be neglected in the frequency regime of GHz since the
viscous dissipation is remarkable and the dimensionless Gol’dberg number is
comparable to the unit. This is the case for the recent experimental works at
GHz [1].

In this work, we revisit the source term of bulk acoustic streaming with two
assumptions: (i) acoustics are rotational, and (ii) the streaming-induced
steady hydrodynamic flow is incompressible. Only weakly nonlinear effects with
harmonic waves are taken into consideration and the shock wave is outside of
scope since it will make the second assumption fail. Both theoretical and
numerical examples are proposed to illustrate the issue of Nyborg’s expression
which should be avoided as shown in section III. More importantly, based on
the peculiar characteristics of bulk streaming at GHz, the nonlinear effect of
acoustic propagation will be considered and a theoretical source term for this
situation is provided in terms of pressure fields of different orders of
harmonics. This work provides a theoretical basis for steady streaming at GHz
with high-order harmonics.

Figure 1:  (Color online) (a) The wavelength, acoustic attenuation length, and
shock distance versus the incident frequency from 0.01 MHz to 1000 GHz (both
axes are in logarithmic coordinates). $B/A=$ 5.1 is taken from Ref. [16] for
water. The vibration velocity to estimate the shock distance is $U_{ac}$ = 1
m/s based on previous experimental results from 0.1 to 10 GHz. The attenuation
length $L_{a}$ is equal to the incident wavelength $\lambda=$ 0.055 $\mu$m at
$f\approx 27.5$ GHz (see the grey dashed line). For a special case at
frequency $f=1.5$ GHz, the wavelength is $\lambda=$ 1 $\mu$m, the attenuation
length $L_{a}=18.4$ $\mu$m, and the shock distance is $L_{s}=$ 67.2 $\mu$m.
The pressure field in the propagation plane with a circular transducer at $f=$
1.5 GHz is shown in (b).

##  II  Characteristic lengths in the frequency regime near 1 GHz

Table 1:  Physical parameters. The fluid medium in this work is water. $f$ is the incident frequency. $A=\rho_{0}c_{0}^{2}$ and B are two classical acoustic coefficients in nonlinear acoustics. Another nonlinear parameter is defined as $\beta=1+B/2A$. $U_{ac}$ is the magnitude of the acoustic velocity perturbation at the transducer surface. Note that a factor of 2 is missing in the definition of $L_{a}$ in Refs. [14] and [15]. Symbol | Physical parameter | Value
---|---|---
$\rho_{0}$ | static density |  1000 kg/m3
$c_{0}$ | acoustic velocity | 1500 m/s
$\omega$ | angular frequency | $2\pi f$
$\mu_{s}$ | dynamic viscosity |  $1.002\times 10^{-3}$ Pa s
$\nu$ | kinematic viscosity | $\mu_{s}/\rho_{0}$
$\mu_{b}$ | bulk viscosity |  $2.8\times 10^{-3}$ Pa s
$b$ | defined coefficients | $4/3+\mu_{b}/\mu_{s}$
$\Gamma$ | Goldberg number | $L_{a}/L_{s}$
$L_{a}$ | acoustic attenuation length | $2\rho_{0}c_{0}^{3}/[\omega^{2}\mu_{s}b]$
$L_{s}$ | acoustic shock distance | $c_{0}^{2}/[\omega\beta U_{ac}]$
Reac | acoustic Reynolds number | $L_{a}/\lambda$
Rehd | hydrodynamic Reynolds number | $\rho_{0}|\mathbf{v_{2}}|L_{c}/\mu_{s}$
$\delta$ | boundary layer thickness | $\sqrt{2\nu/\omega}$
$L_{c}$ | microchannel height |  $60\,\mu$m

Compared with the general frequencies in the medical ultrasonic regime (e.g.,
1-10 MHz), there are two different characteristics in the regime near 1 GHz as
shown in Fig. 1(a): (i) the attenuation length
$L_{a}=2\rho_{0}c_{0}^{3}/[\omega^{2}\mu_{s}b]$ is just 1 order of magnitude
as the wavelength at typical MHz regime, and the acoustic Reynolds number is
$Re_{ac}=L_{a}/\lambda\gg 1$ [15]. $\rho_{0}$ is the static mass density of
the propagation medium, $c_{0}$ is the sound speed, $\omega$ is the angular
frequency, and $b=4/3+\mu_{b}/\mu_{s}$ is the defined coefficient related to
dynamic ($\mu_{s}$) and bulk viscosities ($\mu_{b}$) for convenience. The
related physical parameters and values for the medium of water are listed in
Table 1. At $f=1.5$ GHz in water, the wavelength is $\lambda=1$ $\mu$m and the
attenuation length is $L_{a}=18.4$ $\mu$m [see Fig. 1(b) for the pressure
distribution]. While at $f=1.5$ MHz in water, it has $\lambda=1$ mm and
$L_{a}=1.84\times 10^{4}$ mm. (ii) the attenuation length is comparable with
the shock distance $L_{s}=c_{0}^{2}/[\omega\beta U_{ac}]$ at the typical
vibration velocity of the transducer working at GHz (e.g., $U_{ac}=1$ m/s).
$\beta=1+B/[2A]$ is the nonlinear parameter with two nonlinear acoustic
coefficients $A$ and $B$. At $f=27.5$ GHz as marked with the grey dashed line
in Fig. 1(a), the attenuation length equals the wavelength, i.e.,
$L_{a}=\lambda=0.055$ $\mu$m. That is to say, most of the energy from the
transducer will dissipate in the propagation distance of one wavelength. In
this work, only low GHz transducers for acoustic steaming will be studied and
they have been used most in recent microfluidics experiments in the regime of
GHz. In fact, the shock distance describes the nonlinearity of the acoustic
propagation and the attenuation length is an indicator of the wave dissipation
in the medium. Here, the Goldberg number $\Gamma$ is introduced to provide a
dimensionless measure of the importance of nonlinearity relative to
dissipation with the definition

$\displaystyle\Gamma=\frac{L_{a}}{L_{s}}=\frac{2\rho_{0}c_{0}\beta
U_{ac}}{\omega\mu_{s}b}$ (1)

Note that the attenuation length $L_{a}$ is independent of the activation
velocity $U_{ac}$, while the shock distance $L_{s}$ has a linear relation with
it. In general, there are no high-order harmonics during the acoustic
propagation with $\Gamma\ll 1$ which is the case in Ref. [15]. When $\Gamma$
is around the unit, it is possible to induce the high-order harmonics which
will be discussed in detail in Sec. IV. Note that the shock distance depends
on the $U_{ac}$ and for the case of $U_{ac}=1$ m/s in Fig. 1(a), $\Gamma$
equals 1 (i.e., $L_{s}=L_{a}$) at $f=0.407$ GHz. However, the source term of
bulk streaming with the consideration of high-order acoustic harmonics is not
available and this will be solved in the present work.

##  III  Revisit of source term for bulk streaming

###  III.1 Source term of bulk streaming

Before studying the streaming effect with high-order harmonics, we need to
revisit the source term of bulk streaming. By following the work in Ref. [15],
we briefly recall the source term for the Eckart-type acoustic streaming with
bulk waves. To derive the final formula of the streaming source term, we start
with the constitutive equations including the conservation of mass and
momentum in the fluid medium as

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0$
(2) $\displaystyle\frac{\partial\rho\mathbf{v}}{\partial
t}+\nabla\cdot(\rho\mathbf{v}\otimes\mathbf{v})=-\nabla
p+\mu_{s}\Delta\mathbf{v}+\left(\frac{\mu_{s}}{3}+\mu_{b}\right)\nabla\nabla\cdot\mathbf{v}$
(3)

with $\rho$, $\mathbf{v}$, and $p$ are the density, velocity vector, and
pressure, respectively. $t$ designates the time. The entropy ($s$) balance is
ensured for the system with $\mathrm{d}s=0$ and the state equation is

$\displaystyle p=p(\rho),\text{ with }\left.\frac{\partial
p}{\partial\rho}\right|_{s}=c_{0}^{2}$ (4)

The first-order equation of state is easily obtained by using the Taylor
expansion with $p_{1}=c_{0}^{2}\rho_{1}$, where $p_{1}$, $\rho_{1}$ are the
first-order acoustic pressure and density. In general, compared with the
viscous effect, the thermal effect is negligible because it is proportional to
$\gamma-1$ with $\gamma$ the adiabatic index which is weak in liquids [15].
Since only small perturbation occurs with respect to the hydrostatic
parameters, we can apply the perturbation method and hence decompose the
physical field $\mathcal{X}$ into three parts: hydrostatic $\mathcal{X}_{0}$,
acoustic $\mathcal{X}_{1}$, and hydrodynamic $\mathcal{X}_{2}$. Under the
assumption with the perturbation method, we assume that
$\mathcal{X}_{0}\gg{\mathcal{X}_{1}}\gg{\mathcal{X}_{2}}.$ $\mathcal{X}$ can
be either a scalar or a vector. The time average of acoustic component
$\mathcal{X}_{1}$ is equal to zero, written as
$\left<\mathcal{X}_{1}\right>=0$, while the hydrodynamic part
$\mathcal{X}_{2}$ shows the nonlinear feature with its time average not equal
to zero (i.e., $\left<\mathcal{X}_{2}\right>\neq 0$). Based on the above
assumptions, we can expand the density $\rho$, pressure $p$, and velocity
$\mathbf{v}$ as follows

$\displaystyle\rho=\rho_{0}+{\rho}_{1}+{\rho}_{2},$ (5) $\displaystyle
p=p_{0}+{p}_{1}+{p}_{2}$
$\displaystyle\mathbf{v}=\mathbf{0}+{\mathbf{v}}_{1}+{\mathbf{v}}_{2}.$

Here the particle velocity $\mathbf{v_{0}}$ is zero when not disturbed by
acoustic waves. By expanding the above mass and momentum conservation
equations up to the first- and second-order with some algebraic operations
[see details in Appendix A], the source term of acoustic streaming can be
derived as

$\mathbf{F_{s}}=-\left(\frac{4\mu_{s}}{3}+\mu_{b}\right)\left\langle\frac{\rho_{1}}{\rho_{0}}\Delta\mathbf{v}_{1}\right\rangle,$
(6)

where $\Delta$ is the Laplace operator with
$\Delta\mathbf{v}_{1}=\nabla^{2}\mathbf{v}_{1}=\nabla\cdot\nabla\mathbf{v}_{1}$.
It is noteworthy that this formula excludes the contribution of the gradient
of acoustic Lagrangian, which in fact will not induce the bulk acoustic
streaming [13, 15]. In addition, the gradient of acoustic Lagrangian is very
large compared to the sole source term which is possible to result in large
computational errors. The effect of the acoustic Lagrangian gradient has not
been studied before and will be demonstrated numerically in this work [see
Sec. III.3].

As shown in Appendix B.1 by using the linear wave equation for mono-chromatic
waves, the source term for acoustic streaming in Eq. (6) can be simplified as

$\mathbf{F_{s}}=\left(\frac{4\mu_{s}}{3}+\mu_{b}\right)\frac{\omega^{2}}{\rho_{0}c_{0}^{4}}\langle{p_{1}\mathbf{v}_{1}}\rangle$
(7)

which depends on the time-averaged acoustic intensity
$\langle{\mathbf{I}}\rangle=\langle{p_{1}\mathbf{v}_{1}}\rangle$. It indicates
that once the acoustic field is calculated, the source term to compute the
steady fluid streaming ($\mathbf{v_{2}}$) is available to solve the Stokes
equation derived in Eq. (31)

$-\nabla p^{*}_{2}+\mu_{s}\Delta\mathbf{v}_{2}+\mathbf{F_{s}}=0$ (8)

where $p^{*}_{2}$ is the modified hydrodynamic pressure to remove the
contribution of the average acoustic Lagrangian.

It is noteworthy that the classical source term for acoustic streaming
simulation proposed by Nyborg [11, 12] is widely used although it should be
avoided:

$\mathbf{F_{s}}^{Nb}=-\rho_{0}\nabla\cdot\langle{\mathbf{v_{1}}}\otimes{\mathbf{v_{1}}}\rangle$
(9)

where $\langle\cdot\rangle$ is the time average operator. Recent studies have
shown that this expression is not completely related to acoustic streaming,
which may cause large numerical errors, and the source term can be divided
into two parts through mathematical derivation [15, 14]. In the following, we
propose an analytical example of a one-dimensional (1D) bulk standing wave and
a numerical simulation of a two-dimensional (2D) traveling wave to show the
difference between the streaming field induced by these two source terms,
i.e., based on Eq. (7) and Eq. (9) by Nyborg.

Figure 2:  (color online) The schematic of a 2D acoustic streaming by a planar
transducer. The resonator (black area) excites GHz plane waves in the fluid
(blue) and forms vortices (closed dashed lines) under the activation of body
force (red arrow). In the acoustic simulation, the wall is set as the
impedance boundary condition, and in the flow field simulation, the wall is
set as the no-slip condition. Based on the conservation of fluid inside the
cavity, the Eckart streaming always pushes the fluid up in the center and
rolls back again from the two sides.

###  III.2  One-dimensional bulk standing wave example

In this section, we take the ideal 1D bulk standing waves as an example to
show the difference of the source terms from Nyborg and in Eq. (7). This
simple case will easily show the contradiction of the two source terms for
bulk streaming and can be well understood from the point of view of streaming
physics. The acoustic velocity field and pressure field of the 1D plane
standing waves can be expressed as the addition of two ideal counter-
propagating plane waves:

$\mathbf{v_{1}}=2v_{am}\sin(kz)\cos(\omega t)\mathbf{e_{z}}$ (10a)
$p_{1}=2p_{am}\cos(kz)\sin(\omega t)$ (10b)

where $v_{am}$ and $p_{am}$ are the amplitudes of the acoustic velocity and
pressure, $k$ is the wavenumber, $z$ is the space coordinate with the unit
vector in the propagation direction $\mathbf{e_{z}}$. By substituting them
into the Nyborg source term in Eq. (9), we can get the body force in the
propagation direction as:

$\mathbf{F_{s}}^{Nb}=-2\rho_{0}kv_{am}^{2}\sin(2kz)\left\langle\cos(2\omega
t)+1\right\rangle$ (11)

Obviously, the magnitude of the force $\mathbf{F_{s}}^{Nb}$ does not vanish
since the time dependence $\left\langle\cos(2\omega t)+1\right\rangle$ is not
equal to 0. On the other hand, we can get the following different result if we
use the source term in Eq. (7):

$\mathbf{F_{s}}=\left(\frac{4\mu_{s}}{3}+\mu_{b}\right)\frac{\omega^{2}}{c_{0}^{4}\rho_{0}}\left\langle
p_{1}\mathbf{v}_{1}\right\rangle\propto\left\langle\sin(\omega t)\cos(\omega
t)\right\rangle$ (12)

which vanishes as 0 after the time average procedure in one period. It is
clear that these two results are contradictory.

Figure 3:  (Color online) The flowcharts of the simulation for bulk streaming
(a) when the hydrodynamic Reynolds number is much smaller than 1; (b) when
nonlinear propagation is considered.

As observed from Eq. (7), the body force depends on the average acoustic
intensity which is uniform in space for progressive ideal plane waves. Since
the standing waves can be regarded as the addition of two counter-propagating
plane waves, the body forces of each plane wave should have the same magnitude
while in reversed directions, leading to the null of the total body force for
the streaming effect. This physical explanation agrees with the results of Eq.
(12). Indeed, the source term $\mathbf{F_{s}}^{Nb}$ by Nyborg is not entirely
related to the streaming, and it contains the gradient of acoustic Lagrangian
which contributes to acoustic radiation pressure instead of streaming [13,
15].

###  III.3  Numerical examples of 2D bulk streaming

Since the wavelength at GHz is in the order of nanometers or a few microns, it
often needs small size steps in the computational domain based on the rules of
thumb. Hence, the three-dimensional simulation of bulk streaming problems will
suffer the challenge of a massive amount of numerical computation cost. To
reduce the numerical burden, we build a 2D model as illustrated in Fig. 2, and
this will not hinder our understanding of the acoustic streaming mechanism at
this stage.

Figure 4: (Color online) The streaming simulation of two source terms at 1.5
GHz when the vibration velocity of the transducer surface is $U_{ac}=0.1$m/s.
The radius of the transducer is $a=50$ $\mu$m. The channel height is
$L_{c}=60$ $\mu$m and the width is 200 $\mu$m. The left half is simulated by
using the recast source term [Eq.(7)], and the right half is simulated with
Nyborg’s source term [Eq.(9)]. The uniform colorbar is used. The color
background represents the velocity distribution of the hydrodynamic flow and
the arrows represent the flow directions. Significant differences between the
two source terms are observed to calculate the streaming-induced velocities.
The maximum flow velocity by Nyborg’s expression is 0.02 mm/s, which is much
smaller than that with the present source term in Eq. (7).

The contradiction between the two streaming source terms has been revealed by
the ideal 1D standing wave model in Sec.III.2. This section will further show
the difference of hydrodynamic flow velocities through numerical simulations.
In the following simulations, only the source terms are different [i.e., Eqs.
(7) and (9)]. We use water as the propagation medium with the parameters
listed in Table 1. The excitation frequency is $f=1.5$ GHz and the vibration
velocity of the transducer with a radius 50 $\mu$m is $U_{ac}=$ 0.1 m/s,
leading to the Goldberg number $\Gamma=0.027\ll 1$. Under this condition, the
source term of streaming in Eq. (7) is suitable without the consideration of
nonlinear propagation. The boundary layer thickness (or viscous penetration
depth) is $\delta=\sqrt{2\nu/\omega}=$ 22 nm $\ll L_{c}=60$ $\mu$m with
$L_{c}$ the height of the microchannel. Meanwhile, the microchannel size is
smaller than the shock distance $L_{s}=67.2$ $\mu$m so that the shock waves
will not be accumulated and formed. Hence, the bulk streaming is dominant in
the fluid domain and the Rayleigh boundary streaming can be negligible [17].

The numerical simulations of acoustic streaming in this paper are carried out
by COMSOL Multiphysics 6.0 [18] with the flowchart shown in Fig. 3(a). For the
calculation of the first-order acoustic field, we use the “Thermoacoustic”
Interface. The oscillation velocity boundary is used to replace the effect of
the resonator, and the impedance boundary condition is adopted for the walls.
Then, the steady flow simulations are conducted based on the “creep flow”
interface with the input of the source terms under the circumstance that the
hydrodynamic Reynolds number is $\text{Re}_{hd}=0.072\ll 1$. All the walls
were set with no-slip boundary conditions. Indeed, the simulation with the
“laminar” interface obtains the same streaming results as the “creep flow”
interface case (not shown in the following for brevity).

Figure 5:  (color online) Acoustic pressure based on the theoretical method in
Ref. [19] and frequency spectrum with two different Goldberg numbers. (a) For
$\Gamma=0.1$, the pressure amplitude is $p_{am}=0.6$ MPa with the attenuation
distance $L_{a}=18.4$ um (indicated by the black dotted line) and the shock
distance $L_{s}=168.1$ um. (b) Only the fundamental wave appears along with a
few sidelobes. There are no obvious high-order harmonics in the frequency
spectrum. (c) For $\Gamma=1.8$, the pressure amplitude $p_{am}=10$ MPa with
$L_{a}=18.4$ um and $L_{s}=10.1$ um. (d) Different from the results in (b),
the second-order harmonic wave appears in this case.

In the 2D numerical simulations of bulk streaming, we set the volume forces in
the propagation direction $F_{z}$ and the lateral direction $F_{r}$ with the
flow velocities indicated by the colormap of the background and the directions
indicated by the black arrows in Fig. 4. Because of the symmetry of the
hydrodynamic flows in the microchannel, we show half of the simulation results
and put them together for ease of comparison: the left half is for the source
term in Eq. (7) while the right half is for the source term by Nyborg in Eq.
(9). It is obvious that the streaming velocities of these two terms have a
significant difference at this time. Note that the main differences lie in the
magnitudes of the maximum flow velocities, while the hydrodynamic flow
patterns are similar because of the conservation of fluid flow in the
microchannels. Through this example, it is obvious that there are problems in
using the classical source term by Nyborg for acoustic streaming simulation.
In addition, it will be noteworthy that an empirical expression of the source
term [1, 3] is widely used for the present simulations of GHz streaming which
is limited to the plane wave case and could not predict the exact values of
the hydrodynamic flow velocities.

##  IV  Weakly nonlinear propagation with high-order acoustic harmonics

The source term of the bulk streaming in Eq. (7) is proper under the
assumption of linear acoustic propagation. However, it does not apply if
nonlinear acoustics are considered (e.g., high-order harmonics) for the
streaming phenomenon since the pressure field may include the contribution of
multiple frequencies. To solve the issue, we re-derive a general expression of
the sole source term for acoustic streaming with the multi-frequency acoustic
field in terms of only acoustic pressure $p_{1}$ as follows [see Eq. (38)]

$\mathbf{F_{s}}=\left(\frac{4}{3}\mu_{s}+\mu_{b}\right)\left\langle\frac{{p}_{1}}{c_{0}^{4}\rho_{0}^{2}}\nabla\frac{\partial{p}_{1}}{\partial
t}\right\rangle$ (13)

With the detailed derivation of the source expression given in Appendix B.2.
This formula has the advantage to compute the incompressible hydrodynamic
fluid motion induced by acoustic streaming once the acoustic field can be
calculated.

Recall the Goldberg number as first introduced in Sec. II, which plays an
important role in the viscous process considering nonlinear effects [20]. It
measures the relative importance of nonlinear effects and dissipation effects.
Since only weakly nonlinear propagation is considered in this work and there
are no shock waves as assumed, the following will derive the source terms with
Goldberg numbers within $\Gamma\ll 1$ and $\Gamma\sim 1$ based on the
analytical expressions of harmonic waves from finite-amplitude Gaussian beam
in a fluid by Du and Breazeale [19].

###  IV.1 Acoustic pressure and frequency spectrum at different $\Gamma$

To calculate the acoustic streaming with nonlinear propagation, we need to
solve the pressure field according to Eq.(13). The nonlinear propagation of
acoustics could be calculated with either analytical methods or numerical
simulations. This part will review the analytical expressions of the acoustic
pressure fields of the fundamental and second-order harmonics by Du &
Breazeale for a Gaussian beam from a piston transducer [19]. As Goldberg
points out, the shock wave will not be formed when $\Gamma<1$. This is because
the dissipation effect is significant under the circumstance, and the
harmonics cannot be effectively accumulated. For the piston transducer, we
assume that the vibration velocity of the sound source satisfies the Gaussian
function and take the Gaussian coefficient as a unit for simplicity:

$U(R)=U_{ac}e^{-R^{2}/a^{2}}$ (14)

where $R$ is the radial coordinate, and $a$ is the radius of the piston
transducer. By taking the method from Du and Breazeale, we can expand its
acoustic pressure waveform into Fourier series (see Eqs. (A1) and (A4) in Ref.
[19]):

$p_{1}=\sum_{n=1}^{2}P_{n}\sin\left[n\omega\left(t-\frac{z}{c_{0}}\right)\right]$
(15)

where $P_{n}$ represents the spatial components of the $n_{th}$ harmonic
acoustic pressure, $t$ designates time, and $z$ is the distance from the
transducer surface in the propagation direction. This expression illustrates
the physical mechanism of the nonlinear pressure field in terms of the
addition of different orders of harmonics. For the fundamental component (the
first order)

$P_{1}(f\mid z)=p_{am}e^{-\alpha z}\exp\left(-R^{2}/a^{2}\right)$ (16)

and for the second-order harmonics

$\begin{split}P_{2}(f\mid z)=\frac{p_{am}^{2}k\beta e^{-4\alpha
z}}{4\alpha\rho_{0}c_{0}^{2}}e^{\left(-2R^{2}/a^{2}\right)}\left(e^{2\alpha
z}-1\right)\end{split}$ (17)

where the pressure amplitude at the center of the transducer surface is
$p_{am}=\rho_{0}c_{0}U_{ac}$, and $\alpha=1/L_{a}$ is the attenuation
coefficient. The detailed derivation is briefly organized in Appendix B.3.

An example shows that the nonlinear effect starts around $\Gamma=3$ for a
plane wave by Hamilton and Blackstock, while other scholars select 4.5 as the
threshold for shock formation [21]. To compare the pressure fields with or
without high-order harmonics, two pressure amplitudes are selected to make the
Gol’dber numbers $\Gamma=0.1$ and $\Gamma=1.8$ for the Gaussian beam,
respectively. The radius of the piston transducer is $a=50$ $\mu$m with the
excitation frequency $f=1.5$ GHz. The pressure fields along the propagation
direction $z$ are computed with the analytical method and numerical
simulations based on COMSOL [22]. These two results agree well with each other
and are shown (one blue curve for brevity) in Fig. 5(a) at $p_{am}=0.6$ MPa
and (c) at $p_{am}=10$ MPa. The attenuation distance $L_{a}$ is indicated by
the vertical dashed lines. To illustrate the nonlinear high-order harmonics
when $\Gamma>1$, the frequency spectrum based on fast Fourier transform is
calculated as shown in Fig. 5(b) and (d). It could be observed that there is a
sharp peak at $f=3$ GHz for $\Gamma=1.8$ in the frequency spectrum of Fig.
5(d) which comes from the contribution of the second-order harmonics.

###  IV.2 Source terms of streaming with and without high-order harmonics

Figure 6:  (Color online) The ratios of pressure magnitude of the fundamental
wave $p_{10}/p_{am}$ (red solid line), second harmonic wave $p_{20}/p_{am}$
(blue dotted line), and the body force $F_{z}^{h}/F_{z}$ (black dashed line)
change along the axial direction when the pressure magnitude is $p_{am}=$ 10
MPa and the frequency is 1.5 GHz. $F_{z}^{h}$ is the axial body force
considering the contribution of both the fundamental and the second-order
harmonics, which is defined by Eq. (19). $F_{z}$ is the volume force
considering only the fundamental frequency acoustic pressure, which is
determined by Eq. (18). The fundamental sound pressure decreases with the
increase of axial distance, and the second harmonic sound pressure and body
force increase first and then reduce with the increase of axial distance. The
attenuated energy of the fundamental wave turns into high-order harmonics.

According to the general source term of acoustic streaming in Eq. (13),
explicit expressions with and without high-order harmonics could be obtained
once the pressure fields are known. For the situation with only fundamental
component, by the insertion of Eqs. (15) with (16) into Eq. (13), we can
rewrite the axial body force $F_{z}^{s}$ as (note $n$ is truncated up to 1 for
the fundamental wave case):

$F_{z}^{s}=\frac{\alpha}{\rho_{0}c_{0}^{2}}\left|P_{1}\right|^{2}$ (18)

where $\left|P_{1}\right|^{2}=p_{am}^{2}e^{-2\alpha z}e^{-2R^{2}/a^{2}}$.
Similarly, by substituting Eq. (15) into Eq. (13), the axial body force
$F_{z}^{h}$ considering the second-order harmonics can be expressed as:

$F_{z}^{h}=\frac{\alpha}{\rho_{0}c_{0}^{2}}\sum_{n=1}^{2}\left(n^{2}\left|P_{n}\right|^{2}\right)$
(19)

where the upper right “$h$” represents the body force with the contribution of
the second-order harmonics. Indeed, this source term works for high-orde
harmonics when $n>2$. When only the fundamental wave is considered ($n=1$),
Eq. (19) degenerates into (18). Note that the lateral component of the source
term vanishes after the time average procedures, i.e., $F_{R}=0$ for the cases
with and without high-order harmonics.

To study the contribution of the second-order harmonics on the total pressure
field and body force, the same piston transducer is used with the initial
pressure amplitude $p_{am}=10$ MPa at the excitation frequency $f=1.5$ GHz.
The sizes of the microchannel are the same as those in Fig. (4). The
simulation flowchart is given in Fig. 3(b) with the pressure field computed
based on the theoretical method in Matlab and the induced streaming obtained
by using COMSOL. Note that the hydrodynamic flow speed can reach up to 1 m/s,
leading to $Re_{hd}\gg 1$, the ‘laminar’ interface is applied in this section.
The normalized pressure fields versus the propagation distance $z$ of the
fundamental (red solid line) and second-order (blue dotted line) harmonic
waves are shown in Fig. 6. The surface of the transducer is defined as $z=0$.
The amplitude of the fundamental component decrease versus $z$ because of the
wave absorption. While the pressure amplitude of the second-order harmonics
increases from the transducer surface to a maximum value at $z\approx 6$
$\mu$m, which agrees with the fact that the nonlinearity effect is dominant
over the dissipation close to the source. In addition, the relative axial body
force with the consideration of the second-order harmonics with respect to
only the fundamental wave contribution is plotted versus $z$ with the black
dashed line in Fig. 6. A horizontal grey dashed line is also provided as a
reference. It is shown that the excess of the body force comes from the
contribution of the second-order harmonics.

To further explore the streaming-induced hydrodynamic motion in the
microchannel, we take the same configuration as in Fig. 6 but double the
pressure amplitude $p_{ac}=20$ MPa to enhance the contribution of the second-
order harmonics. The left half part of Fig. 7 shows the streaming patterns
with only the fundamental component, while the right half is for the situation
with both the fundamental and second-order harmonics. The colormap in the
background indicates the velocity amplitude of the hydrodynamic flows with the
arrows giving the flow directions. The maximum hydrodynamic flow velocity
increases more than 20$\%$ if the second harmonic component is taken into
consideration, as shown in Fig. 7.

Figure 7:  (color online) Streaming results without and with the second-order
harmonics at 1.5 GHz. The acoustic pressure is 20 MPa. The left half shows the
acoustic streaming field with only the contribution of the fundamental wave
with the maximum flow speed of 13.47 m/s, it is calculated by Eq.(18). The
right half considers the contribution of the second-order harmonics with the
maximum flow speed of 17.11 m/s, which is based on Eq.(19). The color depth
represents the streaming velocity.

##  V Conclusion and discussion

The theory of the source term for bulk acoustic streaming has been developed
with an emphasis on the nonlinear propagation in the frequency regime of GHz.
A dimensionless number called Goldberg number ($\Gamma$) is introduced to
measure the importance of acoustic nonlinearity relative to dissipation at
GHz. This work makes it possible to compute the streaming-induced microfluidic
motion for the recent GHz devices in a better manner instead of using the
empirical formulas [1, 3]. As shown in the numerical simulation, the
contribution of high-order harmonics could increase the maximum velocity of
hydrodynamic flows up to 20% if only the fundamental component is considered
in linear propagation. In addition, the source term of bulk streaming in
linear acoustic propagation is revisited in the form of the wave damping and
acoustic intensity [15]. More importantly, we propose analytical and numerical
examples to show the contradictory results by the source term of linear
propagation and the classic expression by Nyborg [11, 12]. The term by Nyborg
should be avoided by others (e.g., see Ref. [23]) since it contains the
contribution of acoustic radiation pressure [13, 15].

The alternative method to compute the bulk streaming is the full-model direct
numerical simulation (DNS) [7]. However, it will take more computational cost
and will be challenging to handle the full simulation in the 3D domain. The
computational burden could be relieved with the effective boundary conditions
developed by Bach and Bruus [24] and using the axial symmetry of the
geometrical model [25]. However, attention should be still paid to the
simulation of 3D bulk streaming if there is no axial symmetry with either the
present theory or DNS, for instance, the drop-shaped transducer working at GHz
in a microchannel [3]. It should be noted that this work could be helpful to
design acoustical tweezers for particle trapping with the consideration of
both the acoustic radiation pressure[26, 27] and the streaming-induced drag
force.

###### Acknowledgements.

Z.Gong Thanks for the support from the National Natural Science Foundation of
China (24Z990200542) and Shanghai Jiao Tong University for the startup
funding.

##  Appendix A  Derivation of the sole source of acoustic streaming

###  A.1  First-order perturbation expansion

To derive the source term of acoustic streaming (i.e., sound-induced fluid
motion), we need to bridge the acoustic field (acoustic pressure $p_{1}$ and
hydrodynamic velocity $\mathbf{v}_{1}$) with the fluid velocity
$\mathbf{v}_{2}$. Using the perturbation method, we need to expand the fields
up to the first-order and second-order terms. Firstly, we expand the fields up
to the linear limit with the following perturbation

$\displaystyle\rho=\rho_{0}+\rho_{1}$ (20) $\displaystyle p=p_{0}+p_{1}$
$\displaystyle\mathbf{v}=\mathbf{0}+\mathbf{v}_{1}$

Substituting Eq. (20) into (2) and Eq. (3), we can obtain the first-order
formulas of mass and momentum conservation respectively

$\frac{\partial\rho_{1}}{\partial
t}+\rho_{0}\nabla\cdot\mathbf{v}_{\mathbf{1}}=0$ (21)
$\rho_{0}\frac{\partial\mathbf{v}_{1}}{\partial t}=-\nabla
p_{1}+\mu_{s}\Delta\mathbf{v}_{1}+\left(\frac{\mu_{s}}{3}+\mu_{b}\right)\nabla\nabla\cdot\mathbf{v}_{1}$
(22)

where $\Delta\mathbf{v}_{1}=\nabla\cdot(\nabla\mathbf{v}_{1})$. Since the
first-order acoustic field is defined as irrotational (i.e.,
$\nabla\times\mathbf{v}_{1}=0$), Eq. (22) can be further changed into the
following form:

$\rho_{0}\frac{\partial\mathbf{v}_{1}}{\partial t}=-\nabla
p_{1}+\left(\frac{4}{3}\mu_{s}+\mu_{b}\right)\Delta\mathbf{v_{1}}$ (23)

It should be noted that the derivation of the above equation uses the
mathematical identity of the vector Laplacian:
$\nabla\nabla\cdot\mathbf{v_{1}}=\nabla^{2}\mathbf{v_{1}}+\nabla\times\nabla\times\mathbf{v_{1}}$
[see Eq.(3.70) in Ref. 28].

###  A.2  Second-order perturbation expansion

Since acoustic streaming is a nonlinear phenomenon, we need to expand the
fields up to the second-order as given in Eq. (5). After substitution into the
mass conservation equation, one obtains

$\frac{\partial\rho_{2}}{\partial
t}+\nabla\cdot(\rho_{1}\mathbf{v_{1}})+\rho_{0}\nabla\cdot\mathbf{v_{2}}=0$
(24)

Before expand the momentum conservation equation [Eq. (3)] up to the second-
order, we first use the mass conservation to cancel some items in the momentum
conservation [see Eq. (3)]

$\begin{split}&\left[\frac{\partial\rho}{\partial
t}+\nabla\cdot\left(\rho\mathbf{v}\right)\right]\mathbf{v}+\rho\frac{\partial\mathbf{v}}{\partial
t}+\rho\mathbf{v}\cdot\nabla{\mathbf{v}}\\\ &=-\nabla
p+\mu_{s}\Delta\mathbf{v}+\left(\frac{\mu_{s}}{3}+\mu_{b}\right)\nabla\nabla\cdot\mathbf{v}\end{split}$
(25)

with the terms in the first square brackets vanishing based on the global mass
conservation as given in Eq. (2). Note that the identity
$\nabla\cdot(\rho\mathbf{v}\otimes\mathbf{v})=\rho\mathbf{v}\cdot\nabla(\mathbf{v})+\mathbf{v}\nabla\cdot(\rho\mathbf{v})$
is used here. By insertion of Eq. (5) into (25) and taking the equation up to
the second order:

$\begin{split}&\rho_{0}\frac{\partial\mathbf{v_{2}}}{\partial
t}+\rho_{1}\frac{\partial\mathbf{v_{1}}}{\partial
t}+\rho_{0}\mathbf{v_{1}}\cdot\nabla\mathbf{v_{1}}\\\ &=-\nabla
p_{2}+\mu_{s}\Delta\mathbf{v_{2}}+\left(\frac{\mu_{s}}{3}+\mu_{b}\right)\nabla\nabla\cdot\mathbf{v_{2}}\end{split}$
(26)

Since only the steady acoustic streaming is of interest, the first term of Eq.
(26) vanishes with $\partial\mathbf{v_{2}}/\partial t=0$. Throughout the whole
work, the incompressible flow assumption is hold, so that

$\nabla\cdot\mathbf{v}_{2}=0$ (27)

The last term on the right-hand side in the second-order momentum conservation
equation [Eq.(26)] vanishes. To isolate the hydrodynamic fluid motion from the
acoustic perturbation over time (i.e., $\partial\mathbf{v}_{1}/\partial t$),
we need to recall the first-order momentum conservation in Eq. (23) which is
rewritten as the following for convenience

$\rho_{1}\frac{\partial\mathbf{v_{1}}}{\partial
t}=-\frac{\rho_{1}}{\rho_{0}}\nabla
p_{1}+\frac{\rho_{1}}{\rho_{0}}\left(\frac{4}{3}\mu_{s}+\mu_{b}\right)\Delta\mathbf{v_{1}}$
(28)

Substitute Eq. (28) into Eq. (26), we get

$\begin{split}-\frac{\rho_{1}}{\rho_{0}}\nabla
p_{1}+\frac{\rho_{1}}{\rho_{0}}\left(\frac{4}{3}\mu_{s}+\mu_{b}\right)\Delta\mathbf{v_{1}}+\rho_{0}\mathbf{v_{1}}\cdot\nabla\mathbf{v_{1}}\\\
=-\nabla p_{2}+\mu_{s}\Delta\mathbf{v_{2}}\end{split}$ (29)

If we take the time average of the above Eq. (29) and use the definitions of
kinetic energy density
$\left\langle\mathcal{K}\right\rangle=\left\langle\rho_{0}\mathbf{v_{1}}\cdot\nabla\mathbf{v_{1}}\right\rangle=\left\langle\nabla\left(1/2\rho_{0}\mathbf{v_{1}^{2}}\right)\right\rangle$
and potential energy density
$\left\langle\mathcal{U}\right\rangle=\left\langle{\rho_{1}}/{\rho_{0}}\nabla
p_{1}\right\rangle=\left\langle\nabla\left({p_{1}^{2}}/({2\rho_{0}c_{0}^{2}})\right)\right\rangle$
(note that the first-order equation of state is applied here), the second-
order momentum conservation equation can be expressed as

$\nabla\langle\mathcal{L}\rangle+\left(\frac{4\mu_{s}}{3}+\mu_{b}\right)\left\langle\frac{\rho_{1}}{\rho_{0}}\Delta\mathbf{v}_{1}\right\rangle=\left\langle-\nabla
p_{2}+\mu_{s}\Delta\mathbf{v}_{2}\right\rangle$ (30)

where the time-averaged acoustic Lagrangian is defined as
$\langle\mathcal{L}\rangle=\langle\mathcal{K}-\mathcal{U}\rangle$. As noted
that the acoustic Lagrangian is independent of the viscous effect and hence
not affected by the wave attenuation, while the acoustic streaming is induced
by the transfer of acoustic momentum to the viscous mode through wave
attenuation. That is to say, the averaged acoustic Lagrangian
$\langle\mathcal{L}\rangle$ is not related to the steady acoustic streaming
and can be balanced with a hydrostatic pressure gradient, having $\nabla
p^{*}_{2}=\nabla p_{2}+\langle\mathcal{L}\rangle$. [13, 15] Hence, the second-
order momentum conservation equation can be finally written as

$-\nabla p^{*}_{2}+\mu_{s}\Delta\mathbf{v}_{2}+\mathbf{F_{s}}=0$ (31)

with the sole source for steady acoustic streaming is

$\mathbf{F_{s}}=-\left(\frac{4\mu_{s}}{3}+\mu_{b}\right)\left\langle\frac{\rho_{1}}{\rho_{0}}\Delta\mathbf{v}_{1}\right\rangle$
(32)

We use the hydrodynamic Reynolds number $\mathrm{Re_{hd}}$ to characterize the
relative contribution of the inertia term
$\rho_{0}\left(\mathbf{v}_{2}\cdot\nabla\right)\mathbf{v}_{2}$ and the
viscosity term $\mu_{s}\nabla^{2}\mathbf{v}_{2}$[13]:

$\mathrm{Re_{hd}}=\frac{\rho_{0}{v}_{2}d}{\mu_{s}}$ (33)

where $v_{2}$ is the characteristic streaming velocity of the system, $d$ is
the characteristic length of the system and $d\sim 60\mu m$. The hydrodynamic
Reynolds number is equal to 1 when the streaming velocity is:

${v}_{cr}=\frac{\mu_{s}}{\rho_{0}d}=16.7\mathrm{mm/s}$ (34)

When $\mathrm{Re_{hd}}\gg 1$ $\left(v_{2}\gg v_{cr}\right)$, the inertia term
is dominant, and we use Lighthill theory to simulate the acoustic streaming,
at this point, the governing equation becomes:

$\rho_{0}\left({\mathbf{v_{2}}}\cdot\nabla\right)\mathbf{v_{2}}+\mu_{s}\Delta\mathbf{v_{2}}-\nabla
p_{2}=\mathbf{F}$ (35)

When $\mathrm{Re_{hd}}\ll 1$ $\left({v}_{2}\ll{v}_{cr}\right)$, the viscosity
term is dominant, and the
$\rho_{0}\left(\mathbf{v_{2}}\cdot\nabla\right)\mathbf{v_{2}}$ will disappear,
the governing equation will be restored to (31).

This makes it possible to use commercial software like COMSOL to solve the
acoustic streaming problems with the source term induced by acoustic field.

##  Appendix B  Explicit expressions of source terms

###  B.1  Monochromatic fundamental wave

For most applications using the acoustic streaming effect, the acoustic wave
is considered as a monochromatic wave with the linear wave equation for the
velocity as

$\frac{\partial^{2}\mathbf{v}_{1}}{\partial
t^{2}}-c_{0}^{2}\Delta\mathbf{v}_{1}=0$ (36)

where $\partial^{2}/\partial t^{2}=-\omega^{2}$ for a single-frequency
fundamental harmonic wave with $\omega=2\pi f$ the angular frequency. Hence,
the wave equation can be simplified as
$\Delta\mathbf{v}_{1}=-\omega^{2}\mathbf{v}_{1}/c_{0}^{2}$. Insertion of the
steady wave equation into the source term of Eq. (32) with the combination of
the first-order equation of state, the sole source of acoustic streaming is

$\begin{split}\mathbf{F_{s}}=\left(\frac{4\mu_{s}}{3}+\mu_{b}\right)\frac{\omega^{2}}{c_{0}^{4}\rho_{0}}\left\langle
p_{1}\mathbf{v}_{1}\right\rangle\\\ \end{split}$ (37)

with the average acoustic density
$\left\langle\mathbf{I}\right\rangle=\left\langle
p_{1}\mathbf{v}_{1}\right\rangle$. This formula is also given in Eq. (19) of
Ref. [14].

###  B.2  High-order harmonic waves

When the weak nonlinear effect induces a few high-order harmonics, e.g., up to
the order of $n=2$, the source term for the acoustic streaming given in Eq.
(LABEL:Eq._B2:_source_for_singe_frequency_wave) will not apply which is
limited to the case of single-frequency fundamental wave. Under this
condition, we have to re-derive the source term based on the general
expression of Eq. (32) with only two assumptions: the irrotational acoustic
field ($\nabla\times\mathbf{v}_{1}=0$) and the incompressible streaming fluid
field ($\nabla\cdot\mathbf{v}_{2}=0$). By combining the first-order mass
conservation equation [see Eq. (21)] and the first-order equation of state,
the relation between the acoustic pressure $p_{1}$ and velocity vector
$\mathbf{v}_{1}$ can be derived as $\nabla
p_{1}/(\rho_{0}c_{0}^{2})=\Delta\mathbf{v}_{1}$ with vector identity
$\nabla\nabla\cdot\mathbf{v_{1}}=\nabla^{2}\mathbf{v_{1}}+\nabla\times\nabla\times\mathbf{v_{1}}$.
By insertion of this relation into the general expression of body force
$\mathbf{F_{s}}$ in Eq. (32), one can derive the following source form
including only the acoustic pressure:

$\mathbf{F_{s}}=\left(\frac{4}{3}\mu_{s}+\mu_{b}\right)\left\langle\frac{{p}_{1}}{c_{0}^{4}\rho_{0}^{2}}\nabla\frac{\partial{p}_{1}}{\partial
t}\right\rangle$ (38)

This form is helpful to deal with acoustic streaming problems if the acoustic
field with high-order harmonics is solved with either analytical or numerical
methods.

Consider the acoustic pressure form including low-order harmonics[19]:

$\begin{split}p_{1}&=P_{1}\sin\left[\omega\left(t-\frac{z}{c_{0}}\right)\right]+P_{2}\sin\left[2\omega\left(t-\frac{z}{c_{0}}\right)\right]\\\
&=\sum_{n=1}^{2}P_{n}\sin\left[n\omega\left(t-\frac{z}{c_{0}}\right)\right]\end{split}$
(39)

Substituting the above equation into Equation (38), the following form of body
force can be obtained:

$F_{z}^{h}=\frac{\alpha}{\rho_{0}c_{0}^{2}}\sum_{n=1}^{2}\left\langle
n^{2}P_{n}^{2}\right\rangle$ (40)

The orthogonality of trigonometric functions is used in the derivation of the
above formula, that is $1/T\int_{0}^{T}[sin\left(n\omega
t\right)sin\left(m\omega t\right)]dt$ = 0 when $n\neq m$ and $T=2\pi/\omega$
is the period for the fundamental harmonics.

###  B.3  Pressure fields of high-order harmonic waves by Du & Breazeale

By using the perturbation method to solve the nonlinear wave equation (first
proposed by Kuznetsov [29]), we can obtain the following fundamental
quasilinear solution (see Eq.(4) in see Ref. [19] with $z/z_{0}$$\ll$1 and
$\gamma$$\approx$$\pi$/2):

$P_{1}(f\mid z)=p_{am}\frac{e^{-\alpha
z}}{\sqrt{1+\left(z/z_{0}\right)^{2}}}\exp\left({-\frac{R^{2}/a^{2}}{1+\left(z/z_{0}\right)^{2}}}\right)$
(41)

where $z_{0}=ka^{2}/2$. It is noted that the near-field condition is
applicable ($z\ll z_{0}$) when the channel height is small. In this case, the
diffraction effect is very weak, and therefore the acoustic pressure in Eq.
(41) can be simplified as

$P_{1}(f\mid z)=p_{am}e^{-\alpha z}e^{\left(-R^{2}/a^{2}\right)}$ (42)

Noted that in Blackstock’s study [20], the extra attenuation (EXDB) is very
small (The EXDB refers to the loss beyond the normal small signal attenuation
$e^{-\alpha z}$) when the Goldberg number $\Gamma$ $\sim$ 1\. Under this
condition, low-order harmonics will be generated, while no shock waves will be
formed. With the help of the Hankel transform, $P_{2}$ can be calculated by
following the work of Du & Breazeale [see Eqs. (11) and (12) in Ref. [19]):

$\begin{split}P_{2}(f\mid z)=\frac{p_{am}D_{0}e^{-4\alpha
z}}{4\sqrt{1+(z/z_{0})^{2}}}\exp{\left(-\frac{2R^{2}/a^{2}}{1+(z/z_{0})^{2}}\right)}\\\
\times\sqrt{H_{2}^{2}+F_{2}^{2}}\end{split}$ (43)

with

$H_{2}=\int_{0}^{\sigma}\frac{e^{2\alpha
z_{0}\sigma^{\prime}}}{\sqrt{1+\left(z/z_{0}\right)^{2}}}\cos\left[\tan^{-1}\left(\sigma^{\prime}\right)\right]d\sigma^{\prime}$
(44)

and

$F_{2}=\int_{0}^{\sigma}\frac{e^{2\alpha
z_{0}\sigma^{\prime}}}{\sqrt{1+\left(z/z_{0}\right)^{2}}}\sin\left[\tan^{-1}\left(\sigma^{\prime}\right)\right]d\sigma^{\prime}$
(45)

where $\sigma=z/z_{0}$, and $D_{0}={2z_{0}}/{L_{s}}$. Considering the near-
field condition $z\ll z_{0}$, we can further simplify the second harmonic
acoustic pressure Eq. (43) as:

$P_{2}(f\mid z)=\frac{p_{am}^{2}k\beta e^{-4\alpha
z}}{4\alpha\rho_{0}c_{0}^{2}}\exp\left(-2R^{2}/a^{2}\right)\left(e^{2\alpha
z}-1\right)$ (46)

Note that the solutions for higher harmonics up to the fourth order are also
listed in the Appendix of Ref. [19].

## Reference

  * Cui  _et al._ [2016] W. Cui, H. Zhang, H. Zhang, Y. Yang, M. He, H. Qu, W. Pang, D. Zhang, and X. Duan, Localized ultrahigh frequency acoustic fields induced micro-vortices for submilliseconds microfluidic mixing, Applied Physics Letters 109, 253503 (2016).
  * Cui  _et al._ [2019a] W. Cui, L. Mu, X. Duan, W. Pang, and M. A. Reed, Trapping of sub-100 nm nanoparticles using gigahertz acoustofluidic tweezers for biosensing applications, Nanoscale 11, 14625 (2019a).
  * Yang  _et al._ [2022] Y. Yang, L. Zhang, K. Jin, M. He, W. Wei, X. Chen, Q. Yang, Y. Wang, W. Pang, X. Ren, _et al._ , Self-adaptive virtual microchannel for continuous enrichment and separation of nanoparticles, Science Advances 8, eabn8440 (2022).
  * Shilton  _et al._ [2014] R. J. Shilton, M. Travagliati, F. Beltram, and M. Cecchini, Nanoliter-droplet acoustic streaming via ultra high frequency surface acoustic waves, Advanced materials 26, 4941 (2014).
  * Dentry  _et al._ [2014] M. B. Dentry, L. Y. Yeo, and J. R. Friend, Frequency effects on the scale and behavior of acoustic streaming, Physical Review E 89, 013203 (2014).
  * Cui  _et al._ [2019b] W. Cui, W. Pang, Y. Yang, T. Li, and X. Duan, Theoretical and experimental characterizations of gigahertz acoustic streaming in microscale fluids, Nanotechnology and Precision Engineering 2, 15 (2019b).
  * Steckel [2021] A. G. Steckel,  _Theory and modeling of thin-film actuation in microscale acoustofluidics_ , Ph.D. dissertation, Technical University of Denmark, Kongens Lyngby, Denmark (2021).
  * White and Voltmer [1965] R. M. White and F. W. Voltmer, Direct piezoelectric coupling to surface elastic waves, Applied physics letters 7, 314 (1965).
  * Zhang  _et al._ [2006] H. Zhang, W. Pang, H. Yu, and E. S. Kim, High-tone bulk acoustic resonators on sapphire, crystal quartz, fused silica, and silicon substrates, Journal of applied physics 99, 124911 (2006).
  * Eckart [1948] C. Eckart, Vortices and streams caused by sound waves, Physical review 73, 68 (1948).
  * Nyborg [1953] W. L. Nyborg, Acoustic streaming due to attenuated plane waves, The Journal of the Acoustical Society of America 25, 68 (1953).
  * Nyborg [1965] W. L. Nyborg,  _Acoustic streaming in Physical Acoustics Vol. 2 Pt. B_ (Academic Press, New York, 1965) pp. 265–331.
  * Lighthill [1978] J. Lighthill, Acoustic streaming, Journal of sound and vibration 61, 391 (1978).
  * Baudoin and Thomas [2020] M. Baudoin and J.-L. Thomas, Acoustic tweezers for particle and fluid micromanipulation, Annual Review of Fluid Mechanics 52, 205 (2020).
  * Riaud  _et al._ [2017] A. Riaud, M. Baudoin, O. B. Matar, J.-L. Thomas, and P. Brunet, On the influence of viscosity and caustics on acoustic streaming in sessile droplets: an experimental and a numerical study with a cost-effective method, Journal of Fluid Mechanics 821, 384 (2017).
  * Panfilova  _et al._ [2021] A. Panfilova, R. J. van Sloun, H. Wijkstra, O. A. Sapozhnikov, and M. Mischi, A review on b/a measurement methods with a clinical perspective, The Journal of the Acoustical Society of America 149, 2200 (2021).
  * Vanneste and Bühler [2011] J. Vanneste and O. Bühler, Streaming by leaky surface acoustic waves, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 1779 (2011).
  * [18] Comsol multiphysics 6.0, http://cn.comsol.com/.
  * Du and Breazeale [1986] G. Du and M. A. Breazeale, Harmonic distortion of a finite amplitude gaussian beam in a fluid, Journal of the Acoustical Society of America 80, 212 (1986).
  * Blackstock [1964] D. T. Blackstock, Thermoviscous attenuation of plane, periodic, finite‐amplitude sound waves, Journal of the Acoustical Society of America 36, 534 (1964).
  * Hamilton [2016] M. F. Hamilton, Effective gol’dberg number for diverging waves, The Journal of the Acoustical Society of America 140, 4419 (2016).
  * [22] COMSOL Multiphysics Model Library, Nonlinear Acoustics- Modeling of the 1D Westervelt Equation, 2013.
  * Eisener  _et al._ [2015] J. Eisener, A. Lippert, T. Nowak, C. Cairós, F. Reuter, and R. Mettin, Characterization of acoustic streaming beyond 100 mhz, Physics Procedia 70, 151 (2015).
  * Bach and Bruus [2018] J. S. Bach and H. Bruus, Theory of pressure acoustics with viscous boundary layers and streaming in curved elastic cavities, The Journal of the Acoustical Society of America 144, 766 (2018).
  * Jorgensen [2022] J. H. Jorgensen,  _Theory and modeling of thermoviscous acoustofluidics_ , Ph.D. dissertation, Technical University of Denmark, Kongens Lyngby, Denmark (2022).
  * Silva [2011] G. T. Silva, An expression for the radiation force exerted by an acoustic beam with arbitrary wavefront (l), The Journal of the Acoustical Society of America 130, 3541 (2011).
  * Gong and Baudoin [2021] Z. Gong and M. Baudoin, Equivalence between angular spectrum-based and multipole expansion-based formulas of the acoustic radiation force and torque, The Journal of the Acoustical Society of America 149, 3469 (2021).
  * Arfken  _et al._ [2013] G. B. Arfken, H. J. Weber, and F. E. Harris,  _Mathematical methods for physicists_ , 7th ed. (Academic Press, New York, 2013) pp. 155,756–765.
  * Kuznetsov [1971] V. P. Kuznetsov, Equation of nonlinear acoustics, Sov.phys.acoust (1971).


uniformly. Indeed, $\|E(0)\|\leq C_{\delta}+\max_{i=1,\cdots,N}|Y_{i}|$, where
$C_{\delta}>0$ is the constant introduced in equation C.3, which in turn
implicitly defines $\mathcal{A}_{\delta}$. Now, it follows from the Berry-
Esseen theorem that

$-\frac{1}{\sqrt{n}}+\mathbb{P}(|G|\leq
C_{\delta})\lesssim\mathbb{P}(|f_{\theta(t)}(x)|\leq
C_{\delta})\lesssim\frac{1}{\sqrt{n}}+\mathbb{P}(|G|\leq C_{\delta}),$

with $G$ a Gaussian. Hence, $\mathbb{P}(|G|\leq C_{\delta})\sim
1-e^{-C_{\delta}^{2}}$. Then, equation C.3 is guaranteed to hold if

$-\frac{1}{\sqrt{n}}+(1-e^{-C_{\delta}^{2}})\gtrsim 1-\delta,$

which we can rewrite as

$\delta\gtrsim e^{-C_{\delta}^{2}}+\frac{1}{\sqrt{n}}.$

We conclude the following.

######  Proposition 1.

There are positive constants $C_{1},C_{2}>0$ with the following significance.
Let $\delta_{n}$ and $C_{\delta_{n}}$ be a pair of sequences satisfying

(C.28) $C_{1}\delta_{n}>e^{-C_{\delta_{n}}^{2}}+\frac{1}{\sqrt{n}}.$

Then, there is a set $\mathcal{A}_{\delta_{m}}\subset\mathbb{R}^{4n}$ of
measure $1-\delta_{n}$; such that for any $x$ with $|x|\leq 1$ and
$\theta(0)\in\mathcal{A}_{\delta_{n}}$

$\displaystyle|f_{\theta(t)}(x)-f_{\infty}(x,t)|$ $\displaystyle\leq$
$\displaystyle\left(\max_{i=1,\cdots,N}|K^{\infty}_{xi}|+1\right)\frac{C_{2}\|E(0)\|}{\min\\{\lambda_{\infty},\lambda\\}}\frac{(\log(n))^{1/2+}}{\sqrt{n}}\;,$

for any $n\in\mathbb{N}$. Moreover,

$\|E(0)\|\leq C_{\delta_{n}}+\max_{i=1,\cdots,N}|Y_{i}|.$

######  Example 3.

Consider the case $\delta_{n}=C\frac{(\log n)^{1/2+}}{\sqrt{n}}$ for a
constant $C>0$ and

$C_{\delta_{n}}^{2}=\frac{1}{2}\log n-\Big{(}\frac{1}{2}\Big{)}^{1+}\log\log
n.$

We can check that these satisfy C.28 and so the conclusion of the previous
proposition is valid with

$\displaystyle|f_{\theta(t)}(x)-f_{\infty}(x,t)|$ $\displaystyle\lesssim$
$\displaystyle\frac{(\log(n))^{1+}}{\sqrt{n}}\;,$

## References

  * [1] S. Arora, S. S. Du, W. Hu, Z. Li, R. Salakhutdinov, and R. Wang. On exact computation with an infinitely wide neural net. In Advances in Neural Information Processing Systems, 32, 2018.
  * [2] S. Arora, S. S. Du, W. Hu, Z. Li, and Ruosong Wang. Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks. In 36th International Conference on Machine Learning, ICML 2019.
  * [3] P. L. Bartlett, P. M. Long, G. Lugosi, and A. Tsigler. Benign overfitting in linear regression. In Proceedings of the National Academy of Sciences, 117.48: 30063-30070, 2020.
  * [4] P. Billingsley. Convergence of probability measures. John Wiley & Sons, 2013.
  * [5] V. Bogachev. Weak convergence of measures. American Mathematical Society, 2018.
  * [6] Y. Cao, and and G. Quanquan. Generalization bounds of stochastic gradient descent for wide and deep neural networks. In Advances in neural information processing systems, 32, 2019.
  * [7] L. Carvalho, J. L. Costa, J. Mourão, and G. Oliveira. On the spectral gap of the Neural Tangent Kernel of wide and deep neural networks, in preparation.
  * [8] L. Chizat, E. Oyallon, and F. Bach. On lazy training in differentiable programming. In Advances in neural information processing systems, 32, 2019.
  * [9] S. Du, J. D. Lee, H. Li, L. Wang, and X. Zhai. Gradient descent finds global minima of deep neural networks. In International conference on machine learning, 2019.
  * [10] S. Du, X. Zhai, B. Poczos, A. Singh. Gradient descent provably optimizes over-parameterized neural networks. In International Conference on Learning Representations, 2018.
  * [11] R. Durrett. Probability: Theory and Examples, Cambridge University Press, 2019.
  * [12] N. Radford. Priors for infinite networks, (tech. rep. no. crg-tr-94-1), University of Toronto, 1994\.
  * [13] D. A. Roberts, and S. Yaida. The Principles of Deep Learning Theory: An Effective Theory Approach to Understanding Neural Networks, Cambridge University Press, 2022.
  * [14] Z. Chen, G. Rotskoff, J. Bruna, and E. Vanden-Eijnden. A Dynamical Central Limit Theorem for Shallow Neural Networks, In NeurIPS Proceedings, 2020.
  * [15] Hall, Peter. The bootstrap and Edgeworth expansion. Springer Science and Business Media, 2013.
  * [16] A. Jacot, G. Franck, and C. Hongler. Neural Tangent Kernel: Convergence and Generalization in Neural Networks. In Advances in Neural Information Processing Systems, 2018.
  * [17] B. Ghorbani, S. Mei, T. Misiakiewicz, and A. Montanari. When do neural networks outperform kernel methods?. In Advances in Neural Information Processing Systems, 33, 2020.
  * [18] J. Huang, H.-T. Yau. Dynamics of Deep Neural Networks and Neural Tangent Hierarchy. In Proceedings of the 37th International Conference on Machine Learning, PMLR 119:4542-4551, 2020.
  * [19] G. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang. Physics-informed machine learning. In Nature Reviews Physics, 3.6: 422-440, 2021.
  * [20] J. Lee, Y. Bahri, R. Novak, S. Schoenholz, J. Pennington, J. Sohl-Dickstein. Deep Neural Networks as Gaussian Processes. In International Conference on Learning Representations, 2018.
  * [21] J. Lee, S. Schoenholz, J. Pennington, B. Adlam, L. Xiao, R. Novak, and J. Sohl-Dickstein. Finite versus infinite neural networks: an empirical study. In Advances in Neural Information Processing Systems, 33: 15156-15172, 2020.
  * [22] J. Lee, L. Xiao, S. Schoenholz, Y. Bahri, R. Novak, J. Sohl-Dickstein, and J. Pennington. Wide neural networks of any depth evolve as linear models under gradient descent. In Advances in neural information processing systems, 32 2019.
  * [23] S. Mei, and A. Montanari. The generalization error of random features regression: Precise asymptotics and the double descent curve. Communications on Pure and Applied Mathematics, 75.4: 667-766, 2022.
  * [24] B. Neyshabur, Z. Li, S. Bhojanapalli, Y. LeCun, and N. Srebro. Towards understanding the role of over-parametrization in generalization of neural networks. In International Conference on Learning Representations, 2019.
  * [25] S. Oymak, and M. Soltanolkotabi. Toward moderate overparameterization: Global convergence guarantees for training shallow neural networks. In IEEE Journal on Selected Areas in Information Theory, 1.1: 84-105, 2020.
  * [26] Rivera, Alejandro. Statistical mechanics of Gaussian fields. Diss. Université Grenoble Alpes, 2018.
  * [27] S. Wang, X. Yu, and P. Perdikaris. When and why PINNs fail to train: A neural tangent kernel perspective. In Journal of Computational Physics, 449: 110768, 2022.
  * [28] K. Xu, M. Zhang, J. Li, S. Du, K. Kawarabayashi, and S. Jegelka. How neural networks extrapolate: From feedforward to graph neural networks. In International Conference on Learning Representations, 2021.
  * [29] G. Yang. Scaling limits of wide neural networks with weight sharing: Gaussian process behavior, gradient independence, and neural tangent kernel derivation. Preprint arXiv:1902.04760, 2019.
  * [30] G. Yang, and E. J. Hu. Tensor programs iv: Feature learning in infinite-width neural networks. In International Conference on Machine Learning, 2021.
  * [31] Y. Yamasaki. Measures on infinite dimensional spaces, World Scientific, 1985.