A comparison of the factorization approach to temporal and spatial propagation in the case of some acoustic waves

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.


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. Notably, the impulse-response and causal properties of some of these wave equations have been analyzed by Buckingham [5], and the (perhaps surprising) 'non-causal' 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 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.
A feature of the factorization approach is that it forces us to play close attention to the propagation axis, and to the distinctions between temporal and spatial choices [8]; the main differences being highlighted on figure 1. However, beyond the basic conceptual differences, an additional aspect highlighted here is that for some wave models, the temporal and spatial reference behaviours can be very different. We will see that some terms in a wave equation can 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 2 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 3, acoustic waves in an elastic cylindrical rod in section 4, and the Stokes' equation and its bubbly generalization in section 5. This variety allows a discussion of both common and contrasting features of these directional decompositions. The conclusions are given in section 6.

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 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. A more general discussion of the choices and interpretations possible when choosing a propagation axis is given in [8].
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 equation for the relevant wave field  ( ) g r t , has, on the right hand side (RHS), a term Q which represents some general source term. Typically [5,9] 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 dependence on the density r  ( ) r could be added, setting Q=(1/ρ) ∇ ρ·∇g, and so match a wave equation used for ultrasound propagation [10,11]. 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 [12]; 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 wave property  ( ) g r t , (typically a velocity potential, or a displacement) 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 Figure 1. For temporal propagation (left), initial conditions cover all space at an initial time t i ; the final state at t f also covers all space. For spatial propagation (right), what are normally considered to be boundary conditions take up the role of the initial conditions; at both the starting location x i and the final location of x f , they specify the wave's behaviour over all time. In the unidirectional case, only forward-going wave evolution is treated, so that any reflections-whether from defined interfaces or other changes in the material properties-are neglected. Figures adapted with permission from [8].
non-local character which is allowed for using a convolution (''), with 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 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,13]; but the physical meaning alters with the interchange of the roles of time and space. Factorization methods have a long history [14], but their adaption, application, and adoption to wave propagation of the type and context proposed here has been lacking until more recently [15][16][17].
Since Q can contain any sort of behaviour, equation (2.2) is already useful: e.g. factorization could easily be applied to the nonlinear propagation addressed by Pinton et al [18]. In that case their initial equation (3) can be matched to equation (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. equation (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.

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 figure 1. This is useful for analyzing situations where signals need to be separated from reflections, but requires an explicit modeling of the medium's time-response, 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 equation (2.2) for temporal propagation, . Unfortunately, the frequency dependence of W means that it is not a good reference 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 3 .
If there is no temporal response, then p(t) ≡ p 0 δ(t) and p(ω)=p 0 , and the wave equation becomes Now decompose g into wave components g + and g − that evolve either forward or backward in space, with g=g + + g − . Then split equation (2.5) into two coupled first order wave equations for  3 For example, if restricted to waves near a frequency ω 1 in a narrow bandwidth w w D  1 , we can write w ¢  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 [19], and has been discussed at greater length in [20].
] as we propagate in time. The original second order equation can be recovered [13] 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 equation (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.
, then the wave evolves slowly as it propagates. This 'slow evolution' permits a temporal version of the unidirectional approximation [6,21] to be made 4 , setting g − =0. The forward wave + Note that if Q is an impulse such as a delta function, the slow evolution condition will not hold at that point. If W < 0 2 , then g or g ± do not oscillate as time passes, but instead decay.

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 figure 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 [22,23]. 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 [24,25], even enabling numerical investigations into exotic phenomena by means of artificial dispersions [26]. The angular spectrum approach [27,28] 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 . Next, we need a suitable reference wavevector κ, which is allowed to have a frequency (ω) dependence, but should not have a wavevector (k) dependence. Rewriting equation (2.2) then gives In the converse of the temporal case, where a frequency dependence was inconvenient, here k 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 in the weak spatial dispersion case, approximate 5 it as a frequency response s(ω). If there is no spatial response, then 4 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. 5 For example, if restricted to waves near a wavevector k 1 in a narrow bandwidth D  k 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. Now decompose g into wave components g + and g − that evolve either forward or backward in time: with g=g + + g − . Then split equation (2.13) into two coupled first order wave equations, i.e.
x leading to either wavevector or spatial domain forms ] 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 equation (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 this reference evolution is modified by the additional source terms, either the general term Q, or the diffraction term dependent on ∇ T 2 (or k T 2 , 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.
then the wave evolves slowly as it propagates. This 'slow evolution' permits a unidirectional approximation [6,21] to be made, setting g − =0. The forward waves w 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.

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.

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 [12]. It has a three-dimensional, inhomogeneous form for the velocity potential º 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 wavevector-frequency space, with w º  ( ) g g k, 3.1. 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 = where Ω(k)=ck, and its reference wave speed is c Ω =Ω/k=c. This is essentially the same as the wave equation As already explained, 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 equation (3.5). The forward wave 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 [29]. Fortunately, in this unidirectional wave equation, the two time derivative terms can be combined before re-applying the slow evolution approximations to get

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 2 ). This TDDE contains loss terms dependent on η, but it can be inadvisable to build loss into the reference behaviour of spatially propagated waves [30]. 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).
where κ(ω) is defined as 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 equation (3.1)as in e.g. Equation (2.1).
After defining , follow the same steps as for equations (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  3.11 x In the x, ω domain, equation (3.11) can be rewritten where the different conventions for κ and ω give rise to differing signs in the leading RHS terms (e.g. compare equations (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 T 2 , or the loss term dependent on η 2 . 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 changes slowly as it propagates in space, so that we can make a unidirectional approximation, setting g − =0. Then the condition containing k T 2 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 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 equation (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 ± .

Discussion
Both decompositions of the TDDE have the same (lossless) reference wave speed, i.e. the high frequency speed of sound = = k W 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.

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 [31], e.g. where the Poisson's ratio was negative [32]. 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 [31]) 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. n < | | 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 equation (4.1)Q is a source, driving, or other modification to the wave equation. In wavevector-frequency space, with w º ( ) g g k, , equation (4.1) becomes

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 w c where Ω(k) and new source term ¢ Q k are Now follow the same steps as for equations (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 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.
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 c ¶ = -W -W -+ where V 2+ is just V 2 calculated using only g + , i.e. with g − ≡ 0. This equation (and indeed the bidirectional equations (4.11)) have a RHS without any frequency (time) dependence, and so give the temporal propagation directly.

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 ¢ w Q are defined as Now following the same steps as for equations (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 w  ( ) g x, , becoming Expanding ¢ w Q in equation (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.
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 high order spatial derivatives (or, back in the bidirectional equations (4.19), a wavevector dependence), so that the spatial propagation is not given directly. Since this is not strictly compatible with the computational causality of an explict numerical solution [8,29,33], the alternate choice of temporal propagation seems preferable for this model-although a suitable implicit scheme could instead be used to calculate a solution.

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 3, the reference wave speeds were the same in the lossless case (where η=0), and similar otherwise. ), the wave speed is simply c. In contrast, the high-wavevector limit of W 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 figure 2, along with results for σ>0. We we see that the denominator causes an asymptote at =b k 1 1 2 , such that for σ<0 a stop band extends from =b k 1 1 2 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 s =b k 1 . The alternative case where b 1 is positive is shown on figure 3, where the low-wavevector limit is simply c 2 . In the high-wavevector limit W c 2 tends to b b 2 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. Equations (4.19) or (4.22)) being less straightforward. The reference wave vector κ 2 gives a wave speed 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 w = b c 1 2 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 figure 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.
It is important to note that such differences are caused by the difference between the temporal and spatial environments experienced by the wave field g. Further, as highlighted by figure 1, the starting assumptions are not the same. For example, a case where the reference frequency is positive, but the comparable reference wavevector is imaginary suggests a spatially localized excitation-although the presence of significant source terms may complicate this simple interpretation.

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 [34]. 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 [35]. However, the 3D VWE differs from the common form of the Stokes' equation in that it has spatial operators of  · g and not  g 2 . 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 [35] 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 w º  ( ) g g k, and = = +

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 w  ( ) g k , is now written as where Ω(k) and the new source term Q' are defined as 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 figure 5 the reference wave speed = W W c k 2 2 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 equations (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  Eqn. (5.7) can be transformed so as to apply to  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.
then the wave evolves slowly as it propagates, so a unidirectional approximation can be made, setting g − =0.
The forward waves + Notice that a k-dependence has appeared on the term proportional to Q in equation (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 section 3.

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 equation In both equations (5.13) and (5.14) the reference wavevector κ 0 (ω) is defined as The S/VWE equation (5.14) contains an imaginary part in its reference wavevector κ 0 , which is not always desirable. This can be circumvented by splitting k 0 2 into real and imaginary parts using , and moving κ I to the RHS 6 : The S/VWE can then be re-expressed as although it will also be convenient to merge the κ I term into the Q term as follows Figure 5. The S/VWE reference wave speed c Ω =Ω/k as a function of scaled wavevector k c β. In the bubble-free Stokes' case, b = 0 2 and the wave speed is a constant at = W c c 1 for all k. 6 Note the imaginary κ I part is not squared.
speed the wave becomes evanescent, since k R 2 is negative and so κ R imaginary. Here the directed g ± evolve as specified by κ, but that that reference evolution is modified and coupled by Q, γ, and k T 2 . 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 equation (5.34) will fail at low frequencies; also that for non-zero γ, both the Q and the dynamic viscosity (equation (5.35)) conditions will fail near w b  1 2 2 . Equation (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 Both here in equation (5.37) and in the bidirectional equation (5.33) there is a non-zero γ-dependent source term.

Discussion
For this combined S/VWE model, just as for the ERE in section 4, 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 If γ=0, all describe the same propagation and are in agreement, despite the difference between the high wavevector limit where  ¥  W ( ) c k 0 and the transition from propagation to evanescence in w k ( ) c at b w = 1 2 2 . The behaviour of the frequency-dependent forms c κ are shown in figure 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 b ¹ 0. Since the β 2 term in the VWE has a counterpart in the ERE's equation (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 equations (5.9), or equation (5.11) are integrated forward in time-the solution has no choice but to be causal, as long as the loss term g µ ¶ t is handled appropriately (or more simply, just neglected). 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 equations (5.9) or (5.11), and analyze the temporal history. That behaviour will be dominated by the reference frequency, which is bounded with 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  W ( ) 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 b w g  c c 2 phase 2 3 for large ω; here this feature can be recovered from the complex k w ( ) 0 2 , since it is spatial propagation which allows us access to 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 complex-valued wave speeds, it can be helpful to explicitly exclude loss-like terms when evaluating phase and group velocities [30]; 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.

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 coupled 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