{'fraction_non_alphanumeric': 0.035500282645562466, 'fraction_numerical': 0.01582815149802148, 'mean_word_length': 4.497410399834266, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comparison of fractional wave equations for power law attenuation in ultrasound and elastography

A set of wave equations with fractional loss operators in time and space are analyzed. It is shown that the fractional Szabo equation, the power law wave equation, and the fractional Laplacian wave equation in the causal and non-causal forms all are low frequency approximations of the fractional Kelvin-Voigt wave equation and the more general fractional Zener wave equation. The latter two equations are based on fractional constitutive equations while the former wave equations are ad hoc, heuristic equations. We show that this has consequences for use in modelling and simulation especially for applications that do not satisfy the low frequency approximation, such as shear wave elastography. In such applications the wave equations based on constitutive equations are the preferred ones.


Introduction
Models for ultrasound attenuation based on relaxation losses have a long history, see e. g. Markham et al. (1951); Bhatia (1967). Waag and his coworkers also contributed with their work on multiple relaxation losses (Nachman et al. (1990)). Similar, but much simpler relaxation models are used for attenuation in salt water and air. In salt water, the two important relaxation processes are due to boric acid and magnesium sulphate (Ainslie and McColm (1998)) and in air they are due to nitrogen and oxygen (Bass et al. (1995)). In contrast, in medical ultrasound and elastography, attenuation for both compressional and shear waves often follows a power law which at first sight is very different from a relaxation model: where y usually is between 0 and 2 and ω is angular frequency. The subscript k is used to indicate that this is the imaginary part of the wave number, k, and to distinguish it from the order of the fractional derivative, α, which will be used later.
In such media the number of relaxation processes may in practice be uncountable and it has not been possible to make models for attenuation that are so grounded in physical processes as for salt water and air. Nevertheless, the multiple relaxation model with a few processes can be used to model the power-law attenuation as Waag's group and several others have done (Tabei et al. (2003); Yang and Cleveland (2005)) over a limited frequency range. But, the relaxation parameters lose their clear physical meaning in this case and become just parameters of a mathematical model. Thus on the one hand there is power law attenuation which is experimentally observed in many different complex materials and on the other hand physical models, exemplified by the multiple relaxation model. In the last decade or so advances have been made in providing wave equations that partially bridge this gap. They model power law attenuation through the use of fractional, i.e. noninteger, derivatives. To varying degrees these equations are physics-based, but especially for some viscoelastic polymers there is a good physical foundation, Näsholm and Holm (2013). However, there is seldom as direct a connection as for the relaxation models for air and salt water. Nevertheless, they represent one step on the way to a deeper understanding of the underlying physics. They also provide alternative simulations models for wave propagation, often characterized by being parametrized with a low number of parameters.
For this paper the time fractional derivative is easiest to define in the frequency domain where it is an extension of the Fourier transform of an n'th order derivative: The fractional derivative of arbitrary order can be understood as a generalization where the integer n is replaced with a real number α. The close connection between power laws and fractional derivatives is evident from this definition. The time domain equivalent involves a convolution integral implying that the fractional derivative is non-local and has a powerlaw shaped memory, Podlubny (1999).
The purpose of this paper is to analyze and compare several of the fractional wave equations, to understand their origins, and to relate them to each other. The ultimate objective is to find wave equations which are physically viable. Another more practical objective is to determine which wave equations that are best suited for simulation of medical ultrasound and of shear wave elastography. One of the questions asked is also: What is the most fundamental property to model with a wave equation? Is the objective only to model a power law attenuation, or is it more fundamental to model a medium with viscoelastic properties so that different power law attenuation laws are achieved in different frequency regions?
Due to the fact that it is usually the low frequency region which is of interest, many wave equations model only this region. But there are applications where the high frequency region is important also, one particularly important one is in shear wave elastography. Even when the interest is only in the low frequency solution, such as for ultrasound imaging, then in order for the solution to be physically viable, it is often an advantage that wave equations give correct solutions beyond this region. The key to obtaining such results is the viscoelastic constitutive equation. The paper therefore builds a case for the viewpoint that it is the viscoelastic constitutive equation which is the more fundamental and physical property, rather than the power law characteristics.
It should be noted that low frequencies here means low in comparison to the value of a time constant (τ σ in e.g. Eq. (4)). Thus for e.g. compressional waves in medical ultrasound, low frequencies stretch well beyond the usual range and up into the 100's of MHz or higher. For such an application the power law has an exponent, y, between one and two. On the other hand, it turns out that for shear waves in the human body in elastography, the low frequency limit is at less than 10 Hz (Holm and Sinkus (2010)). Thus elastography operates above the low frequency limit with y less than unity.

Physcially viable solutions
There are some criteria that need to be satisfied by the solutions of wave equations in order for them to be physically viable. The first is that of causality: No output is allowed before there is an input.
The second property is that the loss function cannot rise arbitrarily fast and has to satisfy: This means that the exponent, y, has to be less or equal to one in the high frequency limit. This condition was first given by Weaver and Pao (1981) and more rigorously by Hanyga (2013). It follows from causality, passivity, and linearity and is a consequence of the property that viscoelasticity can be modelled as a sum of many spring-dashpot mechanical models with realistic, i.e. positive, constants, Beris and Edwards (1993).

Fractional wave equations
We will give wave equations for plane waves in a homogenous medium. Fig. 1 illustrates how the wave is modified by transmission through a slab of thickness z. The wave number is in general complex and given by k = β k − iα k where the real part describes dispersion and the imaginary part gives the attenuation. Note that many papers use the opposite sign for the exponent in the travelling wave. This will result in somewhat different wave equations. This is shown in the Appendix A which has been included in order to facilitate comparisons between papers using different sign conventions. Fractional time loss operator, type 1 One of the oldest fractional wave equations is due to Caputo (1967). In the terminology of Wismer (2006) it is where u is the displacement, c 0 is the phase velocity at zero frequency, α 0 = τ σ /(2c 0 ) in Eq. (1), and the condition (ωτ σ ) y−1 1 has to be satisfied. Here we will call it the fractional Kelvin-Voigt wave equation. This equation and all subsequent wave equations in this paper could equally well have been formulated for the particle velocity, v =u, and in the case of compressional waves also the pressure, p.
For low frequencies relative to the relaxation time constant τ σ , α and the power law exponent of Eq. (1) are related by y = α + 1 where 0 < α ≤ 1.
The fractional Zener wave equation derived by Holm and Näsholm (2011) is an extension of Eq. (4) with an additional loss term which contributes to a different behavior for very high frequencies: These equations are here said to have a fractional time loss operator of type 1, because there is a term with a mix of a fractional time derivative and a non-fractional Laplacian in the loss term.

Fractional time loss operator, type 2
Type 2 of the fractional time loss operator is an operator which has only temporal derivatives. The distinction between types 1 and 2 is only to aid in the presentation of the various loss operators. The first type 2 operator is the general Szabo convolution model: where * is the temporal convolution operator. The causal loss operator, L y (t) is given by the power law attenuation that it describes and has different forms for even and odd integer values of y and for the non-integer case, see Szabo (1994Szabo ( , 2004 for details. This wave equation can be written in terms of fractional derivatives (Chen and Holm (2003)) and in the terminology of Zhang et al. (2012) the fractional Szabo equation is: Kelly et al. (2008) added an extra loss term to Eq. (8) in order to model power law losses with exactly the same exponent, y, for all frequencies. Therefore we will call it the power law wave equation: Fractional Laplacian loss operator Chen and Holm (2004) proposed the fractional Laplacian wave equation: It was amended by the addition of an extra loss term by Treeby and Cox (2010): where τ = −2α 0 c y−1 0 and η = 2α 0 c y 0 tan(πy/2). The additional loss term makes the wave equation causal at low frequencies, so a proper name is the causal fractional Laplacian wave equation. Note however that the sign of this last term will depend on the definition of the Fourier transform used in the derivation, as explained in Appendix A.

Convolution loss operator
The most general way to express the loss term in the wave equations is by means of a convolution operator. It can model losses which are different from power laws also. The convolution terminology is taken from Sec. (4.2.2) in Mainardi (2010) and in the creep representation it is: where Ψ(t) = [dJ(t)/dt]/J g is the rate of creep, J(t) is the creep compliance, and J g = J(0 + ) is the glass compliance. The various fractional operators discussed in this paper are such that they result in power law attenuation. Thus by comparison with the fractional Szabo equation, Eq. (8), the convolution with the rate of creep corresponds to the operator: The alternative relaxation representation of Eq. (12) is, Mainardi (2010): where Φ(t) = [dG(t)/dt]/G g is the rate of relaxation, G(t) is the relaxation modulus, and G g = G(0 + ) is the glass modulus. Hence the fractional Kelvin-Voigt wave equation, Eq. (4), can be formulated by using: Towards the end of the paper, we will also formulate the Nachman-Smith-Waag model in the convolution framework.

Wave equations derived from viscoelastic models
The most general constitutive equation considered here is the fractional Zener model, or the fractional standard linear solid. Following Holm and Näsholm (2011) it is: where σ (t) is stress, ε(t) is strain, and the elastic modulus is the equilibrium modulus. In the general case, stress and strain are tensors, but here it is sufficient to represent them as scalars in order to demonstrate the properties of the resulting wave equations. This model is also a special case of a more general convolution model between stress and strain expressed by the creep compliance, J(t), or the relaxation modulus, G(t). This constitutive equation is well founded in material science as a good description of material properties as reviewed in Ch. 3 of Mainardi (2010) and the survey paper Näsholm and Holm (2013). There the conditions τ α σ ≥ τ β ε and α ≤ β are justified also. It is also common to assume the even stricter condition that α = β which will be done in the rest of this paper also. This condition is necessary in order to make the system thermodynamically well-behaved. As an example, the model has been used to model shear wave propagation in human tissue in elastography, Klatt et al. (2007).
Although the constitutive equation gives a good description of materials, our understanding of why this is so is still lacking, i.e. what is it in the underlying microstructure of the materials that gives this behavior? Therefore the fractional parameter α cannot be determined directly from a physical model, but it is closely related to the power law behavior of the attenuation see belowl. One of the advantages of the model over the multiple relaxation model (Nachman et al. (1990)) is parsimony, i. e. that it is described by only four parameters (G e , τ σ , τ ε , and α). This is fewer than the number of parameters required for the multiple relaxation model except for the very simplest of media such as air and salt water. This is usually an advantage in simulations.
When Eq. (16) is combined with conservation of mass and momentum, it leads to the fractional Zener wave equation, Eq. (6), for which the dispersion relation is: The simpler, three-term Kelvin-Voigt constitutive equation is obtained from Eq. (16) by setting τ ε = 0. Examples of fitting of that model to shear wave elastography data can be found in Zhang et al. (2007). The fractional Kelvin-Voigt constitutive equation leads to the wave equation first derived by Caputo,Eq. (4).
It should be noted that there are several ways that this result can be extended for the special case of pressure waves in fluids and liquids. First, the model may include nonlinear acoustics as done in Prieur and Holm (2011);Prieur et al. (2012). Second, in a complete acoustic theory for fluids and liquids a thermal loss term may also be significant. This was analyzed in the former two references as well as summarized in Holm et al. (2013). This paper will build on the asymptotic results previously determined for the low, intermediate and high frequency regions of Eq. (6). For completeness the main results from Holm and Näsholm (2011) are repeated here. The attenuation follows these asymptotes: The phase velocity also has three regimes: In the case of the fractional Kelvin-Voigt model with τ ε = 0, the third regime will disappear for both the attenuation and the phase velocity. In that case one can recognize the wellknown result for the viscous loss model with a non-fractional derivative, α = 1, (Stokes' equation) which gives attenuation proportional to ω 2 for low frequencies and proportional to √ ω for high frequencies.
The following sections builds on a set of approximations derived for the low-frequency regime which are given in Appendix B. They will be used subsequently to find various interesting approximations to the lossy wave equation.

Wave equations with fractional temporal derivative loss operator
An approximation which is valid for low frequencies is obtained by replacing k 2 by (ω/c 0 ) 2 , Eq. (B.8), in the viscous terms of the dispersion equation for the fractional Kelvin-Voigt equation (Eq. (17), with τ ε = 0): This leads to: This is the Szabo wave equation written in terms of fractional derivatives as given in Eq. (8). Thus this is one alternative for a low frequency variant of Eqs. (4) and (6). Because it can be derived directly from the causal Kelvin-Voigt wave equation, it is causal for low frequencies as demonstrated in e.g. Zhang et al. (2012) also. Here it has been derived as a low frequency approximation to Eq. (4) and this is fully consistent with the low frequency approximation which is a central step in Szabo's own derivation. If one for a moment neglects that it is based on a low frequency assumption, then it is straightforward to show that for high frequencies, the attenuation will be proportional to ω 1+α/2 , 0 < α ≤ 1. But at these frequencies, the equation is used for something which violates its assumptions. It is therefore not surprising that this also conflicts with the condition for a physically viable solution given by Eq. (3).

Two loss terms
In the derivation above, the starting point was the Kelvin-Voigt wave equation or its generalization the Zener wave equation. Then low frequencies were assumed according to Eq. (5) so that the time-space substitution k ≈ ω/c 0 could be used. An alternative is to say that the power law attenuation and the accompanying dispersion required from causality is the goal. The objective is then to find a wave equation which satisfies power law attenuation exactly. This is the starting point of Kelly et al. (2008) and it is the same as requiring that the low frequency approximation for the wave number in Eq. (B.5) should be satisfied exactly for all frequencies, not just for low frequencies.
That equation is the same as Eq. (7) in Kelly et al. (2008) when the difference in sign conventions for the Fourier transform have been taken into account (see Appendix A).
Squaring both sides of Eq. (B.5) gives this dispersion relation: − α 2 0 cos 2 (πy/2) (iω) 2y (22) Inverse Fourier transformation in space and time yields Eq. (9) which has an extra loss contribution compared to Eq. (21). It is repeated here: Kelly et al. (2008) say that the last term usually can be neglected since it is proportional to α 2 0 and thus is small compared to the main loss term. Therefore it doesn't make much difference. This is certainly the case for low frequencies according to the condition in Eq. (5), where α 0 is proportional to τ σ .
But for high frequencies, the enforced power law relation of Eq. (B.5) still holds exactly, and does so for however high frequencies one tries for. As the frequency increases above the low frequency condition of Eq. (5), the last loss term will play a larger and larger role. This wave equation therefore models the same power law attenuation for all frequencies. Most other wave equations only model this attenuation for low frequencies.
But, as long as the exponent, y, is larger than 1, the condition for a physically viable solution according to Eq. (3) is violated. This is the limitation of the wave equation of Kelly et al. (2008): At low frequencies the extra loss term can be neglected, and at high frequencies it may make the solution physically unviable. This problem was also illustrated by Johnson and McGough (2011) through numerical simulations.
The paper by Kowar et al. (2011) takes the same approach as Kelly et al. (2008) and tries to model a power law attenuation exactly. As a result their Eq. (36) is similar to Eq. (22). However, Kowar et al. (2011) go further as they have observed that the result is non-causal when y > 1. They then amend their wave equation in order to make it causal. The result, Eq. (61) in Kowar et al. (2011) is a wave equation which after some reorganization can be written on this form: where c ∞ = c 0 (1 + α 1 ) is the asymptote of the phase velocity at infinite frequency. The first four terms are similar to those of the fractional Zener wave equation, Eq. (6). The final Fourier transform term can be approximated for different frequency regions and it can be shown that it will merge into the second and fourth terms. In our view Kowar et al. (2011) give an argument, albeit a rather indirect and involved one, for why the fractional Zener wave equation has the desirable property of physical validity for all frequencies.

Wave equations with a fractional Laplacian loss operator
In this section we go back to Eq. (B.7) and use substitutions based on k ≈ ω/c 0 in the loss terms as in the previous section. The difference is that we only partially substitute time for space. By assuming low frequencies and replacing ω α by ω · (c 0 k) α−1 in the first loss term and by using (c 0 k) α instead of ω α in the second loss term one gets: The equivalent wave equation using the fractional Laplacian introduced by Chen and Holm (2004) is: First let's neglect the last term. As sin(πα/2)τ α σ c α−1 0 = 2α 0 c α 0 = 2α 0 c y−1 0 using Eq. (B.3) and y = α + 1, one then obtains the lossy wave equation, Eq. (10), with a fractional Laplacian attenuation operator first proposed by Chen and Holm (2004). Due to the neglect of the second term, this equation does not satisfy the Kramers-Kronig condition as has been pointed out by several researchers. The fractional Laplacian wave equation of Eq. (10) is therefore a non-causal low frequency approximation to the fractional Kelvin-Voigt wave equation of Eq. (4).
By keeping all terms in Eq. (26), and using the substitutions of the previous paragraph, one gets: This is almost Eq. (11), the equation derived by Treeby and Cox (2010) to be causal. There is however one minor difference and that is the sign of the last term. This is caused by the differences in the sign of the Fourier transform kernel between Treeby and Cox (2010) and this paper, see the Appendix A. This wave equation is therefore a causal low frequency approximation to the fractional Kelvin-Voigt wave equation of Eq. (4) expressed with fractional Laplacians instead of time-fractional operators. Unlike for the Szabo wave equation, it is not possible to find an analytic solution of the fractional Laplacian wave equation beyond the region of validity and into the high frequencies, and in that way show that it violates Eq. (3). The substitution of time with space is entirely due to the implicit low-frequency approximation and the use of Eq. (B.8). Unfortunately we doubt that there is much significance beyond that of the fractional Laplacian, i.e. little or no relationship to the spatial structure of the medium.

Convolution formulation of the multiple relaxation model
Finally we would like to point out some connections between fractional models and the Nachman-Smith-Waag multiple relaxation model (Nachman et al. (1990)). As shown in Section 4.2 of Näsholm and Holm (2013), conservation of mass and conservation of momentum give the general dispersion relation where the frequency-domain generalized compressibility is κ(ω) ≡ ε(ω)/σ (ω) and therefore given directly by a constitutive equation such as Eq. (16). This relation can be inverse Fourier transformed into where κ(t) is the temporal representation of the generalized compressibility. Notice that the formulation of Eq. (29) is related to the wave equation (1) in Verweij et al. (2013) (see also the references therein). These two wave equations are mathematically similar if nonlinear terms are disregarded, and if m(t) in Verweij et al. (2013) is identified as m(t) = ρ 0 c 2 0 κ(t). Related similarity is observed between Eq. (37) and the random effective propagation model Eq. (16) in Garnier and Sølna (2010). Näsholm and Holm (2011) considered the integral generalization of the Nachman-Smith-Waag generalized compressibility, Nachman et al. (1990) where the set of τ ν are the relaxation time constants and κ ν are the corresponding discrete compressibility contributions. For the continuous distribution of compressibility contributions of Eq. (30), κ ν (Ω) denotes the contribution at relaxation frequency Ω. Now combining Eq. (30) with Eq. (28) and taking the inverse Fourier transform gives the convolution form of the multiple relaxation wave equation: where H(t) is the Heaviside step function and the relation c 2 0 = 1/(ρ 0 κ 0 ) is applied. This is a special case of Eq. (12) and now we want to find the convolution kernel, the rate of the creep. Using the arguments of Näsholm (2013), Eq. (32) can be rewritten in the form where L signifies the Laplace transform from the Ω domain to the t domain. Thus the multiple relaxation wave equation can be expressed in the creep representation form of Eq. (12) with This is an equivalent, but often straightforward form of the wave equation than that given in Nachman et al. (1990). In addition, we note that Eq. (33), the convolution representation of the Nachman-Smith-Waag wave equation, is similar to the Szabo model, Eq. (7), so the Szabo kernel can also be found. Other parallels can also be drawn such as the relationship between the fractional Kelvin-Voigt and Zener wave equations. This has been discussed in Näsholm andNäsholm (2013).

Conclusion
Four of the common wave equations with time-and spacefractional loss operators, here denoted as the fractional Szabo equation, the power law wave equation, the fractional Laplacian wave equation, and the causal fractional Laplacian wave equation (Eqs. (8) -(11)), have been analyzed. It has been shown that they can all be derived as low frequency approximations to the fractional Kelvin-Voigt wave equation which again is a special case of the more general fractional Zener wave equation, Eqs (4) and (6).
The latter two wave equations are based on fractional constitutive equations which are fairly well established as descriptions of realistic materials, while the former four wave equations were derived so that they produce the power law solutions which are often observed in measurements. The fractional Kelvin-Voigt and the fractional Zener wave equations therefore possess more desirable physical properties than the others.
The fractional Kelvin-Voigt and Zener wave equations can also be viewed as special cases of wave equations with convolution kernels for loss terms. Moreover, Eq. (33), provides a straightforward convolution form of the Nachman-Smith-Waag wave equation which is linked to the creep representation of Eq. (12). What is so attractive about using fractional derivatives rather than general convolution kernels, is that they provide solutions that fit well with the power laws that are frequently observed.
For use in numerical evaluation where the excitation source and medium of propagation satisfy the low frequency approximation, such as in medical ultrasound, there is a wide choice of equations to use, and one can use either one of the exact or low-frequency approximations given here. The choice may be dictated by numerical considerations which are beyond the scope of this paper.
It should however be pointed out that for simulation of waves that do not satisfy the low frequency approximation, and shear waves in elastography is an important example of that, these low frequency approximations cannot be used. Here either the fractional Kelvin-Voigt or the fractional Zener equations are the only ones which describe the properties of the wave propagation and they are also the best ones for simulations.