Active Nonlinear Acoustic Sensing of an Object with Sum or Difference Frequency Fields

: A number of nonlinear acoustic sensing methods exist or are being developed for diverse areas ranging from oceanic sensing of ecosystems, gas bubbles, and submerged objects to medical sensing of the human body. Our approach is to use primary frequency incident waves to generate second order nonlinear sum or difference frequency ﬁelds that carry information about an object to be sensed. Here we show that in general nonlinear sensing of an object, many complicated and potentially unexpected mechanisms can lead to sum or difference frequency ﬁelds. Some may contain desired information about the object, others may not, even when the intention is simply to probe an object by linear scattering of sum and difference frequency incident waves generated by a parametric array. Practical examples illustrating this in ocean, medical, air and solid earth sensing are given. To demonstrate this, a general and complete second-order theory of nonlinear acoustics in the presence of an object is derived and shown to be consistent with experimental measurements. The total second-order ﬁeld occurs at sum or difference frequencies of the primary ﬁelds and naturally breaks into (A) nonlinear waves generated by wave-wave interactions, and (B) second order waves from scattering of incident wave-wave ﬁelds, boundary advection, and wave-force-induced centroidal motion. Wave-wave interactions are analytically shown to always dominate the total second-order ﬁeld at sufﬁciently large range and carry only primary frequency response information about the object. As range decreases, the dominant mechanism is shown to vary with object size, object composition, and frequencies making it possible for sum or difference frequency response information about the object to be measured from second-order ﬁelds in many practical scenarios. It is also shown by analytic proof that there is no scattering of sound by sound outside the region of compact support intersection of ﬁnite-duration plane waves at sum or difference frequencies, to second-order. Analytic expressions for second-order ﬁelds due to combinations of planar and far-ﬁeld wave-wave interactions are also derived as are conditions for when wave-wave interactions will dominate the second order ﬁeld.


Introduction
A number of nonlinear acoustic sensing methods exist or are being developed for diverse areas ranging from oceanic sensing of ecosystems, gas bubbles, and submerged objects to medical sensing of the human body [1][2][3][4][5][6][7].The approach is to use primary frequency incident waves to generate second order nonlinear sum or difference frequency fields that carry information about an object to be sensed.Here we show that in general nonlinear sensing of an object, many complicated and potentially unexpected mechanisms can lead to sum or difference frequency fields.Some may contain desired information about the object, others may not, even when the intention is simply to probe an object by linear scattering of sum and difference frequency incident waves generated by a parametric array [1,3].Practical examples illustrating this in ocean, medical, air and solid earth sensing are given.To accomplish this, a general and complete second order theory of nonlinear acoustics in the presence of an object is derived.The theory employs consistent asymptotic analysis of the wave equation and time-dependent boundary conditions to second order, following methods developed in fluid dynamics to quantify second order nonlinear waves found in the presence of floating objects [8,9].This second order nonlinear acoustic theory is confirmed by quantitative comparison with classic experimental measurements of Jones and Beyer [1].To our knowledge, no previous nonlinear acoustic theory derived from first principles has been available to explain these measurements.While illustrative examples are computed from exact second order solutions, which may involve integration, analytic solutions and approximate asymptotic expressions are also derived when possible to gain further insight.A detailed discourse on the relationship of this work to previous developments is provided in the Discussion section.
When an object is insonified by two primary frequency incident waves, second order acoustic waves at the sum and difference frequencies of the primary waves arise due to multiple mechanisms.It is found that these generating mechanisms fall into two categories.The first is (A) nonlinear wave-wave interactions of the primary waves, including the interaction between two incident waves (denoted as II), the interaction between two scattered waves (denoted as SS) and the interactions between an incident and scattered wave (denoted as IS and SI).The second is (B) linear scattering of nonlinear waves of the first category due to the presence of the object with radiation from second order boundary advection and wave-force-induced centroidal motion (together denoted as S2).Second order acoustic waves of category (A) contain information about the object's first order response at the primary frequencies but not at the sum or difference frequency.Second order acoustic waves of category (B) contain information about the object's linear response to a sum or difference incident wave but this information is often contaminated by complicated nonlinear effects due to primary scattered wave fields at the object boundary so that it cannot typically be extracted simply by factoring from S2 the II amplitudes in analogy to approaches used in linear scattering theory.
We show analytically that wave-wave interactions dominate the total second order field at sufficiently long ranges from the object, but in some scenarios these ranges may be too large to enable practical measurements of this.At shorter ranges, second order scattering may also dominate depending on physical parameters such as object size, object composition, primary frequencies, and range.We show that centroidal motion of a rigid object due to its interaction with acoustic waves can significantly affect the second order field at shorter ranges.In this case, S2 can be decomposed into scattering from a fixed object at its time-averaged position (denoted as D2), and radiation due to the second order centroidal motion (denoted as R2).We also show that resonance effects from elasto-mechanical effects or gas filled bubbles in liquids can dominate the second order field at shorter ranges.
Following experimental practice, we employ a time domain formulation of finite-duration narrow-band incident waves and provide analytic solutions for the resulting second order fields in the presence of an object.Finite-duration incident waves enable nonlinear field components that contain information about the object to be separated in time and space from those that do not, e.g., the II component, for certain sensing geometries.In such a finite-time duration approach, it is important to demonstrate that all potential effects of sound scattering by sound in the II field are accounted for properly.To quantify the scattering of sound by sound inside and outside a bounded interaction region of primary waves with compact support, we derive analytic expressions for the second-order nonlinear field arising from the interaction of plane waves of arbitrary time dependence.In an early investigation of the scattering of sound by sound, Ingard and Pridmore-Brown [10] employed unphysical collimated primary time-harmonic fields, which led to unphysical results as noted by Westervelt [11].The subsequent time-harmonic plane wave analysis of Westervelt [11] provided no spatial region without interacting primary fields.By analytic proof we show that to second order there is no scattering of sound by sound outside the region of compact support intersection of finite-duration plane waves at sum or difference frequencies in the absence of an object.We also show that non-collinear interaction leads to finite scattering of sound by sound at the primary frequencies only within the region of compact support union through which compact support intersection occurred if the primary waves had zero-frequency spectral components.We rigorously prove that Lamb's collinear [12] and Westervelt's non-collinear [11] time-harmonic solutions at the sum or difference frequency for ideal time-harmonic primary plane wave interactions provide good approximations to the second order field arising from the interaction of narrow-band finite-duration primary plane waves with sufficiently long and smooth envelopes, but only in the region of compact support intersection where the primary fields exist.We show such a time-harmonic approximation can be made in many practical scenarios involving the interaction between scattered waves, as well as scattered and incident waves of narrow band.
Applications of practical sensing in the ocean, human body, atmosphere and solid earth using the second order nonlinear acoustic theory derived here are analyzed and evaluated in roughly twenty case studies, where transitions between at least five potentially dominant mechanisms are explored for each.These examples are illustrated and evaluated both in terms of dimensionless parameters and specific frequencies and lengths.The former makes it possible to interpret the results universally over the broadest range of possible applications by scaling the object size, receiver range, and frequencies as necessary.The range of dimensionless parameters considered are chosen to provide a general characterization of regions where important transitions between dominant mechanism occur.The specific cases discussed include object radii spanning hundreds of micrometers to several meters; primary frequencies ranging from tens of Hertz to tens of mega Hertz; difference frequencies ranging from tens of Hertz to hundreds of kilo Hertz and sum frequencies of tens of mega Hertz.
Nonlinear acoustic sensing of oceanic fish at the difference frequency of a parametric array is studied for two representative species: Atlantic herring (Clupea harengus) and Atlantic cod (Gadus morhua).The former are of the plankton feeding Clupeidae family of keystone species that serve as prey to many other fish, mammals and birds.The latter are of the fish-feeding Gadidae family of top level predators.Both families are found throughout the world's oceans and are important ecologically in circulating biomass between open ocean and coastal waters via feeding to spawning cycles, and as human food sources [13,14].Traditionally, downward directed echosounders in the tens to hundreds of kilo Hertz range have been used for scientific analysis of fish populations.These frequencies, however, are typically well above the swim-bladder resonance frequencies of most oceanic fish species.Difference frequency sensing with a parametric array has been used to potentially investigate swim bladder resonances [7], which are important for wide area sensing with Ocean Acoustic Waveguide Remote Sensing (OAWRS) [15][16][17][18].Here we show that potentially unexpected second order nonlinear mechanisms, different from the intended linear scattering of the parametric array's incident difference frequency field, can mask the difference frequency response of oceanic fish and so make difference frequency sensing of swimbladder resonance difficult or impossible.
Nonlinear acoustic sensing of ocean air bubbles at difference frequencies is also studied.Air bubbles are formed by breaking waves and are ubiquitous in the world's oceans.Quantifying air bubble distributions in the ocean is important to a number of areas including underwater sensing and communication where bubbles cause attenuation, dispersion and randomization of transmitted signals, and climatology where bubbles are a source of critical gas/aerosol exchanges between the ocean and atmosphere [19][20][21][22].Oscillating air bubbles are also a major natural source of ambient noise in the ocean [23].Here it is shown that transitions between dominant nonlinear acoustic mechanisms occur at practical sensing ranges and frequencies and can limit the ability to determine desired bubble parameters.
Applications of nonlinear acoustic difference frequency sensing for medical ultrasound applications [4,5] are investigated for immovable rigid, movable rigid, vacuous, air-filled, and resonating objects.Here, we provide many examples showing that in practical medical applications no single mechanism is always dominant.Instead, a number of mechanisms are found to dominate the received difference frequency field depending on the object composition, size, source and receiver range and orientation, and primary and difference frequencies, which may lead to unexpected consequences for difference frequency sensing.
Nonlinear acoustic difference frequency sensing in the atmosphere for robotic applications, and in the solid earth for detection and sensing of buried structures are also provided, as are opportunities for sensing objects using incident-scattered wave interactions in a variety of media.

Linear and Second Order Nonlinear Formulation
Assuming weak nonlinearity, the total pressure field p(r, t), density field ρ(r, t), fluid particle velocity v(r, t), object boundary displacement ξ(r, t) and total force on the object f(t) due to the acoustic field can be described by [24] where subscripts 0, 1 and 2 indicate zeroth-, first-and second-order quantities, respectively, with {} 0 {} 1 {} 2 .Here p 0 and ρ 0 are ambient pressure and density of the medium without acoustic waves and v 0 = ξ 0 = f 0 = 0.The general problem is sketched in Figure 1, where two primary incident fields p Ia and p Ib are incident on an object, causing two primary scattered fields, p Sa and p Sb .The second order field p 2 generated by wave-wave interactions and second order scattering is measured by a receiver at r R .In practice, the different orders expressed in Equation ( 1) can be isolated by slowly increasing the amplitude of an incident time-harmonic primary field.For very small incident amplitude, only the ambient or zeroth order fields will be measurable.As the incident primary field amplitude is slowly increased, it will be possible to measure the first order quantities in Equation (1) above ambient levels at the primary frequency with no other higher order contributions present.Slowly increasing the primary field amplitude further will eventually lead to the addition of only sum and difference frequency fields that are measurable above the ambient levels that correspond to the second order variables of Equation (1).Increasing the incident amplitude further such that there are more frequencies present than the primary, sum and difference frequency introduces nonlinearities beyond the second order, which are not treated in the present paper.In this case, incident amplitudes should be reduced so that only primary, sum and difference frequencies are present.The second order theory presented here is then valid and can be properly used to explain nonlinear phenomena.

First Order Formulation
The total first order pressure field p 1 satisfies the wave equation with source function q 1 , via where c 0 is the sound speed, ∇ 2 is the Laplacian with respect to the receiver position r and t is time.
The boundary conditions for pressure release, rigid immovable, or rigid movable objects without rotation are pressure release (Dirichlet): rigid immovable (Neumann): rigid movable: where the fluid velocity is related to the pressure by the linearized momentum equation v 1 = −ρ −1 0 ∇p 1 dt, u c1 is velocity of the object's centroidal motion in response to the acoustic field, S is the mean reference boundary position, and n is the surface outward normal vector of the medium.Pressure release and rigid immovable boundaries are limiting cases of a general locally reacting boundary.It is shown in Appendix C that the solution to the first and second order scattering problems for any locally reacting boundary can be constructed from pressure release and rigid immovable scattering solutions.A rigid movable boundary, however, is a special case of a non-locally reacting boundary because u c1 depends on the acoustic field on the whole boundary.
The complete solution for p 1 from Green's theorem is [25,26] where dt 0 dV 0 g(r, t|r 0 , t 0 )q 1 (r 0 , t 0 ), Here, the 0 subscript appears on spatial variables of integration and their gradient; t 0 is a temporal variable of integration for the first two integrals; t init is the time at which the source begins to operate; t + is a time much larger than the period of any harmonic source in operation; p I1 is the known incident wave, p S1 is the scattered wave, which depends on the object boundary condition; and p init1 represents the transient effect associated with initial conditions at t init .The time domain Green function is g(r, t|r 0 , t 0 ) = δ(t − t 0 − R/c 0 )/(4πR) where δ is the Dirac delta function and R = |r − r 0 |.
For pressure release or rigid immovable objects, the scattered wave p S1 can be obtained by solving the integral equation (8) with boundary conditions pressure release: rigid immovable: where For rigid movable objects, let the total scattered wave p S1 be decomposed as where subscript D1 represents first order scattering from the rigid object whose boundary is fixed, and subscript R1 represents first order radiation due to centroidal motion of the rigid object.By definition, p D1 is the same as p S1 for the rigid immovable object.The radiated wave p R1 depends on the u c1 , and it satisfies the boundary condition where v R1 = −ρ −1 0 ∇p R1 dt.The centroidal velocity u c1 is determined by the equation of motion where M is the mass of the object and the integral on the right-hand side of Equation ( 14) represents the total force on the object.The inertial term on the left-hand side of Equation ( 14) can be written as z m u c1 with operator z m = Md/dt.Since the radiated wave is a linear function of the centroidal motion, the force on the body due to radiation can be written as S p R1 ndS = −z r u c1 , where the linear operator z r depends on the object geometry and the medium properties.Then u c1 becomes where the wave-exciting force is defined as which is completely determined by p I1 and p D1 .Once u c1 in Equation ( 13) is determined, p R1 can be obtained with the same procedure as that for p D1 .
The formulation for time-harmonic first-order scattering including centroidal motion of the object is provided in the online supplementary material.

Second Order Formulation
The governing equation for the second order acoustic pressure field p 2 (r, t) in a homogeneous, inviscid medium is [11 where and J 2 = ρ 0 v 1 v 1 is the momentum tensor [24], and A = ρ 0 (∂p/∂ρ)| ρ 0 = ρ 0 c 2 0 and B = ρ 2 0 (∂ 2 p/∂ρ 2 )| ρ 0 are constants from the equation of state.The second order fluid velocity v 2 is related to p 2 via To find the second order boundary condition on the body, we expand the boundary condition on the exact boundary S(t) in a Taylor series with respect to S [8].After substituting the perturbation expansions for p, v and ξ into the series and collecting terms at each order, the second order boundary conditions are obtained as pressure release: rigid immovable: rigid movable: where u c2 is the second order centroidal velocity for rigid objects without rotation, ξ 1 and ξ c1 are first order boundary and centroidal displacements, respectively.The second order boundary condition for a general locally reacting boundary is derived in Appendix C, where it is shown that solution to the second order scattering problem for any locally reacting boundary can also be constructed from pressure release and rigid immovable scattering solutions.
As in the first order problem, the total second order wave field is given by Green theorem, as where dt 0 dV 0 g(r, t|r 0 , t 0 )q 2 (r 0 , t 0 ), Here p I2 is the second order incident wave, which is due to wave-wave interactions of the first order field, p S2 is the second order scattered field determined by the second order boundary condition, and p init2 accounts for the initial conditions.
When the primary wave field consists of incident and scattered waves, we have Substituting p 1 and v 1 from (27) into Equation (18) then into Equation ( 24) leads to a decomposition of p I2 as where in which p I1 , p S1 , v I1 , and v S1 are functions of r 0 and t 0 , ∇ 0 is the gradient operator with respect to r 0 , II represents contributions from incident-incident interactions, SS from scattered-scattered interactions, and IS and SI from incident-scattered interactions.The procedure to obtain p S2 for pressure release or rigid immovable objects is similar to that for p S1 (Section 2.1 and online supplementary material), but the boundary conditions are pressure release: rigid immovable: where v I2 is the second order fluid velocity due to p I2 , which is given by Equation (19) with The pressure component ξ 1 • ∇p 1 describes contributions from wave-boundary interactions due to first order boundary motion.
Following the approach used in first order for rigid movable objects, p S2 is decomposed as where p D2 is the second order scattered wave from a rigid object whose centroid is fixed at second order, i.e., u c2 = 0, and p R2 is the second order radiated wave due to any second order centroidal motion u c2 .
The procedure to obtain p D2 is similar to that for p S2 and p S1 , but the boundary condition is where • n) also describes contributions from wave-boundary interactions due to first order boundary motion.
To obtain p R2 , the procedure is again similar to that for p R1 , but with boundary condition where v R2 = −ρ −1 0 ∇p R2 dt.The centroidal velocity u c2 is determined by the second order equation of motion where the second order pressure component ξ c1 • ∇p 1 arises from Taylor series expansion of the total pressure p at the instantaneous boundary S(t) with respect to S. By defining the second order wave-exciting force as When the primary incident field has two components p I1 (r, t) = { pIa (r, t) + pIb (r, t)}, the primary scattered field is p S1 (r, t) = { pSa (r, t) + pSb (r, t)}, and the total first order field is and v 1 v 1 can be expressed in the same manner.Substituting p 2 1 and v 1 v 1 into q 2 in Equation ( 18) and then Equation (24) The first parenthetical group of terms in Equation ( 42) is due to II interaction, the second to SS interaction, and the last to IS and SI interactions.The II components p II,aa , p II,aa When p I1 has two harmonic components at angular frequencies ω a and ω b , ω a > ω b , the variables p I1 , p S1 , v I1 , v S1 , ξ 1 and ξ c1 also contain these harmonic components, i.e., Then p I2 , p II , p SS , p IS and p SI in Equation ( 28) contain static and double frequency (2ω a and 2ω b ) components due to self-interactions, and sum (ω + = ω a + ω b ) and difference (ω − = ω a − ω b ) frequency components due to cross-interactions.Specifically, the ω ± components for p I2 , p II , p SS , p IS and p SI are where in which the asterisk denotes complex conjugation that applies to the difference frequency case only.Equations ( 45)-( 48) are obtained by substituting p I1 , p S1 , v I1 and v S1 from (43) into Equations ( 29)-( 32) and integrating over t 0 .
The ω ± components of p S2 , p D2 , p R2 , u c2 and f excit 2 can be written as where the complex amplitudes P S2± , P D2± , P R2± , U c2± and F excit 2± can then be determined by further analysis.Specifically, substituting p S2± into Equation ( 25) and integrating over t 0 yields where G ± (r|r 0 ) = e ik ± R /(4πR) and k ± = ω ± /c 0 .The scattered wave P S2± for pressure release or rigid immovable objects can be obtained from Equation ( 50) with boundary conditions pressure release: rigid immovable: where V II± , V SS± , V IS± and V SI± are the second-order complex fluid velocity amplitudes associated with P II± , P SS± , P IS± and P SI± respectively, P a = P Ia + P Sa and P b = P Ib + P Sb are first order complex pressure field amplitudes at frequencies ω a and ω b respectively, and In the pressure release boundary Equation ( 51), the P I2± term leads to linear sum or difference frequency scattering of incident wave-wave fields II, SS, IS and SI, while the last two terms lead to nonlinear radiation from boundary advection that will involve II, SS, IS and SI type components, where SS or SS-type boundary components dominate the amplitude of P S2± for k a , k b 1.In the rigid immovable boundary Equation ( 51), there is only linear sum or difference frequency scattering of incident wave-wave fields.
For rigid movable objects, P D2± can be obtained from Equation (50) with P S2± = P D2± and boundary condition where are first order fluid velocities at frequencies ω a and ω b respectively.The wave-exciting force is The centroidal velocity for a rigid movable object attached to a spring and damper in the direction of motion is where ω n is the undamped natural frequency, ζ is the dimensionless damping ratio, Z m (ω ± ) = −iω ± M is the Fourier transform of z m and Z r (ω ± ) is the Fourier transform of z r , which is sometimes referred to as the radiation impedance [24,27,28] and is determined from its elements by where P R,j (ω), for j = x, y, z are the first order radiated fields from an object due to unit amplitude velocity oscillation at frequency ω in x, y, z directions (online supplementary material).In terms of PR = ( PR,x , PR,y , PR,z ), the second order radiated field P R2± is

Summary
In summary, the total second order field from any object insonified by acoustic waves is generated by wave-wave interactions, including II, SS, IS and SI, and by second order scattering S2, as which, for the sum and difference frequency amplitudes, leads to The interaction components p II , p SS , p IS and p SI are given by Equations ( 29)- (32).The S2 component p S2 is given by Equation ( 25) with boundary conditions ( 33) and ( 34) for pressure release and rigid immovable objects, respectively.The general solution for objects with non-rigid locally reacting boundaries is given in Appendix C in terms of the pressure release and rigid immovable solutions of this section.For rigid movable objects, S2 is decomposed into D2 and R2.The D2 component p D2 is given by Equation ( 25) with boundary condition (36).The R2 component p R2 is given by Equation ( 25) with boundary condition (37), where the centroidal motion is determined from the wave-exciting force by Equation (39).
From the sensing perspective, the II interaction only depends on the primary incident waves and it contains no information about the object.The SS, IS and SI interactions depend on the first order scattered waves p Sa and p Sb , and they contain primary frequency scattering information of the object.S2 contains the object's second order scattering information at the sum or difference frequency.For example for the rigid object, D2 depends on the geometry of the object and R2 depends on the second order motion of the object.
The present theory provides a complete and self-consistent treatment of second order interaction and scattering in the presence of an object, including the effect of second order centroidal motion.More details about the wave-exciting force formulation are provided in Appendix D.

General Approach
Let p I2 of Equation ( 23) and the associated second order velocity v I2 due to wave-wave interactions be further decomposed as [29] where p 2 and v 2 are nonzero only when the primary field exists and can be directly evaluated from first order pressure p 1 and velocity v 1 , as and p 2 satisfies where β = 1 + B/(2A).As in Equation ( 24), solution for p 2 is Finite-duration incident plane waves propagating at directions îa and îb with center frequencies ω a and ω b respectively can be written as where P Ia = P a0 e iω a îa •r/c 0 and P Ib = P b0 e iω b îb •r/c 0 (66) are harmonic plane waves with amplitudes P a0 and P b0 , respectively, and w 1 (t) is a window function with compact support of duration T, which is here assumed to have unit height within the window except at each end where smooth transitions to zero occur over time periods much less than T.
When the object is small compared to the spacial extent of the incident plane wave window, the dominant scattering components are within the narrow band of the incident waves, and the window functions are sufficiently smooth, the primary scattered waves can be approximated as which correspond to harmonic scattered waves P Sa and P Sb at frequencies ω a and ω b , where compact support is preserved by moving windows.Like p I2 in Equation ( 42), p 2 and p 2 have self-interaction and cross-interaction components due to II, SS, IS and SI mechanisms.By substituting Equations ( 65) and (67) into Equation (64) and integrating over t 0 , the cross-interaction components for p 2 are where When the primary incident and scattered waves are purely time harmonic, w II,ab where The cross-interaction components of p 2 are denoted by p II,ab ( * ) , p SS,ab ( * ) , p IS,ab ( * ) and p SI,ab ( * ) .They can be determined by substituting appropriate first-order field products into equation (61).When the primary fields are time harmonic, where P II± , P SS± , P IS± and P SI± are the complex amplitudes at the sum and difference frequency for the cross-interactions, which again can be directly obtained from Equation (61).

Analytic Solutions
While time domain and frequency domain solutions for wave-wave interactions can be directly evaluated numerically via Equations (68)-(71) and Equations ( 77)-(80) respectively, which is the approach we have used in all rigid and pressure release simulations presented in this paper, we have obtained analytic solutions for some cases not considered previously.We also provide approximate analytic expressions for the S2 pressure for small pressure release or gas filled bubbles.

Incident-Incident Interaction of Plane Waves of General Time Dependence
For two collinear time harmonic plane waves given by Equation (66) with îa = îb = îz , P II± was derived by Lamb [12], who showed that it grows linearly with range.For zero second order normal velocity at z s , P II± is Here the complex conjugate and the "−" sign in "±" apply to the difference frequency case only.For two non-collinear time harmonic plane waves given by Equation (66) with îa = îx and îb = îz , Westervelt [11] found that where z = x cos θ + z sin θ and θ is the angle between îx and îz .These time-harmonic results are here generalized to plane waves with arbitrary time dependence to investigate the scattering of sound by sound and also properly model typical experimental scenarios.For two collinear plane waves of arbitrary time dependence p Ia = { pIa (t − z/c 0 )} and p Ib = { pIb (t − z/c 0 )}, the second order fields due to their cross-interaction are (online supplementary material).
where the complex conjugates only apply to p II,ab * .For narrow-band waves with compact support given by Equation (65) with îa = îb = îz , and a sufficiently smooth temporal window w 1 , Equation (84) becomes (online supplementary material) where the iβk ± z term in the square brackets corresponds to p II,ab ( * ) and the other term in the square brackets corresponds to p II,ab ( * ) .It can be seen from Equations ( 84) and (85) that sum and difference frequency components of p II,ab ( * ) only exist within the region of compact support intersection.For narrow-band primary waves with sufficiently smooth windows, the second order field in the intersection corresponds to Lamb's solution with boundary condition p II,ab ( * ) = p II,ab ( * ) at z = z s = 0.For two non-collinear intersecting plane waves of arbitrary time dependence p Ia = { pIa (t − x/c 0 )} and p Ib = { pIb (t − z /c 0 )}, the second order fields due to their cross-interaction are found to be (online supplementary material).
which is the full solution including both p II,ab ( * ) and p II,ab ( * ) components.Sum and difference frequency components of p II,ab ( * ) only exist within the region of compact support intersection, so that no scattering of sound by sound at sum or difference frequencies is found outside the region of compact support intersection to second order.Primary frequency components in p II,ab ( * ) exist within the region of compact support union through which compact support intersection occurred if the primary waves had non-zero zero-frequency components.There is no scattering of sound by sound at any frequency outside the region of compact support union through which compact support intersection occurred.
For narrow-band plane waves with compact support given by Equation ( 65) with îa = îx and îb = îz , where window w 1 is sufficiently long and smooth, Equation ( 86) becomes (online supplementary material) which shows that a time-harmonic approximation involving Westervelt's time harmonic solution (Equation ( 83)) can be made in the region of compact support intersection for narrow-band plane waves with sufficiently long and smooth windows.

Scattered-Scattered Interaction
Let r a and r b be the far field ranges [30] for the primary scattered fields P Sa and P Sb respectively, such that beyond r a and r b , far field approximations [31] apply e ik a r r and P Sb (r where r a,b = l 2 /λ a,b , k a r a , k b r b 1, l is the length scale of the object, λ a,b are the wavelengths of the incident fields, P a0 and P b0 are the amplitudes of the incident fields, S a and S b are the far field scatter functions, and îr = r/r. The integral representation of P SS± , Equation (78), can then be decomposed into where P SS± is due to the SS interaction within r ref± and P SS± is due to SS interaction beyond r ref± , and r ref± satisfies k ± r ref± 1 and r ref± ≥ r far,a , r far,b .The P SS± term for r ≥ r ref± can then be analytically approximated via a stationary phase or spherical wave expansion approach (online supplementary material), yielding where E 1 is the exponential integral [32,33] and for sufficiently small objects r ref± is the object radius.For k ± r → ∞, since the logarithmic function dominates the exponential integrals in Equation (90), P SS± falls off by log(r)/r.The component P SS± falls off by r −1 and P SS± by r −2 .The total second order field due to the SS interaction can be approximated by which is consistent with a general diverging waveform suggested by Westervelt and Radue [34] but without explicit derivation or interpretation.
When the two primary fields radiating from a sphere with radius a are spherically symmetric, S a = S b = constant, r ref± = a and P (1) SS± = 0. Equation (90) reduces exactly to Baxter's solution [33].Compared to Baxter's solution, it is found that Dean's earlier solution [29] has an extra term proportional to the spherical Bessel function j 0 (k ± r) so violates the Sommerfeld radiation condition (online supplementary material).A further approximation made by Dean for the far field k ± r 1 and a small object k ± a 1, however, satisfies the Sommerfeld radiation condition because it contains only the dominant e ik±r log(k ± r)/r term.Jones and Beyer [1,35] developed a heuristic formula based on the radial dependence of Dean's omnidirectional solution with far field and small object approximations, and applied it to model the sum frequency interaction of two scattered fields from a large sphere (k + a = 160), where they found good agreement in angular dependence with data by introducing an arbitrary scaling factor.

Incident-Scattered Interaction
For incident plane waves P Ia and P Ib given by Equation (66) and scattered waves P Sa and P Sb given by Equation (88), the second order fields P IS± and P SI± in the forward directions are analytically derived (online supplementary material) from Equations (79) and (80).At large range, P IS± and P SI± have constant magnitude in range while P IS± and P SI± fall off by r −1 .The total second order fields due to the IS and SI interactions are then Unlike the interaction of collinear plane waves where growth is found along the propagation path, collinearity between planar and spherical wavefronts within an equivalent Fresnel area about the forward direction, together with spreading of the spherical wave, balances out second-order wave growth.For P IS− , the π term in Equation ( 93) is due to the additional contribution from a stationary phase point at range rk b /k a in the forward direction.
The backscatter directions − îa and − îb for the IS and SI interactions respectively lack collinearity and stationary phase contributions so that P IS± and P SI± fall off in range by r −1 as do P IS± and P SI± , and have much lower magnitude than P IS± and P SI± in the forward directions.Detailed derivations for IS and SI interactions are provided in the online supplementary material.

Sound Pressure Level of the IS and SS interactions
The Sound Pressure Level (SPL) of P IS± in the forward direction is determined from Equations ( 92) and (93) as which depends on the SPL of incident waves, the frequencies ω a , ω b , the far field scatter function S b and constants β and A.
At long range and away from the forward direction, P SS± is dominant.The SPL for P SS± for k ± r → ∞ can be determined from Equation (91) as Compared to the SPL of P IS± in Equation ( 95), it can be seen that a higher incident SPL is required for P SS± to compensate for its falloff.

S2 for Small Pressure Release or Gas Filled Bubbles
Consider small pressure release objects with k a a, k b a 1.The righthand side of Equation ( 51), P I2± is dominated by the SS component of ρ 0 v 2 1 (r, t)/2 from Equation (61) and the last two advective terms of Equation ( 51) are dominated by SS-type components, so that where A ± res = 1 and B ± res = 0 for the pressure release case.If the small object is a gas filled bubble, or a fish swim bladder, an additional constraining equation at the boundary is necessary, which makes the dimensionless coefficients A ± res and B ± res dependent on sum, difference and primary frequencies as well as parameters of medium and object such as resonance frequency, damping, and radius as shown in Appendix B via Equations (B28) and (B29) respectively.In these cases, P S2± may attain maxima for probing frequencies at resonance, and approaches the pressure release case for high probing frequencies and zero for low probing frequencies.

Confirmation of Theory with Measurements
The exact second order theory presented here without approximation (Equation ( 59)) is computationally confirmed by comparison with experimental data from Jones and Beyer [1,35] for second-order sum-frequency pressure at range r R = 481.8mm due to a large rigid sphere of radius a = 3.18 mm in water insonified by two perpendicular incident beams at frequencies ω a /2π = 7 MHz and ω b /2π = 5 MHz, such that k a a = 89, k b a = 63.6 and k + a = 152.6 which are all large.Good quantitative agreement is found between theory and measurements with 0.98 correlation across scattered angle and overall mean square difference of less than 0.3 dB (Figure 2) as shown in the online supplementary material.[1,35] and exact theory obtained computationally for the sum frequency second order pressure.Good quantitative agreement is found between measurement and theory with 0.98 correlation across scattered angle and overall mean square difference of less than 0.3 dB.

Transitions Between Dominant Mechanisms in Sum or Difference Frequency Sensing of an Object
At sufficiently long range, the SS, IS and SI interactions always dominate, as seen from the asymptotic behavior given by Equations ( 90) and ( 92)-(94) compared to the inverse radial far field fall-off of S2 dictated by Equation (50).At such long ranges, only the object's response at the primary frequencies can be measured from the second order field, if the field is large enough to be measured.Our analysis here (Figures 3-8) shows that at shorter ranges, however, depending on frequencies, receiver range and the object's scattering properties, the S2 mechanism can also dominate, where the object's response at the sum or difference frequency may be deduced from the second order field at the sum or difference frequency.While it is possible to exploit these complicated transitions between dominant mechanisms for novel remote sensing applications, difficulties also arise in many cases, as we will show.
Examples are computed from the exact governing equations which allow general object shapes and sizes compared to the wavelength, but for simplicity scenarios involving spherical objects are presented.The geometry is shown in Figure E2, which does not include measurements in forward scatter where IS and SI are larger.These forward scattered cases are treated in online material.With the harmonic wave approximation, p II is given by Equation (82) with z s = 0 or (83), p SS , p IS and p SI are determined from Equations (78)-(80) respectively, and p S2 from spherical wave expansions (online supplementary material).The results represent difference frequency wave amplitudes in the time domain between t = r R /c 0 and t = r R /c 0 + T for narrow-band incident fields, as discussed in Appendix E. In all examples a practical object-centered coordinate system is used where k − z s = 0 for the II component in Equation (82).A normalization 2E 0 = |P a0 P b0 |/(ρ 0 c 2 0 ) is used in the examples.To generalize the analysis, dimensionless parameters k 0 a, k 0 r R , k − r R and k − a are used, where k 0 = (ω a + ω b )/(2c 0 ) is the center frequency wavenumber, k − = ω − /c 0 is the difference frequency wavenumber, a is the radius of the object and r R is the receiver range.The values for the parameters are chosen based on the following considerations: (1) k − r R 1, so that the II interaction can be time separated from the other mechanisms, (2) the difference frequency is lower than the primary frequencies, and (3) important transitions between mechanisms can be seen.Illustrative examples are presented and discussed for medical imaging, various ocean sensing applications as well as for sensing in air and the solid earth.Some general observations on the transitions between mechanisms found for rigid and pressure release cases are noted in Appendix F.   p SS dominates for small k − r R and p IS + p SI dominates for large k − r R .The nonlinear parameter is β = 3.6, and the density ratio between the object and medium is 7.6.Lower dimensionless abscissas are for general interpretations.Upper abscissas are for a medical or underwater sensing example as in Figure 3, except the rigid object is movable.Upper abscissas are for a medical sensing or underwater example of perpendicular primary waves with 1 MHz primary and 100 kHz difference frequencies, sphere radius of 0.1 mm at range 1 m from the receiver with sound speed 1500 m/s sound speed.In (a), p SS dominates for large k 0 a, while p IS + p SI and p S2 dominate for small k 0 a.In (b), p SS dominates for large k − a, p IS + p SI and p S2 dominate for small k − a.In (c), p SS is dominant and has the most gradual roll-off with range.In (d), p S2 and p IS + p SI dominate and SS is too small to be shown but appears in Figures S13 or S18.The nonlinear parameter is β = 3.6.
By matching the relevant dimensionless parameters and properly scaling the nonlinear parameter β, results presented in Section 5 can be applied to various sensing scenarios in air, water and solid earth.Depending on the actual dominant mechanism, the response of the object at the primary frequencies or difference frequency can then be deduced from P 2− .It is useful to clearly define the II component of S2, corresponding to the S2 II mechanism, via pressure release: rigid immovable:   The object radius a, receiver range r R and the center frequency are fixed so k 0 a = 2.09 and k 0 r R = 1047.The nonlinear parameter is β = 3.6 and the density ratio between the object and the medium is 2.4.

Medical and Fine-Scale Underwater Sensing
For nonlinear medical ultrasound sensing applications [28,36,37], which are portable to fine-scale high frequency underwater sensing, a transition between nonlinear mechanisms can be seen for rigid (Figures 3-5) and vacuous objects (Figure 6) sensed at the difference frequency in a watery environment.The SS mechanism is dominant (top right point of Figures 3b and 4b) for a 0.5 mm rigid object in a water environment sensed at a range of 1 m in backscatter for 1 MHz primary and 10 kHz difference frequency, and stands above noise with typical incident and received levels.As the object diameter is reduced IS and SI in backscatter become dominant for the immovable collinear case, so that no difference frequency information can be obtained due to lack of a measurable S2 component.When perpendicular waves are incident on the immovable object, the dominant mechanisms are similar to those of the collinear case except S2 is as dominant as the IS and SI as object diameter decreases (Figure 5b).Allowing the rigid object's centroid to move in response to the wave fields for the same parameters enables S2 to become dominant through D2 as the object radius is reduced by two orders of magnitude from 0.5 mm (Figure 4b).The R2 component, due to radiation from object centroidal motion is negligible and unmeasurable for this example when S2 dominates (Figure 4b).If the rigid object was attached to sufficiently underdamped springy material, it is possible that a resonance effect could make R2 and consequently S2 dominant and detectable at resonance even for large object radii, as shown in Figure 8.It may be difficult to extract the difference frequency response at resonance, however, if R2 is not dominated by the response to II as is the case for the example in Figure 8 as shown in Figure S15.
Continuing the medical ultrasound and high frequency underwater analysis with a vacuous bubble of 0.1 mm radius in water, first in backscatter, the S2 and SS contributions are equally dominant with the bubble sensed well above noise at 0.25 m range with a 2.0 MHz primary and 400 kHz difference frequency (Figure 6b, k − a = 0.16) for typical incident levels.This makes it challenging to extract difference frequency response information from the second order field, since this information is not contained in the contaminating SS component.All frequencies in this example are well above the bubble's resonace frequency (k res a ≈ 0.0137 near the water-atmosphere surface), so pressure release theory is accurate.For smaller bubble diameters (k − a < 0.1) the S2 mechanism becomes dominant but the S2 II contribution becomes increasingly negligible (Figure S14) so that a complicated inversion involving primary frequency scattering information is necessary to extract the object's difference frequency response information from S2 especially for k a a and k b a near unity.Dominance of S2 in backscatter occurs when k a a and k b a begin to decrease below unity and SS and SS-type effects at the object boundary begin to dominate S2 which then obeys Equation (97) for k a a and k b a much less than unity.As the difference frequency approaches the air bubble resonance frequency from above, the S2 field (Figure 7b) becomes significantly different from that of the pressure release case (Figure 6b), but the transition from SS to S2 dominance is seen to occur at frequencies above where resonance effects matter and where pressure release theory is sufficient.The S2 pressure peaks at resonance but then eventually decreases to zero as frequency decreases (Figure 7b).The minima in S2 seen in Figure 6a,b occur when SS effects at the boundary become less dominant as k a a and k b a approach unity from below.As the bubble radius increases above 0.1 mm, the SS mechanism transitions to dominance and masks difference frequency responses from the object (Figures 6b and 7b) .For the same 2 MHz primary and 400 kHz difference frequencies and 0.1 mm radius gas filled or pressure release bubble, but with reception in the forward direction at 1 m range, IS and SI are dominant and become increasingly so for larger bubble radii so masking difference frequency information from the object (Figure S19).(Typical primary wave incident levels in medical ultrasound are 240 dB re 1 µPa in water with noise levels in the noise levels in the 50 to 65 dB re 1 µPa range for the difference frequencies considered.) The examples are illustrated and evaluated both in terms of dimensionless parameters and specific frequencies and lengths.The former makes it possible to interpret the results universally over the broadest range of possible applications by scaling the object size, receiver range, and frequencies as necessary.In this way, for example, Figure 3b with 1 MHz primary and 10 kHz difference frequencies at range 1 m from the receiver can be interpreted for lower or higher frequency cases.Consider increasing the frequencies by a factor of fifteen, from 1 MHz to 15 MHz for the primary and 10 kHz to 150 kHz for the difference which are also used in ultrasound, and correspondingly decreasing the receiver range by a factor of fifteen.The nondimensional axes of the same figure show that the object radius where the SS mechanism begins to dominate must then correspondingly be reduced by a factor of fifteen.Similar analysis of scaling up or down frequencies and length scales can be performed readily for all figures with nondimensional axes.
Depending on one's perspective and the physical scenario, the difference frequency responses from an object sought in the medical vibro-acoustics (e.g., [4,5]) may be extractable from some or all S2 field components with varying amount of difficulty, and would certainly exist in an R2 subcomponent of S2 since this corresponds well with the mechanism discussed in Greenleaf et al. [4] and referred to as an "acoustic emission", but would not be available in SS, IS or SI fields.

Sensing Ocean Bubbles
The preceding example of a pressure release bubble in water applies directly to near-surface air bubbles in the ocean of the same size, primary and difference frequencies, and receiver ranges, and can be scaled to other similar cases using the non-dimensional parameters in (Figure 7) in correspondence with Equation (97).Other examples of the dominance of SS IS and SI at ranges of practical interest can be seen in Figure 9a where for a 1 mm bubble with primary frequency 430 kHz and difference frequency 43 kHz, SS will dominate the total second order difference frequency field at 4.88 m range from the object in backscatter, while IS and SI forward will dominate at much shorter ranges in forward scatter (Figure 9b).Scaling up, with a 10 mm bubble at 43 kHz primary and 4.3 kHz difference frequency SS will dominate at 48.8 m range in backscatter with IS and SI dominating at much shorter ranges in forward scatter.In the ocean context, these examples show that for primary and difference frequencies well above the bubble's resonance frequency, wave-wave interactions SS, IS and SI become important and will dominate the difference frequency field from an air bubble even at relatively close sensing ranges.In this regime where wave-wave interactions are important, the total second order field from an aggregate of bubbles is not the superposition of the individual second order bubble fields.In the analysis of air bubbles in water using a primary pump frequency at the air bubble resonance (with wavenumber k b ) and a variable primary image frequency (with wavenumber k a ) moving from to well above resonance (e.g., [2]), Figure 10 shows that wave-wave fields from a bubble are not significant at a receiver in backscatter for many practical sensing ranges and primary image frequencies.The same Figure 10 shows, however, that as k a a approaches unity the wave-wave IS and SI backscatter components can attain similar levels as S2 and become dominant.This transition occurs for smaller k a a as range increases.It also occurs for smaller k a a for receivers in the forward direction (Figure S21).Analytic conditions for when SS, and forward IS and SI are negligible compared to S2 are given in Appendix B. Due to the complicated nature of IS and SI in nonforward directions such conditions are available only through computations such as in Figures 10 and S21.In the other contexts, such as the previous medical and the following ocean fisheries examples, as well as cases where parametric arrays or related approaches are intended to employ II incidence, primary frequencies are typically not used at the frequencies of the desired response from an object otherwise linear acoustic sensing would be employed.For image wave k a a > 1, wave-wave effects dominate at shorter receiver ranges as k a a increases.

Ocean Fisheries Sensing
In ocean fisheries applications, difference frequency sensing has been used in ocean line-transect surveys to potentially probe the near-resonance frequency response of herring schools in the Nordic Seas via a parametric array sonar [7].The linear near resonance response of fish has relevance to linear ecosystem-scale sensing techniques [15,16,38].In the nonlinear work [7], the parametric sonar transmits high intensity primary frequency waves that lead to combined incident primary 15-21 kHz and difference frequency 0.5-6.0kHz fields at the herring shoals, roughly 50-100 m below the research vessel [7].Herring swimbladder equivalent radii are between roughly 0.8-1.0cm for herring in the 100-50 m depth range.Our analysis shows that the SS component, which has no difference frequency information and is essentially a nuisance contamination in this context, begins to make a non-negligible contribution to the total received second order field from any herring in the 50-100 m depth range for difference frequencies above 2 to 3 kHz (Figures 11 and 12).For difference frequencies below 2 kHz S2 is dominant, but not always equal to S2 II as see in Figures 11 and 12 which means that a complicated inversion involving knowledge of the herring scattering response at the primary frequencies and wave-wave boundary effects will be necessary to extract the linear difference frequency response of the herring for difference frequencies this range.This is because k a a and k b a are near unity so that a delicate balance between II, SS, IS and SI boundary values and second order boundary advection may affect the radiated S2 amplitude.It is noteworthy, however, that some of the necessary information for this inversion may be obtained from primary frequency scattered returns that could be received simultaneously from parametric array survey.Similar analysis for cod, (Figure 13) which typically shoal in the 100-150 m range during spawning in Lofoten Norway and have corresponding swimbladder equivalent radii of roughly 4.2 cm where they are neutrally buoyant, shows SS is either dominant or a significant contaminant over the entire difference frequency range from 0.5 to 6 kHz making it impossible or difficult to obtain difference frequency response information from the second order field in this scenario.The relative amplitudes of difference to primary frequency waves radiating from a parametric array may differ from those assumed here, depending on array geometry and other factors [1] which may alter some results.

Sensing Objects From IS Interactions
The IS mechanism can lead to robust sum or difference frequency detection and sensing of objects in air or water over extended ranges since asymptotically IS is constant in range.Perfectly reflecting objects, for example, with circumferences that are on the order of the primary wavelength or larger (|S a | ≥ 1) and incident sound pressure levels on the order of that for a 1 watt monopolar point source at 1 m ( 171 dB re 1 µPa at 1 m in water and 108 dB re 20 µPa at 1 m in air) are seen to have IS levels (Figure 14) that stand above typical ambient noise ambient noise levels (Figure S8).A trip wire employing perpendicular primary beams could be used, for example, so that when an object crosses the beam of a primary field, a sum or difference frequency IS would trigger its presence at a sensor in the forward direction of the primary beam crossed by the object.The advantage is the sensor could be well past the primary shadow of the object where detection at primary frequencies would be difficult since the received primary field would be dominated by the incident primary field.In another example, it is seen from Equations ( 92)-(94) that stealthy long range sensing of an object's primary frequency response via sum and difference frequency measurements is possible.For example, if an object scatters or radiates an outward propagating wave P Sb at angular frequency ω b , a beam P Ia could be transmitted at a different frequency ω a forward scattering past the object to a receiver.By continuing this for angles about the object, it is possible to reconstruct the directivity pattern of P Sb from P 2± , for primary frequency object classification at the difference frequency.

Sensing Objects in Air
In air, the SS mechanism can lead to robust sum or difference frequency detections from objects large compared to the primary wavelength.Consider, for example, nonlinear acoustic sensing of a 4 cm diameter rigid sphere with a primary frequency of 40 kHz and a difference frequency of 4 kHz with bat sonar incident levels of 140 dB re 20 µPa [39] which has relevance to robot sonar.Received levels of at least 53 dB re 20 µPa are expected at 2 m range at the difference frequency and at nearly an order of magnitude greater range at the sum frequency, which are well above noise levels at those frequencies (Figure S8) using Equation ( 90) or (91).As the radius of object increases, so will the detection range.

Sensing in the Solid Earth
Difference frequency sensing of buried structures in the solid earth, such as archaeological remains, is also possible.A 4-m scale rigid structure, for example, can be robustly detected by the S2 mechanism for 50 Hz primary and 10 Hz difference frequency waves at distances greater than 1 km (Figure S7) for standard ground vibrators within tens of meters of the feature, e.g., a 240 dB re 1 µPa incident SPL can be generated from a 150 kN vibrator 30 m away with β = 1000, ρ 0 = 3000 kg/m 3 , c = 3000 m/s.The SPL in air, (a,c), is referenced to 20 µPa and the SPL in water, (b,d), is referenced to 1 µPa.In (a,c) for air, high frequency incident SPL between 120 dB and 140 dB re 20 µPa gives sum frequency SPL between 48 dB to 88 dB re 20 µPa, and difference frequency SPL between 38 dB and 77 dB re 20 µPa respectively, for ω b /ω a = 0.6.Similarly, in (b,d) for water, high frequency incident SPL between 170 dB and 190 dB re 1 µPa gives sum frequency SPL between 48 dB to 90 dB, and difference frequency SPL between 37 dB and 77 dB respectively, for ω b /ω a = 0.6.With moderate to high incident SPLs, these sum and difference frequency second order waves can be measurable at long range.

Discussion
Unique opportunities exist for sensing an object at the sum or difference frequency of primary incident waves.Many of these have been discussed in a variety of contexts, as for example through the use of parametric source arrays and related approaches in ocean fishery, gas bubble and submerged object sensing and medical ultrasonic imaging [1][2][3][4][5][6][7], but in the absence of a complete theory for the second order field in the presence of an object.Here we examine some practical sensing scenarios in water, air and solid earth.In the context of using parametric arrays, for example, to probe an object's sum or difference frequency response, a problem we find is that in many practical scenarios, SS, IS and SI incidence on the object and boundary advection responses from the object contribute more than II incidence to S2, which may not have been expected.This can make accurate extraction of the object's linear response to sum or difference frequency waves from S2 far more complicated than may have been anticipated.The other problem is that the field at the receiver may be dominated by SS, IS or SI, which also may not have been expected, and these contain no such sum or difference frequency linear response information.For parametric arrays and related approaches, these issues arise because the primary fields needed to form the sum or difference frequency incident field (II) typically interact with the object to be sensed in a manner that cannot be ignored when analyzing the second order field.
Inspired by 18th-century accounts of difference frequency waves propagating from musical instruments [1], Lamb developed solutions for the II interaction of collinear time-harmonic primary planewaves [12].Westervelt later solved the time-harmonic non-collinear planewave case [11] and estimated the second order field of an endfire array so inventing the parametric array [3] which is a wave-wave based instrument with mature and well tested theory [1].Then Dean and Baxter obtained time-harmonic solutions for the omnidirectional cylindrical [29] and spherical [29,33] wave cases, and Tjotta et al., investigated II interaction for time-harmonic incident beams [40][41][42].Jones and Beyer measured the second-order acoustic field at the sum frequency caused by an object in the overlap region of two primary incident waves [1,35] and proposed a heuristic formula consistent with an SS-type mechanism based on Dean's solution.Greenleaf et al., proposed measuring an object's difference frequency response from primary wave insonification for applications in medical ultrasound where significant interest in this exists [4,5].Without a complete theoretical formulation for the second order field in the presence of a movable object, it was noted that the interaction of primary scattered waves may dominate the second order field [36,43,44] and mask information about the object's difference frequency response contained in weaker field components.Such a complete theory, as presented here, is necessary to investigate the conditions under which an object's difference frequency response may be detectable and when each of the various second order mechanisms may dominate.To treat the case of a movable object, we define a complete and self-consistent second order acoustic wave-exciting force, and find previous helpful formulations [28,[45][46][47] missed the effects of some second order field components and overcounted others.(An analytic solution for the SS field arising from a rigid immovable sphere was dervied by Silva and Bandeira [48].)Approximate analytic solutions for the second order nonlinear wave fields found in the presence of low-impedance contrast inhomogeneities, that apply to one or many inhomogeneities, are also derived in Appendix A.
Since the SS, IS and SI interactions always dominate the field components that carry object information at sufficiently long range, exact and asymptotic analytic solutions are rigorously obtained for them, which put previous heuristic estimates [1,34,35] of SS into perspective.For spherically symmetric interacting waves, our analytic solution is exact and reduces to Baxter's [33] for this special case, where we show that Baxter's solution satisfies the Sommerfeld radiation condition but Dean's [29] does not.Analytic solutions are derived for the second-order field due to the interaction of a plane and spherical wave.These lead to analytic expressions for IS and SI in the far field of any arbitrarily shaped scatterer.
A number of investigations on the nonlinear acoustics of air bubbles in water are based on leading order approximations near and below resonance frequency where the bubble's circumference is very small compared to the primary field's wavelength [2,49,50].In these approaches SS, IS and SI fields are neglected, and are typically negligible away from the bubble for practical receiver ranges in this low frequency regime.Even in these scenarios, however, we find the leading order effect of SS, IS and SI on S2 through wave-wave boundary incidence can be significant and so are included here as well as the boundary advection effects included in the other approaches [49][50][51], which are of similar order.Our formulation then is consistent to leading order with these previous formulations in this low frequency regime near and below resonance.We find that for primary and difference frequencies well above bubble resonance, however, as the bubble circumference to primary frequency wavelength approaches and exceeds unity, SS, IS and SI wave-wave fields begin to significantly affect and eventually dominate the second order field propagating outward from an air bubble in water at ranges and frequencies important to many practical sensing scenarios and so cannot be neglected.In this case, linear superposition of individual bubble fields in an ensemble, which is valid and sometimes theoretically relied upon for primary frequencies near bubble resonance [49,51], is not valid in this higher frequency regime where wave-wave interaction effects are important.Other approaches [52] heuristically treated bubble and water or seabed mixtures as effective media in a second order wave equation formulation to investigate nonlinear acoustic scattering, but these cannot be applied to or built up from a single scatterer since they lack fundamental boundary conditions at an inhomogeneity and have not been clearly verified experimentally.
Sum and difference frequency methods for detecting and analyzing bubbles in water with one primary wave (the 'pump') near the bubble resonance, and primary wave (the 'image') at much higher frequency (e.g., [2]) also depend on wave-wave interactions away from the bubbles to be negligible.A cautionary note in an introductory paper on this topic "the transducers were also kept at a short distance (less than an inch) from the bubble stream so as to avoid the generation of sum or difference frequencies due to nonlinearity of the medium" [2] refers to possible contamination by wave-wave interactions citing Westervelt [11] and Beyer [1].Here we show that for typical scenarios in this pump-image approach, while SS, IS and SI received fields are significant as incident fields at the bubble boundary, they are usually not significant away from the object in practical implementations of the pump-image approach due to the fact that the pump wavelength is always roughly two orders of magnitude larger than the bubble circumference to pump near resonance, the image wavelength is typically less than the bubble circumference, and sensing ranges are not large.The use of such a long primary pump wavelength, however, may be impractical for other applications involving parametric arrays or related approaches [1,[3][4][5]7] where much of the merit rests on the avoidance of traditional primary sources that transmit at such long wavelengths.

Conclusions
For applications in nonlinear acoustic sensing of objects by sum or difference frequency fields, a general second order theory of nonlinear interaction and scattering of acoustic waves in the presence of an object is derived and confirmed by comparison with experimental measurements.Practical guidance on designing experiments and interpreting measurements to sense objects with second order nonlinear acoustics have been provided, including approaches to determine whether wave-wave interaction effects, second order boundary scattering and advection, object resonance, or radiation from wave-induced object centroidal motion are dominant in a given scenario.We show that the wave-wave fields propagating from the object contain only primary wave frequency information about the object, while the remaining components contain sum and difference frequency information, with the latter typically being more desirable in many applications of nonlinear acoustic sensing.In the latter case, we show that a complicated inversion may be necessary to extract this sum or difference frequency information.
In the use of parametric arrays, for example, we find it is often not possible to simply assume the only incident field on the object is the difference frequency field generated by the parametric array.This is because intense primary frequency waves are needed to generate this incident difference frequency field.We show that in many practical examples of sensing objects these primary fields in fact lead to the dominant difference frequency response from the object by mechanisms other than linear scattering of the incident difference frequency field of the parametric array.
In general, we analytically find that at long ranges from the object, wave-wave interactions always dominate the second order field, but the required ranges in some cases may be too large to enable practical measurements.We find at shorter ranges, any of a number of second order nonlinear mechanisms may dominate, depending on frequencies, object size, object composition, and source and receiver positions with respect to the object.We show that in diverse applications ranging from ocean fishery and bubble to medical sensing, careful analysis using complete and self-consistent nonlinear theory is necessary to determine which mechanism or mechanisms are dominant in a given scenario.
Since finite-duration waveforms are typically used in practical sensing scenarios, we analytically show that there is no scattering of sound by sound outside a region of interacting primary plane waves of compact support.This makes it possible to readily distinguish incident-incident wave interactions from those involving an object by appropriate choice of measurement geometry.
Similarly, eliminating the second order velocity v 2 from Equations (A5) and (A6) yields For the equation of state, it is assumed that the density is a function of pressure only and can be expressed as a Taylor series expansion about the ambient pressure, Taylor series expansion leads to which is the classical wave equation for space containing inhomogeneities [24].Substituting Equation (A11) into Equation (A8) to eliminate δ 2 and further assuming Γ e is stationary in time yields By adding ρ −1 0 ∇ 2 p 2 − κ 0 ∂ 2 p 2 /∂t 2 on both sides and then multiplying by ρ 0 , Equation (A13) can be rewritten as where γ κ = (κ e − κ 0 )/κ 0 and γ ρ = (ρ e0 − ρ 0 )/ρ e0 are the fractional changes in compressibility and density of the medium due to inhomogeneities, κ 0 = ρ −1 0 c −2 0 is the mean compressibility and c 0 is the mean sound speed of the medium.For a homogeneous medium, 2A), and the above equation reduces to Equation (17).Consider time harmonic fields incident on a single target.Let the incident wave consist of two harmonic components, P a and P b with frequencies ω a and ω b .The corresponding Helmholtz equation for the second order pressure at the sum or difference frequency is where The solution to Equation (A15) is which has five components corresponding to the five parenthetical groups in curly brackets in Q inhomo 2± of Equation (A16).Components in the first three sets of curly brackets in Equation (A16) are nonzero in or on the inhomogeneity and describe second order scattering S2.As in the linear case, the changes of compressibility in the medium, k 2 ± γ κ P 2± , scatter as monopoles while the changes of density, ∇ • (γ ρ ∇P 2± ), scatter as dipoles.Components in the last two sets of curly brackets in equation (A16) are nonzero where the primary waves are nonzero and describe nonlinear wave-wave interactions, including II, SS and IS, where only II exists if the medium has no inhomogeneities.

Appendix B. Second Order Pressure Field from a Small Resonant Air Bubble
Consider a small spherical gas bubble in an infinite fluid subject to the action of two collinear harmonic plane incident waves.The incident waves have frequencies ω a and ω b and amplitudes P Sa and P Sb .The bubble has a radius of a at rest.We derive the asymptotic solution for k a a, k b a, k ± a 1.By including the bubble's damping and restoring effects into the equation of motion for a pressure release sphere, we obtain an approximate equation governing the motion of the bubble wall in the radial direction, which takes the first order form ξ1 where the damping constant δ = 4µ/ρ 0 ω 0 a 2 includes dissipation from thermal phenomena, viscosity, acoustic radiation, and surface tension, and ω 0 is the natural frequency of the bubble volume pulsation, which is determined from where σ is surface tension, p 0 is the ambient pressure, and γ is the polytropic constant [2].
For the given incident waves, the first order scattered field from a small pressure release sphere p S1 can be determined from the boundary condition, Equation (10).The amplitude of first-order motion of the bubble wall is obtained from Equation (B18): where the amplification factors α(ω a ) and α(ω a ) are defined by The first order scattered field will be monopolar as given by P Sa (r) = −P a0 α(ω a ) a r e ik a r , P Sb (r) = −P b0 α(ω b ) a r e ik b r .(B22) The results in Equations (B20) and (B22) approach the solution of a pressure release sphere at large frequencies ω a /ω 0 , ω b /ω 0 1.The equation governing the second-order motion of the bubble wall is obtained from Equation ( 19 where f NL denotes the second-order interaction terms associated with surface tension, thermal phenomena and viscous effects, given by For small bubbles with k a a, k b a 1, the second-order incident wave P I2 (Equation ( 44)), which contains II, SS, IS and SI components, is dominated by P SS on the bubble wall In the vicinity of the bubble, v 1 is dominated by the velocity components due to first order scattering, V Sa and V Sb , which are determined by use of Equation (B22) and near the bubble surface take the form The second-order S2 pressure component from a pressure release sphere p S2 can be obtained by solving the linear scattering problem with the boundary condition given by Equation (51).With ξ 2 determined from Equation (B23), the S2 pressure is where the amplification factor A ± res is and the coefficient B ± res is defined by B ± res represents the effect of the second-order interactions associated with surface tension and fluid viscosity, which is found to be insignificant compared to the other effects.Equation (B27) becomes the second-order solution for a pressure release sphere when A ± res =1 and B ± res = 0.If the second-order incident wave P I2 , which contains all II, SS, IS and SI wave-wave fields at the boundary, is omitted in (B23), the solution in Equation (B27) becomes equal to that obtained by Newhouse & Shankar [2].We find, however, in our examples that P I2 is of similar order as the other terms in (B23) and should not be ommitted.
In the far field k ± r 1, P SS can be determined from Equation (91), and P IS and P SI in the forward scatter direction can be determined from Equations (92)-(94).Their ratios to P S2 are found to be ) , (B30) α(ω a )α(ω + ) , (B31) α(ω a )α(ω − ) , (B32) which, when the ratios are small, become conditions on receiver range, probing frequencies, and parameters of the object and medium for which these wave-wave fields at the receiver will be negligible, and when the ratios become large become conditions on when these wave-wave fields at the receiver are dominant.It is noteworthy that no simple expression has been obtained for IS + SI in nonforward directions, so computations of this have been provided in Section 5.
On S, the first and second order waves for a general locally reacting boundary problem are The coefficients α 1 and α 2 are functions of position on S, and they can be determined from the first or second order boundary conditions (C39) and (C40), as ) (C46) The first order scattered wave p S1 from a general locally reacting boundary can be obtained by evaluating Equation ( 8), with Equations (C43) and (C45).The second order scattered wave p S2 from a general locally reacting boundary can be obtained by evaluating Equation (25) with Equations (C44) and (C46).The II, SS, IS and SI overlap regions of compact support are shown in Figure E2.If a receiver is placed in any location where the II overlap region (gray in Figure E2) and the SS overlap region (blue in Figure E2) do not intersect, as in the example shown in Figure E2, it will measure p SS + p IS + p SI + p S2 between t = r R /c 0 and t = r R /c 0 + T with no p II component, where t = 0 occurs when the fronts of the incident waves simultaneously arrive at the object center.The II overlap region entirely passes the receiver before the SS overlap region arrives for the examples shown in Figure E2 if r R > c 0 T/2 for the collinear case, and if r R > c 0 T √ 2/( √ 2 + 1) for the perpendicular case.If a receiver is placed in the forward direction in the collinear case, p II cannot be separated from the other components regardless of the duration of the incident waves, and will mask information about the object carried in p SS , p IS , p SI and p S2 .
Since second order scattering (S2) occurs at the object between t = 0 and t = T when all narrow-band primary waves of compact support overlap, the p S2 component has compact support via p S2 (r, t) ≈ {P S2± e −iω ± t w 2 (t − r/c 0 )} (E53) where P S2± is the time harmonic complex-amplitude of the second order scattered wave of Equation ( 50) and w 2 has compact support, and following w 1 , has duration T with unit height within the window except at the ends where smooth transitions to zero occur over periods much less than T.

Appendix F. Transitions Between Mechanisms for Rigid and Pressure Release Objects
Some general patterns of behavior for various mechanisms and transitions between dominant mechanisms are observed at the difference frequency when the receiver range is not asymptotically large, and the receivers are not in forward scatter where IS and SI are more significant: For rigid objects, it is found that: (1) p SS is dominant for high primary frequencies and large objects, as seen in Figures 3a,b and 4a,b for collinear incidence, and in Figures 5a,b for perpendicular incidence; (2) p SS is most sensitive to the changes of primary frequency and object size, because it is a function of the product of primary scattered fields P Sa P Sb which depend strongly on k a a and k b a; (3) when the primary frequency or object size decreases, the dominant mechanism transitions from SS to D2, IS and SI.Specifically, for movable objects, p D2 due to the wave-boundary interaction can become dominant as seen in Figure 4a,b.For immovable objects, p D2 and p IS + p SI are comparable, as seen in Figures 3a,b and 5a,b; (4) p IS + p SI is dominant at large difference frequency, as seen in Figures 3d, 4d and 5d; and (5) p R2 , which arises from the second order centroidal motion, has very little contribution to the total second order field, as seen in Figure 4.For rigid objects attached to springy material, however, the centroidal motion can be amplified such that p R2 becomes dominant when the difference frequency is close to a resonance frequency, as shown in Figure 8.
For pressure release objects, it is found that (1) p S2 is dominant at small primary frequencies, for small objects or at small difference frequency, as seen in Figure 6a,b,d, and is mainly due to SS and SS-type boundary effects from the first term of p 2 (Equation (61)) and the wave-boundary interaction (Equation ( 51)); (2) p SS is sensitive to changes in object size but not to the primary frequency for small objects, as seen in Figure 6a,b, because P Sa and P Sb do not depend on frequency for small pressure release objects; and (3) p S2 is insensitive to changes in difference frequency when it is small, as seen in Figure 6d, because p 2 and the wave-boundary interaction depend on the primary frequencies, which converge to the center frequency as the difference frequency decreases.
For all rigid, pressure release and rigid-mass-spring-resonating cases, p SS falls off in range more gradually than p IS + p SI in the non-forward direction and p S2 do, as seen in Figures 3c, 4c, 5c and 6c.This is consistent with the asymptotic behavior of each component discussed in Section 3.2.

Figure 1 .
Figure 1.Two primary incident fields p Ia and p Ib at respective angular frequencies ω a and ω b are incident on an object, causing two primary scattered fields, p Sa and p Sb .A receiver at r R measures the second order field p 2 (r R , t).

Figure 2 .
Figure 2. Comparison between Jones and Beyer's experiment[1,35] and exact theory obtained computationally for the sum frequency second order pressure.Good quantitative agreement is found between measurement and theory with 0.98 correlation across scattered angle and overall mean square difference of less than 0.3 dB.

Figure 3 .
Figure 3.A rigid immovable sphere with collinear incident waves propagating at îz .The difference frequency component of p 2 is calculated at r R = −r R îz as a function of (a) center frequency or k 0 a, (b) object radius or k − a, (c) receiver range or k 0 r R and (d) difference frequency or k − r R .The lower dimensionless abscissas are for general interpretations.The upper abscissas are for a medical or underwater sensing example with 1 MHz primary and 10 kHz difference frequencies, 0.1 mm sphere radius at range 1 m from the receiver with 1500 m/s sound speed.In (a), p SS dominates for large k 0 a, while p IS + p SI dominates for small k 0 a.In (b), p SS dominates for large k − a, p IS + p SI dominates for small k − a.In (c), p SS is dominant and has the most gradual roll-off with range.In (d), p SS dominates for small k − r R while p IS + p SI dominates for large k − r R .The nonlinear parameter is β = 3.6.The normalization employs 2E 0 = |P a0 P b0 |/(ρ 0 c 2 0 ).

Figure 4 .
Figure 4.A rigid movable sphere with collinear incident waves propagating at îz .The difference frequency component of p 2 is calculated at r R = −r R îz as a function of (a) center frequency or k 0 a, (b) object radius or k − a, (c) receiver range or k 0 r R and (d) difference frequency or k − r R .In (a), p SS dominates for large k 0 a, while p D2 dominates for small k 0 a.In (b), p SS dominates for large k − a and p D2 dominates for small k − a.In (c), p SS is dominant has the most gradual roll-off with range.In (d), p SS dominates for small k − r R and p IS + p SI dominates for large k − r R .The nonlinear parameter is

Figure 5 .
Figure 5.A rigid immovable sphere with perpendicular incident waves propagating at îx and îz direction, respectively.The difference frequency component of p 2 is calculated at r R = −r R ( îx + îz )/ √ 2 as a function of (a) center frequency or k 0 a, (b) object radius or k − a, (c) receiver range or k 0 r R , and (d) difference frequency or k − r R .Lower dimensionless abscissas are for general interpretations.Upper abscissas are for a medical sensing or underwater example of perpendicular primary waves with 1 MHz primary and 100 kHz difference frequencies, sphere radius of 0.1 mm at range 1 m from the receiver with sound speed 1500 m/s sound speed.In (a), p SS dominates for large k 0 a, while p IS + p SI and p S2 dominate for small k 0 a.In (b), p SS dominates for large k − a, p IS + p SI and p S2 dominate for small k − a.In (c), p SS is dominant and has the most gradual roll-off with range.In (d), p S2 and p IS + p SI dominate and SS is too small to be shown but appears in Figures S13 or S18.The nonlinear parameter is β = 3.6.

Figure 6 .
Figure 6.A pressure release sphere with collinear incident waves propagating in îz .The difference frequency component of p 2 is calculated at r R = −r R îz as a function of (a) center frequency or k 0 a, (b) object radius or k − a, (c) receiver range or k 0 r R , and (d) difference frequency or k − r R .Lower dimensionless abscissas are for general interpretations.Upper black abscissas are for a medical or underwater sensing example of 1 MHz primary and 100 kHz difference frequencies with 0.1 mm sphere radius at 1 m range from the receiver with 1500 m/s sound speed.The upper blue abscissa in (b) is for different medical or underwater sensing example of 2 MHz primary and 400 kHz difference frequencies with 0.1 mm sphere radius at 0.25 m range with 1500 m/s sound speed.In (a), p S2 is dominant.In (b), p S2 dominates for small k 0 a, p SS becomes comparable to p S2 for k 0 a near 1 then transitions to dominance for k 0 a > 1.In (c), p S2 is dominant, but p SS will dominate as range increases.In (d), p S2 dominates for small k − r R , and p SS becomes comparable to p S2 for large k − r R .The nonlinear parameter is β = 3.6.

Figure 7 .
Figure 7. Same as Figure 6 but for an air bubble in water using Equation (97) for p S2 at low k − a.

Figure 8 .
Figure 8.A resonating rigid movable sphere with collinear incident waves propagating in îz direction.P 2− is calculated in r R = −r R îz .In (a), total difference frequency component of p 2 is plotted as a function of difference frequency near resonance for damping ratio ζ = 0.001, 0.01, 1.In (b), p 2 , p SS , p IS + p SI , p D2 and p R2 are shown near resonance for ζ = 0.001.When the object's resonance is excited at the difference frequency, p R2 is amplified and it becomes dominant for small enough damping.The object radius a, receiver range r R and the center frequency are fixed so k 0 a = 2.09 and k 0 r R = 1047.The nonlinear parameter is β = 3.6 and the density ratio between the object and the medium is 2.4.

Figure 9 .
Figure 9. (a) A near surface air bubble in water with collinear incident waves propagating at îz .The difference frequency component of p 2 is calculated in backscatter at r R = −r R îz as a function of difference frequency or k − r R .Lower dimensionless abscissas are for general interpretations.Upper abscissas are for 430 kHz primary and 43 kHz difference frequencies with bubble radius 0.1 mm radius at 4.88 m range from the receiver in water.A transition between mechanisms is seen where p SS is dominant for k − a > 1 while p S2 is dominant for k − a 1 above and near resonance.The nonlinear parameter is β = 3.6.(b) Same as (a) but with p IS + p SI in the forward direction also included for r R = r R îz which is seen to dominate the field over the entire range shown.

Figure 10 .
Figure 10.An acoustically compact near surface air bubble in water will resonate at roughly k b a = 0.014 (Equation (B19)).Primary pump wave frequency is then set to roughly the bubble resonance frequency which is determined by k b a = 0.014.For collinear incident waves propagating at îz , the difference frequency component of p 2 is calculated in backscatter at r R = −r R îz as a function of difference frequency through k − r R .This nondimensional example may be interpreted for a 1 mm radius underwater bubble with 3.34 kHz primary pump wave frequency and image wave frequency varying through k a a = k b a + k − r R (k b a)/(k b r R ) from the vicinity of 0.0167 to 1.0137 for receiver ranges (a) r R = 0.1 m, (b) r R = 1.0 m, (c) r R = 10 m, and (d) r R = 100 m as in the upper black abscissas.The upper blue abscissas are for the corresponding primary image wave frequency for this example calculated by f a = (c/2πa) (k b a + (k − rR)(a/r R )) with sound speed c = 1500 m/s.For image wave k a a well below unity, S2 dominates and a localized peak is seen at difference frequency resonance.For image wave k a a > 1, wave-wave effects dominate at shorter receiver ranges as k a a increases.

Figure 11 .
Figure 11.Equivalent to a herring at 50 m depth modeled as swimbladder with equivalent radius of roughly 1 cm with collinear incident waves of 21 kHz primary and varying difference frequencies propagating at îz using resonant bubble Equation (97) with herring damping.The difference frequency component of p 2 is calculated in backscatter at r R = −r R îz as a function of difference frequency through k − r R or k − a = k − r R k a a/k a r R .Upper abscissa shows the difference frequency for the current example.

Figure 12 .
Figure 12.Same as Figure 11 but for herring at 100 m depth with equivalent bubble radius of 0.80 cm.

Figure 13 .
Figure 13.Equivalent to a cod at 100 m depth modeled as swimbladder with equivalent radius of roughly 4.2 cm with collinear incident waves of 21 kHz primary and varying difference frequencies propagating at îz using resonant bubble equation Equation (97) with cod damping.The difference frequency component of p 2 is calculated in backscatter at r R = −r R îz as a function of difference frequency through k − r R or k − a = k − r R k a a/k a r R .Upper abscissa shows the difference frequency for the current example.

Figure 14 .
Figure 14.Sound Pressure Level (SPL) of P IS± in the forward direction as a function of incident SPL at large range in air and water for ω b /ω a = 0.6, 0.9 and 0.99 assuming P a0 = P b0 and |S a | = |S b | = 1.The SPL in air, (a,c), is referenced to 20 µPa and the SPL in water, (b,d), is referenced to 1 µPa.In (a,c) for air, high frequency incident SPL between 120 dB and 140 dB re 20 µPa gives sum frequency SPL between 48 dB to 88 dB re 20 µPa, and difference frequency SPL between 38 dB and 77 dB re 20 µPa respectively, for ω b /ω a = 0.6.Similarly, in (b,d) for water, high frequency incident SPL between 170 dB and 190 dB re 1 µPa gives sum frequency SPL between 48 dB to 90 dB, and difference frequency SPL between 37 dB and 77 dB respectively, for ω b /ω a = 0.6.With moderate to high incident SPLs, these sum and difference frequency second order waves can be measurable at long range.

Figure E1 .
Figure E1.Second order wave due to the interaction of (a) incident waves of compact support (II), (b) scattered waves (SS) of compact support, (c) incident and scattered waves (IS) of compact support in the forward (c), and backscatter (d) direction by direct numerical integration of the indicated full time domain solutions (solid lines).The second order waves have compact support within which the indicated harmonic approximations are shown to be valid.The top panels show the real pressure and the bottom panels show the pressure amplitude at the difference frequency.Here ω a /(2π) = 500 kHz, ω b /(2π) = 300 kHz, T = 50 µs, c 0 = 1500 m/s and β = 3.6.The pressure is normalized by 2E 0 , where E 0 = |P a0 P b0 |/(2ρ 0 c 2 0 ).In (a), the pressure is calculated at the origin; in (b-d), the pressure is calculated at r R = 0.1 m from the center of a pressure release sphere of radius a = 10 µm and the scattered waves are approximately spherically symmetric with scatter functions S a,b = −k a,b a.

Figure E2 .
Figure E2.Primary incident and scattered waves of compact support and regions of compact support overlap, at different time instances (a-e).The left column is for the collinear incidence case, and the right column is for the perpendicular incidence case.The overlap regions for the incident waves, scattered waves, and incident and scattered wave are volumes in 3-D, and they are projected and shown on a 2-D plane as the gray, blue and green areas, respectively.
, p II,bb and p II,bb * , SS components p SS,aa , p SS,aa * , p SS,bb and p SS,bb * , and IS and SI components p IS,aa , p IS,aa * , p SI,bb and p SI,bb * correspond to self-interactions of each incident wave.The II components p II,ab and p II,ab * , SS components p SS,ab and p SS,ab * , and IS and SI components p IS,ab , p IS,ab * , p SI,ab and p SI,ab * are due to cross-interactions between the p Ia and p Ib incident fields. * = ρ −1 e0 ∂ρ e /∂p| p 0 = ρ −1 e0 c −2 eand Γ e = ∂ 2 ρ e /∂p 2 A10) into Equation (A7) to eliminate δ 1 and assuming the inhomogeneities ρ e0 and κ e are stationary in time yields