Acoustic waves: should they be propagated forward in time, or forward in space?

The evolution of acoustic waves can be evaluated in two ways: either as a temporal, or a spatial propagation. Propagating in space provides the considerable advantage of being able to handle dispersion and propagation across interfaces with remarkable efficiency; but propagating in time is more physical and gives correctly behaved reflections and scattering without effort. Which should be chosen in a given situation, and what compromises might have to be made? Here the natural behaviors of each choice of propagation are compared and contrasted for an ordinary second order wave equation, the time-dependent diffusion wave equation, an elastic rod wave equation, and the Stokes'/ van Wijngaarden's equations, each case illuminating a characteristic feature of the technique. Either choice of propagation axis enables a partitioning the wave equation that gives rise to a directional factorization based on a natural"reference"dispersion relation. The resulting exact coupled bidirectional equations then reduce to a single unidirectional first-order wave equation using a simple"slow evolution"assumption that minimizes effect of subsequent approximations, while allowing a direct term-to-term comparison between exact and approximate theories.

The evolution of acoustic waves can be evaluated in two ways: either as a temporal, or a spatial propagation. Propagating in space provides the considerable advantage of being able to handle dispersion and propagation across interfaces with remarkable efficiency; but propagating in time is more physical and gives correctly behaved reflections and scattering without effort. Which should be chosen in a given situation, and what compromises might have to be made? Here the natural behaviors of each choice of propagation are compared and contrasted for an ordinary second order wave equation, the time-dependent diffusion wave equation, an elastic rod wave equation, and the Stokes'/ van Wijngaarden's equations, each case illuminating a characteristic feature of the technique. Either choice of propagation axis enables a partitioning the wave equation that gives rise to a directional factorization based on a natural "reference" dispersion relation. The resulting exact coupled bidirectional equations then reduce to a single unidirectional first-order wave equation using a simple "slow evolution" assumption that minimizes effect of subsequent approximations, while allowing a direct term-to-term comparison between exact and approximate theories.

I. INTRODUCTION
An important category of acoustic wave models consists of those based on second order wave equations. Because they include second order derivatives in both time and space, they naturally lend themselves to rearrangements designed to focus either on the temporal or the spatial behaviour. In practice, most acoustic wave models have more than just this pair of second order terms, and these extra contributions may (or may not) be suited to a preference for either temporal or spatial analysis. But despite the wide range of different acoustic wave models, a subset of cases is sufficient to cover most important modifications. Here we choose the time-dependent diffusion equation (TDDE) [1], a model for waves in an elastic cylindrical rod [2], and acoustic waves as described by a generalized Stokes' equation [3] which allows for bubbles [4]. For these three situations, harmonic solutions are well known, but transient solutions are harder to find. Recently, the impulse-response and causal properties of some of these wave equations have been analyzed by Buckingham [5], and the (perhaps surprising) "noncausal" claims made therein are briefly addressed.
An advantage of second order wave equations is that once a propagation type is selected -whether into the future (i.e. along time), or along an appropriate spatial axis -they can be decomposed into pairs of explicitly directional and first order acoustic wave equations. The resulting coupled first order equations are not only often easier and faster to solve than the starting point of a second order wave equation, they are also well suited to the case of unidirectional traveling waves. However the unidirectional approximation is not demanded, and the formulation means that any additional approximations tend to be less restrictive [6]. For the purposes of this * Electronic address: Dr.Paul.Kinsler@physics.org article, however, these desirable features are not the core message. Instead, by comparing the results of factorizations aimed at generating temporally propagated wave equations, with those that generate spatially propagated ones, comparisons and contrasts can be made between the natural "reference" behaviors present in each decomposition. We will see that some terms in a wave equation will lend themselves to inclusion in one reference behaviour but not the other. These results then inform us as to how we might choose to trade off the practical efficiency of a spatially propagated picture against the physically accurate temporally propagated one.
Section II provides an overview of the factorization process for both spatial and temporal decompositions, using the time domain diffusion equation as an example. Factorization is then used to analyze the TDDE in section III, acoustic waves in an elastic cylindrical rod in section IV, and the Stokes' equation and its bubbly generalization in section V. This variety allows a discussion of both common and contrasting features of these directional decompositions. The conclusions are given in section VI.

II. DIRECTIONAL DECOMPOSITIONS
To compare the physical representations of waves in temporal or spatially propagated pictures, we use a factorization technique. This makes use of the concept of "reference" or "underlying" wave evolution [6,7]. If these reference behaviors are a close match to the actual wave evolution, subsequent approximations will be less stringent -in particular if we choose to make a unidirectional approximation. Although in some systems an exact match is possible, this is rarely the case if the waves depart from some idealized linear behaviour. Nevertheless, the power of these directional decompositions is in considering such variation from these exactly solvable cases, situations which often require approximation or numerical integration. To "factorize", start by selecting FBACOU Acoustic waves: propagated in time or space? Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ a propagation axis -either time, or along some spatial trajectory -after which the total wave can be decomposed into directional components that evolve forwards or backwards perpendicular to that propagation.
Throughout this article the context and/or arguments determine which domains ( r or k, t or ω) a given instance of some relevant function covers; in addition some symbols (e.g. Q, Ω, κ, c, c Ω , c κ , etc) are reused independently of each other (for different wave equation models) without distinguishing subscripts -this is in order to avoid cluttering the notation. Further, each wave equation given has, on the right hand side (RHS), a term Q which represents some general source term. Typically [5,8] Q is an impulse designed to elicit the primitive response of the system (i.e. Q = Qδ(t)δ(x)), but here it is allowed to be any kind of perturbation, (non)linear modification, or driving term we desire. For example, in the simple wave equations considered in this section, a density ρ( r) could be added, setting Q = (1/ρ)∇ρ · ∇g, and so match a wave equation used for ultrasound propagation [9,10]. Alternatively, adding a loss term to Q with the form η∂ t g would give us the time-dependent diffusion equation (TDDE) [1], which appears in a variety of contexts in physics, including acoustic waves in plasmas or the interstitial gas filling a porous, statistically isotropic, perfectly rigid solid [11]; it also models electromagnetic wave propagation through conductive media and is known as the telegrapher's equation.
In this section, an ordinary second order wave equation for some appropriate property g( r, t) will be used as a test bed on which to demonstrate the factorization process. It contains both a temporal response function p(t) and a spatial one s( r), both have a non-local character which is allowed for using a convolution ("⋆"), with a(u) ⋆ b(u) ≡ a(u ′ )b(u − u ′ )du ′ . Whilst some materials, in some regimes, do indeed have non-local spatial properties, the role of s( r) here is to ensure that the wave equation includes some nontrivial spatial properties. In many situations spatial structure would be introduced using p( r, t) instead of just p(t), while s( r) and its convolution would be absent; however here that would introduce complications beyond the scope of this example. The wave equation is where a non-interacting wave travels with speed c. Upon Fourier transforming into the k, ω domain, where d/dt ≡ ∂ t ↔ −ıω and ∇ ↔ +ı k, with k 2 = k · k, the result is where s( k) is the spatial Fourier transform of s( r), and p(ω) the temporal Fourier transform of p(t).
In both styles of derivation that follow, the steps taken are a good mathematical match to typical spatial propagation procedures [6,12]; but the physical meaning alters with the interchange of the roles of time and space. Factorization methods have a long history [13], but their adaption, application, and adoption to wave propagation of the type and context proposed here has been lacking until more recently [15? , 16].
Since Q can contain any sort of behaviour, eqn. (2.2) is already useful: e.g. factorization could easily be applied to the nonlinear propagation addressed by Pinton et al. [17]. In that case their initial eq. (3) can be matched to eqn. (2.1) above by restricting to the x-axis, using c 2 s(x) = µδ(x), p(t) = ρδ(t), and setting Q to match their RHS nonlinear term. Using directional decomposition, a first order wave equation simpler than e.g. eq.(10) in Pinton et al. can be obtained rapidly with fewer and less restrictive approximations. Of course, some acoustic wave equations reduced down to apply to a single wave property will not have the second order derivatives needed for this factorization scheme. However, such wave equations are already extensively approximated, so there may be scope for factorizing related equations which do contain them.

A. Temporal propagation, spatial decomposition
The most physically motivated factorization is to choose to propagate forward in time, and decompose the system behaviour (waves) into directional components that then evolve either forward or backward in space, as shown in Fig. 1. This is useful for analyzing situations where signals need to be separated from reflections, but requires an explicit modeling of the medium's timeresponse, perhaps involving convolutions or auxiliary differential equations.
To obtain such a temporally propagated, spatially decomposed representation, choose a suitable reference frequency Ω. This is allowed to have a wavevector (k) dependence, but should not have a frequency (ω) dependence. Rewriting eqn. (2.2) for temporal propagation, withΩ( k, ω) 2 = c 2 k 2 s( k)/p(ω). Unfortunately, the frequency dependence ofΩ means that it is not a good ref- Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ erence on which to decompose the propagation. Thus we must either ignore the material's frequency response p(ω), treat it as a correction by incorporating it into Q, or convert its frequency response into a wavevector response p(k) by approximation 1 .
Now decompose g into wave components g + and g − that evolve either forward or backward in space, with g = g + + g − . Then split eqn. (2.5) into two coupled first order wave equations for g ± ( k, t), i.e.
leading to either frequency or time domain forms This last equation tells us how g ± ( k, t) evolves [or if appropriately transformed, g ± ( r, t)] as we propagate in time. The original second order equation can be recovered [12] by substituting one of these into the other. Note that retaining a frequency dependence for Ω would have given rise to time derivatives on the RHS of eqn. (2.9). In addition to mathematical complications, these would also disrupt the otherwise straightforward integration of g ± in time.
If Q is small, i.e. |c 2 Q/p 0 | ≪ 2|Ω 2 g ± |, then the wave evolves slowly as it propagates. This "slow evolution" permits a temporal version of the unidirectional approximation [6,18] to be made 2 , setting g − = 0. The forward 1 For example, if restricted to waves near a frequency ω 1 in a narrow bandwidth ∆ω ≪ ω 1 , we can write ω ≃ c ′ k, where c ′ = ω 1 /k, and replace p(ω) with s(k) -i.e. convert a temporal dispersion into an approximate spatial dispersion. This is what was assumed in [31]. 2 However, when Q is a simple source term, it merely acts to drive both forward and backward waves equally. The two waves are then uncoupled, and no g + dynamics are affected by the neglect of g − , and vice versa. If Q has any dependence on g, e.g. if it were a loss or nonlinear term, then this would not be the case.
wave g + ( k, t) then follows Note that if Q is an impulse such as a delta function, the slow evolution condition will not hold at that point. If Ω 2 < 0, then g or g ± do not oscillate as time passes, but instead decay.

B. Spatial propagation, temporal decomposition
We could also choose to propagate forward in space, and then decompose the system behaviour (waves) into components that evolve either forward or backward in time, as shown in fig. 1. Although spatial propagation can seem non-intuitive, and requires careful handling of reflections or scattering, spatial propagation schemes are nevertheless popular, as in e.g. the KZK equation [19,20]. Their primary advantage is that the full (past and future) time history of a wave is available at whatever point in space the system has reached; in addition transitions across interfaces are also easy to handle. Thus medium parameters such as speed or wavevector can be defined as functions of frequency, so that material dispersion can be easily handled using pseudospectral methods [21,22], even enabling numerical investigations into exotic phenomena by means of artificial dispersions [23]. The angular spectrum approach [24,25] is another pseudospectral method, widely used in underwater and biomedical acoustics.
To obtain a spatially propagated, temporally decomposed representation, we choose x as the propagation axis, so that y, z are the transverse spatial; dimensions, with ∇ 2 T = ∂ 2 y + ∂ 2 z and transverse wavevectors k y and k z with k 2 T = k 2 y + k 2 z . Next, we need a suitable reference wavevector κ, which is allowed to have a frequency (ω) dependence, but should not have a wavevector (k) dependence. Rewriting eqn. (2.2) then gives In the converse of the temporal case, where a frequency dependence was inconvenient, hereκ has an inconvenient wavevector dependence; thus it is not a suitable basis for decomposing the wave evolution. Our choices therefore are to ignore the material's spatial response s( k), incorporate it into Q, or approximate it as a frequency response s(ω) 3 . If there is no spatial FBACOU Acoustic waves: propagated in time or space? Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ response, then s( r) ≡ s 0 δ( r) and s( k) = s 0 so that κ(ω) 2 = ω 2 p(ω)/s 0 c 2 . Eqn. (2.4) can be rearranged to give an expression for g( k, ω) directly, i.e.
Now decompose g into wave components g + and g − that evolve either forward or backward in time: with g = g + + g − . Then split eqn. (2.13) into two coupled first order wave equations, i.e.
leading to either wavevector or spatial domain forms where the different conventions for κ and ω give rise to the differing signs in the leading (reference) RHS terms if eqns. (2.9) and (2.17) are compared. This last equation tells us how to evolve g ± ( r, ω) [or if appropriately transformed, g ± ( r, t)] as we propagate along our chosen spatial axis x. Again, by substituting one of these into the other, the original second order equation can be recovered. Note that any wavevector dependence for κ would have given rise to spatial derivatives on the RHS of eqn. (2.17). In addition to mathematical complications, these would also disrupt the otherwise straightforward integration of g ± along x.
Here the directed g ± evolve according to κ, but that this reference evolution is modified by the additional source terms, either the general term Q, or the diffraction term dependent on ∇ 2 T (or k 2 T , depending on the chosen domain). These source terms not only modulate (e.g. drive, amplify, or attenuate) the wave equations equally, they also couple them together. If the source terms are small, i.e.
can be replaced by p(ω) -i.e. we have converted a spatial dispersion into an approximate temporal dispersion. This approximation allows the spatial dispersion to be added onto the material's temporal dispersion, which is particularly useful when propagating waves along a waveguide.
then the wave evolves slowly as it propagates. This "slow evolution" permits a unidirectional approximation [6,18] to be made, setting g − = 0. The forward waves g + ( r, ω) then follow Again, if Q is an impulse such as a delta function, the slow evolution condition will not hold at that point. If κ 2 < 0, then g or g ± do not oscillate in space, but instead are evanescent.

C. Discussion
There is an interesting tension between these two factorizations. At first sight, the most physical choice of propagation is that in time; however this means at any specific time, a full time-history is unavailable. This means that it is impossible (strictly speaking) to know spectral properties generally taken for granted, e.g. the wave speed c(ω) or the wavevector k(ω). To get frequency-dependent quantities we need either (a) a solution containing a complete time history, as might be obtained analytically for a some problems, or (b) to work in a spatially propagated picture, where a full time history is automatically available. However, the spatially propagated picture differs from our experience of a universe advancing in time.
Further, note that the two factorizations treat diffraction differently -in the temporally propagated case, having the full spatial profile to hand at each time step means that diffraction can be done exactly, even in a unidirectional model, whereas in the spatially propagated (and unidirectional) case a paraxial approximation needs to be applied. The first order wave equations derived above (and below) can be conveniently adapted using simple transformations or restrictions [6]. These are typically applied to unidirectional models, because whilst they make (e.g.) the representation of the forward wave better behaved, the backward wave becomes more problematic.

III. THE TIME-DEPENDENT DIFFUSION EQUATION
The time-dependent diffusion equation (TDDE) [1] is a second order wave equation with a loss term added; it appears in a variety of contexts in physics, including acoustic waves in plasmas or the interstitial gas filling a porous, statistically isotropic, perfectly rigid solid [11]. It has a three-dimensional, inhomogeneous form for the velocity potential g ≡ g( r, t) of [5] Here Q is a source term, such as a driving term or some modification to the wave equation. Next, η is a positive constant that imparts loss, and c is the high frequency speed of sound. In wavevectorfrequency space, with g ≡ g( k, ω) and A. Temporal propagation, spatial decomposition To decompose the TDDE into wave components evolving forwards or backwards in space first choose to propagate forwards in time while utilizing a suitable reference frequency Ω(k), with k = | k|. Eqn. (3.2) for g( k, ω) is now written where Ω(k) = ck, and its reference wave speed is c Ω = Ω/k = c. This is essentially the same as the wave equation eqn. (2.4), but with Q ′ = Q − ıηωg replacing Q. As already explained, k, t domain wave equations can now be given for velocity potentials g + and g − that evolve forward or backward in space, These directed g ± evolve according to Ω, but this reference evolution is modified by both the general term Q and the loss term η. These source terms not only modulate (e.g. drive, amplify, or attenuate) the wave equations equally, they also couple them together. If the source terms are small, i.e.
the wave evolves slowly as it propagates, so that a unidirectional approximation can be made, setting g − = 0; but remain alert to the inconvenient time derivative on the left hand side (LHS) of eqn. (3.5). The forward wave g + ( k, t) then follows Note that here the RHS side has a time derivative. Normally it is preferable for all such terms to be on the LHS, so that the RHS directly and unambiguously describes how the wave evolves. Fortunately, in this unidirectional wave equation, the two time derivative terms can be combined before re-applying the slow evolution approximations to get B. Spatial propagation, temporal decomposition To decompose the TDDE into wave components evolving forwards or backwards in time, first choose to propagate forwards along a spatial axis while utilizing a suitable reference wavevector κ(ω). To do this we need to select a primary propagation direction, here chosen to be along the x axis, so that k y and k z are the transverse wavevectors (with k 2 T = k 2 y + k 2 z ). This TDDE contains loss terms dependent on η, but it can be inadvisable to build loss into the reference behaviour of spatially propagated waves [26]. Therefore it is best to allow for the possibility that η might be retained on either the LHS (in which case η 1 = γ and η 2 = 0) or the RHS (in which case η 2 = η and η 1 = 0). Eqn. (3.2) for g( k, ω) is now written where κ(ω) is defined as so that the lossless (η 1 = 0) reference wave speed c κ = ω/κ is constant at c. If complex valued κ or c κ are acceptable, then (3.10) If a physical justification could be imagined, the spatial decomposition used here would allow the parameters c, η to have a dependence on ω, although the appropriately matching time dependence (i.e. convolutions over a temporal history) would need to be present in eqn. (3.1) -as in e.g. eqn. (2.1).
After defining Q ′ = Q + k 2 T g − ıη 2 ωg, follow the same steps as for eqns. (2.13) to (2.15), and decompose g into velocity potentials g + and g − that evolve forward or backward in time, with g = g + + g − . The two coupled first order wave equations for g ± ( k, t) are then In the x, ω domain, eqn. (3.11) can be rewritten where the different conventions for κ and ω give rise to differing signs in the leading RHS terms (e.g. compare eqns. (3.4) and (3.12)).
Here the directed g ± evolve according to κ, but that this reference evolution is modified by the additional source terms, either the general term Q, the diffraction term dependent on k 2 T , or the loss term dependent on η 2 . These source terms not only modulate (e.g. drive, FBACOU Acoustic waves: propagated in time or space? Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ amplify, or attenuate) the wave equations equally, they also couple them together. If the source terms are small, i.e.
|Q| , ηω g + + g − , k 2 T g + + g − ≪ 2 κ 2 g ± , (3.13) the wave changes slowly as it propagates in space, so that we can make a unidirectional approximation, setting g − = 0. Then the condition containing k 2 T is making a paraxial approximation, which is appropriate where the wave propagates primarily in a narrow beam oriented along some particular direction. The forward wave g + (x, k y , k z , ω) then follows (3.14) Here the choice of whether to put the loss (η) dependent term into the reference wavevector κ or not is not necessarily so important if the loss is small, since it will not break the unidirectional approximation. Further, just as for the time-propagated eqn. (3.6), again there is a factor of ω (or if transformed into the time domain, a time derivative ∂ t ) applied to the η term on the RHS. Now, however, because we are propagating forward in space, it can be easily calculated using the known ω (or t) dependence of g ± .

C. Discussion
Both decompositions of the TDDE have the same (lossless) reference wave speed, i.e. the high frequency speed of sound c = c Ω = c κ . However, they treat diffraction differently -in the temporally propagated case, having the full spatial profile to hand at each time step means that diffraction can be done exactly, even in a unidirectional model, whereas in the spatially propagated (and unidirectional) case a paraxial approximation needs to be applied.

IV. AN ELASTIC ROD WAVE EQUATION
Another acoustic system to consider is waves traveling along an infinite, isotropic and elastic cylindrical rod of radius R. Following Murnaghan's free energy model, Porubov has derived a wave equation governing propagation of the solitary waves along such a rod [2]. There is no impediment in this model against the rod having "auxetic" parameters [27], e.g. where the Poisson's ratio was negative [28]. This "elastic rod equation" (ERE) describes the displacement g ≡ g(x, t) with a second order wave equation of the form where g is the longitudinal displacement in the rod, and (following [27]) the coefficients read: Here β is the nonlinear coefficient, E and ν are Young's modulus and Poisson's ratio respectively, l, m, n specify Murnaghan's modulus, and ρ 0 denotes the density.
Poisson's ratio is typically rather small (i.e. |ν| < 1), in which case b 1 will be negative. In contrast b 2 can cover a wide range of values, especially if auxetic materials are considered, but usually ν, E > 0, so that b 2 < 0. In eqn. (4.1) Q is a source, driving, or other modification to the wave equation. In wavevector-frequency space, with g ≡ g(k, ω), eqn. (4.1) becomes

A. Temporal propagation, spatial decomposition
To decompose the ERE into wave components evolving forwards or backwards in space, choose to propagate forwards in time. This relies on a suitable reference frequency Ω(k), which results from a careful partitioning of the terms in the wave equation (4.4). The ERE for g(k, ω) is now written as where Ω(k) and new source term Q ′ k are (4.8) Now follow the same steps as for eqns. (2.5) to (2.7), and decompose g into velocity potentials g + and g − that evolve forward or backward in space, with g = g + + g − . The coupled wave equations for g ± (k, t) are then Eqn. (4.9) can be transformed so as to apply to g ± (k, t), becoming Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ Here the directed displacement waves g ± evolve as specified by Ω, but that this reference evolution is modified and coupled together by Q and the nonlinear term χV 2 . Note that V 2 needs to be re-expressed in terms of the sum of g ± , a procedure best done in the x, t domain. If the source terms are small, i.e.
|Q| , χk 2 V 2 ≪ 2 Ω 2 g ± 1 + b 1 k 2 , (4.12) then the wave evolves slowly as it propagates, so we can follow the prescription in [6] and make a unidirectional approximation, setting g − = 0. The forward waves g + (x, t) then follow where V 2+ is just V 2 calculated using only g + , i.e. with g − ≡ 0. This equation (and indeed the bidirectional eqns. (4.11)) have a RHS without any frequency (time) dependence, and so give the temporal propagation directly.

B. Spatial propagation, temporal decomposition
To decompose the ERE into wave components evolving forwards or backwards in time, choose to propagate forwards in space. This relies on a suitable reference wavevector κ(ω), which results from a careful partitioning of the terms in the wave equation (4.4). The ERE for g(k, ω) is now written as Here κ(ω) and the new source term Q ′ ω are defined as Now following the same steps as for eqns. (2.13) to (2.15), we decompose g into velocity potentials g + and g − that evolve forward or backward in space, with g = g + + g − . The coupled wave equations for g ± (k, t) are Eqn. (4.18) can be transformed so as to apply to g ± (x, ω), becoming Expanding Q ′ ω in eqn. (4.19), and taking care to transform its k dependence correctly, gives us Again, the directed g ± evolve according κ, with this reference evolution being modified and coupled together by Q, the nonlinear term χ, and the high-order dispersive term b 2 . If the source terms are small, i.e.
|Q| , χ∂ 2 then the wave evolves slowly as it propagates, so we can make a unidirectional approximation, setting g − = 0. The presence of the spatial derivatives in the latter two these conditions means that more care needs to be taken when judging whether or not they are satisfied. If the conditions can be relied upon to hold, the forward waves g + (x, ω) then follow This equation has a RHS containing spatial derivatives (or, back in the bidirectional eqns. (4.19), a wavevector dependence), so that the spatial propagation is not given directly. Consequently, for this model, the alternate choice of temporal propagation seems preferable.

C. Discussion
The nature of the ERE wave equation means that either choice of propagations, based either on Ω(k) or κ(ω), provide characteristically different evolutions, so that there is no perfect way of connecting the two pictures in all circumstances. In contrast, for the decompositions of the TDDE discussed in section III, the reference wave speeds were the same in the lossless case (where η = 0), and similar otherwise.
For a temporally propagated ERE wave, the reference wave speed (squared) depends on Ω(k), and is . If c 2 Ω < 0, as in the case for negative σ and large enough k, waves no longer propagate and instead decay. At low enough wavevector (i.e. where b 1 k 2 , σb 1 k 2 ≪ 1), the wave speed is simply c. In contrast, the high-wavevector limit of c 2 Ω tends to σc 2 , which, depending on ν may indicate either propagation or loss, i.e. a pass band or stop band.
The common case where b 1 and σ are both negative is shown on fig. 2, along with results for σ > 0. We we see that the denominator causes an asymptote at b 1 k 2 = −1, such that for σ < 0 a stop band extends from b 1 k 2 = −1 upwards; but for 0 < σ < 1 the stop band has a finite extent, and high-k propagating waves again exist. For σ > 1 the numerator switches sign first as k increases, and the stop band exist between σb 1 k 1 = −1 and b 1 k 2 = −1. The alternative case where b 1 is positive is shown on fig. 3, where the low-wavevector limit is simply c 2 . In the The temporally propagated ERE reference wave speed c 2 Ω as a function of wavevector k; where both b1, σ < 0, for selected σ. If c 2 Ω < 0, the waves decay (dashed lines) instead of propagating. The temporally propagated ERE reference wave speed c 2 Ω , in the b1 > 0 case, for selected σ. If c 2 Ω < 0, the waves do not propagate, but instead decay (dashed lines).
high-wavevector limit c 2 Ω tends to b 2 /b 1 , meaning that for σ < 0 (the usual case) only low-k waves propagate.
For the spatially propagated ERE wave, the reference wave speed has a much simpler behaviour; although this simplicity is counteracted by the corresponding wave equations (e.g. eqns. (4.19) or (4.22)) being less straightforward. The reference wave vector κ 2 gives a wave speed (4.24) In the low-frequency limit this is simply c 2 , but that the other limit is more complicated. For a high enough frequency (or large enough b 1 , so that b 1 ω 2 = c 2 ) the reference wave speed vanishes. For higher frequencies, the wave becomes evanescent, since κ 2 is negative and hence κ is imaginary. These are shown on fig. 4. Between the two choices of propagation direction (i.e. either of time or of space) we can see that the partitioning of terms leading to the reference frequency or wavevector are not the same. The b 2 term only appears in the reference behaviour for the choice of temporal propagation, but the b 1 term appears in both. Lastly, despite the differing sign of the b 1 term (or rather, because of it), to first order in k or ω, either reference behaviour gives the same dispersion.

V. THE STOKES' AND THE VAN WIJNGAARDEN'S EQUATIONS
The propagation of an acoustic velocity field is commonly described by the Stokes' equation [3], which can even be be extended to describe propagation of acoustic waves in an isothermal, viscous, bubbly liquid [29]. This extension was first done for one dimension by van Wijngaarden [4], hence the "van Wijngaarden's equation" (VWE), but has also been generalized to three dimensions [30]. However, the 3D VWE differs from the common form of the Stokes' equation in that it has spatial operators of ∇∇ · g and not ∇ 2 g. Here, to allow for a compact description that includes the 3D Stokes' equation, but remains compatible with the 1D VWE, a hybrid Stokes'/ VWE equation (S/VWE) for the velocity potential g ≡ g( r, t) is introduced. It is Here Q is a source term, such as a driving term or some modification to the wave equation. One could even reconcile this equation perfectly with the 3D VWE form [30] if the differences from the ∇ 2 → ∇∇· substitution were incorporated in Q, and any side effects properly considered. The equilibrium bubble radius scales as β, despite it having dimensions of time. With β = 0, this becomes the ordinary Stokes' equation, where the full 3D behaviour of the wave equation is valid. The parameter γ is proportional to the dynamic viscosity of the (bubbly) mixture. In wavevector-frequency space, with g ≡ g( k, ω) and k 2 = k · k = k 2 x + k 2 y + k 2 z , eqn. (5.1) becomes

A. Temporal propagation, spatial decomposition
To decompose the S/VWE into wave components evolving forwards or backwards in space, choose to propagating forwards in time. This relies on a suitable reference frequency Ω(k), which results from a careful partitioning of the terms in the wave equation (5.2) to give an Ω with only a k dependence. The S/VWE for g( k, ω) is now written as where Ω(k) and the new source term Q ′ are defined as Ω(k) 2 = c 2 k 2 1 + c 2 β 2 k 2 , (5.5) Here Ω tends to 1/β in the high-wavevector limit -or is constant at c 2 k 2 in the Stokes' case where β = 0. Note that as can be seen on fig. 5 the reference wave speed c 2 Ω = Ω 2 /k 2 in the same high-wavevector limit decreases towards zero as 1/kβ, and that the low-wavevector limit is simply c. Thus all components of the wavevector spectrum have an ordinary oscillatory (and non-lossy) reference evolution. Then, as for eqns. (2.5) to (2.7), decompose g into velocity potentials g + and g − that evolve forward or backward in space, with g = g + + g − . The two coupled first order wave equations for g ± ( k, t) are Eqn. (5.7) can be transformed so as to apply to g ± ( k, t), and becomes Here the directed g ± evolves as specified by Ω, but that that reference evolution is modified and coupled by Q and the dynamic viscosity γ. If the source terms are small, i.e. Ω k 2 Q , |γΩ∂ t (g + + g − )| ≪ 2 |Ωg ± | , (5.10) then the wave evolves slowly as it propagates, so a unidirectional approximation can be made, setting g − = 0. The forward waves g + ( k, t) then follow Notice that a k-dependence has appeared on term proportional to Q in eqn. (5.9) and (5.11), and also that its usual dependence on the reference frequency Ω has been canceled out. Although an inconvenient RHS time derivative appears, applied to the γ-dependent term, both time derivatives could be combined to give a directly solvable propagation, just as for the TDDE model in sec. III.

B. Spatial propagation, temporal decomposition
To decompose the S/VWE into wave components evolving forwards or backwards in time, choose to propagate forwards in space. This uses a suitable reference reference wavevector κ(ω), resulting from a careful partitioning of the terms in the wave equation (5.2) to give an κ with only a ω dependence. The S/VWE eqn. (5.2) for g( k, ω) is now written (5.13) and if the propagation axis lies along x, this becomes In both eqns. (5.13) and (5.14) the reference wavevector κ 0 (ω) is defined as The S/VWE eqn. (5.14) contains an imaginary part in its reference wavevector κ 0 , which is not always desirable. This can be circumvented by splitting κ 2 0 into real and imaginary parts using κ 2 0 = κ 2 R + ıκ I , and moving κ I to FBACOU Acoustic waves: propagated in time or space?
Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ the RHS 4 : The S/VWE can then be re-expressed as (5.18) although it will also be convenient to merge the κ I term into the Q term as follows If the propagation axis is along x, the RHS of the wave equation (still for g( k, ω)) can be written in either of two forms, namely Let us proceed without deciding yet whether κ 0 and Q (as in eqn. (5.14)) or κ R and Q R (as in eqn. (5.21)) is the best choice. To achieve this, a subscript-free κ and Q ′ are used to stand in for whichever pair of {κ 0 , Q} or {κ R , Q R } is convenient, noting that the RHS is given by either of If a physical justification could be imagined, this spatial decomposition would allow the parameters c, β, γ to have a dependence on ω, although the appropriately matching time dependence would need to be present in eqn. (5.1). Now following the same steps as for eqns. (2.13) to (2.15), decompose g into velocity potentials g + and g − that evolve forward or backward in time, with g = g + + g − . The coupled wave equations for g ± ( k, ω) are 4 Note the imaginary κ I part is not squared Eqn. (5.26) can be transformed so as to apply to g ± (x, k y , k z , ω), and becomes The two possible choices of κ are now addressed separately.

The bare κ0 form
On the basis of κ 2 0 , the reference wave speed in the lowfrequency limit is simply c, but the other limit is more complicated. For γ = 0, reference wave speed vanishes at ω 2 = 1/β 2 , and above this the wave becomes evanescent, since κ 2 0 is negative and κ 0 imaginary. Choosing κ = κ 0 and Q ′ = Q + k 2 T g, means that eqn. (5.27) for g ± (x, k y , k z , ω) becomes Again, the directed g ± evolve as specified by κ 0 , but that reference evolution is modified and coupled by Q and k 2 T . If the source terms are small, i.e.
a unidirectional approximation [6] can be made, with g − = 0; noting also that the k 2 T condition demands paraxial propagation. The forward waves g + (x, k y , k z , ω) then follow

The dressed κR form
On the basis of κ 2 R , the reference wave speed in the lowfrequency limit is simply c, but the other limit is more complicated. For zero dynamic viscosity (i.e. γ = 0), it matches the bare behaviour, although for finite γ the zero-speed transition is brought down down to a lower frequency. Again, above this vanishing reference wave speed the wave becomes evanescent, since κ 2 R is negative and so κ R imaginary. With κ = κ R and Q ′ = Q R + k 2 T g, eqn. (5.27) can be transformed to apply to FBACOU Acoustic waves: propagated in time or space?
Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ g ± (x, k y , k z , ω), becoming Here the directed g ± evolve as specified by κ, but that that reference evolution is modified and coupled by Q, γ, and k 2 T . If the source terms are small, i.e.
then the wave evolves slowly as it propagates. Note that the constraint on the Q term in eqn. (5.34) will fail at low frequencies; also that for non-zero γ, both the Q and the dynamic viscosity (eqn. (5.35)) conditions will fail near ω 2 β 2 ≃ 1. Eqn. (5.36) is equivalent to demanding paraxial propagation. If these conditions, and hence the unidirectional approximation hold, it is reasonable to set g − = 0. The forward waves g + (x, k y , k z , ω) then follow Both here in eqn. (5.37) and in the bidirectional eqn. (5.33) there is a non-zero γ-dependent source term.

C. Discussion
For this combined S/VWE model, just as for the ERE in section IV, the wave equation is such that each choice of propagation direction results in reference evolutions that differ both conceptually and in practice. This is evident from the reference wave speeds, which are wavevector limit where c Ω (k → ∞) → 0 and the transition from propagation to evanescence in c κ (ω) at β 2 ω 2 = 1. The behaviour of the frequency-dependent forms c κ are shown in fig. 6.
Buckingham [5] takes issue with the physical properties of analytic solutions of the VWE equation, notably its pressure Green's function; this has non-physical high frequency response when β = 0. Since the β 2 term in the VWE has a counterpart in the ERE's eqn. (4.1)) (i.e. the b 1 term), the same remarks could also apply to that model. However, if the (physically correct) temporal propagation is chosen, and either eqns. (5.9), or eqn. (5.11) are integrated forward in time -the solution has no choice but to be causal. The only spectral information available at a given time t depends on the wavevector k, not the frequency ω.
To understand the frequency spectrum, we must obtain a solution to the temporally propagated eqns. (5.9) or (5.11), and analyze the temporal history. That behaviour will be dominated by the reference frequency, which is bounded with 0 ≤ Ω(k) < 1/β. An impulsive Q, (e.g. one proportional to a delta function in time), does indeed force onto g(x, t) temporal frequency components above the maximum reference frequency Ω max = 1/β, but the natural behaviour of the wave equation itself will not. It is then the spatial profile of the impulse that reveals how wave disturbances will evolve forward or backward in space, at finite & bounded speeds c Ω (k) ≤ c. The step discontinuity reported by Buckingham [5] is not a problem with the wave equation, but an image of the impulsive source term driving the velocity potential.
Another issue relates to the case of non-zero γ, where that loss has been incorporated in the reference behaviour. In the discussion in subsection III.D of Buckingham [5], the frequency dependent phase speed diverged as c phase ≃ 2cβ 2 ω 3 /γ for large ω; here this feature can be recovered from the complex κ 2 0 (ω), since it is spatial FBACOU Acoustic waves: propagated in time or space?
Dr.Paul.Kinsler@physics.org http://www.kinsler.org/physics/ propagation which gives us frequency space properties. However, if loss is excluded from κ 0 (and hence c κ0 ), or is zero, or if the dressed κ R (c κR ) is used, then no divergence is seen. Instead the phase speed is pure imaginary, corresponding to evanescent (and not propagating) waves. Perhaps most importantly, the wave speed c Ω (k) relevant to temporal propagation exhibits no such anomaly.
Indeed, given the difficulty of interpreting complexvalued wave speeds, it can be helpful to explicitly exclude loss-like terms when evaluating phase and group velocities [26]; i.e. only ever calculate real valued wave velocities. Directional decompositions are then ideal, since they extract easily understood reference behaviors from complicated wave equations.

VI. CONCLUSION
Temporal and spatial propagation schemes for some typical acoustic wave equations have been compared and contrasted. The comparison enables a better judgment to be made as to which scheme is more practical in a given circumstance -e.g. if reflections are unimportant, a spatially propagated scheme is advantageous due to its more efficient handling of dispersion. Alternatively, for the ERE, temporal propagation incorporates two parameters (b 1 , b 2 ) into its reference Ω(k), but spatial propagation incorporates only one (b 1 ) into κ(ω). This suggests that for the ERE that temporal propagation is the natural choice, although if material dispersion were present the judgment might well become less simple.
Factorization methods mean that such comparisons can be made in a very transparent way -structurally similar temporally and spatially propagated wave equations can be compared term by term; just as the exact cou-pled bidirectional wave equations can be compared with the approximate unidirectional form that results from the single (and physically motivated) "slow evolution" assumption. This enables both quantitative and qualitative judgments to be made as to the significance of approximations, and/or the effect of any "source" terms that perturb or modify the free propagation. Thus these factorized wave equations have practical advantages for systems that are affected by driving terms, additional material dispersion, or nonlinearity -i.e. effects that make a numerical simulation the most practical way of finding a solution.
The potential for two competing ways of analyzing a situation can also add clarity to debate on specific properties of particular acoustic wave equation. As an example, the apparently non-physical response function for the VWE equation [5] is revealed as an image of the driving impulse, and not of the evolving wave profile. That profile is primarily controlled by c Ω (k), and not c κ (ω), since to be physically correct the waves must be propagated in time.
In summary, directional decompositions have been used to highlight the distinction between the physically accurate choice of temporal propagation and the often more convenient spatial propagation of waves. Some typical wave equations, containing terms for loss and with high-order derivatives have been analyzed, and the handling and consequences of these discussed.