Kinetic and potential energy densities of fields predicted by linear acoustics are first derived from exact nonlinear expressions of the energies. [The acoustic Poynting theorem (a theorem that relates power and intensity) is then derived, which is written in terms of the energy densities of the linear acoustic fields. This derivation is not actually needed in the calculation of radiation force.] The acoustic radiation stress tensor is derived from the conservation of momentum and mass. Diagonal elements of the matrix representation of this tensor are related to the mean excess pressure, which can be written in terms of energy densities. The acoustic radiation force is calculated by appealing to the conservation of momentum and mass, which is written as a single differential equation in terms of the acoustic radiation stress tensor.
The derivation of the potential and kinetic energy densities below are based on my class notes from Prof. Hamilton's Acoustics I course. These energy densities are quadratic quantities calculated by taking the product of two field variables that solve the linear acoustic wave equation.
The exact potential and kinetic energy densities of sound are \begin{align} E_T &= \tfrac{1}{2} \rho v^2\,, \label{eq:acoust:ET} \\ E_V &= -\frac{1}{V_0} \int_{V_0}^{V} p\,dV\,, \label{eq:acoust:EV} \end{align} where \(p = P-P_0\) is the linear acoustic field variable, and where \(v^2 = \vec{v} \cdot \vec{v} = |\vec{v}|^2\). Equation \eqref{eq:acoust:ET} is linearized by writing \(\rho = \rho_0 + \rho'\), where \(\rho_0\) is a constant and \(\rho'\) is the acoustic (perturbation) quantity that is a function of space and time: \[E_T = \tfrac{1}{2} (\rho_0 + \rho')v^2\,.\] Neglecting the cubic term \(\rho'u^2\) yields \begin{align}\label{eq:kinetic} \boxed{E_T = \tfrac{1}{2}\rho_0v^2\,.} \end{align}
Linearizing Eq. \eqref{eq:acoust:EV} is more challenging, because \(dV\) is not a convenient differential. It is desired to instead integrate over the acoustic pressure \(p\). First, \(dV\) is written in terms of \(d\rho'\), and then \(\rho'\) is related to the pressure through the linearized equation of state, as is now shown.
Volume is related to density through the relation \(M = \rho V\), where \(M\) is mass. Writing \(V = M/\rho\) allows for the calculation of the derivative \[\frac{dV}{d\rho} = - \frac{M}{\rho^2}\,.\] Solving for \(dV\) yields \begin{align} dV &= - \frac{M}{\rho^2}\,d\rho \notag\\ &= -\frac{V}{\rho}d\rho\,,\label{eq:acoust:dV} \end{align} where in Eq. \eqref{eq:acoust:dV} the fact that \(M = \rho V\) has again been invoked. Noting that \(d\rho = d(\rho_0 + \rho') = d\rho'\) and \(V = V_0 + V'\), where \(V'\) is the acoustic volume, yields \begin{align} dV &= -(V_0 + V')(\rho_0 + \rho')^{-1} d\rho'\,,\label{eq:acoust:dV:1} \end{align} where \(1/\rho\) has been written as \((\rho+\rho')^{-1}\). Noting that \begin{align*} (\rho_0+\rho')^{-1} &= \frac{(1 + \rho'/\rho_0)^{-1}}{\rho_0} \simeq \frac{1 - \rho'/\rho_0}{\rho_0}\,, \end{align*} Eq. \eqref{eq:acoust:dV:1} becomes \begin{align} dV &= -(V_0 + V')(1 - \rho'/\rho_0)\frac{d\rho'}{\rho_0}\notag\\ &= -\frac{V_0}{\rho_0} d\rho' + \text{higher-order terms}\,,\label{eq:acoust:dV:2} \end{align} where in Eq. \eqref{eq:acoust:dV:2} only the linear term has been retained. Thus Eq. \eqref{eq:acoust:EV} becomes \begin{align} E_V = \frac{1}{\rho_0} \int_{0}^{\rho'} p\, d\tilde{\rho}'\,, \label{eq:acoust:EV:1} \end{align} where the integration is over \(\tilde{\rho}'\) so that the upper limit of integration is distinct from the integration variable. Equation \eqref{eq:acoust:EV:1} can be related to the acoustic pressure by appealing to the linearized equation of state \[\rho' = \frac{1}{c_0^2}p\,,\quad d\rho' = \frac{1}{c_0^2}dp,\] which is derived by Taylor expanding \(P(\rho)\) as \(P_0 + (\partial P/\partial \rho)_0 (\rho-\rho_0)\) and identifying \((\partial P/\partial \rho)_0\) as \(c_0^2\). Equation \eqref{eq:acoust:EV:1} becomes \begin{align}\label{eq:potential} \boxed{E_V = \frac{1}{\rho_0c_0^2} \int_0^p \tilde{p}\, d\tilde{p} = \frac{p^2}{2\rho_0c_0^2 }} \end{align}
The following derivation is based on Sec. 1-11 of Ref. [5], in which it is called the "acoustic energy corollary." The theorem is not needed in deriving the radiation force, which follows from the conservation of momentum.
To derive the acoustic Poynting theorem, start with the linearized conservation of momentum, \begin{align}\label{eq:acoust:poynt:mom} \gradient p + \rho_0 \frac{\partial \vec{v}}{\partial t} = 0\,. \end{align} Dotting Eq. \eqref{eq:acoust:poynt:mom} into \(\vec{v}\) yields \begin{align}\label{eq:acoust:poynt:mom:2} \vec{v} \cdot \gradient p + \rho_0 \vec{v} \cdot \frac{\partial \vec{v}}{\partial t} = 0\,. \end{align} Noting that \(\vec{v} \cdot \partial \vec{v}/\partial t = \tfrac{1}{2} \partial v^2/\partial t\), Eq. \eqref{eq:acoust:poynt:mom:2} becomes \begin{align}\label{eq:acoust:poynt:mom:3} \vec{v} \cdot \gradient p + \tfrac{1}{2} \rho_0\frac{\partial v^2}{\partial t} = 0\,. \end{align} Noting that \(\divergence(p\vec{v}) = \vec{v} \cdot \gradient p + p \divergence \vec{v}\), and identifying the instantaneous Poynting vector \(\vec{S} = p\vec{v}\), allows for Eq. \eqref{eq:acoust:poynt:mom:3} to be written as \begin{align}\label{eq:acoust:poynt:mom:4} \divergence \vec{S} - p\divergence\vec{v} + \tfrac{1}{2}\rho_0 \frac{\partial v^2}{\partial t} = 0\,. \end{align} The third term in Eq. \eqref{eq:acoust:poynt:mom:4} is simply the time derivative of Eq. \eqref{eq:kinetic}, i.e., \(\partial E_T/\partial t\). Meanwhile, the second term of Eq. \eqref{eq:acoust:poynt:mom:4} resembles the quantity that appears in the linearized conservation of mass equation, \({\partial \rho'}/{\partial t} + \rho_0 \divergence \vec{v} = 0\), which when combined with the linearized state equation \(p/c_0^2 = \rho'\), yields \begin{align}\label{eq:acoust:poynt:mass} \frac{1}{c_0^2}\frac{\partial p}{\partial t} + \rho_0 \divergence\vec{v} = 0\,. \end{align} Multiplying Eq. \eqref{eq:acoust:poynt:mass} by \(p/\rho_0\) and rearranging yields \begin{align}\label{eq:acoust:poynt:mass:2} -p\divergence\vec{v} = \frac{p}{\rho_0c_0^2}\frac{\partial p}{\partial t} = \frac{1}{2\rho_0c_0^2}\frac{\partial p^2}{\partial t} = \frac{\partial E_V}{\partial t}\,, \end{align} where the last equality holds by Eq. \eqref{eq:potential}. Inserting Eq. \eqref{eq:acoust:poynt:mass:2} into Eq. \eqref{eq:acoust:poynt:mom:4} yields \begin{align}\label{eq:acoust:poynt:thm} \boxed{\divergence \vec{S} + \frac{\partial}{\partial t}(E_V + E_T) = 0\,.} \end{align}
This section is based on Sec. 2.1 of Ref. [6]. There, the derivation of the acoustic radiation stress tensor is performed in index notation, while it is done symbolically here.
Neglecting terms related to viscocity, the exact conservation of momentum requires that \begin{align}\label{eq:mom:exact} \rho\frac{\partial \vec{v}}{\partial t} + \rho\vec{v} \cdot \gradient \vec{v} = -\gradient P\,. \end{align} Meanwhile, the exact continuity equation states that \begin{align}\label{eq:cont:exact} \frac{\partial \rho}{\partial t} + \divergence (\rho \vec{v}) = 0\,. \end{align} Multiplying Eq. \eqref{eq:cont:exact} by \(\vec{v}\) yields \begin{align}\label{eq:cont:exact:1} \vec{v}\frac{\partial \rho}{\partial t} + \vec{v} \divergence (\rho \vec{v}) = \vec{0}\,. \end{align} Adding Eqs. \eqref{eq:mom:exact} and \eqref{eq:cont:exact:1} yields \[\rho\frac{\partial \vec{v}}{\partial t} + \rho\vec{v} \cdot \gradient \vec{v} + \vec{v}\frac{\partial \rho}{\partial t} + \vec{v} \divergence (\rho \vec{v}) = -\gradient P\,.\] Rearranging terms yields \[\rho\frac{\partial \vec{v}}{\partial t} + \vec{v}\frac{\partial \rho}{\partial t} + \rho\vec{v} \cdot \gradient \vec{v} + \vec{v} \divergence (\rho \vec{v}) = -\gradient P\,.\] Using the product rules \(\rho\partial \vec{v}/\partial t + \vec{v}\partial \rho/\partial t = \partial(\rho\vec{v})/\partial t\) and \(\rho\vec{v} \cdot \gradient \vec{v} + \vec{v} \divergence (\rho \vec{v}) = \divergence (\rho \vec{v}\otimes \vec{v})\) (see Item 1 of "Useful identities") yields \begin{align}\label{eq:take-avg} \frac{\partial (\rho\vec{v})}{\partial t} + \divergence (\rho \vec{v}\otimes \vec{v}) = -\gradient P\,. \end{align} It is now assumed that \(\rho\), \(\vec{v}\), and \(P\) in Eq. \eqref{eq:take-avg} are time-harmonic. The time average of Eq. \eqref{eq:take-avg} is taken. By Item 2 of the useful identities, \(\langle {\partial (\rho\vec{v})}/{\partial t}\rangle = 0\), leaving \begin{align}\label{eq:mom:2} \langle\divergence (\rho \vec{v}\otimes \vec{v})\rangle = -\langle\gradient P\rangle\,. \end{align} Since the time averages amount to integrals over time, the spatial differential operators can be moved outside the time averages: \begin{align}\label{eq:mom:3} -\gradient \langle P\rangle -\divergence \left\langle\rho \vec{v}\otimes \vec{v}\right\rangle = \vec{0}\,. \end{align} Note by Item 3 of the useful identities that \(\dyad{I} \cdot \gradient \langle P\rangle = \divergence(\dyad{I} \langle P\rangle ) - \langle P\rangle \divergence \dyad{I}\), where \(\dyad{I}\) is the identity tensor. Since \(\divergence \dyad{I} = \vec{0}\), one can equivalently write \( \dyad{I} \cdot \gradient \langle P\rangle = \divergence(\dyad{I} \langle P\rangle)\). Equation \eqref{eq:mom:3} then becomes \begin{align}\label{eq:mom:4} \divergence (-\dyad{I} \langle P\rangle -\left\langle\rho \vec{v}\otimes \vec{v}\right\rangle) = \vec{0}\,. \end{align} Noting that \(\langle P\rangle = \langle P-P_0\rangle\) and retaining terms up to quadratic order in Eq. \eqref{eq:mom:4} yields \begin{align}\label{eq:mom} \divergence (-\dyad{I} \langle P-P_0\rangle - \rho_0\left\langle \vec{v}\otimes \vec{v}\right\rangle) = \vec{0}\,. \end{align} The quantity in parentheses is identified as the acoustic radiation stress tensor: \begin{align}\label{eq:stress} \boxed{\dyad{S} = - \dyad{I} \langle P-P_0 \rangle - \rho_0 \left\langle \vec{v}\otimes \vec{v}\right\rangle\,.} \end{align} In orthonormal index notation, the \(ij^\mathrm{th}\) component of Eq. \eqref{eq:stress} is \begin{align}\label{eq:stress:ind} {S_{ij} = - \delta_{ij} \langle P-P_0 \rangle - \rho_0 \left\langle v_i v_j\right\rangle\,,} \end{align} recovering Eq. (7) of Ch. 6, Sec. 2.1.2 of Ref. [6]. Note the distinction between the stress tensor, "\(\dyad{S}\)" (sans-serif) and the Poynting vector "\(\vec{S}\)" (serif). In terms of the acoustic radiation stress tensor, Eq. \eqref{eq:mom} becomes \begin{align}\label{eq:mom:stress} \boxed{\divergence \dyad{S} = \vec{0}\,.} \end{align}
Attention is now turned to the quantity \(\langle P - P_0 \rangle\), referred to as the "mean excess pressure," that appears in Eq. \eqref{eq:stress}. The derivation of this quantity is based on Sec. 2.1.3 of Ref. [6] (though in a slightly different order, and with the details worked out).
It is desired to calculate the time average of the quantity \(P-P_0\), which is reminiscent of the Taylor expansion of the pressure in the enthalpy per unit mass \(w\), where \(dw = T ds + dP/\rho\), which reduces to \begin{align}\label{eq:enthalpy} dw = \frac{dP}{\rho} \end{align} for an isentropic process (\(ds = 0\)): \[P = P_0 + \left(\frac{\partial P}{\partial w}\right)_{s,0}w + \frac{1}{2}\left(\frac{\partial^2 P}{\partial w^2 }\right)_{s,0}w^2 + \dots\,.\] The motivation for expanding in the enthalpy (rather another thermodynamic variable, like density) will be seen presently. Since by Eq. \eqref{eq:enthalpy} \((\partial P/\partial w)_{s,0} = dP/dw = \rho\), the Taylor expansion becomes \begin{align} P - P_0 &= \rho w + \frac{1}{2}\left(\frac{\partial\rho}{\partial w} \right)_{s,0}w^2 + \dots\,.\label{eq:taylor:simp} \end{align} Since \[\frac{\partial \rho}{\partial w} = \frac{\partial \rho}{\partial P}\frac{\partial P}{\partial w} = \frac{1}{c^2} \frac{dP}{dw} = \frac{\rho}{c^2}\,,\] Eq. \eqref{eq:taylor:simp} becomes \begin{align} P - P_0 &= \rho_0 w + \frac{\rho_0}{2c_0^2} w^2 + \dots\,.\label{eq:taylor:simp:2} \end{align} Next, to evaluate \(w\) and \(w^2\) in Eq. \eqref{eq:taylor:simp:2}, Eq. \eqref{eq:enthalpy} is divided by \(d\vec{r}\), yielding \begin{align}\label{eq:enthalpy:2} \gradient w = \frac{\gradient P}{\rho}\,. \end{align} Equation \eqref{eq:enthalpy:2} is equated to the left-hand side of the momentum equation given by Eq. \eqref{eq:mom:exact} divided by \(\rho\): \begin{align}\label{eq:ex:mom} -\gradient w = -\frac{\gradient P}{\rho} = \frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \gradient \vec{v}\,. \end{align} Noting that \(w\) appears in Eq. \eqref{eq:taylor:simp}, it is desired to integrate Eq. \eqref{eq:ex:mom} over volume so as to remove the gradient operator. To facilitate this manipulation, the velocity potential \(\phi\) is introduced, where \(\vec{v} = \gradient \phi\). No generality is lost because in the absence of viscocity terms in Eq. \eqref{eq:mom:exact}, the velocity field is irrotational (\(\curl \vec{v}=\vec{0}\)), and thus \(\vec{v}\) is equal to the gradient of a scalar potential \(\phi\). In terms of \(\phi\), Eq. \eqref{eq:ex:mom} becomes \begin{align}\label{eq:ex:mom:1} -\gradient w =\gradient \frac{\partial \phi}{\partial t} + \gradient \phi \cdot \gradient \gradient \phi\,. \end{align} Noting by the product rule that \begin{align*} \gradient \phi \cdot \gradient\gradient \phi = \tfrac{1}{2}\gradient (\gradient \phi \cdot \gradient \phi) = \tfrac{1}{2}\gradient(|\gradient \phi|^2) \,, \end{align*} Eq. \eqref{eq:ex:mom:1} becomes \begin{align}\label{eq:ex:mom:2} \gradient w = -\gradient \left[\frac{\partial \phi}{\partial t} - \tfrac{1}{2} |\gradient \phi|^2 \right]\,. \end{align} Equation \eqref{eq:ex:mom:2} can be integrated over volume, introducing an arbitrary constant \(C'\): \begin{align}\label{eq:ex:mom:3} w = -\frac{\partial \phi}{\partial t} - \tfrac{1}{2} |\gradient \phi|^2 + C'\,. \end{align} The motivation for expanding \(P\) in \(w\) can now be appreciated: it was possible to use the momentum equation to express \(w\) in terms of \(\phi\), its gradient, and its time derivative. Combining Eq. \eqref{eq:ex:mom:3} with Eq. \eqref{eq:taylor:simp:2} yields \begin{align} P - P_0 &= \rho_0\left(-\frac{\partial \phi}{\partial t} - \tfrac{1}{2} |\gradient \phi|^2 + C' \right) + \frac{\rho_0}{2c_0^2} \left( -\frac{\partial \phi}{\partial t} - \tfrac{1}{2} |\gradient \phi|^2 + C'\right)^2 + \dots\,.\label{eq:taylor:simp:3} \end{align} Equation \eqref{eq:taylor:simp:3} matches Eq. (11) of Ref. [6]. Taking the time average of Eq. \eqref{eq:taylor:simp:3} eliminates terms that are of odd powers, namely \(\partial \phi/\partial t\). According to Ref. [6], \(C'\) is a quadratic quantity, and thus it has a nonzero time average \(\langle C'\rangle \equiv \rho_0 C\). Eliminating quantities of quartic order and higher yields, for the time average of Eq. \eqref{eq:taylor:simp:3}, \begin{align} \langle P - P_0\rangle &=\frac{\rho_0}{2c_0^2} \left\langle\left(\frac{\partial \phi}{\partial t} \right)^2\right\rangle - \tfrac{1}{2} \rho_0 \langle |\gradient \phi|^2\rangle + \rho_0C \,.\label{eq:taylor:simp:4} \end{align} \(\vec{v}\) is now reinstated for \(\gradient \phi\). To relate \(\partial \phi/\partial t\) to the traditional acoustic field variables \(p\) and \(v\), it suffices to consider the linearized version of Eq. \eqref{eq:mom:exact}, since this quantity is squared in Eq. \eqref{eq:taylor:simp:4}: \begin{align*} \gradient p = -\rho_0 \frac{\partial \vec{v}}{\partial t} = -\rho_0 \gradient \frac{\partial \phi}{\partial t}\,. \end{align*} Integrating over volume, noting that the integration constant is arbitrary (since the value of the potential \(\phi\) can be set to \(0\) anywhere without changing the values of \(p\) and \(v\)), and solving for \(\partial \phi/\partial t\) yields the relation \({\partial \phi}/{\partial t} = -p/ \rho_0\). Thus Eq. \eqref{eq:taylor:simp:4} becomes \begin{align} \boxed{\langle P - P_0\rangle =\frac{1}{2\rho_0 c_0^2} \langle p^2\rangle - \tfrac{1}{2} \rho_0 \langle v^2\rangle + \rho_0C \,,} \label{eq:excess} \end{align} matching Eq. (13) of Ref. [6]. Combining Eq. \eqref{eq:excess} with Eqs. \eqref{eq:kinetic} and \eqref{eq:potential} allows for the mean excess pressure to be defined in terms of the Lagrangian density \(L = E_T - E_V\): \begin{align} \langle P - P_0\rangle = - \langle L \rangle + \rho_0C \,. \label{eq:excess:L} \end{align} Equation Eq. \eqref{eq:excess:L} recovers Eq. (4) of Ref. [9].
The term "Lagrangian" above should not be confused with the Lagrangian coordinates used in continuum mechanics. Confusingly, there is a way of writing Eq. \eqref{eq:excess} in said Lagrangian coordinates—see Sec. 2.1.4 of Ref. [6]. Meanwhile, Eq. \eqref{eq:excess} as written "is Eulerian because it is evaluated at a fixed point in space" [6]. As such, \(P\) in Eq. \eqref{eq:excess} is sometimes superscripted with an "\(\mathrm{E}\)" to emphasize that the quantity is Eulerian. In these notes, only the Eulerian quantity is used, so the superscript is suppressed.
Equation \eqref{eq:mom:stress} has units of stress per unit length, which equals force per unit volume. Integrating over volume therefore results in quantity with dimensions of force. To calculate the radiation force, first consider the spherical surface \(A\) immediately surrounding an object (enclosing a corresponding volume \(V\)), and a much larger spherical surface \(A_0\) (enclosing a corresponding volume \(V_0\) and concentric with \(A\)), as shown below:
The approach to calculate the radiation force on the object enclosed by \(A\) is attributed to Ref. [10] (which is in turn based on Ref. [9]). Reference [6] (Sec. 3.1.2) summarizes the procedure in those works as follows, where the notation has been modified to fit the present conventions:
By integrating \(\divergence \dyad{S} = 0\) over the space between [the surface] \(A\) and a much larger spherical surface \(A_0\) concentric with the sphere, and using Gauss's theorem, we obtain a surface integral of \(\dyad{S}\cdot \vec{e}_n\) over \(A' = A + A_0\), equated to zero, where \(\vec{e}_n\) is the normal unit vector on the surface \(A'\) pointing away from the enclosed space... [T]he integral over \(-\dyad{S} \cdot \vec{e}_n\) over \(A\) is the force acting on the sphere...
Equation \eqref{eq:mom:stress} is integrated over the volume within surface \(A_0\) but outside surface \(A\): \begin{equation}\label{eq:mom:stress:int:1} \int_{V_0} \divergence \dyad{S}\,dV_0 - \int_V \divergence \dyad{S}\,dV = \vec{0}\,. \end{equation} Upon invoking the divergence theorem, the first term in Eq. \eqref{eq:mom:stress:int:1} becomes \(\oint_{A_0} \dyad{S}\cdot \vec{e}_r\, dA_0\), while the second term becomes \(-\oint_A \dyad{S}\cdot \vec{e}_r\, dA\), where \(\vec{e}_r\) is the radial unit vector in spherical coordinates: \begin{equation}\label{eq:mom:stress:int:2} \oint_{A_0} \dyad{S}\cdot \vec{e}_r\,dA_0 -\oint_A \dyad{S}\cdot \vec{e}_r\, dA =\vec{0}\,. \end{equation} Since it is desired for the radiation force on the object enclosed by \(A\) to be defined with respect to the inward normal, \(\vec{e}_n\) is identified as \(-\vec{e}_r\) on \(A\), while it \(\vec{e}_n\) is identified as \(\vec{e}_r\) on \(A_0\), as shown in the figure above: \begin{equation}\label{eq:mom:stress:int:3} \oint_{A_0} \dyad{S}\cdot \vec{e}_n\,dA_0 + \oint_A \dyad{S}\cdot \vec{e}_n\, dA =\vec{0}\,. \end{equation} Solving Eq. \eqref{eq:mom:stress:int:3} for \(-\oint_A \dyad{S}\cdot \vec{e}_n\, dA\) results in the radiation force on the object enclosed by \(A\): \begin{equation}\label{eq:mom:stress:int} {\vec{F} = -\oint_A \dyad{S}\cdot \vec{e}_n\, dA = \oint_{A_0} \dyad{S}\cdot \vec{e}_n\,dA_0 \,.} \end{equation} The first equality of Eq. \eqref{eq:mom:stress:int}, when expressed in index notation, recovers the first equation (not numbered) in Ref. [7]. When expressed in terms of Eq. \eqref{eq:stress:ind} (the indicial form of the acoustic radiation stress tensor), Eq. \eqref{eq:mom:stress:int} recovers Eq. (82) of Ref. [6]: \begin{equation}\label{eq:mom:stress:int:ind} F_i = \oint_A [\delta_{ij} \langle P-P_0 \rangle + \rho_0 \left\langle v_i v_j\right\rangle] \, {n}_j\, dA = -\oint_{A_0} [\delta_{ij} \langle P-P_0 \rangle + \rho_0 \left\langle v_i v_j\right\rangle]\, {n}_j\,dA_0 \,. \end{equation} [Note that typo in Eq. (82) of Ref. [6]: "\(F_z\)" should read "\(F_i\)." The force becomes \(F_z\) in Eq. (84), in which azimuthal symmetry has been assumed in Eq. (83).] Symbolically, Eq. \eqref{eq:mom:stress:int:ind} reads \begin{align} \vec{F} = -\oint_{A_0} \left(\dyad{I} \langle P-P_0 \rangle + \rho_0 \left\langle \vec{v}\otimes \vec{v}\right\rangle \right)\cdot \vec{e}_n\,dA_0 \,, \label{eq:mom:stress:int:symb} \end{align} where Eq. \eqref{eq:stress} has been used. Equation \eqref{eq:mom:stress:int:symb} can be written in terms of the acoustic fields \(p\) and \(\vec{v}\) by using Eq. \eqref{eq:excess}, resulting in \begin{align} \boxed{\vec{F} = -\oint_{A_0} \left[\dyad{I} \left(\frac{1}{2\rho_0 c_0^2} \langle p^2\rangle - \tfrac{1}{2} \rho_0 \langle v^2\rangle \right) + \rho_0 \left\langle \vec{v}\otimes \vec{v}\right\rangle \right]\cdot \vec{e}_n\,dA_0\,,}\label{eq:int:F:acous} \end{align} where \(C\) has been set to zero to ensure "zero perturbation at infinity" {See Ref. [6], p. 182}. When reduced to index notation and expressed in terms of velocity potential, Eq. \eqref{eq:int:F:acous} recovers the equation at the top of the second column on the first page of Ref. [7]. Equation \eqref{eq:int:F:acous} reveals the utility of the surface integral over \(A_0\) in the second equality of Eq. \eqref{eq:mom:stress:int}: since the surface \(A_0\) can be arbitrarily large, the calculation of the radiation force \(\vec{F}\) amounts to solving the far-field linear scattering problem (which provides \(p\) and \(\vec{v}\)), taking time averages of \(p^2\), \(v^2\), and \(\vec{v} \otimes \vec{v}\), and integrating the results over \(A_0\). Indeed, Gor'kov writes [7],
[F]or the calculation of the average force correct up to terms of the second order in the velocity, it is sufficient to find the solution of the linear scattering problem.Equation \eqref{eq:int:F:acous} will be evaluated for (1) traveling wave fields and (2) standing wave fields.
The following approach to calculating radiation force on objects subjected to planar traveling waves is credited to Westervelt [9, 10].
It is convenient to write the mean excess pressure in Eq. \eqref{eq:int:F:acous} (which is the quantity in parentheses) in terms of the Lagrangian density \(L = E_T - E_V\), as in Eq. \eqref{eq:excess:L}: \begin{align} \vec{F} = \oint_{A_0} \left[\dyad{I} \langle L\rangle - \rho_0 \left\langle \vec{v}\otimes \vec{v}\right\rangle\right]\cdot \vec{e}_n\,dA_0\,.\label{eq:int:F:L} \end{align} Equation \eqref{eq:int:F:L} is equivalent to Eq. (4) of Ref. [10]. For both traveling and standing wave configurations, the pressure and particle velocity fields are written in terms of the incident fields (subscripted "\(i\)") and scattered fields (subscripted "\(s\)"), \begin{align} p &= p_i + p_s \label{eq:scat:total:p}\\ \vec{v} &= \vec{v}_i + \vec{v}_s\,, \label{eq:scat:total:v} \end{align} and thus the Lagrangian density can be expanded as \(L = L_i + L_s + L_{is}\), where \begin{align} L_i &= \tfrac{1}{2}\rho_0 v_i^2 - p_i^2/2\rho_0c_0^2 \,, \label{eq:L:i}\\ L_s &= \tfrac{1}{2}\rho_0 v_s^2 - p_s^2/2\rho_0c_0^2 \,, \label{eq:L:s}\\ L_{is} &= \rho_0 \vec{v}_i \cdot \vec{v}_s - p_i p_s/\rho_0c_0^2 = \rho_0 v_i v_s \cos\psi - p_i p_s/\rho_0c_0^2 \,, \label{eq:L:is} \end{align} where \(\psi\) is the angle between \(\vec{v}_i\) and \(\vec{v}_s\). (Westervelt calls this angle \(\theta\), but on this page, \(\theta\) is reserved for the spherical polar angle). Equations \eqref{eq:L:i} and \eqref{eq:L:s} are the "contributions to the flux of momentum from the incident and scattered waves respectively" [10]. Westervelt calls Eq. \eqref{eq:L:is} the "interaction Lagrangian" since it contains the cross terms of the squares of Eq. \eqref{eq:scat:total:p} and \eqref{eq:scat:total:v}. Similarly, the quantity \(\vec{v} \otimes \vec{v}\) that appears in the second term in the integrand of Eq. \eqref{eq:int:F:L} can be expanded as \begin{align} \vec{v}\otimes \vec{v} = \vec{v}_i\otimes \vec{v}_i + \vec{v}_s\otimes \vec{v}_s + \vec{v}_i\otimes \vec{v}_s + \vec{v}_s\otimes \vec{v}_i\,. \label{eq:vv} \end{align}
In the absence of an object from which sound is scattered, Eqs. \eqref{eq:scat:total:p} and \eqref{eq:scat:total:v} reduce to \(p = p_i\) and \(\vec{v} = \vec{v}_i\), and \(\vec{F} = \vec{0}\). Thus Eqs. \eqref{eq:L:s} and \eqref{eq:L:is} vanish, and \(L = L_i\). Similarly, Eq. \eqref{eq:vv} becomes \(\vec{v} \otimes \vec{v} = \vec{v}_i \otimes \vec{v}_i\). Equation \eqref{eq:int:F:L} becomes \begin{align} \oint_{A_0} \left[\dyad{I} \langle L_i\rangle - \rho_0 \left\langle \vec{v}_i\otimes \vec{v}_i\right\rangle \right]\cdot \vec{e}_n\,dA_0 = \vec{0}\,, \label{eq:no-force} \end{align} which according to Westervelt "is simply a statement of the fact that the unperturbed incident wave loses no momentum in passing through the control surface [\(A_0\)]" [10]. Combining Eqs. \eqref{eq:L:i}-\eqref{eq:no-force}, Eq. \eqref{eq:int:F:L} becomes \begin{align} \vec{F} = \oint_{A_0} \left[\dyad{I}\langle L_s + L_{is} \rangle - \rho_0 \langle \vec{v}_s\otimes \vec{v}_s + \vec{v}_i\otimes \vec{v}_s + \vec{v}_s\otimes \vec{v}_i\rangle\right] \cdot \vec{e}_r \, dA_0\,,\label{eq:int:F:L:1} \end{align} where it has been noted that \(\vec{e}_n\) equals \(\vec{e}_r\) on the surface \(A_0\) (See the figure at the beginning of this section). Next, since the scattered wave considered is in the far field (\(r/a \gg 1\)), \begin{equation}\label{eq:vs} \vec{v}_s = p_s \vec{e}_r /\rho_0c_0\,, \end{equation} where \(\vec{k}_s = k\vec{e}_r\) is the scattered wave vector, which is to say that the scattered wave is locally planar in the far field. Thus \(\langle L_s\rangle = 0\), because the kinetic and potential energy densities are equal in plane progressive waves. Making use of the definition of the outer product (See the introduction to the useful identities), Eq. \eqref{eq:vs}, and the fact that \(\vec{v}_i \cdot \vec{e}_r = v_i\cos\psi\), the outer products in Eq. \eqref{eq:int:F:L:1} become \begin{align} (\vec{v}_s\otimes \vec{v}_s) \cdot \vec{e}_r\, dA_0 &= (\vec{v}_s \cdot \vec{e}_r \, dA_0) \vec{v}_s = p_s v_s \, \vec{e}_r\, dA_0/\rho_0c_0 \label{eq:outer:simp:1}\\ (\vec{v}_i\otimes \vec{v}_s) \cdot \vec{e}_r\, dA_0 &= (\vec{v}_s \cdot \vec{e}_r\, dA_0) \vec{v}_i = p_s \vec{v}_i \, dA_0/\rho_0c_0 \label{eq:outer:simp:2}\\ (\vec{v}_s\otimes \vec{v}_i) \cdot \vec{e}_r\, dA_0 &= (\vec{v}_i \cdot \vec{e}_r\, dA_0) \vec{v}_s = p_s v_i \cos\psi\, \vec{e}_r\, dA_0/\rho_0c_0\,. \label{eq:outer:simp:3} \end{align} By similar manipulations, Eq. \eqref{eq:L:is} becomes \begin{equation} L_{is} \label{eq:L:is:2} = \frac{1}{c_0} (p_s v_i \cos\psi - p_i v_s) \,. \end{equation} Using Eqs. and \eqref{eq:outer:simp:1}-\eqref{eq:L:is:2}, along with the fact that \(\dyad{I} \cdot \vec{e}_r = \vec{e}_r\), Eq. \eqref{eq:int:F:L:1} becomes \begin{align*} \vec{F} &= \frac{1}{c_0}\oint_{A_0} \left\langle p_s v_i \vec{e}_r\cos\psi_0 - p_i v_s \vec{e}_r - p_s v_s \vec{e}_r - p_s \vec{v}_i - p_sv_i \vec{e}_r \cos\psi_0\right\rangle \, dA_0\,, \end{align*} where the subscript "\(0\)" has been included on \(\psi\) to indicate that it is an integration variable corresponding to the differential \(dA_0\). The first and last terms in the integrand cancel, resulting in \begin{align} \boxed{\vec{F} = -\frac{1}{c_0}\oint_{A_0} \left\langle p_s \vec{v}_i + p_i v_s \vec{e}_r + p_s v_s \vec{e}_r \right\rangle \, dA_0 \,,\label{eq:int:F:Westervelt}} \end{align} where the terms in Eq. \eqref{eq:int:F:Westervelt} have been rearranged to match the ordering of terms in Eq. (12) of Ref. [10].
Equation \eqref{eq:int:F:Westervelt} can be dotted into \(\vec{e}_i = \vec{k}_i/k\) to determine the component of the force in the direction of the incident wave. Noting that \(\vec{v}_i \cdot \vec{e}_i = v_i\) and \(\vec{e}_r \cdot \vec{e}_i = \cos\psi\), this component is \(F_\parallel = -c_0^{-1}\oint_{A_0} \left\langle p_s v_i + p_i v_s\cos\psi_0 + p_s v_s \cos\psi_0\right\rangle \, dA_0\). Westervelt prefers an alternate form, however. Since the scattered and incident fields are considered to be planar, one can write \(p_i \vec{e}_i = \vec{v}_i \rho_0c_0\) and \(p_s \vec{e}_r = \vec{v}_s \rho_0c_0\). Thus \begin{align} p_s\vec{v}_i = v_s \rho_0c_0 p_i\vec{e}_i/\rho_0c_0 = p_iv_s \vec{e}_i \label{eq:squares:switch:1}\\ p_i \vec{v}_s = v_i\rho_0c_0 p_s \vec{e}_r/\rho_0c_0 = p_s v_i\vec{e}_r \label{eq:squares:switch:2} \end{align} Inserting Eqs. \eqref{eq:squares:switch:1} and \eqref{eq:squares:switch:2} into Eq. \eqref{eq:int:F:Westervelt} yields \begin{align} \vec{F} = -\frac{1}{c_0}\oint_{A_0} \left\langle p_iv_s \vec{e}_i + p_s v_i\vec{e}_r + p_s v_s \vec{e}_r \right\rangle \, dA_0 \,,\label{eq:int:F:Westervelt:switch} \end{align} Dotting Eq. \eqref{eq:int:F:Westervelt:switch} into the incident wave unit vector \(\vec{e}_i\) gives the force in the direction of the incident wave, \begin{align} \boxed{F_\parallel = -\frac{1}{c_0}\oint_{A_0} \left\langle p_i v_s + p_s v_i \cos\psi_0 \right\rangle \, dA_0 -\frac{1}{c_0}\oint_{A_0} \left\langle p_s v_s \right\rangle \cos\psi_0 \, dA_0 \,,}\label{eq:int:F:Westervelt:parallel} \end{align} recovering Eq. (13) of Ref. [10].
The quantity \(\langle p_s v_s\rangle\) in the second integral of Eq. \eqref{eq:int:F:Westervelt:parallel} is identified as the time-averaged Poynting vector of the scattered wave. Since \(v_s = p_s/\rho_0c_0\), this quantity equals \(\langle p_s^2 \rangle/\rho_0c_0\), which can be recast in terms of energy density by noting that this quantity equals \(2c_0 \langle E_{s,V}\rangle\) from Eq. \eqref{eq:potential}. Since \(\langle E_{V}\rangle = \langle E_{T}\rangle = \langle E\rangle/2\) in plane progressive waves, \(\langle p_s v_s\rangle = c_0 \langle E_s\rangle\). Normalizing \(\langle E_s\rangle\) to \(\langle E_i\rangle = \langle p_i^2 \rangle/ \rho_0c_0^2 = p_0^2/2\rho_0c_0^2\), where \(p_0\) is the pressure amplitude of the incident plane wave, the second integral of Eq. \eqref{eq:int:F:Westervelt:parallel} can be written as \(-\langle E_i\rangle \oint \gamma \cos\psi_0\, dA_0\), where \(\gamma \equiv \langle E_s\rangle/\langle E_i\rangle\) is the "the magnitude of the mean asymptotic scattered intensity per unit incident intensity" [10].
Meanwhile, Westervelt identifies \(-c_0^{-1} \oint \langle p_i u_s + p_s u_i \cos\psi_0 \rangle dA_0 \) as \(\langle E_i\rangle (\pi_s + \pi_a)\) {See Eq. (14) of Ref. [10]}, where \(\pi_a\) and \(\pi_s\) are "the cross sections for absorption and scattering respectively" [10]. The argument for this identification is based on Lamb's discussion of cross terms in the product \((p_i + p_s)(\vec{v}_i + \vec{v}_s)\). Lamb writes that the surface integral over the time average of these cross terms in the far field equals "the total rate at which energy is withdrawn from the primary waves, in consequence of the presence of the obstacle" [14].
Similarly, the force in the direction perpendicular to \(\vec{v}_i\) can be calculated by dotting Eq. \eqref{eq:int:F:Westervelt} into a unit vector \(\vec{e}_m\) that is perpendicular to \(\vec{v}_i\). Since \(\vec{v}_i \cdot \vec{e}_m = 0\), the first term in Eq. \eqref{eq:int:F:Westervelt} is zero, leaving \begin{align} F_\perp &= -\frac{1}{c_0}\oint_{A_0} \left\langle p_i v_s \vec{e}_r\cdot \vec{e}_m + p_s v_s \vec{e}_r\cdot \vec{e}_m \right\rangle \, dA_0 \,,\label{eq:int:F:Westervelt:perp} \end{align} Westervelt claims in Ref. [10] that the surface integral over the first term in the integrand of Eq. \eqref{eq:int:F:Westervelt:perp} evaluates to zero. His argument is based on replacing the scattering object with a simple sources with source strength \(q\). Westervelt writes (The notation has been adjusted to match that of the present discussion) [10],
Consider the scattering object to be replaced by a volume distribution of simple sources. Choose the phase and amplitude of these sources in such a way that the field radiated by the sources is identical (this identity need only occur asymptotically at a large distance from the source distribution) with the scattered field of the object. This is equivalent to replacing the scattering object by a suitable collection of multipoles[,] although there is here no restriction that all the multipoles be located at the same point.
The force on the scattering object will be equal to the force on the source distribution by virtue of the identity of the far field in each case. The force on the sources is equal to the force arising from their radiation (equal to \(-\rho_0 \oint d\vec{A}_0\cdot \langle \vec{v}_s \otimes \vec{v}_s\rangle\)) plus the force arising from the interaction of the unperturbed incident wave with the sources. The latter force is the volume integral, over the source distribution, of an interaction force density equal to \(\rho_0 \langle q \vec{v}_i\rangle\), where \(\rho_0\) is the density of the medium, \(q\) is the instantaneous source strength density (that is, the volume of fluid generated per unit time per unit volume), and \(\vec{v}_i\) is the particle velocity of the incident wave. It is clear that, since the interaction force density is in the direction of the incident wave, the total interaction force can have no component at right angles to the incident wave. This proves that \(c_0^{-1} \vec{e}_m \cdot \oint \langle p_i v_s \rangle \vec{e}_rdA_0\) vanishes and that the perpendicular force is given simply by \begin{align} F_\perp = -c_0^{-1} \oint_{A_0} \langle p_s v_s\rangle \vec{e}_r\cdot \vec{e}_m \, dA_0\,. \end{align}
Meanwhile, Westervelt regards \(\vec{e}_r \cdot \vec{e}_m\) as \(\cos\beta\), where [10]
\(beta\) is the angle formed by \(\vec{e}_m\) and...the mean scattered intensity vector. Should the object scatter but a single beam at the angle \(\theta\) with respect to the incident wave, \(\cos\beta = \sin\psi\) in the direction of [the mean scattered intensity vector], and the expression given in 1951 [9] yields the correct answer.
By "the expression given in 1951," I think Westervelt is referring to the imaginary part of Eq. (13) of Ref. [9], but he does not point to which specific equation of the 1951 paper he is referring to. Westervelt elaborates [10],
The force at right angles to the incident beam is in a direction which depends on the nature of the scatterer. Its component in any particular direction is clearly given by \(-E_I \oint \gamma \cos\beta \, dA_0\), where \(\beta\) is the angle between the particular direction and \(\boldsymbol{\gamma}\). An obvious error was made in 1951 where this force was stated to be \(-E_I \oint \gamma \sin\psi \, dA_0\) which, as can be seen from Fig. 1, is true only if \(\gamma = 0\) everywhere except when \(\cos\beta = \sin\psi\); such a condition can for example occur when the object scatters but one specular component.
Westervelt's distinction between \(\beta\) and \(\psi\) (Westervelt's \(\theta\)) is not clear. Does Westervelt implicitly consider \(\theta\) to be the spherical polar angle? (There is nothing in his papers that indicate this. Here too, the argument is written independent of coordinate systems, and to emphasize this, \(\psi\) was used to denote the angle between the incident and scattered wave vectors, while \(\theta\) is reserved for the spherical polar angle). If Westervelt tacitly assumed that \(\theta\) is the spherical polar angle, then \(\vec{e}_r \cdot \vec{e}_m = \sin\psi = \cos\beta\), which is generally not equal to \(\sin\theta\), as is now seen:
Suppose the incident wave is oriented along the Cartesian unit vector \(\vec{e}_z\). Only if \(\vec{e}_m\) is taken to be the cylindrical radial unit vector \(\vec{e}_s = \cos\phi \vec{e}_x + \sin\phi \vec{e}_y\) (which is perpendicular to \(\vec{e}_z\)) does \(\vec{e}_r \cdot \vec{e}_s = \sin\theta\). Indeed, the spherical radial unit vector in terms of the Cartesian unit vectors is \(\vec{e}_r = \sin\theta \cos\phi \vec{e}_x + \sin\theta \sin\phi \vec{e}_y + \cos\theta \vec{e}_z\), the quantity \(\vec{e}_r \cdot \vec{e}_s\) is \begin{align} \vec{e}_r \cdot \vec{e}_s &= \sin\theta \cos^2\phi + \sin\theta \sin^2\phi = \sin\theta\,. \end{align} Meanwhile, if \(\vec{e}_m = \vec{e}_x\) (which is also perpendicular to \(\vec{e}_z\)), then \[\vec{e}_r\cdot \vec{e}_x = \sin\theta \cos^2\phi \neq \sin\theta\,.\] Is an argument like this along the lines of what Westervelt means?
Under construction
In the derivation of the acoustic radiation stress tensor, nonlinear terms in the momentum \eqref{eq:mom:exact} and continuity \eqref{eq:cont:exact} equations must be retained to obtain a nonzero stress tensor. Since radiation force is given by the surface integral of the acoustic radiation stress tensor [See Eq. \eqref{eq:mom:stress:int}], the radiation force is \(0\) if only linear terms in the momentum equation are retained. Thus radiation force is a nonlinear effect.
Another argument utilizes the superposition principle, one of the two mathematical criteria of a linear function: \[f(a+b) = f(a) + f(b)\,.\] Radiation forces do not superpose, as can be seen by considering two plane waves with pressure amplitudes \(p_1\) and \(p_2\). If \(p_2 = 0\), the force on the object can be calculated by Eq. \eqref{eq:mom:stress:int} to be \(F_1\), and if \(p_1 = 0\), the force on the object can similarly be calculated to be \(F_2\). The radiation force \(F_{12}\) on an object insonified by both \(p_1\) and \(p_2\) is not \(F_1 + F_2\). Instead, the radiation stress tensor given by Eq. \eqref{eq:stress} must be re-calculated, and it will contain cross terms proportional to \(p_1p_2\) (and \(v_1v_2\), the corresponding velocity fields). Radiation force is a nonlinear effect because the superposition principle, a hallmark of linear phenomena, does not apply.
While admitting the above arguments, one might still hesitate to categorize radiation force as "nonlinear acoustics," a term that implies that the governing wave equation is nonlinear (For example, see the introduction of the Wikipedia entry for "Nonlinear Acoustics"). The acoustic radiation stress tensor \eqref{eq:stress} and mean excess pressure \eqref{eq:excess} are formed by the squares of quantities that are solutions of the linear wave equation, and thus the calculation of radiation force does not involve solving any nonlinear wave equations.
One might argue that linear quantities are squared when solving nonlinear wave equations by the perturbation techniques—see Sec. 2 of Ref. [15]. But these squares are then substituted into the nonlinear wave equations, which are in turn solved for higher-order solutions. Radiation force, in contrast, involves no "higher-order solutions" for the acoustic fields.
For another perspective, recall from linear acoustics that the linear wave equation (which admits that the field variable \(p\) is much smaller than the ambient pressure \(p_0\)) yields solutions that can be characterized by a nonzero energy, intensity, and power {See Sec. 1-E-3 of Ref. [12]}. These quantities are all proportional to \(p^2\) and are therefore of the same order as the mean excess pressure, suggesting that radiation pressure/force is non-zero in linear acoustics. For radiation force to vanish in linear acoustics, \(p\) must be sufficiently small such that \(p^2 = \Order(0)\), in which case the energy, intensity, and power would also vanish. So, if one admits the solutions of the linear acoustic wave equation carry energy, then one must also admit radiation force.
The same is true about the small angular deflection \(\theta\) of a pendulum: from Eq. (5), \(F\simeq mg\langle\theta^2\rangle/2\) predicts zero radiation force if \(\theta^2 = \Order(0)\). If we admit that the pendulum possesses potential and kinetic energy (which are proportional to \(\theta^2\)), then we must also admit the existence of radiation force, because all of these quantities are of the same order. When we make the small-angle approximation \(\theta \ll 1\), we mean to say that \(\theta\) is sufficiently small such that \(\sin\theta \simeq \theta\) (leading to a linear equation of motion), but that \(\theta^2 \neq 0\), such that the pendulum can have a nonzero energy. Note that \(\sin\theta \simeq \theta\) is consistent with \(\theta^2 \neq 0\), because the next-highest order term in \(\sin\theta\) is \(-\theta^3/6\), which is neglected.
On another (possibly related) note, Beyer writes [5],
[A]lthough the phenomenon is primarily one of nonlinear acoustics, it can be observed down to the lowest sound intensities under certain conditions. Thus, the Rayleigh radiation pressure vanishes for the linear case, but the usually measured Langevin pressure does not.