Tractable calculation of the Green's tensor for shear wave propagation in an incompressible, transversely isotropic material.

Assessing material properties from observations of shear wave propagation following an acoustic radiation force impulse (ARFI) excitation is difficult in anisotropic materials because of the complex relations among the propagation direction, shear wave polarizations, and material symmetries. In this paper, we describe a method to calculate shear wave signals using Green's tensor methods in an incompressible, transversely isotropic (TI) material characterized by three material parameters. The Green's tensor is written as the sum of an analytic expression for the SH propagation mode, and an integral expression for the SV propagation mode that can be evaluated by interpolation within precomputed integral functions with an efficiency comparable to the evaluation of a closed-form expression. By using parametrized integral functions, the number of required numerical integrations is reduced by a factor of 102-109 depending on the specific problem under consideration. Results are presented for the case of a point source positioned at the origin and a tall Gaussian source similar to an ARFI excitation. For an experimental configuration with a tilted material symmetry axis, results show that shear wave signals exhibit structures that are sufficiently complex to allow measurement of all three material parameters that characterize an incompressible, TI material.

Introduction

Elastic properties of materials can be measured by observing shear wave propagation following an acoustic radiation force impulse (ARFI) excitation and relating the propagation speed to a model of the material. For example, linear, elastic, homogeneous, and isotropic materials can be characterized using two elasticity constants such as the Lamé parameters λ and μ. For nearly incompressible materials, including many soft biological tissues, these parameters differ by a factor on the order of 106, and the corresponding difference in longitudinal and shear wave speeds is on the order of 103. Typically, ultrasonic tracking methods observe only shear wave propagation and calculate the shear modulus μ from the wave speed c and material density ρ by the relation (Graff 1991, Lai et al 2010)

The relation between wave speed and material properties is more complicated in anisotropic materials because of the complex dependence of the velocity relative to the propagation direction, material symmetries, and wave polarizations. For example, in a linear, elastic, transversely isotropic (TI) material, a symmetry axis exists and the material can be characterized by five elasticity constants (Lai et al 2010). Skeletal muscle is an example of a TI material with the symmetry axis defined by the direction of the muscle fibers (Gennisson et al 2010, Wang et al 2013). Stretched or compressed polyvinyl alcohol phantoms also have TI symmetry with the symmetry axis determined by the axis of deformation (Gennisson et al 2007, Chatelin et al 2014, Urban et al 2015). Measurements of shear wave propagation in these materials are typically performed using the experimental geometry shown in figure 1(a) (see, also, figure 2 of Chatelin et al (2014)) and observe propagation in the horizontal (Z = 0) plane in the direction at an angle θ relative to the symmetry axis. These measurements observe the SH propagation mode with displacements in the direction so that the wave polarization is perpendicular to the Z = 0 plane. The phase velocity v (θ) for this propagation mode is given by the relation (Wang et al 2013) where μ and μ are shear moduli for wave propagation in the longitudinal and transverse directions relative to the material symmetry axis. The corresponding group velocity V(θ) has an elliptical shape (Wang et al 2013),

Figure 1.

(a) Experimental configuration commonly used for investigations of shear wave propagation in an incompressible, TI material (see, also, figure 2 of Chatelin et al 2014). The left side shows a sketch of a linear ultrasound transducer and transversely isotropic material with the material symmetry indicated by the gray lines respresenting, for example, skeletal muscle fibers. The transducer can rotate about a vertical axis to observe shear waves for a range of propagation directions. The experimental X, Y, Z coordinate system shows the ARFI excitation force along the Z axis and the material symmetry axis and propagation direction in the X − Y (Z = 0) plane. Polarization vectors for the SH (slow shear), SV (fast shear), and P (longitudinal) propagation modes are defined relative to the − plane shown in gray. Ultrasonic tracking measures the Z component of the shear wave displacement signal and is sensitive only to the SH propagation mode. (b) A more complicated experimental configuration in which the material symmetry axis and propagation direction are not restricted to the X − Y plane. Measurements of the Z component of shear wave displacement are sensitive to both the SH and SV propagation modes.

Experimental observations of shear wave propagation typically measure the group shear wave speed and determine μ and μ from the axes of the ellipse.

The relation between shear wave speed and propagation direction is even more complicated for experimental configurations such as figure 1(b) where the material symmetry axis is not oriented perpendicular to the excitation axis, or where waves are tracked in three dimensions including axial positions above and below the Z = 0 plane. For example, using finite element simulations, Rouze et al (2013) observed complicated shear wave signals for an excitation and tracking configuration with a symmetry axis tilted at an angle of 45° relative to the horizontal plane. The advantage of these geometries is that they allow experimental measurements of both the SH and SV propagation modes and are sensitive to more material properties than can be measured from just the SH propagation mode. However, analysis of these shear wave signals is difficult because the group propagation speeds do not reduce to simple expressions such as (3). Instead, this analysis is typically performed by calculating the shear wave signals using either finite element simulations (Palmeri et al 2005, Rouze et al 2013) or Green’s tensor calculations (Bercoff 2004, Chatelin 2015) and comparing the calculated and measured signals.

In this paper, we use Green’s tensor calculations to model shear wave propagation in an incompressible, TI material. This process calculates the shear wave signal (, t) at an observation position and time t by dividing the excitation force (, t) into spatial and temporal source voxels s and summing the contributions from each voxel. The Green’s tensor G ( − , t − t) gives the relative contribution to the signal from each source voxel. Often, these calculations are performed in the temporal frequency domain with angular frequency ω = 2π f so that, assuming the force can be factored into a spatial function () and temporal window W(t), components u(, ω) of the shear wave signal can be written as where is the relative position between the source and observation positions, the indices i and n take on the values 1, 2, and 3 corresponding to the x, y, and z vector components, and summation over the component index n of the force is implied.

One difficulty with Green’s tensor calculations of shear wave signals using (4) is that, for typical calculations, G(, ω) must be evaluated for 108 − 1015 combinations of r, r, and ω. Thus, these calculations have been used primarily for linear, isotropic materials (Bercoff 2004, Rouze 2018) where the Green’s tensor is given by a closed-form expression (Aki and Richards 2002, Kausel 2006). Unfortunately, for more complicated materials, closed-form expressions for the Green’s tensor are known only for a few special-case materials (Vavryčuk 2001, 2007), and are not known for the case of a TI material, or an incompressible TI material. Instead, the Green’s tensor can be evaluated using integral expressions (Willis 1980, Wang and Achenbach 1995), including the specific case of a TI material (Gridin 2000). However, these expressions must be evaluated using numerical integration, and because the integrand in these expressions can oscillate rapidly, it is necessary to sample the integrand on a dense mesh. Thus, these calculations are too slow to be practical for applications involving a large number of combinations of r, r, and ω.

In this paper, we describe a tractable method to perform integral calculations of the Green’s tensor for an incompressible, TI material. This process allows the Green’s tensor to be expressed as the sum of two terms, each given by the product of multiplicative factors and a relatively simple, parametrized integral. Numerical integration can be used to precompute the integral factors. Then, when calculating the shear wave signal using the sum (4), interpolation can be used to evaluate the integral factors efficiently. This procedure allows the Green’s tensor to be evaluated with an efficiency comparable to the evaluation of a closed-form expression. Furthermore, the number of numerical integrations required to precompute the integral factors is typically reduced by a factor of 102 − 109 depending on the number of combinations of r, r, and ω that would be required if numerical calculations were used to evaluate the Green’s tensor for each combination. This reduction in computational complexity is the key factor that allows tractable calculation of shear wave displacements using Green’s tensor methods.

Background

Incompressible, transversely isotropic (TI) materials

In Sections 2 and 3, we use a Cartesian x, y, z coordinate system such as shown in figure 2 to model the material, wave propagation, and Green’s tensor. This coordinate system is distinct from the X, Y, Z coordinate system used in figure 1 to describe the experimental configuration. The calculation of shear wave signals for a specific experimental configuration is described in Section 4.2.

Figure 2.

Coordinate system with x, y, and z axes used for the analysis of the material, wave propagation, and Green’s tensor in Sections 2 and 3. This coordinate system is distinct from the experimental XYZ coordinate system in figure 1. The relation between the xyz and XYZ coordinate systems is described in Section 4.2. (a) The xyz coordinate system is oriented so that the z axis is aligned with the material symmetry axis , and the position vector is specified by the angles θ and ϕ. To simplify the analysis in Section 3.1, the axes are rotated so that ϕ = 0. (b) Wave propagation is in the direction with wave vector = k at an angle θ relative to the material symmetry axis . The polarization vectors are defined relative to the − plane shown in gray with the SH polarization perpendicular to the − plane, and the SV and P polarizations in the − plane. Explicit expressions for polarization vectors are given in terms of θ and ϕ in (24).

In the limit of small displacements, the stress-strain relationship in an anistropic material is linear and can be described by a generalized Hooke’s law as where σ and ε are components of the time-dependent stress and strain tensors, respectively, and c are the components of a time-independent, fourth-order stiffness tensor. Each index can assume the value 1, 2, or 3, and summation over repeated indices is implied. Symmetries of the stress and strain tensors and the existence of a strain energy allow the stiffness tensor to be expressed in terms of 21 independent elements (Lai et al 2010). For the case of a TI material with symmetry axis , rotation and reflection symmetries allow the stiffness tensor to be expressed in terms of five independent elements. By orienting the coordinate system so that , the stress strain relation (5) can be expressed using Voigt notation as a matrix product with the stiffness matrix C (Lai et al 2010), where missing elements are zero.

The stress-strain relation (6) can also be written in terms of the compliance matrix S = C−1 using Young’s moduli E and E, shear moduli μ and μ, and Poisson’s ratios ν and ν where the longitudinal (L) and transverse (T) directions are defined relative to the material symmetry axis (Lai et al 2010), with the relation

Explicit relations for the elements C11, C13, C33, C44, and C66 in (6) can be expressed in terms of E, E, μ, μ, ν, and μ through the relation S = C−1,

For the case of an incompressible TI material, the fractional volume change, or dilation, of an infinitesimal volume subjected to stresses is zero. The dilation is given by the trace of the strain tensor (Lai et al 2010) so that, using (7)

Both terms in this expression must be zero, and the Poisson ratios for an incompressible TI material are given by

Thus, three material parameters are required to characterize an incompressible, TI material. In this study, we use the parameters μ, μ, and E/E as used previously by Rouze et al (2013). However, with the relations (8) and (11), these parameters are not unique, and any three independent combinations of the parameters μ, μ, E, and E can be used. In particular, muscle strength, characterized by the longitudinal Young’s modulus E, can be determined from the three material parameters.

In addition, using the Poisson ratios (11) in expressions (9) indicates that the C11, C33, and C13 elements of the stiffness matrix diverge in the limit of an incompressible, TI material. However, differences between pairs of these elements remain finite,

As indicated in figure 1, wave propagation in a TI material can be described in terms of the P (longitudinal), SH (slow shear), and SV (fast shear) propagation modes identified by their polarization relative to the propagation direction and material symmetry axis (Tsvankin 2012, Carcione 2015). Polarization of the P mode is quasi-longitudinal and corresponds to the acoustic wave. The SH propagation mode corresponds to shear wave motion with a purely transverse polarization perpendicular to the to the − plane. The SV propagation mode has shear wave propagation with a quasi-transverse polarization in the − plane. Note that Rouze et al (2013) refer to these modes as the QL, PT, and QT modes, respectively. For plane wave propagation in a direction oriented at an angle θ relative to the symmetry axis, the phase velocity v of the SH propagation mode is given by (2). In the limit of an incompressible TI material, the polarization of quasi-longitudinal P mode becomes purely longitudinal, and the propagation velocty v diverges. Similarly, for an incompressible TI material, the polarization of the quasi-transverse SV mode is purely transverse, and the propagation velocty v is given by (Chadwick 1993, Papazoglou et al 2006, Rouze et al 2013)

Thus, all three parameters required to characterize an incompressible TI material can be measured by observing SH and SV shear wave propagation in the material.

Green’s tensors for constrained TI materials

Chatelin et al (2015) analyzed their measurements of shear wave propagation in a transversely isotropic material using the “Anisotropy III” special-case Green’s tensor from Vavryčuk (2001) and Vavryčuk (2007). The stiffness matrix for this case is given by expression (25) of Vavryčuk (2001) and also in the text just above expression (B3) of Vavryčuk (2007). This expression indicates that the stiffness matrix for the special-case Green’s tensor is similar to the stiffness matrix (6) for a general TI material with two additional assumptions,

The explicit expression for the Green’s tensor in this special case is given by expression (B3) of Vavryčuk (2007) and expresion (5) of Chatelin et al (2015) and will not be reproduced here. Note that the expressions from Vavryčuk (2007) and Chatelin et al (2015) are presented for the case of viscoelastic materials with complex, frequency-dependent moduli. We do not consider viscoelastic materials in this paper and will assume the moduli are real. However, because our development is presented in the temporal frequency domain, it can easily be extended to the case of viscoelastic materials by considering complex, frequency-dependent moduli.

The two constraints (14) on elements of the stiffness matrix imply that three parameters characterize the TI material for this special case. Chatelin et al (2015) analyzed their experimental measurements using the acoustic wave speed and the two wave speeds for shear wave propagation along and across the material symmetry axis. These parameters were sufficient to analyze the wave speeds for shear wave propagation in their experimental setup where the material symmetry axis was oriented perpendicular to the excitation axis, and only the SH propagation mode was observed, see figure 1(a) or figure 2 of Chatelin et al (2014). However, this analysis is not sensitive to all three material properties that are needed to characterize shear wave propagation in an incompressible, TI material because these measurements are not sensitive to the ratio E/E in the phase velocity for the SV propagation mode in (13). Thus, in the incompressible limit, the Green’s tensor used by Chatelin et al (2015) does not fully characterize an incompressible, TI material specified by three independent parameters.

In addition to the special-case Green’s tensors described by Vavryčuk (2007), several studies (Payton 1975, 1983, Chadwick and Norris 1990, Burridge et al 1993) have investigated the Green’s tensor for a TI material described by a stiffness matrix with four independent parameters and the constraint

For the case of an incompressible, TI material considered in this study, this relation can be simplified by substituting the relations (9) for each term C and evaluating the result using the Poisson ratios (11) for an incompressible, TI material. The result of this process demonstates that the constraint (15) is equivalent to the relation

Thus, the ratio E/E is not an independent material parameter, and the expression for the Green’s tensor based on the relation (15) cannot fully describe shear wave propagation in an incompressible, TI material.

Green’s tensor for an incompressible, TI material

Integral expression for the Green’s tensor

Assuming a stress-strain relation of the form (5), the equation of motion for the Green’s tensor G(, t) for the i component of displacement due to a spatial and temporal Dirac delta distribution in the x direction is given by (Kim et al 1994, Červený 2001) where δ is the Kronecker delta symbol.

In this paper we define the sign convention for the Fourier transform so that the forward transform from the coordinate to frequency domain is given by the relation and the inverse transform is given by

This sign convention agrees with the MATLAB convention for the forward and reverse transforms. Also note that we use the symbol i to denote as in (18) and (19), and as an index as in (5). The specific meaning should be clear from the context.

Calculating the 4-dimensional Fourier transform of (17) with spatial frequency and temporal frequency ω = 2πf gives where the matrix L has elements L = ck − ρω2δ. This expression can be solved for G (, ω) by multiplying by L−1 and calculating the inverse Fourier transform,

The inverse matrix L−1 in (21) can be expressed as the sum of three terms using the eigenvectors and eigenvalues λ of L,

The eigenvectors in (22) correspond to the polarization vectors for the SH, SV, and P propagation modes. As shown in figure 2(b), for plane wave propagation with wave vector = , the propagation direction is given by and the polarization vectors are given by

The eigenvalues λ of L are related to the phase velocities of the SH, SV, and P propagation modes. For the case of an incompressible, TI material, the phase velocity v diverges, and λ and λ are given by with v and v given by (2) and (13), respectively. Here, we neglect the P propagation mode in (21) and (22) due to its large speed and only consider the Green’s tensor describing shear wave motion. Then, from (21) and (22), G(, ω) can be written as the sum of two terms corresponding to the SH and SV propagation modes,

As shown in the Appendix, (, ω) in (26) can be evaluated in closed-form using tabulated integrals. For indices i = 1, 2 and n = 1, 2, where , , , and the sign of k0 has been chosen to give outgoing waves when combined with the e+ time dependence. For i = 3 or n = 3, in (24) and . After accounting for the sign convention used in the Fourier transform in (18) and (19), this result agrees with expression (4.14) from Gridin (2000) that was obtained using a different approach.

We have not found a closed-form expression similar to (27) for the SV propagation mode in (26). Instead, we can obtain an expression for (, ω) suitable for numerical integration by writing the expression in spherical coordinates, where = k from (23), components of are shown in figure 2(a), and

Following the procedure described by Kim et al (1994), this expression can be simplified by extending the k integration to the range −∞ ≤ k ≤ ∞ and performing the θ integration over the range 0 ≤ θ ≤ π/2. The same result can be obtained using the θ integration range π/2 ≤ θ ≤ π, or as one-half the result using the full range 0 ≤ θ ≤ π. Then, using the full θ range and evaluating the k integration using equation (14) of Kim et al (1994), (, ω) is given by the sum of two terms,

In this expression, the integral in the first term is over the surface of the unit sphere. The delta function in the second term reduces the integration range to a line integral over a great circle of the unit sphere in a plane perpendicular to . Both of these integrals have a finite range and can be evaluated using numerical integration. After accounting for the sign convention used in the Fourier transform in (18) and (19), the result (31) is equivalent to the results from Gridin (2000), Willis (1980), and Wang and Achenbach (1995) that were obtained using forward and inverse Radon transformations instead of the Fourier transforms used here. Expression (31) is also equivalent to equation (15) of Kim et al (1994) which involves integrals over a hemisphere of the unit sphere.

Efficient calculation of the Green’s tensor

Using (26), components of the Green’s tensor are expressed in terms of seven parameters; four dynamic variables r, θ, ϕ, and ω, and three variables which characterize the material, μ, μ, and E/E. Here, we assume that (, ω) is evaluated using the closed-form expression (27). For the evaluation of (, ω), only the product μE/E appears in the phase velocity v in (13), and six parameters are needed to evaluate the integral expression (31). To reduce this number further, we can orient the coordinate system so that ϕ = 0. Also, the phase velocity v from (13) can be written as where is the phase velocity along the symmetry axis at θ = 0, and

After inserting this expression in (31), the first integral can be parametrized using the phase α = rω/v0 so that (, ω) can be written as where the surface and line integrals are given by and These integrals can be precomputed for materials with specific values of Δ using a dense grid of values of α and θ. Then, when summing over combinations of r, r, and ω in (4), (α, θ, Δ) and (θ, Δ) can be evaluated efficiently by interpolation, and these values can easily be combined with the multiplicative factors in (34). With this procedure, the Green’s tensor can be calculated with an efficiency comparable to the evaluation of a closed-form expression.

Methods

Precomputation of the surface and line integrals

The surface integral (35) was evaluated by tabulating the integrand on a two dimensional mesh and using the trapezoidal rule to evaluate the integrals for both θ and ϕ numerically. The mesh used a step size of 0.005 rad for both Δθ and Δϕ. This step size was selected after varying the size and comparing the calculated Green’s tensor to exact results, see section 5.1. The surface integral was evaluated using phase angles α in the range 0 ≤ α ≤ αmax with a step size Δα = 0.2 rad, and angles θ in the range 0° ≤ θ ≤ 180° with a step size of Δθ = 0.5°. For negative values of α, the surface integral was evaluated using the relations Re [(−α, θ, Δ)] = Re [(α, θ, Δ)] and Im [(−α, θ, Δ)] = −Im [(α, θ, Δ)]. The value of αmax was set to 400 rad based on the maximum values of ω and r and the minimum value of v0 for the materials and propagation geometries considered.

The line integral (36) was evaluated using the procedure described by Kim et al (1994). For this approach, the k axis is rotated so that it is oriented in the direction with the great circle of the unit sphere in a plane perpendicular to . The integration is performed as a one dimensional integral using a single variable ϕ around the rotated k axis. However, when performing this integration, the phase velocity v and polarization vectors in the integrand of (36) are evaluated relative to the (unrotated) axis. This process was performed by defining a set of points {x, y, z} distributed around the unit circle in the z = 0 plane, and then rotating these points to the positions {, , } using a rotation matrix determined from the angles θ and ϕ with the specific value ϕ = 0 from section 3.2,

These coordinates were used to calculate the polar angle and azimuthal angle for each of the rotated points, and these angles were used to evaluate the phase velocity v and polarization vectors in the integrand of (36). The integration was performed using the trapezoidal rule for the same values of θ as used to evaluate the surface integral.

Results are presented in section 5 for materials with μ = 25 kPa, μ = 9 kPa, and E/E = 0.16, 0.36, and 0.64. These are the same materials considered by Rouze et al (2013). All calculations were performed on a Linux cluster with an average CPU speed of 2.6 GHz using Matlab (The MathWorks, Natick, MA). For each material, and for the range of variables α and θ described above, the (non-optimized) computation time for the surface and line integral data was roughly 128 minutes.

Calculation of shear wave signals

Shear wave signals were calculated by dividing the excitation force into voxels and performing the Green’s tensor sum (4) over the source positions for each desired observation position . In this sum, the force (r), relative position = − r, and shear wave displacement (, ω) are defined relative to an experimental XYZ coordinate as shown in figure 1. However, the Green’s tensors (, ω) and (, ω) have been calculated using the xyz coordinate system shown in figure 2. To transform between these two coordinate systems, we define a rotation matrix R with components R = · so that, for example, the vector can be written as

Then, shear wave displacements in the XYZ coordinate system can be calculated using the Green’s tensor sum (4) as

The rotation matrix in (38) and (39) can be calculated by expressing , , and in terms of the X, Y, Z coordinates as follows. First, in figure 2(a), the z axis is aligned with the material symmetry axis so that = . Also from figure 2(a), with the xyz coordinate system oriented so that ϕ = 0 as in Section 3.2, the y axis is perpendicular to the − plane and . Finally, . Then, elements R of the rotation matrix can be calculated using the relation R = · with each unit vector expressed relative to the XYZ coordinate system. For the sums in (39), the symmetries of (, ω) and (, ω) with respect to the indices i and n, and the choice ϕ = 0 from section 3.2 imply that it is only necessary to calculate the (i, n) = (1, 1) and (2, 2) components of (, ω), and the (i, n) = (1, 1), (2, 2), (3, 3), and (1, 3) components (, ω). In addition, for the results presented in section 5, the ARFI force is assumed to be directed along the Z axis, and only the Z component of the shear wave displacement is measured. Thus, the component sums in (39) can be simplified to include only these components.

Components of the Green’s tensors in (39) were calculated using (27) for (, ω) and (34) for (, ω) using the surface and line integral functions (35) and (36), respectively. These integral functions were precomputed for the specific materials described in the last paragraph of section 4.1, and the functions were evaluated for specific combinations of , and ω in the sum (39) by interpolation at the coordinates α = ωr/v0 and θ within these functions. For the line integral, a linear interpolation was used to evaluate the function for a specific value of θ. For the surface integral, the interpolation was performed in two steps. First, a linear interpolation was used for the θ variable to obtain a function of α, and then spline interpolation was performed in this function for specific values of α. For this last step, the spline interpolation was performed once to find the piecewise interpolation polynomial, and then this polynomial was evaluated for all required values of ω.

The excitation window function W(t) was assumed to be a rectangular function with duration T = 200 μs so that the Fourier transform W(ω) in (39) is given by

Calculations were carried out using frequencies in the range −10 kHz ≤ f ≤ 10 kHz with a step size Δf = 20 Hz corresponding to a temporal step size of 50 μs. After calculation of u(, ω) using (39), the time-dependent shear wave signal u(, t) was calculated using a discrete, inverse Fourier transform.

Results are presented in section 5 for the case of two force functions, a point source located at the origin, and a tall Gaussian source similar to an ARFI excitation described by the relation with σ = 0.4 mm, σ = 0.4 mm, and σ = 7.5 mm. For the case of the tall Gaussian source, the force function was divided into 0.33×0.33×0.33 mm3 voxels, and the Green’s tensor sum (39) was truncated to include only the source voxels with an amplitude greater than 5% of the maximum force (Rouze et al 2018).

Results are presented in section 5 for the two possible experimental configurations shown in figure 1 with the symmetry axis oriented in a plane perpendicular to the excitation axis as in figure 1(a), and with the symmetry axis tilted at an angle of 45° in the X − Z plane as in figure 1(b). These results show shear wave signals in the X = 0, Y = 0, and Z = 0 planes for the case of the three materials described in the last paragraph of section 4.1. These planes are divided into 160 × 160 pixels, each with a dimension of 0.25 × 0.25 mm for a total of 25,600 observation positions in each plane. For the case of the tall Gaussian source, the (non-optimized) computation time for the shear wave signals in each plane was roughly 16 hr.

Results

Validation of the numerical integration procedure

The numerical integration procedure used to calculate the Green’s tensor for the SV propagation mode was validated by applying exactly the same procedure to the case of the SH propagation mode and comparing the results with the analytic, closed-form expression (27) for (, ω). The surface and line integral functions in (35) and (36) were calculated using the same values of α and θ as used for the SV calculations, with the only changes being the replacement of the SV polarization vectors by the SH polarization vectors from (24), and the replacement of v from (32) and (33) with v given by (2),

Interpolation within the tabulated integral functions was performed using exactly the same procedure described in section 4.2 for the calculation of shear wave signals.

Figure 3 shows a comparison of the analytic (true) results (black line) using (27) with results from the numerical calculation (dashed red line) over the frequency range 0 ≤ f ≤ 5 kHz for the specific case of the position r = 10 mm and angle θ = 45°. Comparisons are presented for the real (top row) and imaginary (bottom row) parts of the (, ω), (, ω), and (, ω) components. The signals have MKS units of mN−1 or Pa−1 m−1. These results show nearly perfect agreement between the numeric and analytic signals. Differences between these results can be quantified using the RMS percentage difference ΔRMS defined as

Figure 3.

Validation of the numerical integration procedure by comparison of the analytic results (black lines) calculated using (27) with results from the numerical calculations (dashed red lines) for the real (top row) and imaginary (bottom row) parts of the (, ω), (, ω), and (, ω) components of the Green’s tensor for the SH propagation mode. The signals have MKS units of mN−1 or Pa−1 m−1. Results are shown for a material with μ = 25 kPa and μ = 9 kPa for the specific position r = 10 mm and angle θ = 45°. Near perfect agreement is observed for the (, ω) and (, ω) components. For the (, ω) component, the analytic value is zero, and any nonzero signal from the numerical calculation is the result of numerical inaccuracies from the use of finite step sizes in the numerical integrations (35) and (36), or from roundoff errors.

This calculation gives ΔRMS = 0.392% and ΔRMS = 0.012% for the real and imaginary parts of (, ω), respectively, and ΔRMS = 0.008% and ΔRMS = 0.002% for the real and imaginary parts of (, ω), respectively. For the case of (, ω), the choice ϕ = 0 from section 3.2 requires the analytic result (27) to be zero, and any nonzero signal calculated numerically is the result of numerical inaccuracies that result from the use of finite step sizes for the angles θ and ϕ in the numerical integrations (35) and (36), or from roundoff errors. The results shown in figure 3 indicate that these effects give very small errors on the order of 10−15 of the calculated signals.

Shear wave signals calculated using a point source excitation

Figure 4 shows an example of shear wave signals in the X = 0, Y = 0, and Z = 0 planes calculated using a point source excitation located at the origin. The three planes overlap and are shown with semi-transparency in the center, and are also duplicated without transparency at positions displaced from the center. The excitation configuration is shown in figure 1(a) with the material symmetry axis oriented along the X axis and the excitation axis in the direction. The images show the component of displacement at time t = 2.2 ms. These signals are shown for a material with μ = 25 kPa, μ = 9 kPa, and E/E = 0.16, and an excitation duration T = 200 μs. A movie showing the time evolution of these signals is included with supplementary material associated with this paper.

Figure 4.

Views of shear wave signals in the X = 0, Y = 0, and Z = 0 planes at time t = 2.2 ms from a point source excitation located at the origin. The planes overlap and are shown with semi-transparency in the center, and are also duplicated without transparency at positions displaced from the center. The experimental configuration is shown in figure 1(a) with the material symmetry axis positioned along the X axis and the excitation force directed along the direction. The component of displacement is shown. Signals are shown for the case of a 200 μs excitation duration in a material with μ = 25 kPa, μ = 9 kPa, and E/E = 0.16. A movie showing the time evolution of these signals is included with supplementary material associated with this paper.

The shear wave signals in the Z = 0 plane of figure 4 are determined only by the SH propagation mode and show the expected elliptical shape from (3) with semi-major and semi-minor axes determined from the μ = 25 kPa and μ = 9 kPa shear moduli. However, for shear wave signals observed at axial positions above and below the Z = 0 plane, the − plane in figure 1 is tilted, and both the SH and SV propagation modes contribute to the shear wave signals.

Figure 5 shows shear wave signals in the X = 0, Y = 0, and Z = 0 planes as shown in figure 4 for experimental configurations with the symmetry axis oriented (a) along the X axis as in figure 1(a), and (b) tilted at an angle of 45° relative to the X axis as in figure 1(b). As in figure 4, these signals are shown for the case of a point source excitation located at the origin at time t = 2.2 ms. Signals are shown for materials with μ = 25 kPa, μ = 9 kPa, and E/E = 0.16, 0.36, and 0.64, and an excitation duration T = 200 μs. The top row of figure 5 shows the same signals as in figure 4. Movies showing the time evolution of these signals are included with supplementary material associated with this paper.

Figure 5.

Normalized shear wave displacement signals in the X = 0, Y = 0, and Z = 0 planes (see figure 4) from a point source excitation located at the origin for cases with the material symmetry axis oriented (a) along the X axis as in figure 1(a), and (b) tilted at an angle of 45° relative to the X axis in the − plane as in figure 1(b). The excitation force is oriented in the direction, and the component of displacement is shown. Signals are shown at a time t = 2.2 ms for materials with μ = 25 kPa, μ = 9 kPa, and E/E = 0.16, 0.36, and 0.64. The top row of signals are the same as the signals in figure 4. Movies showing the time evolution of these signals are included with supplementary material associated with this paper.

As with figure 4, the signals in the Z = 0 plane in figure 5(a) with the symmetry axis directed along the X axis show the expected elliptical shape from (3) with the same semi-major and semi-minor axes. However, we note that for the case with E/E = 0.64, structure from the SV propagation mode contributes to the signals and distorts the elliptical shape near Y = 0. For the case of the symmetry axis tilted at an angle of 45° relative to the X axis in figure 5(b), both the SH and SV propagation modes contribute to the signals and give more complicated shapes as compared to figure 5(a).

Shear wave signals calculated using a tall Gaussian excitation

Figure 6 shows shear wave signals in the X = 0, Y = 0, and Z = 0 planes as shown in figure 4 for the case of a tall Gaussian excitation source and experimental configurations with the symmetry axis oriented (a) along the X axis as in figure 1(a), and (b) tilted at an angle of 45° relative to the X axis as in figure 1(b). The size of the source was determined using σ = 0.4 mm, σ = 0.4 mm, and σ = 7.5 mm in (41). As described in section 4.2, the source was divided into 0.33 × 0.33 × 0.33 mm3 voxels, and the Green’s tensor sum (39) was truncated to include only the source voxels with an amplitude greater than 5% of the maximum force. As in figures 4 and 5, the signals in figure 6 are shown at time t = 2.2 ms for materials with μ = 25 kPa, μ = 9 kPa, and E/E = 0.16, 0.36, and 0.64, and an excitation duration T = 200 μs. Movies showing the time evolution of these signals are included with supplementary material associated with this paper.

Figure 6.

Normalized shear wave displacement signals in the X = 0, Y = 0, and Z = 0 planes (see figure 4) from a tall Gaussian source for cases with the material symmetry axis oriented (a) along the X axis as in figure 1(a), and (b) tilted at an angle of 45° relative to the X axis in the − plane as in figure 1(b). The excitation force is oriented in the direction, and the component of displacement is shown. The source had dimensions σ = 0.4 mm, σ = 0.4 mm, and σ = 7.5 mm in (41) with an excitation duration of 200 μs. Signals are shown at a time t = 2.2 ms for materials with μ = 25 kPa, μ = 9 kPa, and E/E = 0.16, 0.36, and 0.64. Movies showing the time evolution of these signals are included with supplementary material associated with this paper.

The shear wave signals obtained using the tall Gaussian excitation in figure 6 show two key differences compared to the signals from the point source excitation in figures 4 and 5. First, the width of the shear wave signals in the Z = 0 plane from the tall source is greater than from the point source due to contributions from source points at positions above and below the Z = 0 plane that arrive in the Z = 0 plane at times later than contributions from a point source at the origin. Second, the Z-dependent structure seen in the point source calculations with the tilted excitation in figure 5(b) is washed out in figure 6 due to the tall source. Nevertheless, structure remains in the Z = 0 plane of figure 6(b) and can be used to distinguish among the signals obtained with different values of E/E.

Discussion

In this paper, we have described a method for the tractable calculation of the Green’s tensor for shear wave propagation in an incompressible, TI material, and have presented sample shear wave signals for the case of a point source, and a tall Gaussian source similar to an ARFI excitation. The motivation for this work is based on two key factors. First, previous applications of the Green’s tensor methods in TI materials have been based on special-case solutions of the equation of motion where closed-form expressions were available so that the sum over source voxels in (4) or (39) could be performed efficiently. However, these solutions only depend on the shear moduli μ and μ, and not on the ratio of Young’s moduli E/E, and thus, are not sensitive to all three material parameters needed to characterize an incompressible, TI material. In the current work, we have described a procedure to calculate the Green’s tensor and shear wave signals that is sensitive to all three of these material parameters. Thus, these techniques can be used to analyze shear wave signals in nearly incompressible TI materials such as skeletal muscle, and thereby allow the measurement of material properties such as the longitudinal Young’s modulus E which characterizes muscle strength.

The second key factor described in this work is a tractable method to calculate the Green’s tensor. Integral expressions for G(, ω) have been described previously, but have been used sparingly, because of the large computational effort to evaluate the integral numerically due, in part, to the rapid oscillation of the integrand in (28). For the case of an incompressible, TI material considered here, it is possible to obtain relatively simple expressions for the phase velocity v in (13) and polarization vectors in (24). Then, the integral expression (34) for (, ω) can be written in terms of integral functions (α, θ, Δ) and (θ, Δ) that are parametrized in terms of a small number of variables. After precomputing these functions, components of the Green’s tensor can be evaluated by interpolation within these functions and combining the results with multiplicative factors as in (34). This interpolation can be performed with an efficiency that is comparable to the evaluation of a closed-form expression.

Of course, this procedure also requires the numerical computation of the integral functions (α, θ, Δ) and (θ, Δ). However, because these functions are parametrized, they only need to be calculated for a relatively small number of parameter values. For the calculations reported here, the integral functions are calculated for 7.2 × 105 combinations of α and θ. Without this parametrization, the numerical integration would need to be performed for each combination of , , and ω used in the Green’s tensor sum (4) or (39). For example, a typical ARFI excitation force might be divided into 103 − 104 voxels and calculations might be performed using 103 − 104 frequencies. Also, the number of observation points might vary from 102 positions along a single axis, to 107 positions throughout a volume. Then, it would be necessary to perform numerical integrations for 108 − 1015 combinations of , , and ω depending on the specific problem under consideration. Thus, by using parametrized integral functions, the number of required numerical integrations is reduced by a factor of 102 − 109. This reduction in computational complexity is the key reason that allows the tractable calculation of shear wave displacement signal using the Green’s tensor procedure described herein.

Finally, we note that only elastic materials have been considered in this study. However, soft biological tissues are typically viscoelastic (Fung 1993), and it is common to model shear wave propagation in these materials using viscoelastic materials (Bercoff 2004, Gennisson 2010, Chatelin 2015). Characterization of viscoelastic materials is often done by working in the temporal frequency domain and using a stress-strain model that includes a loss mechanism dependent on the viscosity of the material. A common material model is the Voigt model with a frequency dependent shear modulus μ(ω) = μ0 + iωη where μ0 is the material stiffness and η is the viscosity. Other material models include the fractional Kelvin model (Zhang et al 2007) and models with attenuation expressed as a power-law function of frequency (Waters et al 2000, Rouze et al 2018). Because the Green’s tensor model considered in this study has been presented in the temporal frequency domain, it would be straightforward to extend the model to include viscoelastic materials with frequency-dependent material properties.

Conclusion

This paper describes a method for the tractable calculation of the Green’s tensor for shear wave propagation in an incompressible, transversely isotropic (TI) material following an acoustic radiation force impulse (ARFI) excitation. Wave propagation is modeled using SH and SV propagation modes defined by the wave polarization relative to the plane determined by the propagation direction and the material symmetry axis. Components of the Green’s tensor are written as the sum of an analytic expression for the SH propagation mode, and an expression involving numerically calculated integral functions for the SV propagation mode. These integral functions are parametrized in terms of two dynamic parameters and can be precomputed for each material under consideration. Then, interpolation within these functions allows the Green’s tensor to be calculated with an efficiency comparable to the evaluation of a closed-form expression. The numerical integration procedure is validated by applying it to the SH propagation mode and comparing the calculations to the known analytic results. Parametrization of the integral functions reduces the computational complexity by a factor on the order of 102 − 109 depending on the specific problem under consideration, and thereby makes the Green’s tensor approach tractable. Results are presented for the case of a point source positioned at the origin and a tall Gaussian source similar to an ARFI excitation. For an experimental configuration with the material symmetry axis tilted relative to the excitation axis, results show that shear wave signals exhibit structure that are sufficiently complex to allow measurement of all three material parameters that characterize an incompressible, TI material.