é...
; :21! L.
4.; 25w
e. C.
3.1,.353
s .
. a
. m»?
h. .r.
.a , “w.
i
:
x. .. .
.. 3.3.
$-
nugn s
unﬂnvuxs ".‘L‘
dankngf .
£131.;
Q "a... :- Any-5.2.x f’!‘
:Isbul‘ 51:32:513‘ 0x!» .6
yo!>.lt.l3 |.a.,l. Ill! )5!ka
.. .1 .ﬁrﬁfaﬁaﬂ
..(‘...oxf.~o.1¢
:. 1.1:! 0.7.0, )I
II Ilﬁr
i,‘4¥:r§!”ll
2.! i. 3'! 4.? .3415!
l. ...‘>I$tc,’z 9t
1‘ .. .58.:2."avﬂ.uv
(a 01?. 11’! all} 5..
{$1. (2... ..v,« I...
I! .1112
17v. Iii.
«v.51! kl. In
l 173!!!»
aa‘mﬁizéaag
2% a..:.::.n¢l «3
.‘ .Y .4
‘3» J4 .V .2...
.1..a.y.. . . .v.(... ...
LIBRARY
«@5213 Michigan §tate
21 ”3 University
This is to certify that the
dissertation entitled
TRANSIENT ULTRASONIC FIELDS IN POWER LAW
ATI'ENUATION MEDIA
presented by
James F. Kelly
has been accepted towards fulﬁllment
of the requirements for the
PhD. degree in Electrical Engineen‘gq
Mo V. I 3r 2008
Date
MSU is an Afﬁrmative Action/Equal Opportunity Employer
PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.
DATE DUE DATE DUE DATE DUE
5/08 K:lProj/Acc&Pres/ClRC/DateDue.indd
TRANSIENT ULTRASONIC FIELDS IN POWER LAW
ATTENUATION MEDIA
By
James F. Kelly
A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Electrical Engineering
2008
ABSTRACT
TRANSIENT ULTRASONIC FIELDS IN POWER LAW
ATTENUATION MEDIA
By
James F. Kelly
Ultrasonic waves in biological media experience frequency dependent attenuation.
Extensive measurement of the attenuation coefﬁcient in human and mammalian tissue
in the ultrasonic range has revealed a power law dependence on frequency, where
the power law exponent ranges between 0 and 2. For most tissue, the power law
exponent ranges between 1 and 1.7, which cannot be explained by classical theories for
ultrasonic absorption, such as thermo—viscosity or molecular relaxation. The purpose
of this thesis is threefold: 1) to understand the analytical structure of transient ﬁelds
in power law media, 2) to provide a possible description of the physical mechanism
responsible for power law attenuation in biological media, and 3) to develop analytical
models for transient, three-dimensional sound beams in power law media.
Insight into general dissipative media is gained by studying the approximations
available in viscous media. The Stokes wave equation is considered in the time do-
main, and an asymptotic, causal Green’s function is rigorously derived and veriﬁed
using the material impulse response function (MIRF) approach. A lossy impulse
response for the Stokes wave equation is derived for calculations of transient ﬁelds
generated by ﬁnite aperture radiators. Expressions for the uniform circular aper-
ture (in both the nearﬁeld and the farﬁeld), the uniform rectangular aperture in the
nearﬁeld, and the spherical shell in the nearﬁeld are then derived . Power-law media
is then studied using fractional partial differential equations (FPDEs), which add loss
to the wave equation with a time-fractional or space-fractional derivative. A FPDE
is derived that exactly describes power law attenuation, and analytical, time-domain
Green’s functions in power law media are derived for exponents between 0 and 2. To
construct solutions, stable law probability distributions are utilized, which are widely
used in the study of anomalous diffusion and in the study of fractal media. For ex-
ponents strictly less than 1, the Green’s functions are causal, while for exponents
greater than or equal than 1, the Green’s functions are noncausal.
To address the lack of causality, alternate power law FPDES based on fractional
spatial operators are considered: the Chen-Helm wave equation and a spatially disper-
sive wave equation. Green’s functions are derived for both equations, yielding causal
solutions for all applicable power law exponents. The Chen-Holm equation is shown
to be non-dispersive, while the spatially dispersive wave equation supports a phase
velocity predicted by. the Kramers-Kronig relations. To address the physical basis for
FPDEs, a fractal ladder network is proposed as a model for the stress-strain relation-
ship in tissue. This constitutive equation is based on a lumped-parameter inﬁnite-
ladder topology involving alternating springs and dashpots to capture the viscoelastic
and self-similar properties of tissue. This ladder network yields a stress-strain con-
stitutive equation involving a time-fractional derivative. The Caputo—Wismer FPDE
is derived from this constitutive equation. Finally, the impulse response derived for
viscous media is generalized to power-law media. Expressions for ﬁnite apertures
are then derived in dispersive media, thus forming the basis for ultrasonic image
simulation in biological media.
Dedicated to my aunt, Cmdr. Judith E. Epstein, M.D. (Naval
Medical Research Center, Bethesda, MD), and grandmother Mary C.
Sinclair (1914 - 2004).
iv
ACKNOWLEDGMENTS
Firstly, the author thanks his advisor, Robert J. McGough (Department of Elec-
trical and Computer Engineering), for his consistent support during the past four
years. Secondly, the author thanks his doctoral committee for providing guidance
and insightful discussion: Shanker Balasubramaniam (Department of Electrical and
Computer Engineering), Leo Kempel (Department of Electrical and Computer En-
gineering), and Brian Feeny (Department of Mechanical Engineering). The author
also acknowledges the following people for additional discussion, guidance, software,
and illustrations: John P. Nolan (Department of Mathematics and Statistics, Ameri-
can University), Mark M. Meerschaert (Department of Department of Statistics and
Probability, Michigan State University), Thomas L. Szabo (Department of Biomedi-
cal Engineering, Boston University), Stephen W. Wheatcraft (Department of Geolog-
ical Sciences and Engineering, University of Nevada, Reno), Geoff Robinson (CSIRO
Mathematical and Information Sciences), Liyong Wu, (GE Research), Jason Biel
(Department of Electrical and Computer Engineering, Michigan State University),
T. Douglas Mast (Department of Biomedical Engineering, University of Cincinnati),
and Amy Albin (Department of Zoology, Michigan State University). Finally, the
author acknowledges the Michigan State University Graduate School for providing a
dissertation completion fellowship during the summer semester of 2008.
TABLE OF CONTENTS
LIST OF FIGURES .............................
1 Introduction ................................
1.1 Ultrasound Absorption Mechanisms ...................
1.1.1 Thermoviscous Model ......................
1.1.2 Molecular Relaxation .......................
1.1.3 Micro-Heterogeneity and Viscoelasticity ............
1.2 Phenomenological Power Law Models ..................
1.2.1 Power Law and Szabo Wave Equation .............
1.2.2 Chen-Holm and Spatially Dispersive Wave Equations .....
1.2.3 Fractal Ladder Models and the Caputo—Wismer Wave Equation
1.3 Fourier Tiansform Convention ......................
1.4 Content of Thesis .............................
2 Stokes Wave Equation ..........................
2.1 Introduction ................................
2.2 Formulation ................................
2.3 Green’s Function Solution ........................
2.3.1 Causal Approximation (n = 0) ..................
2.3.2 First Order Term (n = 1) ....................
2.4 Green’s Function Evaluation .......................
2.5 Green’s Flinction Decomposition .....................
3 Lossy Impulse Response in Viscous Media ..............
3.1 Lossy Impulse Response: Theory ....................
3.1.1 Impulse Response in Lossless Media ...............
3.1.2 Impulse Response in Viscous Media ...............
3.2 Uniform Circular Piston: Nearﬁeld ...................
3.2.1 Lossless Media ..........................
3.2.2 Viscous Media ...........................
3.2.3 Numerical Results ........................
3.3 Uniform Circular Piston: Farﬁeld ....................
3.3.1 Lossless Media ..........................
3.3.2 Viscous Media ...........................
3.3.3 Numerical Results ........................
3.4 Rectangular Sources: Nearﬁeld .....................
3.4.1 Lossless Media ..........................
3.4.2 Viscous Media ...........................
CDCDQOTCﬁuD-ODH
10
11
11
13
13
14
15
17
18
19
21
25
26
26
27
3.4.3 Numerical Results ........................
3.5 Spherical Shell ...............................
3.5.1 Lossless Media ..........................
3.5.2 Viscous Media ...........................
3.5.3 Numerical Results ........................
3.6 Discussion .................................
3.6.1 Implementation Issues ......................
3.6.2 Numerical Evaluation and Aliasing . _ ..............
3.6.3 Power Law Media .........................
Power Law and Szabo Wave Equations ................
4.1 Introduction ................................
4.2 FPDE Formulations of the Szabo and Power Law Wave Equations . .
4.2.1 The Szabo Wave Equation (0 < y < 1 and l < y < 2) .....
4.2.2 The Szabo Wave Equation (y = 1) ...............
4.2.3 Frequency-Independent (y = 0) and Viscous (y = 2) Media . .
4.2.4 Power Law Wave Equation ....................
4.2.5 Linear Attenuation Media (y = 1) ................
4.3 3D Green’s Function ...........................
4.3.1 Sub-Linear Power Law Media (0 < y < 1) ...........
4.3.2 Super-Linear Power Law Media (1 < y S 2) ..........
4.3.3 Linear Power Law Media (y = 1) ................
4.4 2D Green’s Functions ...........................
4.5 Stochastic Interpretation .........................
4.6 Numerical Results .............................
4.7 Discussion .................................
4.7.1 Causality and the Super-Linear Attenuation Paradox .....
4.7.2 Physical Interpretation and Application to Ultrasonic Imaging
4.7.3 Comparison to Frequency Domain Models ...........
Chen-Helm and Spatially Dispersive Wave Equations .......
5.1 Introduction ................................
5.2 Chen-Holm Equation ...........................
5.2.1 Formulation ............................
5.2.2 Dispersion Relationship .....................
5.2.3 Causality .............................
5.3 Spatially Dispersive Wave Equation ...................
5 . 3. 1 Formulation ............................
5.3.2 Dispersion Relationship .....................
5.3.3 Relationship to Szabo Wave Equation ..............
5.3.4 Causality .............................
vii
39
40
40
41
43
43
43
44
58
58
59
59
60
61
63
65
68
70
73
76
78
79
80
83
83
84
85
91
92
92
93
94
95
99
5.4 3D Green’s Emotion: Chen-Helm Equation .............. 100
5.4.1 General case (0 < y < 1 and 1 < y S 2) ............. 100
5.4.2 Symmetric Stable Distributions ................. 101
5.4.3 Leading Term (n = 0) ...................... 103
5.4.4 Higher-Order (n 2 1) Terms ................... 104
5.4.5 Linear Attenuation Media (y = 1) ................ 105
5.5 3D Green’s Function: Spatially Dispersive Wave Equation ...... 106
5.5.1 Leading Term (n = 0) ...................... 108
5.5.2 Relationship to the Power Law Green’s Function ........ 109
5.6 1D and 2D Green’s Functions ...................... 109
5.6.1 1D Green’s Functions ....................... 109
5.6.2 2D Green’s Functions ....... _ ................ 111
5.7 Numerical Results ............................. 112
5.7.1 Chen-Holm Wave Equation ................... 112
5.7.2 Spatially Dispersive Wave Equation ............... 113
5.8 Discussion ................................. 117
5.8.1 Stochastic Model for the Chen-Holm Wave Equation ..... 117
5.8.2 Stochastic Interpretation of the Spatially Dispersive Wave
Equation .............................. 119
Fractal Ladder Models and the Caputo-Wismer Equation . 122
6.1 Introduction ................................ 122
6.2 Fractal Ladder Model and Fractional Constitutive Equation ..... 123
6.2.1 Biological Motivation ....................... 124
6.2.2 Stress-Strain Relationships ....... ' ............. 127
6.2.3 Fractal Ladder Model ....................... 129
6.2.4 Recursive Fractal Ladder Models ................ 132
6.2.5 Constitutive Equation ...................... 135
6.3 Derivation of the Caputo—Wismer Equation ............... 137
6.3.1 Homogeneous Media ....................... 137
6.3.2 Inhomogeneous Medium ..................... 139
6.3.3 Simpliﬁed Equations ....................... 141
6.4 Qualitative Properties of the Caputo—Wismer Equation ........ 141
6.4.1 Low Frequency Limit ....................... 141
6.4.2 High Frequency Limit and Causality .............. 142
6.4.3 Relationship to the Szabo Wave Equation ........... 144
6.5 Discussion ................................. 144
6.5.1 Bio-Mechanical Interpretation .................. 144
6.5.2 Fractal Geometry Interpretation ................. 146
6.5.3 Stochastic Interpretation ..................... 148
6.5.4 Nonlinear Media ......................... 149
7 Lossy Impulse Response in Power Law Dispersive Media ..... 150
7.1 Introduction ................................ 150
7.2 Green’s Function Decomposition: Power Law Dispersive Media . . . . 152
7.2.1 Power Law Impulse Response .................. 153
7.3 Circular Piston: Nearﬁeld ........................ 153
7.3.1 Numerical Results ........................ 155
7.4 Uniform Circular Piston: Farﬁeld .................... 160
7.4.1 Numerical Results ........................ 160
7.5 Rectangular Aperture ........................... 162
7.5.1 Numerical Results ........................ 164
7.6 Spherical Shell ............................... 164
7.6.1 Numerical Results ........................ 166
7.7 Discussion ................................. 166
7.7.1 Stationary Approximation .................... 166
7.7.2 Physical Interpretation ...................... 169
7.7.3 Complex Apertures and B-Mode Imaging ............ 170
8 Conclusion ................................. 171
8.1 Summary ................................. 171
8.2 Relationships between FPDE Models .................. 175
APPENDICES ................................ 176
A Hactional Derivatives .......................... 176
Al Riemann-Liouville Fractional Derivatives ................ 176
A2 Fractional Laplacian ........................... 177
A3 Skew Fractional Laplacian ........................ 177
B Stable Distributions ........................... 179
8.1 Deﬁnitions ................................. 179
B.l.1 Maximally Skewed Stable Distributions ............. 180
8.1.2 Symmetric Stable Distributions ................. 181
8.2 Cumulative Distribution Functions (CDFs) ............... 181
C Special Functions ............................. 182
CI The Fox-H Function ........................... 182
C2 The Wright Function ........................... 183
1.1
2.1
2.2
3.1
3.2
3.3
LIST OF FIGURES
Measured attenuation coefﬁcients for kidney, liver, and tendon taken
from [1]. The data was ﬁt to Eq. (1.1), yielding power-law exponents
of y = 1.09, 1.14, and 0.76 for kidney, liver, and tendon, respectively.
Comparison of the asymptotic form of the Green’s function for the
Stokes wave equation in Eq. (2.20) and the reference Green’s function
computed numerically via the MIRF approach. The Green’s function
is shown for a) R = 0.1 mm and 'y = 0.01 us, b) R = 0.1 mm and
7 = 0.001 as, c) R = 1.0 mm and '7 = 0.01 ps, and d) R = 1.0 mm
and 'y = 0.001 as ..............................
RMS error for the asymptotic Green’s function relative to the exact
MIRF result. RMS error is displayed on a log-log plot for a range
of viscous relaxation times '7 at four different observation distances:
R = 0.01 mm, R = 0.1 mm, R = 1 mm, and R = 10 mm. The error
increases with respect to 'y and decreases with respect to R. .....
Piston and coordinate geometry. The piston, centered at the origin and
surrounded by an inﬁnite rigid baffle in the z = 0 plane, has radius a
and radial coordinate 0'. The radial and axial observation coordinates
are denoted by a: and z, respectively. The distance between origin and
observer is given by 1' = V3132 + y:22 and the angle between the radial
and axial coordinates is 6. ........................
On-axis (a: = 0 mm) velocity potential for a circular piston of radius
a = 10 mm in a viscous medium. The piston is excited by a Hanning-
weighted toneburst with center frequency f0 = 6.0 MHz and pulse
length W = 0.5 us for 4 combinations of axial distance 2 and relaxation
time '7: a) z = 0.1mm and 'y = 0.01ps, b) 2 = 0.1mm and ’7 = 0.001
ps, c) z = 1 mm and 'y = 0.01 as, d) z = 1 mm and 7 = 0.001 ps.
The velocity potential for each combination of z and ’y is computed via
both the lossy impulse response approach and a MIRF approach for
veriﬁcation. ................................
Snapshots of the lossless impulse response generated by a circular pis-
ton with radius a = 10 mm at t = 22 ps and t = 34 ,us are displayed
in panels a) and b). The constant amplitude component within the
paraxial region lml < (1 represents the direct wave. The remaining
component represents the edge wave generated by the discontinuity in
particle velocity at cc = a. ........................
23
24
46
3.4
3.5
3.6
3.7
3.8
3.9
3.10
Snapshots of the lossy impulse response generated by a circular piston
with radius a = 10 mm with 7 = 0.01/1.8 at fort = 22 ps and t = 34 us
are displayed in panels a) and b). Unlike the lossless impulse response
depicted in Fig. 3.3, the direct wave is attenuated due to viscous dif-
fusion. The edge wave also experiences additional attenuation relative
to Fig. 3.3. ................................
Normalized velocity potential ﬁeld produced by a circular piston of ra-
dius a = 10 mm excited by a Hanning-weighted toneburst in a viscous
medium with relaxation time 7 = 0.001 ps. Snapshots of the normal-
ized velocity potential for t = 22 us and t = 34 us are displayed in
panels a) and b), respectively. ......................
Lossy impulse response for a circular piston (a = 5 mm) and 7 = 0.001
,us. In panel a), the nearﬁeld impulse response and the farﬁeld impulse
response were evaluated at r = 50 mm both on—axis (0 = 0) and off—
axis (0 = 1r/ 12 and d) = 0). Panel b) shows the nearﬁeld and farfield
impulse responses evaluated at r = 200 mm on-axis (6 = 0) and off-axis
(6 = 7r/12 and a = 0). ..........................
Coordinate axis used in the derivation of the impulse response for a
rectangular source. A rectangular piston of half-width a and half-
height b is excited by a uniform particle velocity. The radiator is sur-
rounded by an inﬁnite rigid baffle in the z = 0 plane. Reprinted from
[2]. .....................................
Lossy nearﬁeld impulse response for a rectangular piston of half-width
a = 2.5 mm and half-height b = 10 mm evaluated in lossy media for 7
= 0.001, 0.0001, and 0.00001 ps. The impulse response is evaluated at
a)z=100mmandb) z=200mm ....................
Point spread function for a focused, 64-element linear array with half-
width a = 0.3 mm, half-height b = 2.4 mm, kerf of 0.15 mm, focused
on—axis at 80 mm. The point-spread function is measured in the focal
plane at 241 lateral locations in a) lossless media and b) viscous media
with 7 = 0.0001 us .............................
Normalized pulse echo envelope for a focused, 64-element linear array
with half-width a = 0.3 mm, half-height b = 2.4 mm, kerf of 0.15 mm,
focused on-axis at 80 mm. The point-spread function is calculated at
the focal point in lossless media and viscous media with three different
relaxation times. .............................
48
49
50
51
51
52
3.11
3.12
3.13
3.14
4.1
4.2
4.3
4.4
A spherical shell with radius a and radius of curvature R, surrounded
by an inﬁnite rigid baffle. The origin is located at the geometric focus
of the shell, the radial coordinate denoted by 7‘ = V :32 + y2, and the
axial coordinate is indicated by z .....................
Lossy impulse response generated by a bafﬂed spherical shell with ra-
dius a = 10 mm and radius of curvature R = 70 mm in a viscous
medium with viscous relaxation time 7 = 0.001 as. Eq. (3.35) is
evaluated on an 80 by 80 grid at times a) t = 38 ps and b) t = 46 ps.
Velocity potential generated by a spherical shell with radius a = 8 mm
and radius of curvature R = 50 mm in a lossless, water-like medium,
and a lossy, liver-like medium. Velocity potentials are displayed on
axisata) z=—20mmandb)z=0mm. ...............
Velocity potential generated by a spherical shell with radius a = 8 mm
and radius of curvature R = 50 mm in a lossless, water-like medium,
and a lossy, liver-like medium. Velocity potentials are displayed off-axis
atz=0mmanda):r=4mmandb):v=8mm. ...........
Plots of stable PDFs fy(t) for four sub-linear attenuation cases: y =
1/3, 1/2, 2/3, and 9/10. Eqns. (4.44) are evaluated for y = 1/3,
1/2, and 2/ 3, while the y = 9/10 case is evaluated with the STABLE
toolbox [3]. ................................
Plots of stable PDFs fy(t) for four super-linear attenuation cases: 3; =
11/10, 3/2, 19/10, and 2. Eqns. (4.50) are evaluated for y = 3/2 and 2,
while the y = 11/10 and 19/ 10 cases are evaluated with the STABLE
toolbox [3]. ................................
Comparison of the 3D Green’s function using a) the MIRF approach
and b) the analytical Green’s function approach given by Eq. (4.38).
Green’s functions are displayed at three different depths with y = 1.5,
co = 0.15 cm/ps, and a0 = 0.1151 Np/MHz1'5/cm. ..........
Snapshots of the 3D power law Green’s function for y = 0.5, 1.5, and
2.0 for em = 0.05 mm’lMHz_y. Snapshots of the Green’s function are
shown for t = 20, 30, 40, and 50 us ....................
xii
54
55
56
57
76
87
88
4.5
4.6
5.1
5.2
5.3
5.4
5.5
5.6
Velocity potential produced by a point source excited by a broad-
band pulse deﬁned by Eq. (4.68) in two power-law media: 1) a
fat-like medium with y = 1.5, CO = 1.432 mm/ps, and a0 = 0.086
Np/MHzl's/cm and 2) a liver-like medium with y = 1.139, c0 = 1.569
mm/ps, and a0 = 0.0459 Np/MHz1'139/cm. The velocity potential is
displayed at radial distances 10 mm, 25 mm, 50 mm, and 100 mm. . .
The 2D power law Green’s function given by Eq. (4.59) and lossless
Green’s function given by Eq. (4.62) at a) p = 10 mm and b) p = 50
mm. In each panel, (10 = 0.05 mm'lMHz—y for the sub-linear case
y = 2/ 3 and the super-linear case y = 3/ 2. ...............
Attenuation coefficient and phase velocity associated with the Chen-
Holm wave equation, the spatially-dispersive wave equation, and the
power law model for breast fat (y = 1.5). ................
Attenuation coefﬁcient and phase velocitiy associated with the Chen-
Holm wave equation, the spatially-dispersive wave equation, and the
power law model for liver (y = 1.139) ...................
Plots of symmetric stable PDFs wy(t) for four cases: y = 0.5 , 1
(Cauchy) , 1.5 (Holtsmark), and 2 (Gaussian) evaluated with the STA-
BLE toolbox [3]. .............................
Snapshots of the 3D approximate Chen-Holm Green’s function in Eq.
(5.27) for y = 1.0, 1.5, and 2.0 for do = 0.05 Np/mm/Msz. Snapshots
of the Green’s function are shown for t = 20, 30, 40, and 50 ps.
Velocity potential produced by a point source excited by a broadband
pulse deﬁned by Eq. (3.11) in two media governed by the Chen-Holm
equation: 1) a fat-like medium with y = 1.5, CO = 1.432 mm/ps, and
a0 = 0.086 Np/MHzl's/cm and 2) a liver-like medium with y = 1.139,
co = 1.569 mm/ps, and a0 = 0.0459 Np/MHz1'139/cm. The velocity
potential is displayed at radial distances 10 mm, 25 mm,’ 50 mm, and
100 mm ...................................
Comparison between the spatially dispersive Green’s function given
by Eq. (5.47) and the power law Green’s function given by Eq.
(5.49) in a fat-like medium (3; = 1.5) with attenuation coefﬁcient 0.086
Np/cm/MHzl'5. The two Green’s functions are evaluated as a function
of time t at observation points a) R = 1 mm and b) R = 10 mm.
xiii
89
90
97
98
106
114
115
120
5.7
6.1
6.2
6.3
6.4
6.5
Comparison between the spatially dispersive Green’s function given by
Eq. (5.47) and the power law Green’s function given by Eq. (5.49)
in a liver—like medium (y = 1.139) with attenuation coefﬁcient 0.0459
Np/cm/MHzl'139. The two Green’s functions are evaluated as a func-
tion of time t at observation points a) R = 1 mm and b) R = 10
Schematic showing tissue structure at three different spatial scales (tis-
sue, cellular, and sub-cellular). The ﬁrst panel (A ~ 200 pm) dis-
plays an ensemble of mammalian cells, each bounded by an elastic
membrane, shown suspended in viscoelastic ECM. The second panel
(A/ 10 ~ 20 pm) displays an individual cell at a higher level of mag-
niﬁcation. The third panel (A/ 40 ~ 5 pm) displays the cell nucleus,
consisting of a double membrane, a fluid—like nucleoplasm, and an elas-
tic nucleolus in the interior, which contains chromatin [4]. Although
the speciﬁc biological structures vary at each successive spatial scale,
the essential features are the same: ﬂuid substrates containing elastic
compartments. This self—similar pattern forms the basis for the fractal
structure shown in Fig. 6.2. .......................
Layered fractal model for biological tissue based on the schematic
shown in Fig. 6.1. The ﬁrst panel displays an inﬁnite number of
thin elastic membranes with Young’s modulus E alternating with vis-
cous compartments that have coefficients of viscosity 17. The second
panel zooms in on the ﬁrst panel, thus showing the self-similar layered
structure. .................................
Fractal ladder model for tissue micro-structure. The continuum model
depicted in Fig. 6.2 is described with springs with Young’s modulus
E and coefﬁcients of viscosity 1) ......................
Recursive fractal ladder model for tissue micro—structure which gen-
eralizes the ladder shown in Fig. 6.3. The dashpots in the simple
ladder model are replaced with springpots with a modulus K 37, where
each of the springpots are built from recursive arrangements of ladders,
springs, and dashpots. ..........................
Spring-dashpot network constructed in a Sierpinski gasket network.
This fractal exhibits self-similarity at ﬁve levels of magniﬁcation. Each
121
126
130
131
node corresponds to a dashpot, while each edge corresponds to a spring. 148
xiv
7.1
7.2
7.3
7.4
7.5
Comparison of the power law impulse response with the Field II
model. A baffled circular piston with radius a = 10 mm radiates
in linear attenuation media (y = 1) with an attenuation slope 0.5
dB/MHz/cm. Impulse responses are shown at a) (:r,y, z) = (0,0, 50)
mm, b) (m,y,z) = (10,0, 50) mm, c) (:L‘,y,z) = (0,0,100) mm, and
(any, z) = (10,0,100) mm. ........................
Evaluation of the power law impulse response generated by a circular
piston with radius a = 10 mm with attenuation slope do = 0.086
Np/mm/Msz, zero-frequency sound speed c0 = 1.5 mm/ps, and
power law exponents y = 1.139, 1.5, and 2. Panel a) shows the on-axis
impulse response for z = 50 mm, while panel b) displays the result for
z = 100 mm. ...............................
Velocity potential generated by a circular piston of radius a = 10 mm
in power law media: 1) a fat-like medium with y = 1.5, CO = 1.432
mm/ps, and a0 = 0.086 Np/MHzl°5/cm and 2) a liver-like medium
with y = 1.139, c0 = 1.569 mm/ps, and a0 = 0.0459 Np/MHz1-139/cm.
The input pulse is deﬁned by Eq. (4.68). The velocity potential is
displayed on-axis at a) z = 50 mm and b) 2 = 100 mm .........
Lossy impulse response for a circular piston (a = 10 mm) in a
power law medium with exponent y = 1.5, attenuation slope 0.086
Np/cm/MHzl'5, and phase velocity c0 = 1.5 mm/ps. In panel a), the
nearﬁeld impulse response and the farﬁeld impulse response are eval-
uated at r = 50 mm both on-axis (0 = O) and off-axis (6 = 7r/12 and
d) = 0). Panel b) shows the nearﬁeld and farﬁeld impulse responses
evaluated at r = 200 mm on-axis (t9 = O) and off-axis (6 = 1r/ 12 and
43:0) ....................................
Lossy nearﬁeld impulse response for a rectangular piston with half—
width a = 2.5 mm and half-height b = 10 mm evaluated in power
law media with attenuation slope 00 = 0.086 Np/mm/Msz, zero-
frequency sound speed c0 = 1.5 mm/ps, and power law exponents y =
1.139, 1.5, and 2. Panel a) shows the on—axis impulse response for z =
50 mm, while panel b) displays the results for z = 100 mm .......
XV
157
159
161
163
165
7.6 Lossy impulse response generated by a baffled spherical shell with ra-
dius a = 10 mm and radius of curvature R = 70 mm in power law me-
dia with attenuation slope ao = 0.086 Np/mm/Msz, zero-frequency
sound speed c0 = 1.5 mm/ps, and power law exponents y = 1.139, 1.5,
and 2. Panel a) shows the on—axis impulse response for z = —20 mm,
while panel b) displays the result for z = -1 mm ............. 167
8.1 Diagram showing the relationships between the spatially dispersive,
Caputo-Wismer, Chen-Helm, Szabo, Stokes, and Blackstock wave
equations studied in this thesis. ..................... 175
xvi
CHAPTER 1
Introduction
Ultrasonic waves in biological media experience frequency-dependent attenuation.
This frequency-dependent attenuation has ramiﬁcations in ultrasonic imaging. For
instance, the absorption of ultrasonic energy in tissue places constraints on the max-
imum depth that can be imaged at a particular frmuency. Since the lateral and
axial resolution are frequency-dependent, the spatial resolution of B-mode systems
is affected due to absorption of ultrasound. In addition, frequency dependent at-
tenuation implies dispersion by the Kramers-Kronig relations [5—7], which produces
a frequency-dependent phase velocity. As most forms of ultrasound imaging uti-
lize a broadband pulse, each frequency component of the propagating pulse travels
at a different phase velocity, which yields depth-dependent pulse distortion. This
combination of frequency-dependent attenuation and dispersion degrades the spa-
tial resolution of B-mode images, hence limiting the diagnostic utility of traditional
ultrasonography.
Although attenuation is a limiting factor in B-mode imaging, the attenuation of
ultrasonic energy can yield additional diagnostic information that is not captured in
traditional B-mode scans. For instance, in the ﬁeld of tissue characterization, the
attenuation coefficient has been correlated with the pathological state of liver [8, 9]
and breast [10] tissue; hence, knowledge of the attenuation coefficient may be used to
discriminate healthy tissue from diseased tissue. Several studies have also proposed
methods to extract the attenuation coefﬁcient of tissue from backseattered ultrasonic
echoes [11] or via acoustic macroscopes [10]. Since the attenuation coefficient of soft
tissue is 1) largely independent of sound-speed and density variations captured in
B—mode and 2) displays a larger differential variation than sound-speed and density,
accurate measurement of the attenuation coefficient may provide additional diagnostic
information for clinicians. For these reasons, a rigorous understanding of both the
effect of frequency-dependent attenuation and the underlying mechanism responsible
for the observed absorption, is desirable.
Extensive measurement of the attenuation coefﬁcient in human and mammalian
tissue in the ultrasonic range has revealed a power law dependence on frequency
[1, 10, 12]. For most tissue, the power law exponent typically ranges between 1
and 1.7 [10, 13], although exponents slightly less than one have also been reported
[11]. Although most measurements conﬁrm that the attenuation coefﬁcient in tissue
is governed by a power law, there is no consensus on the explanation why attenu—
ation should follow a power law. Classical theories for ultrasonic absorption, such
as thermo—viscosity [14] and Biot’s porous media theories [15] predict a frequency-
squared dependence in the low-frequency limit, while classical relaxation predicts an
attenuation coefficient with a resonant peak at the relaxation frequency of the ma-
terial (i.e. anomalous dispersion). However, neither of these behaviors is observed
in biological tissue. Multiple-relaxation mechanism models [16] predict power law
behavior over a narrow frequency band by empirically choosing the proper weights
and relaxation frequencies, yet these fail to explain power law behavior over large
frequency bands. Therefore, the power law behavior of the attenuation coefﬁcient
over a large bandwidth in biological media cannot be explained from an underlying
physical principle at the present time.
Finally, our mathematical understanding of the consequences of power law atten-
uation in time-domain acoustics is lacking. Mathematically, the power law frequency
dependence of the attenuation coefficient cannot be modeled with standard dissipa-
tive partial differential equations with integer-order derivatives such as the Stokes
wave equation [17] or the telegrapher’s equation [18]. Rather, integro—differential
equations, or fractional partial differential equations (FPDEs), have been proposed
to capture this power law frequency dependence. These FPDE include modes by
Szabo [19], Chen-Holm [20], and Caputo-Wismer [21]. Unlike their simpler integer-
ordered counterparts, FPDE cannot be solved by classical methods of mathematical
physics discussed in canonical texts such as [18]. Therefore, analytical solutions for
general power law media, generated either by point-sources or ﬁnite apertures, have
not been reported in the literature.
In light of our incomplete knowledge of the behavior of ultrasound in power law
media, the purpose of this thesis is threefold: l) to understand the analytical structure
of transient ﬁelds in power law media, 2) to provide a possible description of the
physical mechanism responsible for power law attenuation in biological media, and 3)
to develop analytical models for transient, three-dimensional sound beams in power
law media. Although many FPDE models exist for power law media, the Green’s
function solutions to these equations have neither been report or analyzed. Therefore,
this thesis develops Green’s functions for the recently proposed phenomenological
power law models discussed above. In the process, the physical basis for power law
attenuation is explored. By utilizing these Green’s functions, analytical models for
transient sound beams in power law media are derived. These models may then be
applied to the simulation of ultrasonic phased arrays in dispersive media.
1.1 Ultrasound Absorption Mechanisms
To orient the reader, this chapter surveys three basic mechanisms for ultrasonic at-
tenuation: 1) thermo—viscosity, 2) molecular relaxation, and 3) micro-heterogeneity.
1.1.1 Thermoviscous Model
The underlying physics of both viscous diffusion and thermal conduction are well
documented. All real ﬂuids (except superﬂuids) possess an internal friction called
viscosity, which results from the diffusion of momentum between ﬂuid particles with
different velocities [14]. Since momentum diffusion is a localized near velocity gradi-
ents, a viscous ﬂuid is memoryless and only weakly dispersive. As shown in standard
textbooks [14, 22], viscous dissipation of sonic energy results in an attenuation coef-
ﬁcient possessing frequency-squared dependence.
Besides viscous diffusion, thermal conduction contributes to absorption in simple
materials like water or air. Since acoustic compression is not perfectly reversible (or
adiabatic), acoustic energy may be converted into thermal energy via heat conduction.
Molecules in compressed regions have higher temperatures than molecules in rareﬁed
regions, causing higher temperature molecules to diffuse to lower temperature regions,
thus decreasing the potential energy stored in the sound wave. Like viscous diffusion,
thermal conductivity is quadratic with respect to frequency and weakly dispersive.
The combination of viscosity and thermal conduction is known as the thermo-
viscous theory, which has been analyzed extensively in both linear [14] and non-linear
acoustics [23]. Viscous diffusion and thermal conduction play a dominant role in
materials with a simple molecular structure, such a monoatomic gases. For media
with more complicated structure, such as soft tissue, viscosity and thermal conduction
play a reduced role. The viscous model for acoustic loss in linear media, which is
derived by linearizing the well-known Navier-Stokes equations for a Newtonian ﬂuid,
results in the Stokes wave equation, which is discussed [in detail in Chapter 2. For
nonlinear media, various approximate models include the Westervelt equation, the
Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation, and the well-known 1D Burger’s
equation [23].
l.1 .2 Molecular Relaxation
The dual processes of thermal conduction and viscosity transform the organized lon-
gitudinal motion of ultrasonic waves into disorganized translational motion, or heat,
thus increasing the entropy of the medium. In polyatomic molecules, energy may
be transformed into additional degrees of freedom, such as rotational, vibrational, or
chemical modes. These additional degrees of freedom are determined by the molecular
structure of the molecule and the thermodynamic state of the system (e.g. temper-
ature, pressure, etc.). During each compression cycle, there will be a net transfer
of energy from the sound wave to each of these internal degrees of freedom, thereby
attenuating the amplitude of the wave [14].
Each degree of freedom is characterized by a relaxation time associated with the
internal energy state. If period of acoustic oscillation is on the same order as the
relaxation time, translational energy will be transferred into internal energy. Poly-
atomic molecules and mixtures may involve multiple relaxation processes, each with
its own characteristic relaxation time. Theoretically, these relaxation times may be
computed from statistical mechanics via a population function. In practice, however,
these relaxation times are ﬁtted to experimental data. Although relaxation models
may be ﬁtted to observed power law attenuation data, these models do not predict a
power law over a wide band of frequencies. Examples of multiple relaxation models
include the time-convolutional model described in [16] and a nonlinear model based
on the KZK equation [24].
1.1.3 Micro-Heterogeneity and Viscoelasticity
Rayleigh’s theory of sub-wavelength scattering, or diffusive scattering, assumes an
ensemble of scatters with a characteristic diameter. As is well known, diffusive scat-
tering predicts a quartic dependence on frequency [25], which conﬂicts with observed
attenuation and backseatter data. To resolve this conﬂict, scattering by heteroge-
neous structures at multiple spatial scales may be considered. Unlike water and air,
mammalian tissue is heterogeneous on scales as small as 1 pm (micron), which con-
tributes to absorption and scattering, as well as larger structures extending upwards
of an acoustic wavelength. Class O scattering, which occurs on the scale of microns,
is associated with the dissolved macromolecules present in tissue. This class of scat-
tering manifests itself as absorption and dispersion on the macro-scale (~ 1 mm)
that is resolved by ultrasound [25]. Class 1 scattering occurs on the level of cells
and ensembles of cells (~ 10 - 100 ,um), and may be measured via the backscattering
coefﬁcient. To explain the sub-quartic dependence of Class 1 scattering in tissue,
Ref. [26] considered fractal arrangements of cells to explain the observed power law
dependence of the backscatter coefﬁcient.
In addition, tissue possesses both solid-like and ﬂuid-like properties, which also
contributes to absorption. Tissue consists of water, protein, lipids, and minerals
[25] arranged in a hierarchical structure (tissues, cells, organelles, etc.). In addition,
individual tissue is highly heterogeneous and composed of over a hundred distinct
cell types. These tissues consists of aggregates of cells suspended by an ﬂuid-like
extra-cellular matrix (ECM) [4]. The ECM is often modeled as an aqueous solution
of viscoelastic polymers, which possess both solid and ﬂuid-like properties. This
combination of micro-heterogeneity and Viscoelasticity may be responsible for the
observed power law behavior of tissue and is explored in Chapter 6.
1.2 Phenomenological Power Law Models
To describe the power law dependence of the attenuation coefﬁcient, attenuation in
the ultrasonic range is assumed to arise from multiple mechanisms: thermo-viscosity,
molecular relaxation, and micro-heterogeneity. The cumulative effect of these mech-
anisms in most biological tissue may be approximated by a power law [13] of the
form
afw) = aolwly (1-1)
over a band of frequencies, where w is angular frequency, a(w) is the frequency-
dependent attenuation coefﬁcient, on is the attenuation slope, and y is the power
law exponent. Measured attenuation coeﬂicients of soft tissue typically have linear
or greater than linear dependence on frequency [1, 10, 13]. For example, breast fat
has an exponent y = 1.5, while most tissue ranges from y = 1 to y = 1.5. Stratiﬁed
tissue, such as muscle or tendon, may have an exponent slightly less than 1. For
non-biological materials, such as castor oil or other organic solutions, power laws are
also observed, although the exponent y is typically closer to 2 [6]. Power-laws are also
found in geophysics (e.g. dispersion by micro—heterogeneous porous earth [27]) and in
electromagnetism (e.g. Cole-Cole relaxation of polar molecules [28]). In underwater
acoustics applications, shallow water waveguides with sandy sediments have power
law exponents 1.6 S y S 2.0 [29].
Figure 1.1 shows published attenuation coefﬁcients as a function of frequency
[1]. Attenuation coefﬁcients, which include the combined effects of absorption and
micro-heterogeneous scattering, of kidney, liver, and tendon were measured via a
transient thermoelectric technique at frequencies ranging from 0.5 MHz to 7.0 MHz
at 37°C. Using least-squares regression, the data was ﬁtted to Eq. (1.1), yielding
power-law exponents y =1.09, 1.14, and 0.76 respectively. Thus, the kidney and liver
considered has a slightly greater than linear dependence of frequency. To account
for the attenuation coeﬂicient’s power-law dependence exhibited in Figure 1.1, many
phenomenological tissue models have been proposed. These models may be divided
into frequency-domain formulations, which utilize transfer functions and/ or ﬁlters,
and time—domain formulations, which utilize partial differential equations (PDE),
Frequency-domain formulations, based on the Kramer-Kronig relations [5] be-
tween attenuation and dispersion, have been proposed for linear frequency depen-
U1
X measured data (Kidney)
—power-law fit (Kidney)
E
9 ;
g 4-. + measured data (Liver) ........... ....................... _
8' - - - power-law fit (Liver) :
5 *- measured data (Tendon) ,
Tue; 3 .. III-l power—law fit (Tendon) ....... ...................... .7
.6 . “"‘
E \,\¢‘*
0) 1 .m ‘
8 2 ......................................... ....... ;‘.;\.‘.O‘. ............................. _
c *i“
.9 ‘9‘“ .
a ‘.\*
3 1 (e
C
G)
3:
(U
O
4
frequency (MHz)
Figure 1.1. Measured attenuation coefﬁcients for kidney, liver, and tendon taken from [1].
The data was ﬁt to Eq. (1.1), yielding power-law exponents of y = 1.09, 1.14, and 0.76 for
kidney, liver, and tendon, respectively.
dence [30] and power—law dependence [31—35]. These models evaluate a frequency-
dependent attenuation coefﬁcient and phase velocity over a bandwidth of frequencies,
and then invert the spectra into the time-domain, usually with the aid of fast Fourier
transforms (FFTs). Although these are ﬂexible, frequency-domain methods provide
neither analytical solutions nor physical insight into the problem.
Dispersion has been incorporated within time-domain models via FPDE, which
add loss to the wave equation with a time-fractional derivative [19, 27, 36] or a space-
fractional derivative [20, 37]. These fractional operators are deﬁned in Appendix A.
Unlike loss models which utilize integer-order derivatives [38, 39], these FPDE support
power-law attenuation coefficients for a range of non-integer exponents y. These
equations are typically solved numerically via ﬁnite-difference time-domain (FDTD)
[40, 41] or ﬁnite element method (FEM) [21] analysis. All of FPDE considered in this
thesis possess the form
1 62p
V2 —-—-—
1) 03 6t2
+ L(p) = 0 (1.2)
where p is acoustic pressure, CO is a reference speed of sound, and L(p) is a linear op-
erator involving time-fractional and/ or space-fractional derivatives. Following Szabo
[19], L(p) is termed a loss operator. The three major FPDE explored in this thesis
are the 1) power law and Szabo wave equations, 2) the Chen-Holm and spatially dis-
persive wave equations, and 3) the Caputo—Wismer equation, each which is discussed
below.
1.2.1 Power Law and Szabo Wave Equation
In particular, the Szabo equation [19] interpolates between the telegrapher’s equation,
which approximates frequency-independent attenuation, and the Blackstock equation
[42], which approximates frequency-squared dissipation, allowing a large range of
lossy behavior to be encapsulated within one equation with two loss parameters:
an attenuation slope a0 and a power-law exponent y. The Szabo wave equation was
originally presented as an integro-differential equation for fractional power law media.
Subsequently, Ref. [43] expressed the Szabo wave equation in terms of fractional
derivatives, thereby allowing the machinery of fractional calculus to be utilized. In
Chapter 4, the Szabo wave equation, and a generalization, called a power law wave
equation, are considered in detail. In particular, Green’s functions to the power law
wave equation are constructed in homogeneous, 3D media.
1.2.2 Chen-Holm and Spatially Dispersive Wave Equations
Recently, alternate FPDE models for power law attenuation have been reported,
including the model by Chen-Holm [20]. The Chen-Holm model utilizes a space-
fractional derivative of order y to support power law attenuation with exponent y in
the zero-frequency limit, which is desirable for biomedical applications. The Chen-
Holm models build on the theory of fractional diffusion [44] and random walk models
[45] by utilizing a fractional Laplacian operator to model absorption. Like the Szabo
equation, the Chen-Holm interpolates between frequency-independent attenuation
(y = 0) and frequency-squared attenuation (y = 2). The Szabo wave equation utilizes
higher-order temporal derivatives of order y + 1 > 2, which require three initial
conditions: 1) p(r,0), 2) p(r, 0), and 3) p'(r,0). Since the third initial condition may
not be available, Chen and Holm incorporated dispersive loss via a symmetric, space
fractional derivative [20] which only requires the ﬁrst two initial conditions. To explore
the ramiﬁcations of the space-fractional derivative, the Chen-Holm equation and a
natural generalization, called the spatially dispersive wave equation, are explored in
Chapter 5.
1.2.3 Fractal Ladder Models and the Caputo—Wismer Wave Equation
Chapters 4 and 5 concentrate on an analytical description of power law attenuation
using the probabilistic tools of stable distributions as solutions. Can the underlying
physics of ultrasonic absorption be described geometrically? Since the pioneering
work of Mandelbrot [46], the concepts of fractal self-similarity have provided an inge-
nious framework for modeling both the physical and the natural world. For instance,
fractals are used to model the branching of air passages in the lung and to model
the folding of cell membranes [47]. Mandelbrot noted the close connection between
fractional calculus, stable distributions, and fractals in his seminal monograph [46].
Therefore, fractals may provide a geometric description of the underlying physics. To
make the connection between fractional calculus and fractals, self-similar networks of
springs and dashpots are introduced to model the dissipative process. Using these
networks, a recently proposed FPDE by Caputo-Wismer [21, 36] is derived. The
Caputo-Wismer FPDE utilizes a fractional time-operator of order y — 1, which yields
a power law attenuation coefficient in the low frequency limit, thus providing an ab
10
initio explanation for the observed power law dependence.
1.3 Fourier Transform Convention
Fourier transforms, and, to a lesser extent, Laplace transforms, are used extensively
throughout this thesis. To maintain consistency with the bulk of the literature de-
voted to acoustic dispersion, the i-convention for \/—1 is used. Temporal Fourier
transforms are denoted by a hat and deﬁned via the following forward/ inverse pair:
A m -
1,0(w) = / 1,!)(t)eu"t dt (1.3a)
—OO
1 °° . -
at) = — / wee-WW (Lab)
27r _00
Spatial Fourier transforms are denoted by an upper-case letter and deﬁned via
in.) = wee-1"” dnr (1.4a)
Rn
_ 1 ik~r
we) _ —(2,,)n [Rn ‘Il(k)e d"k (1.4b)
where 7?." denotes n-dimensional Euclidean space. In the bulk of this thesis, Eq.
(1.4) is specialized to n = 3.
1.4 Content of Thesis
This thesis consists of three major sections. Chapter 2 and 3 concentrate on the
Special case of time-domain diffraction in viscous media, which is easier to treat
mathematically. Chapters 4, 5, and 6 form the core of the thesis, wherein the three
FPDE models discussed in Section 1.2 are explored in depth. Finally, the tools
developed in Chapters 4, 5, and 6 are applied to time-domain diffraction in Chapter
7, thus generalizing the results of Chapters 2 and 3.
Since viscous media is a special case of power law media, insight into general dis-
sipative processes is gained by studying the special case of viscous media. Therefore,
11
the transient theory of viscous media governed by Stokes wave equation, including an
asymptotic Green’s function, is developed in Chapter 2. To understand the behavior
of transient ﬁelds produced by ﬁnite apertures, the classical impulse response theory
[48] is extend to viscous media in Chapter 3. Impulse responses of circular, rectan-
gular, and spherical shell radiators in viscous media are then computed, allowing the
coupled effects’of diffraction and viscous loss to be studied in the time-domain.
Chapters 4 and 5 concentrate on the analytical structure of power law transient
ﬁelds. In order to extend the lossy impulse response to the general power law case,
Green’s functions need to be calculated. In Chapter 4, a power law FPDE is derived
and 3D Green’s functions in power law media are calculated in homogeneous media
use the machinery of stable distributions [49]. Chapter 5 introduces alternate power
law models based on fractional spatial derivatives. Green’s functions to these space
fractional FPDEs are derived and compared to the time-fractional models studied
in Chapter 4. The FPDEs considered in Chapters 4 and 5, unlike the Stokes wave
equation, are phenomenological and not based on an underlying physical principle.
To shed light on the physical mechanisms responsible for power law attenuation, a
fractal spring-dashpot model in considered in Chapter 6. Using this fractal model,
a constitutive equation is deduced, allowing a previously proposed FPDE, known
as the Caputo-Wismer equation [21, 36], to be derived in both homogeneous and
inhomogeneous media.
Chapter 7 generalizes the lossy impulse response developed in Chapter 3 to lossy
media using the Green’s function theory developed in Chapters 4-6. Impulse responses
of circular, rectangular, and spherical shell radiators in power law media are computed
in terms of stable distributions and compared to their viscous counterparts in Chapter
3. Chapter 8 summarizes the original results of this thesis.
l2
CHAPTER 2
Stokes Wave Equation
2.1 Introduction
Viscous loss and thermal conduction, which possess a quadratic dependence on fre-
quency, are modeled by the Stokes wave equation [14]. As discussed in Chapter 1,
viscosity and thermal conduction are the dominant loss mechanisms in simple ﬂuids
such as monoatomic gases and water. For more complicated polyatomic materi-
als, such as polymers, or for materials with heterogeneous micro-structure, such as
tissue, viscosity and thermal conductivity do not accurately describe the observed
attenuation coefficient. Although viscosity fails to describe the mechanism for loss
in biological media, the Stokes wave equation is considered in this chapter for the
following reasons: 1) the basic physics are well understood, allowing a rigorous anal-
ysis, 2) an approximate Green’s function is available in closed-form, and 3) insight
into general dissipative media is gained by studying the approximations available in
viscous media, which is a special case of power law media.
Time-domain studies of wave propagation in viscous and other lossy media have
been primarily limited to one-dimensional plane wave propagation. For example,
the classic study by Blackstock [42] approximated the viscous term in the Stokes
wave equation by a third-derivative, leading to both a noncausal PDE and approx-
imate solution. Wave prOpagation in viscous media has also been modeled by the
13
telegrapher’s equation [50], leading to causal solutions with sharp wavefronts similar
to electromagnetic propagation in conductive media. Recent studies derived causal
Green’s function solutions to the Stokes wave equation in the form of contour integrals
[51], inﬁnite series [52], and closed-form approximations [17]. Although the combined
effects of diffraction and loss have been thoroughly studied in the frequency domain
both analytically [53] and numerically [54], most time-domain studies have been ana-
lytical. For example, the combined effects of diffraction, dissipation, and nonlinearity
were studied numerically in Ref. [55] under the parabolic approximation.
In this chapter, an approximate Green’s function is derived as a shifted and scaled
Gaussian distribution. This Green’s function is numerically compared to a reference
solution and shown to provide an excellent approximation for 1) observation distances
sufﬁciently far from the source and 2) moderate values of the attenuation coefﬁcient.
The Green’s function is then decomposed into diffraction and loss components, where
the diffraction component satisﬁes the lossless wave equation and the loss component
satisﬁes the diffusion equation.
2.2 Formulation
This chapter considers the Stokes wave equation
1 62p
2 a 2 _
V p+7—V p—gw—
6t 0 (2.1)
where Co is the thermodynamic speed of sound and 7 is the relaxation time of the
medium. The relaxation time is related to the coefficient of shear viscosity p by the
following relation
_ a
’7 —' :3ng (2'2)
where p is the density of the homogeneous medium. Physically, 7 measures the time
necessary to restore equilibrium to the translational degrees of freedom in the ﬂuid
following a small thermodynamic disturbance. The velocity potential *(r, t) and the
14
particle velocity u(r, t) also satisfy Eq. (2.1). If 7 = 0, equilibrium is restored im-
mediately, and Eq. (2.1) reduces to the lossless wave equation. The loss operator in
the Stokes wave equation, given by the third term in Eq. (2.1), is a singular pertur-
bation which transforms the lossless, second-order hyperbolic wave equation into a
third-order parabolic equation [56]. The dispersion relation between wavenumber k
and angular frequency w are computed via a space-time Fourier transform applied to
Eq. (2.1), which yields (see Eqns. (58) and (59) in Ref. [52]):
kw = w y/l—iw 2.3
() co ’_—__51+(w) 7 ( )
thereby yielding an attenuation coefﬁcient that is approximately proportional to w2
for w7 << 1. The effects of thermal conductivity may also be included by modifying
the constant 7 [57].
2.3 Green’s Function Solution
Unfortunately, no closed-form analytical Green’s function is known to exist for the
Stokes waves equation. However, approximate Green’s functions may be derived,
which are useful for analytical evaluations. The ﬁrst approach, developed by Black-
stock [42], approximates the viscous term in Eq. (2.1) by a third-order time-derivative.
This approximate Blackstock equation is then solved via frequency-domain methods,
yielding a non-causal solution. A second approach works directly with the Stokes
wave equation and derive a causal approximation [17, 56].
In this section, a causal Green’s function to the Stokes wave equation is derived
as an inﬁnite series. 1 The ﬁrst term in this series reduces to a previously derived
solution (see [56] or Appendix A in Ref. [59]). Higher order terms may then be
evaluated in terms of Hermite polynomials if necessary.
1Similar analytical techniques are utilized in the solution of fractional ordinary differential equa-
tions (FODE) with constant coefﬁcients [58]. This series expansion technique will be generalized to
FPDE in Chapter 5.
15
The transient Green’s function to the Stokes wave equation g(r, t) satisﬁes
2 “929 78v29= —5 6R
Vg—goa—tg‘l' ”YE (t)( ) (2-4)
where R = [r — r'| is the distance between observer r and source r’. Fourier trans-
forming Eq. (2.4) from the space-time domain (R, t) to the spectral-frequency domain
(k, 3) via Eq. (1.4a) and the Laplace transform produces
2 32
—k ———7sk2 G=—1 (2.5)
Co
where k = [k] is the magnitude of the spatial frequency vector k. Solving for G and
completing the square yields
63
c“; k, =
( S) (s + 7c8/2k2)2 + c3192 — 72/41:4
(2.6)
Rewriting Eq. (2.6) yields
6(2) 72/463194 —1
C(k’ 3): (s + 7c0/2k2)2 + ogre? 1 — (s + 7/2k2)2 + c3192 (2'7)
Expanding the second factor as a geometric series yields
00
s) = Z Gn(k, s) (2.8)
n=0
where
c 2k2 2n 1
G u(k, s)— _ cg 7 0 M (2.9)
2 [(3 + 7031.2 /2)2 + ogre]
The inverse Laplace transform of Eq. (2.9) is computed via [60]
1 _
—l _ n+1
.C { [(s + (1)2 + A2]n+l}— u(t)2nn :t A ”jn (tA) (2.10)
where jn(z) is the spherical Bessel function of the ﬁrst kind of order n and u(t) is the
Heaviside step function. Letting A = cok and a = 7c(2,k2 / 2 yields
71634-371163?! ')’C2th
72 tn+l exp (_7C0
Gn(k, t) = u(t)
8%!
> 971030,“) (2:11)
16
The 3D Green’s function is recovered via a three-fold inverse Fourier transform deﬁned
by Eq. (1.4b) in spherical coordinates (k, 9),, (bk), yielding
oo 11' 21r
g(R, t) = —l—/ f / G(k, t) exp (ikR cos 6k) sin 6kk2 dcpk dﬂk dk (2.12)
8W3 o o 0
where 9k = cos—1(kz / k) and ¢k = tan-1(ky/kx). The (bk and the 6k integrals are
evaluated, producing
1
27r2R
Inserting Eq. (2.8) into Eq. (2.13) yields
g(R, t) = [000 ksin (kR) G(k, t) dlc (2.13)
00
g(RJ) = 2912020, (2-14)
n=0
where
7211 2+3n n+1 00 3n+l 7%leth
Rt = t——— t k —- ° kR ° kt all:
at ,> u<>8.,,.,,2Rc0 /0 exp sm< >2n> 1, the incoming wave is
negligible, yielding
u _ 2
g(R, t) z 41r—R\(/t—_2)7rt—7 exp (—(i—§R%@—) (2.20)
Eq. (2.20), which is valid for small viscosity (7 -+ 0) and large times (t -—2 00),
supports inﬁnite propagation speed, yet is strictly causal. As the following section
shows, Eq. (2.20) is an excellent, causal approximation except in the extreme nearﬁeld
or highly viscous media.
2.3.2 First Order Term (n = 1)
Higher order terms in the series expansion Eq. (2.14) may be computed if necessary.
In this subsection, the term 91 (R, t) is computed explicitly in terms of Hermite poly-
nomials. If higher-order terms are required, the technique presented below may be
applied. To evaluate Eq. (2.15) for n = 1, an identity is ﬁrst required.
First, let
_ exp (—a:2/(4a))
_ 2.
which may also be expressed as an inverse cosine transform
1 0° _ak2
72(3) = '7; e cos(k:r) dk (2.22)
0
Differentiating both sides of Eq. (2.22) m times and appealing to the deﬁnition of
the Hermite polynomial Hn(:z:)
Hna) = (—1)"e$ — (ta-$2) (2.23)
yields the identity
1 0" —ak2 m mrr (-1)m :1: —:r2
- k —) dk: H __ _
”/0- e k cos(:z:+ 2 2m+1\/7—r(\/a)m+l m 2J5 exp 4a
(2.24)
To evaluate gl(R, t), we need
, sinz cosz
31(Z)= 2 — (2.25)
z z
Inserting Eq. (2.25) into Eq. (2.15) with n = 1, applying the sine and cosine addition
formulae, and utilizing Eq. (2.24) with a = 7cgt/ 2 and :c = cot :l: R yields
(7/ (2600214
4m 0 [91a(R, t) + 91:,(R, 0] (2.26)
91(R’ t) :
with A = (Mt/2, where
—[(—)(—)(%e)w2(t—%e)l
(2.27 )
W = —— [H3 W, w (as) «3 (ex—co) w as)
(2.27b)
111200 = % eXP (-;\ (2.27c)
Higher-order terms gn(R, t) for n 2 2 may be calculated by applying the expansions
of jn(z) and applying the same procedure.
2.4 Green’s Function Evaluation
The error incurred by utilizing the asymptotic Green’s function given by Eq. (2.20)
is quantitatively analyzed in this section. A reference solution to the Stokes wave
equation is evaluated using the material impulse response function (MIRF) approach
outlined in Ref. [31], and the result is compared to Eq. (2.20). The wavenumber
k(w) is calculated via the reference dispersion relation given by Eq. (2.3), allowing
19
the material transfer function (MTF) to be calculated in the frequency domain. As
speciﬁed by Eq. (10) in Ref. [31], the reference Green’s function, or MIRF, is then
recovered by inverse Fourier transforming the MTF, or frequency-domain Green’s
function, via the fast Fourier transform (FFT)
g(R t) = f—1 (8176003) ’ (2.28)
’ 41rR
Since g(R, t) is necessarily real, conjugate symmetry of the MTF is enforced in the
frequency domain. The relative L2 error, deﬁned via
00 _ re 2d
error: \/f-w|¢(t) (1) fan t (2.29)
W320 new)? dt
is calculated with the MIRF approach used as the reference.
Figure 2.1 summarizes the results of this comparison for various combinations of
the observation distance R and the viscous relaxation time 7. Large values of 7 were
chosen to determine the range of acceptable values for this approximation. Panel a),
which shows the numerical reference and approximate Green’s functions using R = 1
mm and 7 = 0.01ps, displays a signiﬁcant disparity between the reference and asymp-
totic Green’s functions when 7 is relatively large and R is relatively small. Panel b),
which is evaluated with R = 1 mm and 7 = 0.01ps, shows a much smaller difference
between the asymptotic and reference Green’s functions. Panels c) and d) show the
Green’s functions for R = 1 mm and 7 = 0.01ps and 7 = 0.001ps, respectively; these
panels show very close agreement between the reference and asymptotic Green’s func-
tions. The asymptotic and reference Green’s functions are virtually identical in panel
(I) of Fig. 2.1, implying that Eq. (2.20) is an excellent approximation for R > 1 mm
and 7 < .001ps.
A quantitative error analysis comparing the asymptotic Green’s function to the
MIRF result is displayed in Figure 2.2. Error as a function of relaxation time 7 is
displayed on a log-log plot for four different observation distances: R = 0.01 mm,
20
R = 0.1 mm, R = 1 mm, and R = 10 mm. As expected, the error increases as
7 increases and decreases as R increases. For a 1% error threshold and a given
observation distance R in mm, a maidmum acceptable relaxation time is determined
from Fig. 2.2 by ﬁxing R and identifying the range of 7 satisfying the 1% error
threshold, yielding the relation 7 S 2 x 10-4R. By ﬁxing 7 and letting R vary, the
minimum acceptable distance is determined by R Z 5 x 1037.
The viscous relaxation time is related to the attenuation coeﬂicient a at a ﬁxed
frequency f via the relation [57] 7 = COO/(2W2f2). For instance, soft tissue at f = 3
MHz with an attenuation coefﬁcient of 01 = 0.345 cm‘1 and sound speed on = 1.5
mm/ps corresponds to 7 = 3 x 10‘4 us. Applying this error analysis, R must be
chosen larger than 1.5 mm to maintain an error less than 1% in this media.
2.5 Green’s Function Decomposition
Mathematically, Eq. (2.20) consists of a composition of a) the Green’s function to
the 3D wave equation and b) the Green’s function to the 1D diffusion equation. This
observation motivates decomposing Eq. (2.20) into individual diffraction and loss
factors via
g(R,t) =1: [Mt—221260)] [u(t) ﬂimexp 0%)] dt’, (2.30)
where the sifting property of the Dirac delta function applied to Eq. (2.30) yields Eq.
(2.20). Identifying the ﬁrst bracketed expression as a diffraction component and the
second bracketed expression as a loss component via
_ 5(t - R/co)
gD(R, t) — 47rR . (2.318.)
(t t') - u(t)—1—-— ex ——t’—2- (2 31b)
91’ ’ _ (/27rt7 p 2t7 ' '
allows Eq. (2.30) to be expressed as
00
g(R, t) = gut. t’) @9002, t) E / gD(R,t — t’)gt(t,t’) dt’ (2.32)
—00
21
where the convolution denoted by (8 is performed with respect to the t’ variable. In
Eq. (2.32), the subscript “D” refers to diffraction, whereas the subscript “L” refers
to loss. Eq. (2.31a) is responsible for the diffraction of the radiating aperture. Eq.
(2.31b) is a Gaussian function centered at t’ = 0 having width, or standard deviation,
J57. In Eq. (2.31b), the t’ variable is interpreted as “propagation time” (fast),
whereas t is interpreted as the “diffusion time” (slow). In fact, Eq. (2.31b) is the
causal Green’s function for the 1D diffusion equation given by [61]
2
(g - gatig) gut. t’) = 66(t') (233)
where cot’ plays the role of the spatial variable. Thus, Eq. (2.32) is interpreted as
the non-stationary convolution of the 1D Green’s function for the diffusion equation
with the Green’s function of the 3D wave equation. As shown in the next chapter, the
decomposition embodied in Eq. (2.32) facilitates the generalization of the impulse
response for ﬁnite apertures in viscous media.
22
—reference I T . -reference 1
12 ‘ - - -approximation 35 ' - - -approxlmation
1o- 30’
g g 25
a 8’ s
I.L I.L 20
6.
g g 15' | 1
0 4. , (5 l
10* l
2' l 5' J k
c . 4 . ...' ----- c . .
0 0.05 0.1 0.15 0.2 0 0.05 0.1 0.15 0.2
time (microseconds) time (microseconds)
(a) 7 = 0.01 pa and R = 0.1 mm. (b) 7 = 0.001 ps and R = 0.1 mm.
r - ﬂ 1.2 ' ' ' M
—relerence _'°f°'°"_°°
0.35 a - uappmxlmaﬁon " ' 'appi’OXImatlon
1 r .
0.3
5 0.3- l
g 0.25' g
E. 0-2' l l: 0.6r
$0.15. 3
(9 (5 0.4' 1
0.1 ' '
60 0:5 1 1:5 2 60 0:5 i f; 2
time (microseconds) time (microseconds)
(c) 7 = 0.01 is and R = 1.0 mm. (d) 7 = 0.001 us and R = 1.0 mm.
Figure 2.1. Comparison of the asymptotic form of the Green’s function for the Stokes
wave equation in Eq. (2.20) and the reference Green’s function computed numerically via
the MIRF approach. The Green’s function is shown for a) R = 0.1 mm and 7 = 0.01 as, b)
R = 0.1 mm and 7 = 0.001 [1.3, c) R = 1.0 mm and 7 = 0.01 us, and d) R = 1.0 mm and
7 = 0.001 as. '
23
: —R=0.01 mm ;
'---R=0.1 mm '
nun-R =1 mm
RMS error
relaxation time 7 (us)
Figure 2.2. RMS error for the asymptotic Green’s function relative to the exact MIRF
result. RMS error is displayed on a log-log plot for a range of viscous relaxation times 7 at
four diﬂ'erent observation distances: R = 0.01 mm, R = 0.1 m, R = 1 mm, and R = 10
mm. The error increases with respect to 7 and decreases with respect to R.
24
CHAPTER 3
Lossy Impulse Response in Viscous Media
In this chapter, the approximate Green’s function derived in Chapter 2 is applied to
diffraction by ﬁnite apertures in viscous media. To solve these problems, the impulse
response method, which is widely used to solve transient problems in acoustics, is
generalized to viscous media by invoking the Green’s function decomposition given
by Eq. (2.32). The tools developed in this chapter will form the basis for modeling
time-domain diffraction in power law media in Chapter 7.
In order to understand the interplay between diffraction and frequency-dependent
attenuation, expressions that describe time-domain impulse response for simple pis-
ton geometries are needed. Although impulse response expressions exist for circular
[62, 63] and rectangular [48] apertures in lossless, homogeneous media, no analyti-
cal expressions for the transient ﬁeld radiated by ﬁnite planar apertures in a viscous
medium have been published previously. An analytical expression for the on-axis
velocity potential produced by a focused spherical shell was recently derived in Ref.
[64], although ﬁeld points off-axis were not considered. In particular, the effect of
viscous loss on the impulse response of bafﬂed circular and rectangular pistons will
be addressed in this chapter.
25
3.1 Lossy Impulse Response: Theory
3.1.1 Impulse Response in Lossless Media
The impulse response method, developed independently by Oberhettinger [62] and
Chadwick-Tupholme [65], is the de facto standard for computing transient acous-
tic ﬁelds produced by ﬁnite apertures [66]. The (lossless) impulse response theory
assumes a homogeneous, linear acoustic medium with a ﬁnite radiating aperture S
surrounded by an inﬁnite, rigid bafﬂe. The normal velocity is prescribed on the aper—
ture S via uz(r,0) = q(r)v(t) where v(t) is the time-domain input pulse, q(r) is the
apodization which incorporates spatial variation in the surface velocity, and r’ lies in
the plane 2 = 0. Thus, a Neumann boundary condition is imposed via
6(1)
5; 2:0 = q(r)v(t) (3-1)
where **(r,t) is the velocity potential. By solving the lossless 3D acoustic wave
equation via the Helmholtz integral equation, the velocity potential is written
@(r, t) = u(t) =1: h(r, t) (3.2)
where “*” denotes temporal convolution and the impulse response is expressed by the
spatial convolution
h(r, t) = /Sq(r’)6(t ;:;::,Il/CO) dzr'. (3.3)
For simple geometries (circle, rectangle, spherical shell, etc.) and apodization func-
tions q(r’), Eq. (3.3) may be computed in closed form. The velocity potential,
pressure, and transmit-receive pressure may then be computed via fast temporal con—
volutions.
In the remainder of this chapter, a uniform surface velocity across the face of
the aperture S (q(r) = 1) is assumed in order to develop relatively simple analyti-
cal expressions. Since the focus here is understanding the basic features of coupled
26
diffraction and loss, apodized radiators with non-uniform surface velocities (e. g. sim-
ply supported pistons) [67-69] are not considered in this chapter.
3.1.2 Impulse Response in Viscous Media
The impulse response function h L(r, t) in the half space is deﬁned by integrating Eq.
(2.20) over the surface of the radiating aperture and multiplying by two to account
for the inﬁnite, planar baffle in the z = 0 plane, yielding
hL(r, t) = 2/Sg(r — r',t) dr’ (3.4)
For example, Eq. (2.32) allows the impulse response to be analytically evaluated for
ﬁnite apertures. Integrating the decomposed Green’s function over the aperture S
and multiplying by two yields
hL(r,t) = gL(t, t’) ® [S W d2rI. (3.5)
The second factor in Eq. (3.5) is identical to the standard impulse response for a
uniform aperture given by Eq. (3.3), which has been calculated previously for a
circular source. Therefore, the lossy impulse response function is expressed as a non-
stationary convolution of the standard impulse response with a loss function, where
the convolution is taken over t’:
hL(r. t) = gut. t’) s h a, Eq. (3.7) is continuous
29
everywhere, yet not differentiable at both the leading and trailing edges.
3.2.2 Viscous Media
The nearﬁeld solution given by Eq. (3.7) is an exact solution to the transient wave
equation assuming a uniform piston in an inﬁnite rigid bafﬂe. The lossy impulse
response is computed by inserting Eq. (3.7) into Eq. (3.6) and evaluating the unit
step functions in Eq. (3.7). The integration over t’ is evaluated in terms of the error
function [72] erf(z), yielding
60a 7' xcostb—a t-T1 t—Tz
h t = — t f — — f d .
L($’ 2’ ) 27r u( )_/0 x2 + a2 — 2axcosw [er (V2t7) er ((/2t7)] lb
On axis (:1: = O), the lossy impulse response reduces to
_ c0u(t) t-z/c t—V22+a2/co
hL(0,z,t) — 2 [erf( m ) — erf( \/2t_7 )]. (3.10)
As in the lossless case, the ﬁrst and second delays in Eq. (3.10) correspond to the
closest point and the furthest point on the radiating piston, respectively. The “slow”
time scale is embodied in Eqns. (3.9) and (3.10) by the denominator of the erf
functions V277 whereas the “fast” time scale is embodied in the numerator. As
Eqns. (3.9) and (3.10) show, the lossy impulse response possesses an inﬁnitely long
tail that decays like a Gaussian.
3.2.3 Numerical Results
In this section, the lossy impulse response expressions for the circular piston in the
nearﬁeld are numerically evaluated and physically interpreted. The velocity potential
\Il(r, t) resulting from the temporal convolution of a transient excitation with the
lossy impulse response is also computed and discussed. All single-integral expressions
are numerically integrated using Gauss-Legendre quadrature [73]. In the following
30
velocity potential computations, the Hanning—weighted toneburst given by
v(t) = -12— [1 — cos (27rt/W)] rect (57) sin (27rf0t) (3.11)
with center frequency f0 = 6.0 MHz and duration W = 0.5 as is applied, where
rect(t/W) = u(t) — u(t — W).
Since the lossy impulse response solutions utilize an approximate transient Green’s
function to the Stokes wave equation, a numerical comparison is made to a reference
frequency-domain solution. The numerical reference is computed by computing the
Fourier-transform of the velocity potential **(r, w) = 6(w)h(r, k(w)) using the disper-
sion relations k(w) given by Eq. 2.3. This product is computed over the effective
bandwidth of 6(a)) and then numerically inverse Fourier transformed using an inverse
FFT. For a uniform circular piston of radius a, the on—axis (r = 0) transfer function
h(0, z; w) exists in closed-form, thereby obviating the need for numerical integration.
Figure 3.2 displays the on-axis velocity potential for a circular piston of radius»
a = 10 mm in a viscous medium for the same combinations of z and 7 utilized in Fig
2.1. The velocity potential is obtained for a) z = 0.1 mm and 7 = 0.01 ps, b) 2 =
0.1mm and 7 = 0.001ps,c) z = 1mm and 7 = 0.01ps, and d) z = 1mm and 7 =
0.001 ps. The error is computed between the lossy impulse response and the MIRF
method using Eq. (2.29) using the MIRF ﬁeld as reference. The resulting errors for
the four cases are a) 26.7 %, b) 3.26 %, c) 6.89 %, and d) 2.71 %, respectively. Thus,
larger errors are observed for points closer to the piston and large relaxation times
due to the approximation in the Green’s function given by Eq. (2.20). For 2 > 10
mm and 7 5 0.001 ps, the relative L2 error is less than or equal to 1%. Since the
effect of viscous dissipation is negligible in the extreme nearﬁeld, the lossy impulse
response yields an accurate solution that captures the combined effects of diffraction
and viscous dissipation. Finally, Eq. (3.6) accounts for the small differences in the
arrival time of each attenuated spherical wave emitted by the radiating aperture via
the non-stationary convolution.
31
Figures 3.3 and 3.4 display the nearﬁeld impulse response for a circular piston
with radius a = 10 mm for lossless media and viscous media, respectively. Fig. 3.3
shows two snapshots of the lossless impulse response (7 = 0) on a 2D computational
grid extending from cc = —40 mm to :c = 40 mm in the lateral direction and z = 20
to z = 60 mm in the axial direction. Snapshots of the impulse response are shown
at two instances in time: t = 22 ,us and t = 34 us. As shown in Fig. 3.3, the field
within the paraxial region (le S 10 mm) in the lossless case is unattenuated, while the
corresponding region in Fig. 3.4 experiences attenuation and stretching due to viscous
diffusion. For increasing relaxation times, the ﬁeld becomes progressively closer to
that generated by an omni-directional source as the directivity of the aperture is
weakened far from the source.
The nearﬁeld velocity potential ﬁeld generated by a piston with radius a = 10
mm within a viscous medium with relaxation time 7 = .001 us is displayed in Figure
3.5 at two successive snapshots in time. At t =22 ,us, distinct direct and edge waves
are evident, with little apparent viscous spreading. As time increases, the edge and
direct waves becomes less distinct, while spreading and attenuation due to viscous
loss become more pronounced. Unlike the lossless case, the direct wave in Fig. 3.5
experiences a signiﬁcant decrease in amplitude. The physical interpretation of the
uniformly excited circular piston discussed in Ref. [74] is also applicable to transient
propagation in viscous media. Physically, the ﬁrst term in Eq. (3.9) corresponds to
the edge wave generated by the discontinuity at the of the piston, whereas the second
term corresponds to the direct wave due to the bulk motion of the piston. In the
lossless case given by Eq. (3.7), the direct wave contribution localized within the
paraxial region of the radiator maintains a constant amplitude; however, for the lossy
wave, the direct wave decays as 2 increases. As Fig. 3.3 shows, the direct and edge
waves are clearly discernible in the lossless case, while Fig. 3.4 shows the smooth
transition between edge and direct waves in the presence of loss.
32
3.3 Uniform Circular Piston: Farﬁeld
Although Eq. (3.9) is valid for all ﬁeld points, a simpler expression exists in the
farﬁeld, where the effects of attenuation are more pronounced. Assuming the deﬁ-
nition of the farﬁeld utilized by Morse and Ingard [75], the observation distance is
assumed to be much larger than the radius of the aperture (7' >> a) . This deﬁnition
of the farﬁeld is distinguished from the classical farﬁeld distance r > a2/A, which is
only valid over a ﬁnite frequency band.
3.3. 1 Lossless Media
A single-integral representation of the farﬁeld lossless impulse response function is
derived from the well-known frequency-domain expression. In the farﬁeld region,
the impulse response is expressed in spherical coordinates r = ng + 22 and 6 =
tan-1(x/z) via an inverse temporal Fourier transform (Eq. (1.4b)):
A 1 oo eik‘r _
h(r,0,t) = f—1 [h(r,9,w)] = / D(0;w)e_wtdw (3.12)
g _00 21W
where D(6;w) is the far-ﬁeld pattern and k = w/c is the (lossless) wavenumber. Eq.
(3.16) is computed by inverse Fourier transformng the classical frequency-domain
result. The farﬁeld pattern is given by
2J (ka sin 6)
2 1
7m ka sin 0
Damn) = , (3.13)
where J1(z) is the Bessel function of the ﬁrst kind of order 1. To invert the Bessel
function in Eq. (3.12), an integral representation is utilized [76]:
in
Jn(z) = — f” cos ml) 6.” C051!) cit/2. (3.14)
0
7r
Inserting Eqns. (3.14) with n = 1 and (3.13) into Eq. (3.12) yields
eiwr/co
—iw
h(r,6,t) = CO“ r1 [
1r
- —ikasin6costb
7rrsin6 f0 COSt/JG dtD . (3.15)
33
Noting the spectral integration and applying the shifting properties of the Fourier
transform yields the ﬁnal result
h(r, 6, t) = f“ costb u(t — r/co + asin dcos rb/co) dip, (3.16)
ﬁrmnd 0
which is valid for 7‘ >> a and 6 > 0. Note that Eq. (3.16) has support on the
time interval [r/co — asin 0/c0,r/c0 + asin 0/c0], which follows from the fact that
f; coszp d2!) = 0. Thus, the temporal duration of the impulse response increases as 0
increases, which physically corresponds to the piston motion becoming blurred [75].
Unlike the nearﬁeld impulse response [48], Eq. (3.16) is symmetric about t = r/co.
Eq. (3.16) is evaluated in closed form as
Mr, 19, 7'):
«rand
[u(T — r/CO + asin B/co) — u(T — r/co — asin 6/c0)] \[_(r—
(3. 17)
which has the form of a semi-ellipse centered at r/q). Unlike the exact impulse
response [48], which is valid for observation points in both the nearﬁeld and the
farﬁeld, Eq. (3.17) is symmetric about t = r/co. On axis (0 = O), the lossless impulse
response reduces to
a2 ’
h(z,0,t) = -2—26(t — z/co), (3.18)
which is the limiting case of Eq. (3.7) as 2 —* 00.
3.3.2 Viscous Media
The farﬁeld viscous lossy impulse response is now computed from Eqns. (3.16) and
(3.18). For the off-axis case, Eq. (3.16) is substituted into Eq. (3.6). Interchanging
the order of integration and evaluating the unit-step function in the integrand yields
ca T’(‘l/)) tl2
h L(r, 6’ ,=t) 7rrsin0\/u27tr_t7_/ﬂ foo xep(- 2t )—cosz,bdt' dip, (3.19)
where
T,(1/)) = t — r/co + asin 6 cos w/co. (3.20)
34
a2 816:2
7-)2
. 2
The integration over t’ is evaluated usmg 72; ffoo 6‘2 dz = l + erf(a:). Due to
the coszp term in Eq. (3.19), the constant term vanishes, yielding a single integral
expression
_ coau(t) /” t—r/c0+asin6cosz,b/c0
hL(r,6,t) —— _2nr sin6 0 erf m cost/)dw (3.21)
The on-axis case is computed with Eq. (3.18), yielding
a2 _ z 2
hL(Z, O, t) = U(t)m exp ("LQ—tT/yﬂ), (3.22)
which can also be derived by letting z —+ 00 in Eq. (3.10). Unlike the lossless case,
the lossy impulse response does not have compact support in time. Since both the
Gaussian and error functions have inﬁnite support, the lossy impulse response defined
by Eqns. (3.21) and (3.22) are also non-zero for all positive time values. In physical
terms, this semi-inﬁnite region of support implies wave components traveling with
phase speeds varying from zero to inﬁnity, but the inﬁnite phase speeds are inﬁnitely
attenuated [17], so the result is causal.
3.3.3 Numerical Results
Figure 3.6 displays both the nearﬁeld and farﬁeld lossy impulse responses for a circular
.piston (a = 5 mm). Panel a) evaluates the lossy impulse response at r = 50 mm both
on-axis (6 = 0) and off-axis (6 = 7r/ 12 and d = 0). In this region, the farﬁeld
approximation differs signiﬁcantly from the nearﬁeld solution. Panel b) evaluates
the impulse response at r = 200 mm both on—axis (6 = O) and off-axis (6 = 7r/ 12
and d = 0). In this region, the farﬁeld approximation agrees with the nearﬁeld
solution, indicating that the simpliﬁed farﬁeld expressions accurately represent the
lossy impulse response at distances far from the source.
Numerical evaluations of the lossy impulse response show distinct behavior within
the nearﬁeld, the transition region, and the farﬁeld. For the values of 7 evaluated
in Fig. 3.6, the diffraction component dominates in the nearﬁeld, with only a small
35
amount of smoothing in the vicinity of the head and tail of the response. For obser-
vation points in the transition region between the near and farﬁelds, the effects of
diffraction and loss are both apparent. In this transition region, the impulse response
is both smooth and asymmetric. Finally, in the farﬁeld, the effects of viscous loss
predominate, yielding increasingly symmetric and broad responses with a reduced
directivity.
3.4 Rectangular Sources: Nearﬁeld
In ultrasound imaging applications, sources typically consist of large phased arrays of
rectangular or cylindrically curved strips. Since rectangular radiators are not axisym—
metric, the ﬁelds produced by these apertures depend on all three spatial Cartesian
coordinates (:17, y, z) or spherical coordinates (r,6, (1)).
3.4.1 Lossless Media
In this section, a single-integral representation for the nearﬁeld impulse response of a
rectangular radiator is derived. In order to apply a previous continuous-wave result
[2], the rectangle is decomposed into 4 sub-rectangles, as displayed in Figure 3.7. As
shown in [2], the Fourier transform of the impulse response (or transfer function)
hs,[(z;w) due to a sub-rectangle of width 3 and length l is given by
l eilcvrzgﬂriﬁsQ _ eikz s eisz§+a§+l§ _ eikz
do +1]
0
A pc
h z, w = , .9 do
3’“ ) 2mm [0 02 + 32 02 + 12
(3.23)
The total transfer function h(:r,y, z;w) is given by the superposition of the contri-
bution from each sub-rectangle. The resulting expression is an exact solution to the
nearﬁeld diffraction problem and analytically equivalent to the Rayleigh-Sommerfeld
integral.
36
The time-domain impulse response from each sub-rectangle is recovered via ana-
lytically inverse Fourier transforming Eq. (3.23), which yields
2.3,,(2, t) = % (13,,(2, t) + 11,3(2, t)) (3.24)
for a single rectangle, where
l
u(t — Ts) — u(t — z/c)
13,1(2, t) = 3/0 02 + 32 do, (3.25)
and the delay T3 is speciﬁed by
T3 = \/22 + 02 + 32/0. (3.26)
The impulse response is then written as
4
Way, 2,15) = Z ihs,,z,(z, t), (3.27)
i=1
where the lengths s,- and widths l,- and signs are determined from the geometry in
Fig. 3.7. In the case of y = 0, the total number of integrals in Eq. (3.27) reduces
from 8 to 4 due to symmetry.
3.4.2 Viscous Media
The lossy impulse response is computed by inserting Eq. (3.24) into Eq. (3.6) and
evaluating the unit step functions in each term of Eq. (3.25). The term I 3,1(2, t) in
Eq. (3.24) is replaced by a lossy term L3,1(z, t) given by
l t—Ts t—z/
L3,)(z,t)=u(t)s/O 1 [/ gL(t,t’)dt'—/ cogL(t,t’)dt’] do. (3.28)
s2 + 02 _oo -00
Note that the lossless speed of sound c has been replaced by the adiabatic speed of
sound co. The integrations over t’ are evaluated in terms of the error function [72]
erf(z), yielding
L.,z 0) were not considered. In this section, an analytical
nearﬁeld solution is provided for all ﬁeld points. The solution is derived for a spherical
shell centered at (0,0, —R) and surrounded by an inﬁnite, rigid baffle. As shown in
[77], the transient solution for a spherical shell with radius a and radius of curvature
R is given by
2a
h(r. 2. t) = f0“ f0 1’ N.,R(r. «0) [ac — n) — u(t — 72)] «w (331)
where the kernel function Na,R(r, 2,1/1) is given by
r'cosr,D\/R2 — a2 + za
N r, z, = 3.32
a’R( 7(1) R27‘2 + 2am cos w R2 — a2 — a2r2 cos2 1,0 + z2a2 ( )
and the delays 7',- are given by
1'1 = \/R2 + r2 + 22 — 2dr cosz/J + 22\/ R2 - a2/c0 (3.33a)
7'2 = (R + sign(z)\/ 7'2 + 22)/co, (3.33b)
Eq. (3.31) is valid at all observer locations except the geometric focus (0,0,0), where
the impulse response is given by
h(O, 0, t) = 2 (R — 6R2 — a2) 6(t — R/co). (3.34)
3.5.2 Viscous Media
Convolving Eq. (3.31) with gL(t,t’) via Eq. (3.6) yields
hL(r,z, t) = u(t)a n for Na,R(r,z,z/)) [erf 6%) — erf (5%)] dd), (3.35)
40
where erf(z) is the error function [72]. Eq. (3.35) is valid for all points in the acoustic
half—space excluding (0,0,0). The impulse response at the focus is computed by
inserting Eq. (3.34) into Eq. (3.6). The ensuing integration yields
hL(O,O, t) = u(t) i (R — V R2 — (12) exp (W) (3.36)
1rt7 2t7
Unlike the lossless case, the lossy impulse response is not inﬁnite at the geometric
focus.
As shown in [64], the on-axis impulse response for the spherical shell is described
by a closed-form expression (no numerical integration is required). Evaluating Eq.
(3.35) for r = 0 and z 7E 0 yields I
erf (t — \/R2 + 22 + zzJRﬂ—a/co) _ erf(t — (R+z)/c0)
Rco
hL(0, z, t) = U(t)—z /_2t,7 m
(3.37)
which corresponds to Eq. (6) in [64] multiplied by a factor of twol. The delay in the
ﬁrst term of Eq. (3.37) corresponds to the greatest distance from the observation point
and the shell, whereas the second delay term corresponds to the shortest distance from
the observer to the shell.
3.5.3 Numerical Results
Figure 3.12 displays the lossy impulse response generated by a baffled, spherical piston
with radius a = 10 mm and radius of curvature R = 70 mm in a viscous medium
with 7 = 0.001 ,us and c0 = 1.5 mm/ps. Eq. (3.35) is evaluated on an 80 by 80
grid at times t = 38 us and t = 46 us. Eq. (3.35) is evaluated with Gauss-Legendre
quadrature at each point in space-time. Note that the focus is signiﬁcantly degraded
due to the absorption of sonic energy by the medium. In particular, the peak value
of the lossy impulse response decreases as 7 —> 00.
To demonstrate the effect of loss on pulse propagation, velocity potentials in a
water-like medium and a liver-like medium are evaluated by convolving the impulse
1The baffled boundary condition was not utilized in [64].
41
response h L(r, t) with an input pulse v(t). In the following simulations, u(t) is given
by a ﬁve-cycle Hamming-weighted toneburst in Eq. (3.11) with center frequency f0 = 5
MHz. For liver, the attenuation coefﬁcient a given in [13] is converted via the relation
7 = coo/(27r2f2), yielding c0 = 1.58 mm/us and 7 = 7x 10“5 us. For water, c0 =
1.54 mm/ps and 7 = 10'6 us. Thus, the viscous relaxation time in liver is almost
two orders of magnitude larger than that in water. The velocity potential is evaluated
on-axis at a) z = -20 mm and b) 2 = 0 mm by performing FFT based convolutions,
yielding the waveforms displayed in Figure 3.13. As shown, the effect of loss has
little effect in water, allowing the pulse to form a geometric focus near z = 0 mm.
In the liver-like medium, however, loss noticeably impacts the focusing of the pulsed
waveform. In the pre-focal regions, shown in a), the pulse has been attenuated and
slightly stretched due to the absorption of sonic energy by the medium. In the focal
region shown in b), the amplitude of the waveform relative to water is diminished by
about a factor of 5.
Off-axis propagation is displayed in Figure 3.14. The velocity potential is shown
in the focal plane 2 = 0 mm at lateral positions a) a: = a/2 = 4 mm and b) a: =
a = 8 mm in both water-like and liver-like media. Relative to the on-axis velocity
potential shown in Fig. 3.13, the duration of the off-axis velocity potential is longer
due to differences in phase observed across the ﬁnite extent of the radiating shell. As
expected, this destructive interference increases as the observation point moves for
a: = 4 mm to :c = 8 mm. In the liver-like media, there is additional broadening of the
waveform due to viscous absorption, as well as a down-shift in the center frequency.
42
3.6 Discussion
3.6.1 Implementation Issues
Unlike the MIRF method and other schemes that synthesize ﬁelds in the frequency
domain, the lossy impulse response solution presented here is computed directly in
the time-domain without inverse FFTs. Numerical inverse FFT’s pose several diffi-
culties, including 1) additional computational burden and 2) possible numerical error
due to undersampling in the frequency domain. Since the bandwidth of the MIRF
decreases as a function of distance from the source, the Nyquist sampling frequency is
a function of distance. Therefore, an efﬁcient MIRF implementation requires multi—
rate sampling. The lossy impulse response method eliminates these problems by
replacing numerical inverse Fourier transforms that utilize a compact time-domain
expression with a constant sampling rate for all observation points. In addition, com-
puting snapshots of the impulse response at one particular point in time is particularly
straightforward with the lossy impulse response method presented in Eq. (3.9).
However, the major beneﬁt of the closed-form expressions derived in this chapter
is intuition for the behavior of the lossy impulse response. The error function is
essentially a smoothed unit step function. As the relaxation time 7 increases, the
error function becomes progressively smoother, indicating that the bandwidth of the
pulsed ﬁeld decreases. Hence, the basic physics of transient propagation in viscous
media are contained in the analytical lossy impulse response expressions.
3.6.2 Numerical Evaluation and Aliasing
As discussed in Ref. [78], evaluation of the lossless impulse response requires artiﬁ-
cially high sampling rates in order to accommodate the discontinuities in the deriva-
tive of the impulse response, which result in large bandwidths. For circular and
rectangular pistons, these discontinuities are magniﬁed on-axis and in the farﬁeld,
43
where the lossless impulse response is represented by a short-duration rect function
in the nearﬁeld and a scaled delta function in the farﬁeld. However, these numerical
difﬁculties are caused, in part, by an inaccurate physical model which assumes zero
loss. In reality, some loss is always present, effectively acting as a low-pass ﬁlter and
removing the high-frequency components of the impulse response. This ﬁltering ef-
fect is reﬂected in the lossy impulse response expressions, which are inﬁnitely smooth,
implying a Fourier transform that decays faster than 1 / f n for any n > 1 where f rep-
resents the frequency [79]. In contrast, the Fourier transform of the lossless impulse
response, due to discontinuities on-axis, decays as 1 / f . Due to this rapid decay in the
frequency domain, the lossy impulse response is bandlimited whereas the lossless re-
sponse is hand unlimited. The numerical difﬁculties arising from the band unlimited
lossless impulse response are thus eliminated by the inclusion of a loss mechanism,
which effectively bandlimits the impulse response. Since both the diffraction and loss
components are evaluated simultaneously in Eq. (3.3), sampling errors are reduced
with the lossless impulse response.
3.6.3 Power Law Media
Since the Stokes wave equation serves as a ﬁrst approximation for loss in biological
tissue, more accurate models, such as power-law media [19, 33—35], also demand con-
sideration. In power-law media, phase velocity is an increasing function of frequency,
resulting in an asymmetric loss function function 91,. Transient ﬁelds in power law
media are considered in Chapters 4, 5, and 6. In Chapter 7, loss functions for power-
law media will be derived, allowing the lossy impulse responses generated by various
piston sources to be computed in tissue-like media. The general power-law case will
utilize the same machinery as the viscous case with an altered loss function. There-
fore, the main beneﬁts of the lossy impulse response approach (no inverse FFT’s, lack
of aliasing, etc.) should apply in the more general power-law setting. Since all calcula-
44
tions are performed directly in the time-domain, the lossy impulse response approach
should provide a more intuitive solution to transient, dispersive problems. Unlike the
frequency-squared case, the general loss function for power law media cannot be ap-
proximated via a simple expressions such as a Gaussian. Thus, the derivation of the
impulse response for ﬁnite aperture radiators in power law media becomes more com-
plicated, necessitating more advanced mathematics. These mathematical subtleties
are discussed in detail in the following three chapters.
45'
' - - -iossy impulse response " ' 'lOSSY lmPUl8° WSW
2f ’9: coal .
‘ 0.0125> ‘
E 5 0.015»
E a
3 o 3 ° ‘
8' 8 —0.015~
E -0.01 25 1 g
, E -0.03i
‘°'°25’ . . . . . ‘ M45) . . . .
0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5
time (microseconds) time (microseconds)
(a) 7 = 0.01 ps and z = 0.1 mm. (b) 7 = 0.001 as and z = 0.1 mm.
x10'a . 1 .
4.5~ —MIRF reference ~ —M|FlF reference
. 0.025»
- - -lossy impulse response - - -lossy impulse response
a 3 ‘ g
E E 0.0125"
é 5
‘8 '6
E ‘3? °
8 8
g g —0.0125
> >
-0.025
0 0:5 5 1:5 2 2:5 3 0.5 i 125 2 25
time (microseconds) time (microseconds)
(c) 7 = 0.01 us and z = 1.0 mm. (d) 7 = 0.001 as and z = 1.0 mm.
Figure 3.2. On-axis (:1: = 0 mm) velocity potential for a circular piston of radius a = 10
mm in a viscous medium. The piston is excited by a Hanning-weighted toneburst with
center frequency f0 = 6.0 MHz and pulse length W = 0.5 us for 4 combinations of axial
distance 2 and relaxation time 7: a) z = 0.1 mm and 7 = 0.01 us, b) 2 = 0.1 mm and 7 =
0.001 as, c) z = 1 mm and 7 = 0.01 as, d) z = 1 mm and 7 = 0.001 [.18. The velocity
potential for each combination of z and 7 is computed via both the lossy impulse response
approach and a MIRF approach for veriﬁcation.
46
_A
U"
impulse response (mm / us)
.0
0|
8o
40
-2
axial distance (mm) 20 '40 lateral distance (mm)
(a) t = 22 us.
’g 1.5
E
V 1
3
g 0.5
0
ﬂ
3' 0
4O
—-2
axial distance (mm) 20 ‘40 lateral distance (mm)
(b) t = 34 as.
Figure 3.3. Snapshots of the lossless impulse response generated by a circular piston with
radius a = 10 mm at t = 22 us and t = 34 [1.5 are displayed in panels a) and b). The constant
amplitude component within the paraxial region |:r:| < a represents the direct wave. The
remaining component represents the edge wave generated by the discontinuity in particle
velocity at :1: = a.
47
impulse response (mm / us)
o
in
40
axial distance (mm) 2° ‘40 lateral distance (mm)
(a) t = 22 us.
—L
01
f
d
I
0",
.
-“
i"
n'.
\\\\\\\\
V
impulse response (mm / us)
.0
01
80
,,,,,,,
’’’’’
20
axial distance (mm) 20 ‘40 lateral distance (mm)
(b) t = 34 us.
Figure 3.4. Snapshots of the lossy impulse response generated by a circular piston with
radius a. = 10 mm with 7 = 0.01 as at for t = 22 us and t = 34 us are displayed in panels
a) and b). Unlike the lossless impulse response depicted in Fig. 3.3, the direct wave is
attenuated due to viscous diffusion. The edge wave also experiences additional attenuation
relative to Fig. 3.3.
48
.0
or
f
. ,~-,,,._ I‘m/141,"
. ."”.'.'. 00"IO'IO’I ’l’l’”, 3.le,! f. ”
pressure (nonnaiized)
.c'>
2' 9
I
.5
1
X (mm) 20 2 (mm)
(a) t = 22 ps.
'2": 1914/99, I,’
pressure (normalized)
x (mm) 20
2 (mm)
(b) t = 34 ,us.
Figure 3.5. Normalized velocity potential ﬁeld produced by a circular piston of radius a =
10 mm excited by a Hanning—weighted toneburst in a viscous medium with relaxation time
7 = 0.001 us. Snapshots of the normalized velocity potential for t = 22 us and t = 34 us
are displayed in panels a) and b), respectively.
49
. —on-axis (near)
- - - on—axis (far)
'---- off-axis (near)
off-axis (far)
.0
01
.9
.p.
.0
ro
impulse response (mm / us)
9 o
.1 co
up
[0
OJ
0.)
35
(a) r = 50 mm.
A —on-axis (near)
2'1 0.06- - - -on-axis (far)
E “off-axis (near)
g off—axis (far)
31104»
C
8
U)
‘2
g 0.02.
3
O.
.E
55>
1 132 133 134 135 136 137
time (us)
(b) r = 200 mm.
Figure 3.6. Lossy impulse response for a circular piston (a = 5 mm) and 7 = 0.001 us.
In panel a), the nearﬁeld impulse response and the farﬁeld impulse response were evaluated
at r = 50 mm both on-axis (6 = 0) and off-axis (6 = 1r/12 and ()5 = 0). Panel b) shows
the nearﬁeld and farﬁeld impulse responses evaluated at r = 200 mm on-axis (6 = 0) and
off-axis (6 = 1r/12 and (b = 0). «
50
I 2 I
:4
2b . Z
Figure 3.7. Coordinate axis used in the derivation of the impulse response for a rectangular
source. A rectangular piston of half-width a and half-height b is excited by a uniform particle
velocity. The radiator is surrounded by an inﬁnite rigid bafﬂe in the z = 0 plane. Reprinted
from [2].
0'6? —1I 10-5
" '1: 10—4
A A 0.5' 1.1-17.1”
g $0.4
' 8
g, gas]
. "e
3 3 0.2- ' '
a a. .‘ 1
.§ .5 . L
‘ 0.1 .' --._ 1 1
10“.? Q~~~l
1 ‘ ‘ Q l , ,.¢"’ 0 ‘. ..'e.,
o - G . ----- - ___.. - - - 2.11....—
66 66.5 67 67.5 68 132 132.5 133 133.5 134 134.5
time (us) time (us)
(a) z = 100 mm. (b) 2 = 200 mm.
Figure 3.8. Lossy nearﬁeld impulse response for a rectangular piston of half-width a = 2.5
mm and half-height b = 10 mm evaluated in lossy media for 7 = 0.001, 0.0001, and 0.00001
,us. The impulse response is evaluated at a) z = 100 mm and b) 2 = 200 mm.
51
-6 -4 -2 0 2 4 6
lateral distance (mm)
(a)7=0-
-2 0 2
lateral distance (mm)
(b) 7 = 0.0001 [1.3.
Figure 3.9. Point spread function for a focused, 64-elernent linear array with half-width
a = 0.3 mm, half-height b = 2.4 mm, kerf of 0.15 mm, focused on-axis at 80 mm. The
point-spread function is measured in the focal plane at 241 lateral locations in a) lossless
media and b) viscous media with 7 = 0.0001 us.
52
O
I
— lossless
-10 "-7: 5.6—5 -
'-'-"y=1e-4
—20
my: 2e—4
normalized pulse echo envelope (dB)
-30 ...................................
-40
—50
—60 ..........................................................................
'7906 105.5 107 102.5 103
time (us)
Figure 3.10. Normalized pulse echo envelope for a focused, 64-element linear array with
half-width a = 0.3 mm, half-height b = 2.4 mm, kerf of 0.15 mm, focused on-axis at 80
mm. The point-spread function is calculated at the focal point in lossless media and viscous
media with three different relaxation times.
53
Figure 3.11. A spherical shell with radius a and radius of curvature R, surrounded by
an inﬁnite rigid baffle. The origin is located at the geometric focus of the shell, the radial
coordinate denoted by r = V32 + 312, and the axial coordinate is indicated by z.
54
(D
N
l
W
14]
l
impulse response (mm/us)
(I; IMIIQ/
20
X (mm) _ 2 (mm)
(a) t= 38113.
(D
I] V
"479’ IW’I’
I
I [Mil/[Ill
.2. ... ”ﬁll/I’ll”) l
I.
I
I
I), 1,,
impulse response (mm/us)
(b) t: 46 us.
Figure 3.12. Lossy impulse response generated by a baffled spherical shell with radius a
= 10 mm and radius of curvature R = 70 mm in a viscous medium with viscous relaxation
time 7 = 0.001 us. Eq. (3.35) is evaluated on an 80 by 80 grid at times a) t = 38 p5 and
b) t = 46 us.
55
—water
1 ’ - - - liver '
.0
01
velocity potential (mm/us)
.c'>
U1 0
l
A
I
10.5 20 20.5
time (u s)
.4
(O
(a) z = -20 mm.
'3’
E
E,
E
E - - I I
93
8
.1?!
§
0
> 0 4
31.5 32 32.5 33 33.5
time (u s)
(b) 2 = 0 mm.
Figure 3.13. Velocity potential generated by a spherical shell with radius a = 8 mm and
radius of curvature R = 50 mm in a lossless, water-like medium, and a lossy, liver-like
medium. Velocity potentials are displayed on axis at a) z = —20 mm and b) 2 = 0 mm.
56
0.05 ‘
—water
A ---liver
‘él
E
1.5.
E.
E I
93
O
D.
92‘
8
d)
>
"0'05 31.5 32 32.5 33 33.5 34
time (u s)
(a) (2:,2) = (4,0) mm
x 10-3
A 6' —water
a 4_ "'llvel’
E
s 2_
g a I“
8, Q :5": 0
5—2
a -4'
>
-6-
31 32 33 34 35
time (u s)
(b) (33,2) = (8,0) mm.
Figure 3.14. Velocity potential generated by a spherical shell with radius a = 8 mm and
radius of curvature R = 50 mm in a lossless, water-like medium, and a lossy, liver-like
medium. Velocity potentials are displayed off-axis at z = 0 mm and a) :1: = 4 mm and b)
a: = 8 mm.
57
CHAPTER 4
Power Law and Szabo Wave Equations
4.1 Introduction
Although the power law model for the attenuation coefﬁcient in Eq. (1.1 is widely
accepted, few closed-form solutions that describe the effects of dispersion in the time-
domain have been reported. To date, approximate closed form Green’s function
solutions have only been previously provided for the special cases of y = 0 (frequency-
independent media) and y = 2 (viscous media). In order to construct Green’s func-
tions for general power-law exponents y, stable law probability distributions, which
have previously been used to study fractional diffusion [80], are required. As shown
below, these stable distributions facilitate analytical descriptions of power-law atten-
uation directly in the time-domain.
The purpose of this chapter is to derive analytical, time-domain 3D Green’s func-
tions in power law media for the intermediate cases of 0 < y < 2 in terms of stable
probability densities. These Green’s functions are exact solutions for power law me-
dia and approximate Green’s functions to the Szabo wave equation. In Section 4.2,
the Szabo wave equation is formulated as a FPDE, as well as a more general power
law wave equation. In Section 4.3, 3D Green’s functions to be obtained for fractional
power law media in terms of the Fox H -function [81] and the Wright function [80].
These functions have previously been used in the study of fractional diffusion [82],
58
fractional relaxation [83], and fractional advection—dispersion [37]. In the special cases
of y = 0, 1/3, 1 / 2, 2/3, 3/ 2, and 2, the Green’s function is expressed in terms of stan-
dard functions (Dirac delta, Airy, Lévy, hypergeometric, and Gaussian functions).
In the remaining cases, asymptotic and series expansions are provided. 2D Green’s
functions are also derived in power law media.
4.2 FPDE Formulations of the Szabo and Power Law Wave
Equations
4.2.1 The Szabo Wave Equation (0 < y < 1 and 1 < y < 2)
The Szabo wave equation [19] approximates power-law media with an attenuation
coefﬁcient given by Eq. (1.1). In this section, the Szabo equation is expressed in
terms of fractional derivatives, thereby allowing the machinery of fractional calculus
to be utilized. Linear attenuation media (y = 1) must be treated separately from
sub-linear (0 < y < 1) and super-linear (1 < y < 2) attenuation media.
The Szabo equation for acoustic pressure in the regime of 0 < y < 1 and l < y < 2
is written as (of. B7 in [19]):
13% 4h - t (t’)
v2 _ __ _ _._121 ___E_ I:
p c3612 Co ,/_00 lt—t’|y+2 dt 0 (4'1)
where Co is a reference speed of sound and hm; = —aOI‘(y + 2) cos [(y +1)y/2]/7r.
The reference sound speed co in Eq. (4.1) is taken as the high-frequency limit coo for
the sub-linear (0 < y < 1) case, while c0 assumes the low-frequency limit c; for the
super-linear( 1 < y < 2) case. Noting that (t — t') > 0 in the integrand and applying
trigonometric identities yields
dt’ :0 (4.2)
V2 _ i921? _ 200 P(y + 2) )sin (7ry)
1’ (96%? 2 jf;
0 C0 cos(7ry/ )
(t - pt') )3!”
The term in brackets in Eq. (4.2) is now identiﬁed as the hypersingular representation
of the Riemann-Liouville fractional derivative given by Eq. (A.1) in Appendix A.
59
where y > 0 and f (t) possesses ceil(y) integer-ordered derivatives. Calculating the
(y + 1)-th derivative and applying the reﬂection formula for Gamma functions
I’(z)I’(1 — z) = sing”) (4.3)
yields the bracketed term in Eq. (4.2). Thus, Eq. (4.2) is equivalent to
2,22: 2.. 3...,
cg 3t2 cocos(7ry/2) at?“
= 0, (4.4)
where the third term accounts for dispersive loss which is valid for 0 < y < 1 and
1 < y < 2. In the special cases of y = 0 and y = 2, the fractional derivative
term reduces to a ﬁrst and third temporal derivative, respectively. For fractional y,
the third term is deﬁned by the Riemann-Lionville fractional derivative deﬁned in
Appendix A. For y = 1, Eq. (4.4) is invalid since cos 7ry/ 2 -—> 0 as y —-> 1.
For fractional y, Eq. (4.4) is a FPDE that approximates lossy, dispersive media
satisfying a power-law. As discussed in Ref. [19], Eq. (4.4) is valid for 010 << 1
in the a) high—frequency limit for y < 1 and b) low-frequency limit for y > 1. As
shown below, a more general model for power law media is obtained by including a
higher-order temporal derivative. As in the case of the Stokes wave equation, the
particle velocity and velocity potential also satisfy Eq. (4.4). From the deﬁnition
of the fractional derivative given in Appendix A, Eq. (4.4) is also recognized as a
singular integro-differential equation, where the fractional derivative term depends on
the past history of pressure p(r, t).
4.2.2 The Szabo Wave Equation (y = 1)
In the special case of linear attenuation media (y = 1), the Szabo wave equation is
expressed as (see Eq. (36) in [19] with ho = 2010/ rt)
162 ,t 8 t
V2p(r,t)—— [9(1‘ )_ 00/ p(ritl) dtl=0. (45)
c? 3t2 7rc1 _00 (t _ t’)3
60
Like Eq. (4.1), Eq. (4.5) is a singular integro-differential equation. Unlike the
fractional cases 0 < y < 1 and 1 < y < 2, Eq. (4.5) cannot be expressed as a FPDE
and must be treated as a special case.
4.2.3 Frequency-Independent (y = 0) and Viscous (y = 2) Media
In the special cases of frequency-independent media (y = 0) and viscous media
(y = 2), the Szabo wave equation reduces to well-known integer-order PDEs which
are solved via standard methods. The y = 0 case corresponds to the telegrapher’s
equation
2 1 321) 24.03;; _
p ego 6t2 coo 8t —
which models one—dimensional damped string motion and electromagnetic wave prep-
0, (4.6)
agation in conductive media. The reference frequency coo is the phase speed in the
limit of inﬁnitely high frequency, and do is the attenuation coefﬁcient for w = l. The
3D Green’s function of Eq. (4.6) is [84]
— 2 I1 00600 t — R /
g(R, t) = €_OOR6(t R/COO)+u (t — R/Coo) LOCK) e—OOCoot ( 00),
47rR 47f (10600 W. _ R /coo
(4.7)
where 11(2) is the modiﬁed Bessel function of order 1. Eq. (4.7) consists of two terms:
an exponentially attenuated spherical wave and a dispersive wake. For small values
of the attenuation coefﬁcient cm, the dispersive wake is negligible. If 00 is small, the
asymptotic approximation 11(2) z z/ 2 is utilized, yielding
_ 6t—R 2
g(R.t) a: 6 “0R ( /C°°) +u(t — R/coo) 00306—3030.. (4.8)
47rR 87r
Thus, the dispersive wake is of order 043, yielding the approximation
g(R, t) z e—0036(t 1514600). (4.9)
Eq. (4.9) is an attenuated lossless Green’s function, thus indicating the lack of dis-
persion of the frequency-independent case for 00 << 1.
61
In the viscous case (y = 2), the Szabo equation reduces to the Blackstock equation
1 02p 2003319
V2p—ggt'2—‘F—Cz—‘6? =0, (4.10)
2
which models acoustic wave propagation in viscous media [42] for small 00 at low
frequencies. Due to the third-order temporal derivative in Eq. (4.10), Green’s func-
tion solutions to the Blackstock equation are noncausal. To demonstrate the loss of
causality, an impulse point-source is applied to the right—hand side of Eq. (4.10),
followed by a temporal Laplace transform and a spatial Fourier transform, yielding
. 1
CU” 8’ = — Hue, s)
(4.11)
where H(k, s) = 2aO/czs3 — sz/cg — k2. Since 1701:, 0) = —k2 < 0 and H(k,s) -> oo,
G(k, s) has at least one pole in the right hand plane. The inverse Laplace transform,
deﬁned by
1 .
k t = — Std 4.12
G( ,) ,m./La(k,s>e s < >
where L is the Bromwich contour, is now shown to be nonzero for t < 0. Consider
the closed contour C = L + Coo, where Coo is a semi-circle oriented clockwise in the
right-half plane of radius R —) 00. By Jordan’s lemma, the integral of G(k, s)eS’ for
t < 0 along Coo vanishes. Application of Cauchy’s theorem yields
1/CG(k,s) )e-“ds = ——/LG(,)ks )estals+-2—L G(k, s)eStds
2_7rz' 27ml Coo
Residue G(k,s)e"t = —./LG( 11:, s)e eStds
27ri
Due to the right half-plane pole, the residue is nonzero. Hence, G(k, t) > 0 to t < 0.
Therefore, the solution to Eq. (4.10) is noncausal.
By applying a plane-wave approximation to the third term in the Stokes wave
equation given by Eq. (2.1), Eq. (4.10) is obtained. The reference frequency c; is the
phase speed in the limit of zero frequency. Since the spatial bandwidth of a waveﬁeld
is large near the source, Eq. (4.10) fails to capture the high-frequency content in the
62
extreme near-ﬁeld of a radiator [17, 42]. The 3D Green’s function of Eq. (4.10) is
calculated via standard Fourier transform techniques, yielding
N 1 {—1 (t — R/cz)2
g(R,t) ~ 47rR 47rRa0 exp ( 4Ra0 ) ’ (413)
which is an approximate solution to Eq. (4.10). Similar to the asymptotic solution
to the Stokes wave equation [56] given by Eq. (2.20), Eq. (4.13) predicts Gaussian
spreading of a spherical wave as the ﬁeld radiates away from the source. However, Eq.
(4.13) is noncausal, which means that Eq. (4.13) speciﬁes non-zero ﬁeld values fort <
0. However, for (10 small and R sufﬁciently large, Eq. (4.13) and the causal Green’s
function to the Stokes wave equation (see Eq. (2.20) are virtually indistinguishable,
which is demonstrated numerically in Ref. [17] for the 1D case and in Chapter 2 for
the 3D case.
By utilizing the loss operator deﬁned in Eq. (A.2), the Szabo wave equation
interpolates between the telegrapher’s equation and the Blackstock equation. The
telegrapher’s equation is causal and hyperbolic, whereas the Blackstock equation is
noncausal and parabolic. Mathematically, Eq. (4.4) is a second-order, hyperbolic
equation for y < 1. For 3; > 1, however, Eq. (4.4) becomes an order y + 1 parabolic
equation where high-frequency components propagate with inﬁnite velocity. As y —+
1‘ in Eq. (4.4), coo -> co, and as y —> 1+ in Eq. (4.4), oz —> 00, which agrees with
the transition between hyperbolic and parabolic behavior at y = 1. Thus, the Szabo
wave equation interpolates between these two different behaviors as the power-law
exponent y varies from O to 2.
4.2.4 Power Law Wave Equation
As discussed in Ref. [19], Eq. (4.4) is an approximate model for wave propagation
in a power law medium for small values of (10. To facilitate the analytical solution of
wave propagation in power law media, a F PDE is derived that exactly describes power
law attenuation for arbitrary do. The Szabo wave equation is an approximation to
63
this exact power law wave equation. This exact FPDE is then solved for a point-
source in space-time, thereby yielding the desired Green’s function solution. Since
both temporal and spatial Fourier transforms are utilized extensively throughout the
derivation, the Fourier transform given by Eq. (3) in Ref. [19] is utilized here. In
order to derive an analytical time-domain Green’s function solution, a power law
dispersion relationship relating wavenumber k and angular frequency w is required.
This dispersion relationship should yield 1) a power law attenuation coefﬁcient and
2) a frequency—dependent phase speed. The power law dispersion relationship that
satisﬁes these requirements is
w a0(-i)y+1wy
k = — — 4.14
CD cos(7ry/2) ( )
for w 2 O and k(—w) = k*(w) to ensure real solutions. The imaginary part of Eq.
(4.14) yields the power law attenuation coefﬁcient given by Eq. (1.1). The real part
produces
1 _ Re k(w)
_ 1 7ry y_1
C(w)_ w —Co+aotan(2)|w| . (4.15)
Letting c1 denote the phase speed at wo = 1 yields the expression
ﬁ = ;11- + 001331163!) (Iw|y’"1 — 1) , ~ (4-16)
which is valid for all y aé 1. Eq. (4.16) corresponds to the phase velocities computed
via the Kramers-Kronig relations [7, 34] and the time-causal theory [6]. In addition,
Eq. (4.15) is in close agreement with the experimental dispersion data presented
in Refs. [5—7, 32]. Thus, Eq. (4.4) is consistent with the power law attenuation
and dispersion that is predicted by the Kramers-Kronig relationships and supported
by experimental measurements. Note that for O S y g 1, Eq. (4.15) predicts a
monotonically decreasing speed of sound with respect to frequency , while for 1 <
y < 2, the speed of sound increases monotonically with respect to frequency. Finally,
the phase velocities for all y is normal since there are no rapid changes in either the
64
attenuation coefﬁcient or phase velocities over any band of frequencies, as opposed to
relaxation-based dispersion observed in acoustics and electromagnetics [28].
In order to invert the wavenumber—frequency relationship into space-time, the
four-fold Fourier transform 13(k,w) of the time-domain pressure p(r, t) is multiplied
by the dispersion relationship. Squaring both sides of Eq. (4.14) and multiplying by
15(k, w) yields
—2 93— 20:0 —iw +1————(ﬁ—-——iw2 A w - l
[k +% comm/2% ’y cos2 1 since
cos(7ry/ 2) -—> 0 in the denominator. To circumvent this problem, a convolutional
wave equation is derived that exactly models linear power law media. To derive this
wave equation, a dispersion relationship at a reference frequency wo = 1 is considered.
The phase velocity is computed by taking the limit of Eq. (4.16) as y —> 1. Without
65
loss of generality, let y approach one from the right. Deﬁning y = 1 + e and letting
e —> 0, then
1 1 6
— = — + 00 tan(7r(1+ e)/2) (le — 1) (4.20)
6(0)) Cl
_ 1 a IwIe — 1
— c1 0tan(7re/2)
1
—» --+-2ﬂ’-(lw|‘—1)
c1 we
where the limit tanz —> z is used in the last line. To ﬁnish the computation, the
following limit is invoked [85]
'7 _
w 1 —) ln(w) (4.21)
7
yielding
1 1 2 a0
= — U . 2
—c(w) Cl 7r lnlwl (42 )
Eq. (4.22) is identical to the Hilbert dispersive model derived in Ref. [30] and the
Kramers-Kronig relation given by Eq. (13) in Ref. [7]. Combining Eqns. (1.1) and
(4.22) yields a dispersion relationship between the wavenumber k (spatial frequency)
and the temporal frequency w:
. w 2
k(w) = zaolwl + — — —ozowln [ca]. (4.23)
CI 1r
A convolutional 3D wave equation that satisﬁes Eq. (4.23) is now derived. In
order to invert the wavenumber-frequency relationship into Space-time, the four-fold
Fourier transform 13(k,w) of the time-domain pressure p(r,t) is multiplied by the
dispersion relationship. Squaring both sides of Eq. (4.23) and multiplying by 13(k, 0))
yields
2 (122 2010(—z'w) .. 2 .2 .
—k + C—2 — _—Cl—L(w) — 00L (4.)) P(k,Ld) = O (4.24)
1
where
A 2'
L(w) = M + £41111le (4.25)
66
Note that L(-w) = L" (4.2), which is required to ensure real solutions. Performing an
inverse spatial Fourier transform on Eq. (4.24) produces the Helmholtz equation
V2p(r,w) + [LU—22 — ga—(i-lﬂﬂw) - 03122010] p(r,w) = 0. (4.26)
‘31
Performing an inverse temporal Fourier transform using Eq. (4.26) with the aid of
the convolution theorem yields the wave equation
2 1' a
V2p(r, t) — 21,8 2],”) + 261°; [L(t) * p(r,t)] — a8L(t) =1: L(t) * p(r, t) = o (4.27)
where
L(t) = r1 [L(w)] = 27“]? (4.28)
is a generalized function (or distribution) [86]. Eq. (4.27) is hyper-singular since the
convolution integrals in the third and fourth terms diverge. To mitigate this problem,
Eq. (4.28) is integrated twice, yielding
2
L(t) = —%%M(t) (4.29)
where M (t) = u(t) In It]. Note that M (t) is weakly singular at the origin and hence
integrable. Inserting Eq. (4.29) into Eq. (4.27) yields
_ 1 3219(1‘, t) 400 53
64
V219“, 15) 2772—— Tel? [M(t) * P(r,t)l - 065,“: [M(t) * M(t) * 19(1‘, tll = 0-
(4.30)
Eq. (4.30) is a convergent, fourth order integro-differential equation.
The Szabo wave equation [19] is an approximation to the convolutional wave
equation when do is small. Assuming 010 << 1 and neglecting the 0% term in Eq.
(4.27) yields
162120, t) 200 6
of at2 c1 6t
Applying a temporal differentiation to L(t) yields
_4u(t)
7rt3 '
V2p(r, t) — [L(t) =1: p(r, 12)] = o. (4.31)
L'(t) =
(4.32)
Finally, expressing the convolution in Eq. (4.31) as an integration over the entire
history of p(r, t) yields Eq. (4.5).
67
4.3 3D Green’s Function
Time-domain 3D Green’s functions for power law media are derived in this section.
These Green’s functions are exact solutions to the power law wave equation in Eq.
(4.19) and approximate solutions to the Szabo wave equation in Eq. (4.4). Applying
a point-source at time t = 0 to Eq. (4.19) and considering the Green’s function g(R, t)
yields
2 y+l 2 2y
1 8 g 2040 3 9 “0 6 9 ._ —6(R)6(t) (4.33)
2 _ __ __ _ _
V 9 c3 (91:2 c()cos(7ry/2)6t3/+1 cos2(7ry/2) 3t2y
where R is the relative displacement between the source and the observer and R =
IRI. Applying a temporal Fourier transform to Eq. (4.33) yields
v29 + k2(w)§ = —6(R) (4.34)
where h(w) is given by Eq. (4.14). The well-known Green’s function for Eq. (4.34)
is given by 'k( )R
62 O)
47rR '
g(RM) = (4.35)
Inserting Eq. (4.14) into Eq. (4.35) and grouping terms yields
ﬁ(R,w) = [exP(:j’,ﬁ/C°’] [exp («04:0wa — Haney/244142540] (4.36)
where the ﬁrst factor solves the lossless Helmholtz equation. The time-domain Green’s
function for y aé 1 is recovered by applying an inverse Fourier transform and the
convolution theorem, yielding
gum) = Flame)» (437)
= [____"(t 452%] 4 r1 [exp (—a0R(|w|y — itan(7ry/2)w|w|y_1))]
which is expressed as a temporal convolution
g(R, t) = 900%) * 9L(R, 73). (4-38)
68
Within the framework of linear time-invariant (LTI) ﬁlters, Eq. (4.38) combines the
effects of dispersive loss and diffraction as a cascade of two system functions. The
ﬁrst factor in Eq. (4.38) is the Green’s function for the lossless wave equation given
by Eq. (2.31a), which is responsible for the effect of propagation and diffraction. The
second term is a loss function deﬁned via
gL(R, t) = 3:4 [exp (—aoR(|w|y — itan(7ry/2)w|w|y_l))] . (4.39)
Interestingly, the inverse Fourier transform deﬁned by Eq. (4.39) may be calculated in
closed form using the machinery of stable distributions [87]. The expression within the
square brackets in Eq. (4.39) is identiﬁed as the characteristic function of a stable
distribution with index y, center 0, skewness 1, and scale (Rag)1/ 3’. The inverse
transform given by Eq. (4.39) is then expressed concisely as
1 ~ t
9L(R: t) = mfg; (W) (4-40)
where the function fy(t), which is independent of R and cm, is the maximally skewed
stable distribution of index y with scale 1. An expression for fy is given by Eq. (B3)
in Appendix B by setting the skew parameter 6 = 1. This family of functions has
been studied extensively [49, 88, 89]. Software is available for the computation of
stable probability density functions (PDFs) and cumulative distribution functions [3]
based on formulas developed in Ref. [90]. In the following analysis, an alternate
parameterization fy(t) is deﬁned via Eq. (B4) in Appendix B in order to utilize
existing expressions for the stable density fy(t). In the case of both sub-linear power
law media (y < 1) and super-linear power law media (y > 1), fy(t) possesses the
following analytical properties:
1. fy(t) 2 0 for all t.
2. fy(t) is inﬁnitely continuously differentiable (or smooth).
3. fy(t) is unimodal.
69
4. ffo00 fy(t) dt = 1.
Although these properties are common to all maximally-skewed stable distribu-
tions, the behavior of fy(t) differs for y < 1 and y > 1. Properties of these func-
tions are studied in the following two subsections for sub-linear (y < 1), super-linear
(y > 1), and linear power law attenuation media.
Eq. (4.38) is an approximate Green’s function, valid for small 00, for the Szabo
wave equation given by Eq. (4.4). However, Eq. (4.38) is an exact solution to
the power law FPDE given by Eq. (4.19) for linear, homogeneous media since the
imaginary part of Eq. (4.14) exactly describes the power law behavior that the Szabo
wave equation approximates. Thus, the analytical Green’s function given by Eq.
(4.38) is not restricted to small 00.
4.3.1 Sub-Linear Power Law Media (0 < y < 1)
For the general case of 0 < y < 1, fy(t) was ﬁrst considered in Ref. [91] and later in
Ref. [92] as an inverse Laplace transform. Rewriting Eq. (B.4) as a inverse Laplace
transform yields
fy(t) = £‘1[e‘sy} = i- f1. 6%“ d3, (4.41)
27ml
where L is a contour in the left half plane. The notation is Eq. (4.41) was previously
utilized in Ref. [93]. The following analytic properties of fy(t) for 0 < y < 1 are
known [88, 89]
l. fy(t) has support on [0,00).
2. Ast—» oo, fy(t) ~t—y—1.
The inverse Laplace transform deﬁned by Eq. (4.41) is expressed in terms of the
Fox H -function [81] using formulas presented in Ref. [83]. The Fox H -function, which
is deﬁned in terms of a contour integral by Eq. (C.1), generalizes many of the special
70
functions, such as Bessel, hypergeometric, and Mittag-Leffler [81]. Using Eq. (62) in
Ref. [83], the fractional exponential possesses the representation
1 1,0 (1 1)
ex —sy— — —’—H1 3 4.42
and the inverse Laplace transform is then calculated using Eq. (30) in Ref. [83],
yielding
H,10 1 (0 1)
t _ __ 4.43
fy-UW11(I(01/y)) ( )
Eq. (4.43) possesses a simpler form for y = 0, 1/3, 1/2, and 2/3. For y = 0, 1/3, 1/2,
and 2/ 3, Eq. (4.43) possesses the following exact representations:
M0 = 5(t), (4.44a)
f1Mt):(_3__':1:4()t1)/3A"[(—1—'3,g)1/3]’(44%)
f1/2(t)= ut)( 2 [113/2 xI)-( 41), anda (4.44c)
f2/3(t)=U(t)-3}Wexr> (— ,,—,-,)2F («é—é; 12112), (4.444)
where Ai(z) is the Airy function and 2F0(a, b; ; z) is the conﬂuent hypergeometric
function of the second kind.
For the general case, an asymptotic expression, valid for t << 1, is derived using
the asymptotic form [81] of H11,’?(z). Alternatively, the inverse Laplace transform
may be approximated via the method of steepest descents [93], yielding
A B
fy(t) ~ u(t); exp (—t—#-) . (4.45)
The coefﬁcients A, B, u and A depend only on y and are given by
yl/(2-2y)
= —, (4.46a)
2”(1 - y)
B = yl/(l-y)1—_y (4.46b)
y
71
, and (4.46c)
=__, 4.
(4 1_y ( 46d)
Eqns. (4.38), (4.45), and (4.46) provide a closed form, asymptotic approximation to
the time-domain 3D Green’s function for sub-linear power law media (y < 1). These
equations exactly satisfy Eq. (4.440) at all times t for y : 1 / 2 and are valid for y 75 1 / 2
when t < 1. For y z 1 / 2, Eq. (4.45) also provides an excellent approximation [94]
for all values of t. For t > 0, fy is inﬁnitely continuously differentiable (or smooth);
moreover, fén)(t) = 0 as t -—> 0+ for all n 2 0, implying strong causality [17, 93].
However, for y 75 1 / 2 and t >> 1, the asymptotic solution given by Eq. (4.45) decays
exponentially for large times, while fy(t) decays algebraically. Therefore, for large t,
the ﬁrst term in the asymptotic series for Eq. (4.43) is utilized, yielding [88]
~ i ( /2)I‘( +1)
fy(t)z‘°’“”ym,+,y ,
(4.47)
which indicates a slowly decaying tail.
To show the behavior of fy(t) for 0 < y < 1, the PDFs for the sub-linear attenua-
tion cases y = 1 / 3, 1 / 2, 2/ 3, and 9/10 are displayed in Figure 4.1. Since the ordinate
fy(t) is a PDF, the abscissa t is unitless. The expressions given by Eq. (4.44) are
evaluated for y = 1 / 3, 1 / 2, and 2/ 3, while the STABLE toolbox 1 is utilized [3] for
y = 9/ 10. In all of these cases, fy(t) is identically zero for t < 0 and smooth. As y
increases from 1 / 3 to 9/ 10, the main plume spreads out, while the asymptotic decay
accelerates from t’4/3 to t‘19/10 due to the behavior predicted by Eq. (4.47).
As y -—> 0, fy(t) —+ 6(t), while as y -—> 1’, fy(t) —> 6(t — 1). Thus, there is a delay
incorporated into the stable distribution which accounts for the slower phase velocity
in a dispersive medium. This delay is apparent in the asymptotic representation given
by Eq. (4.45). For t very small, the argument of the exponent in Eq. (4.45) is very
large, causing fy(t) to approach zero. The delay td of fy(t) can be approximated by
llittzp : //academic2 . american . edu/ " jpnolan/ stable/ stable . html
72
2.5 . .
—y=1/3
: ---y=1/2
2. 51 ""“'v=2/3 .
!! '-'-'y=9/10
[I
A i-
5 n!
~..>1.5- ii 1
u. 2'
D I !
O. i |
2 ! !
.0 1- | I .1
.59 ",i I
”’ -.i -.
1 I
‘ I': i
0.51- ‘é'g 1 -
.- 1).
92 0 A ‘ 2 4 6
t(unitless)
Figure 4.1. Plots of stable PDFs fy(t) for four sub-linear attenuation cases: y = 1/3, 1/2,
2/3, and 9/10. Eqns. (4.44) are evaluated for y = 1/3, 1/2, and 2/3, while the y = 9/10
case is evaluated with the STABLE toolbox [3].
setting the exponential argument of Eq. (4.45) equal to unity, yielding td z 81/“.
Letting y vary from zero to one, td is also seen to increase smoothly from zero to one.
4.3.2 Super-Linear Power Law Media (1 < y S 2)
In super-linear power law media, the qualitative behavior of the loss function gL(t, t’)
differs from that observed in sub-linear power law media. The maximally skewed
stable distributions possess the additional analytical propertias for 1 < y S 2:
1. fy(t) has support on (—oo,oo).
2. As t —> oo, fy(t) ~ t-3’_1 for y < 2.
3. As t —> —oo, fy(t) decays with exponential order.
73
Hence, for y > 1, fy(t) is non-zero for all times t. Since the 3D Green’s function
is simply a delayed and scaled version of fy(t), the power law Green’s function g(R, t)
is non-zero for all t and is therefore noncausal in super-linear power law media for
all values of y such that 1 < y S 2. Since fy(t) decays with exponential order as
t —) —00, the solution is very small relative to the peak value for t << R/co. This
property, which was discussed by Szabo [19] for the quadratic case (y = 2), where
fy(t) is Gaussian, is also true for all 1 < y S 2 according to the stable law properties
listed above.
The stable distribution fy(t) is evaluated in terms of the Fox H -function (see Eq.
(2.13) in Ref. [95]) via
fy(t)=%H,1;{’(—t[ (1 ‘(t/ﬁsl/y) ), (4.48)
where 1 < y S 2. Eq. (4.48) may be simpliﬁed in terms of the Wright function
d>(a, b; 2:), yielding
1 1 1
fy(t) — ;¢(—§,1—g,t) , (4.49)
which is a known solution of the time—fractional diffusion problem [96]. For the
special cases of y = 3 / 2 and y = 2, Eq. (4.49) is expressed in terms of the conﬂuent
hypergeometric function of the second kind and a Gaussian, respectively:
5722 0 616): 1‘3 1 >
f3/2(t) = 4:7“ M3 25 t ) (4.50a)
_ exP l _1.. 27 -
3E71rt ZFO (6’ 6’ ’ 473) I” < O
eXP(—t2/4)
t = —. 4.50b
Note that the y = 2 case in Eq. (4.50b) yields Eq. (4.13), which is an approximate
solution to the Blackstock equation [42] in Eq. (4.10).
Unlike the solution described earlier for sub-linear power law attenuation (y < l),
fy(t) is non-zero for all t. Therefore, three separate cases are considered: early-time
(t < R/co), late-time (t >> R/co), and times near the wavefront. To obtain an
74
early-time estimate, Eq. (4.48) is expanded in a Taylor series [95]. If t << R/co, then
tr -—2 00, allowing an asymptotic representation of Eq. (4.48) to be utilized [95]:
fy(t) z Altluexp (—B|t|") (4.51)
where
A y—l/(2y_2) (4 52 )
= —, . a
V2401 - 1)
B = y“y/(y‘1)(y — 1). (4.5210)
_ 2 - v
V — 2y _ 2, and (4.52c)
_ y
p —— —y _ 1. (4.52d)
For y = 2, Eq. (4.51) reduces to f2 (t) for all t, yielding the non-causal Blackstock so-
lution. Eq. (4.51) has the form of a “stretched exponential,” which is a characteristic
of anomalous diffusion [97]. In other words, the early-time response of super-linear
power law media is characterized by fractional diffusive behavior. For large t, the
asymptotic formula given by Eq. (4.47) is utilized.
Since Eq. (4.51) is non-zero for all t, power law Green’s functions are noncausal
for 1 < y S 2. However, the rapid exponential decay for t << R/co makes the
solution very small for large negative times (t << —1). Thus, the non-causal nature
of the Green’s function for 1 < y S 2 is not evident in numerical evaluations for
small do. From a physical perspective, the phase velocity becomes unbounded as
frequency increases, implying that high-frequency components propagate inﬁnitely
fast. However, as noted in Ref. [17] for the special case of y = 2, these high-frequency
components experience high absorption and attenuate over a short distance, implying
the Green’s function g(R, t) << 1 for t << R/co.
To show the behavior of fy(t) for 1 < y S 2, the stable PDFs for the super-linear
power law attenuation cases y = 11/ 10, 3/ 2, 19/ 10, and 2 are displayed in Figure
4.2. The expressions given in Eq. (4.50) are evaluated for y = 3/ 2 and 2, while the
75
A
O)
.4
—y=11/10
14_ ” ---y=3/2 _
~-~~y=19/10
1.2- ""'Y=2 ~
3:; 1- .
u.
E 08- -
33
3 0.6- .
(D
O.4~ 4
02* .
Gun-All! 1 s
.4 6 8
2
t (unitless)
Figure 4.2. Plots of stable PDFs fy(t) for four super-linear attenuation cases: y = 11 / 10,
3/2, 19/10, and 2. Eqns. (4.50) are evaluated for y = 3/2 and 2, while the y = 11/10 and
19/ 10 cases are evaluated with the STABLE toolbox [3].
STABLE toolbox [3] is utilized for y =11 / 10 and 19/ 10. In all cases, fy(t) is non-zero
for all t and smooth. For y = 11/ 10, fy(t) is skewed to the right and decays very
rapidly for t < —1. As y increases from 11/10 to 2, the PDF becomes increasingly
symmetric, eventually converging to a Gaussian. Although y = 11/ 10, 3/ 2, and 19/ 10
decay like t—y‘l, y = 2 decays much more rapidly.
4.3.3 Linear Power Law Media (y = 1)
In this section, Green’s functions to the convolutional wave equation in Eq. (4.30)
are derived following the same procedure used in Section 4.3. Applying an impulsive
76
point source at the origin to Eq. (4.30) produces
1 029(r,t) 20:02
ET Kat [1;(1) * g(r, 1)] — 0.3m) * L(t) =1: g(r,t) = —6(t)6(r).
(4.53)
V2g(ri t) _
Applying a temporal Fourier transform to Eq. (4.53) yields Eq. (4.34) where k(w) is
given by Eq. (4.23). The 3D time-domain Green’s function is computed by inserting
Eq. (4.23) into the frequency-domain Green’s function given by Eq. (4.35) and
assigning the delay term 02/ c1 to §D(R, w). Inverse Fourier-transforming this product
of diffraction and loss factors yields Eq. (4.38) with a loss function
gL(R, t) = .7-“1 {exp [—a0R|w| (1+ 2isgn(:)ln |w|)] }, (4.54)
where the signum function sgn(w) is introduced to enforce conjugate symmetry. Eq.
(4.54) is recognized as a stable distribution with index 1, center 0, skewness l, and
scale 00R using the parameterization in Ref. [87]. An expression for f1(t) is given
by Eq. (B2) in Appendix B by setting 6 = 1. Thus, the inverse Fourier transform is
expressed as
gL(R.t)=#f1( ’ ) (4.55)
00R 00R
where f1(t) is the inverse Fourier transform of Eq. (4.54) with 00R = 1. Carrying
out the transform by integrating over negative and positive frequencies and expressing
the cosine in terms of complex exponentials yields
f1(t) = i/Ooo exp(—w) cos [wt + gwln w] dw. (4.56)
Although f1(t) cannot be expressed in terms of the Fox or Wright functions, this
expression is computable using the STABLE toolbox [3]. Moreover, f1 (t) has the same
analytical properties as the super-linear power law solutions shown earlier, including
smoothness and exponential decay as t —> —00. Since fy(t) is nonzero for all t, the
power law Green’s function is noncausal for y = 1. However, for do small, the Green’s
function is very small for negative times. For t >> 1, the asymptotic behavior of f1 (t)
is given by Eq. (4.47), which estimates the tail of the Green’s function.
77
4.4 2D Green’s Functions
The two-dimensional Green’s function g(p, t), which models radiators with inﬁnite
extent in the z-dimension, may also be calculated via the machinery of stable dis-
tributions. The 2D frequency domain Green’s function may be represented as a
superposition of plane waves via
. (1)
40.4) = ”’0 [WP] (4.57)
= % [000 exp (ik(w)pcosh (b) dip, (4.58)
where H61) (2) is the Hankel function of the ﬁrst kind of order zero and p = \/ :02 + y2
is distance in 2D. Inserting the dispersion relation given by Eq. (4.14) into Eq. (4.58)
and performing the inverse Fourier transform yields
1 0° 1 t—pcoshz/J/c)
g(p’b’ 27/0 (ﬁop008h¢)’/yfy((HopcoshWI/y d’b (4'59)
Direct substitution of t = p/c into Eq. (4.59) shows that g(p, t) is smooth at the
“arrival time,” unlike the lossless case. To gain insight into Eq. (4.59), the integrand
is decomposed into a loss function convolved with a plane wave:
1 0° 1 t
9(10,t)--/027r (ﬂopcosh ¢)1/yfy ( (ﬂopcosh ¢)1/y)5(t-pcosh¢/C) 61¢ (4 60)
Noting that the loss function is a slowly varying function of cosh 4p, while the delta
function is a rapidly varying function allows the approximation cosh 1/) z 1 in the loss
function. Evaluating the resulting integral in closed form yields
t
where gD(R, t) is the lossless 2D Green’s function given by
90(p,t) = 24% (4.62)
Thus, the 2D Green’s function is approximated as the lossless Green’s function con-
volved with a loss function. Unlike the 3D case, this decomposition is not exact
78
due to the approximation costh x 1. For t z p/c, the loss function in Eq. (4.61)
smoothes the discontinuity due to the ﬁrst arrival of sound. However, fort >> p/c, the
convolution in Eq. (4.61) is dominated by the lossless diffraction term, yielding long-
term behavior that is similar to the lossless Green’s function. This idea is explored
numerically below.
In the sub-linear case, fy(t) is zero for t < 0, thus yielding a ﬁnite range of
integration
u (t —- p/coo) /°°Sh_1(C°Ot/p) 1 (t — pcosh (D/coo)
dtb.
27r 0
gm” = (ﬂopcosW/y ” (ﬁopcosh 4)”?!
(4.63)
Although the range of integration is inﬁnite in the super-linear case, Eq. (4.51) indi-
cates that fy(t) decays rapidly for t < 0, allowing the upper-limit to be approximated
as ﬁnite. Since fy(t) is smooth, g(p, t) must also be smooth (i.e. C°°).
4.5 Stochastic Interpretation
Since the Green’s function for the power law wave equation involves a PDF, a natural
question arises: can power law attenuation be explained by an underlying stochastic
process? To explore this question, consider the velocity potential produced by a point
source with velocity potential v(t):
**(R, t) = g(R, t) at v(t) (4.64)
Utilizing Eq. (4.38) and the deﬁnition of convolution, Eq. (4.64) is written as
@(R, t) = 17% I: gL(R, t' — R/co)v(t — t') dt’ (4.65)
Since the loss function 91, is expressed as a PDF, the integral in Eq. (4.65) is recog-
nized as an expected value
**(R, t) = 171—RE [u(t — T)] , (4.66)
79
where the random variable T is distributed by the shifted and scaled stable distribu-
tion
_ 1 ~ t- 55/60
9L(R,t — R/CO) — ny (W) . (4.67)
Physically, the random process T(R) represents time delays that arise from the inher-
ent heterogeneity of the medium. The process T(R) consists of two components: 1)
a bulk delay R/co and 2) additional, random delays due to the micro—heterogeneity
of the intervening medium. Since ﬂ, is skewed to the right, there is a much larger
probability that the outwardly propagating spherical wave encounters delay larger
than R/co.
4.6 Numerical Results
To numerically verify the analytical results presented above, the analytical Green’s
function was compared to the material impulse response function (MIRF) approach
developed in Ref. [31]. In Fig. 4.3a), the 3D Green’s function is computed in a
power-law medium with parameters y = 1.5, co = 0.15 cm/as, and a0 = 0.1151
Np/MHz1'5/cm using the analytical formula in Eq. (4.38). In Fig. 4.3b), the inverse
Fourier transform deﬁned by Eq. (4.39) was calculated numerically with the MIRF
approach via an inverse FFT. The two panels demonstrate excellent agreement. For
power law exponents y = 1 / 3, 1 / 2, 2/ 3, 3/ 2, and 2, the function fy(t) is represented
in terms of the Airy, hypergeometric, exponential, and Gaussian functions. Eqns.
(4.44) and (4.50) are evaluated numerically using the GNU Scientiﬁc Library (GSL)
[98]. For other values of y, the function fy(t) is represented as a maximally-skewed
stable distribution computed using the software package STABLE [3], which is a
general-purpose software for analyzing stable distributions. The results of STABLE
have been independently veriﬁed using codes developed by G. Robinson [99] and a
combination of asymptotic expressions, power series, and inverse power series. All
80
demonstrate excellent numerical agreement. Transient ﬁelds are then evaluated using
FFT-based convolutions.
Figure 4.4 shows the 3D power law Green’s function for y = 1 / 2, 3/ 2, and 2 for
do = 0.05 mm"1MHz"y. The Green’s functions are shown as functions of the radial
coordinate R for times t = 20, 30, 40, and 50 us in order to emphasize the spatial
distribution of acoustic energy in power law media. The qualitative properties of
sub-linear power law media (y = 0.5), super-linear power law media (y = 1.5), and
viscous media (y = 2) are displayed in these plots. In these ﬁgures, the impulsive
excitation has been stretched and smoothed by the ﬁltering effect of the medium. The
Green’s function for sub—linear power law media (y = 0.5) has a sharp wavefront at
R = cot, where the ﬁeld is identically zero to the right of this wavefront. Neither the
super-linear power law attenuation Green’s function nor the viscous Green’s function
have this sharp wavefront. In the y = 1.5 case, most of the energy is located in the
slowly decaying tail of the Green’s function, whereas in the viscous case (y = 2),
energy is symmetrically distributed about the location R = cot.
As time evolves in Figure 4.4, the peak-amplitudes of all three of the Green’s
functions decrease, thus accounting for the transfer of acoustic energy into random
thermal energy within the intervening medium. In the y = 0.5 and y = 1.5 cases,
the duration of the Green’s function increases, which is a consequence of the heavy
tail behavior predicted by the asymptotic expression in Eq. (4.47). In the viscous
case (y = 2), spreading of the main plume occurs, although the Green’s function is
more localized relative to the fractional cases due to the rapid decay of the Green’s
function for large times. The absence of a slowly-decaying tail is a consequence of the
minimal dispersion associated with viscous and thermo-viscous loss.
To explore the effect of frequency dependent attenuation on broadband pulse
propagation, the velocity potential due to a point source at the origin was computed
81
by convolving Eq. (4.38) with the pulse [100, 101]
v(t) = A0t3 exp (—ﬁt) sin (27rf0t)rect (17;) . (4.68)
The following simulations utilize the center frequency f0 = 2.5 MHz, pulse length
W = 1.2 as, damping factor 6 = 9.3750 [is—1, and a normalization factor of A0 =
616.3 mm/ps. Two examples of velocity potential calculations are shown in two
different power law media using parameters from Ref. [13]: 1) a fat-like medium
with y = 1.5, CO = 1.432 mm/as, and a0 = 0.086 Np/MHzl°5/cm, and 2) a liver-like
medium with y = 1.139, c0 = 1.569 mm/ps, and do = 0.0459 Np/MHz1-139/cm.
Figure 4.5 displays the velocity potential as a function of time t at radial distances
10 mm, 25 mm, 50 mm, and 100 mm to show the effect of dispersion on the pulse
shape. As R increases in Fig. 4.5, several qualitative features become evident. First,
the pulse experiences rapid attenuation due to the combined effects of loss and spher-
ical spreading. Second, the pulse changes shape due to the more rapid attenuation of
the high-frequency components. In other words, the spectrum of the pulse experiences
a frequency down-shift. There are signiﬁcant differences between the fat-like medium
and the liver-like medium. The fat-like medium is more dissipative; therefore, greater
attenuation and greater frequency down-shifts are observed for large depths (R = 100
mm). Although the liver-like medium is increasingly attenuated as the propagation
distance increases, there is little change in the pulse shape. The pulse distortion in
the fat-like medium indicates that the effects of dispersion become signiﬁcant for large
depths, which is consistent with the conclusion reached in Ref. [34].
Figure 4.6 displays the 2D Szabo Green’s function given by Eq. (4.59) and the
lossless Green’s function given by Eq. (4.62) at distances a) p = 10 mm and b) p = 50
mm. In each panel, 010 = 0.05 mm‘lMHz'y for the sub-linear case y = 2/ 3 and the
super-linear case y = 1.5. Several properties are notable. Unlike the lossless Green’s
function given by Eq. (4.62), the lossy Green’s function is continuous at the arrival
time p/c. While the sub-linear (y = 2/3) Green’s function has a sharp wavefront
82
at p/c, the super-linear (y = 3/ 2) Green’s function is non-zero to the left of this
arrival time. Unlike the 3D case, the long-term behavior of the lossy Green’s function
approaches the lossless case. This limiting behavior is predicted by the approximate
decomposition given by Eq. (4.61), which is dominated by the slow decay of Eq.
(4.62) for t > p/c. That is, the geometric tail in Eq. (4.62) decays like t—l, which is
slower than the decay of the stable PDF fy(t), which behaves like t—i"l for t large.
4.7 Discussion
4.7.1 Causality and the Super-Linear Attenuation Paradox
In Ref. [19], the “Quadratic Loss Paradox” associated with the 1D solution of the
Blackstock equation given by Eq. (4.10) is presented. In the discussion of Sec. VII
of Ref. [19], Szabo notes that the 1D solution is non-causal, yet provides an excellent
approximation to the thermoviscous wave equation for small values of e = ozoc0|w|3=’_1
where w is the upper frequency limit. Szabo’s analysis is applied to three dimensions
and super-linear power law media (y > 1).
By letting R = R//\, where A is the wavelength associated with w, the ratio of the
Green’s function at time t = 0 to the peak value is given by
_ a 1/
n = fy( R/co/( OR) 3’) (4'69)
maxfy(z)
f, (ml—lax)
maxfy(z)
where x = [sec(1ry/2)|1/y/(c0a(l)/y). Letting the normalized distance R = 1 in Eq.
(43) yields the upper limit 2010g10 n < -136 dB for 1 < y g 2. However, values of 1] are
typically much smaller than this upper limit. For instance, for y = 1.5, 2010g10 17 < -
360 dB. Similar computations for other y 2 1 demonstrate that for observation points
only one wavelength from the radiating source, Eq. (4.38) is effectively causal for all
0 1). Firstly, the Green’s function
is not causal. Secondly, the fractional derivative operator has order y + 1 > 2,
which requires three initial conditions: 1) p(r,0), 2) p(r,0), and 3) p(r,0). Since
the third initial condition is not physically meaningful, Chen and Holm modiﬁed the
Szabo wave equation given by Eq. (4.4) by incorporating dissipation via a symmetric,
space fractional derivative [20]. Unlike the Reimann-Louville fractional derivative,
which utilizes a convolution on (—oo, t), the fractional Laplacian involves a symmetric
convolution over either the entire real line 'R in a 1D setting, or the real Euclidean
space ’R.” in a general n-D setting. The resulting Chen-Holm wave equation is second-
order and local in time and order y in space, thus only requiring two initial conditions.
Unlike the Szabo and power law wave equations, which are nonlocal in time and local
in space, the Chen-Holm formulation is local in time yet non-local in space. Like the
Szabo equation, the Chen-Holm model assumes a small attenuation parameter do.
In this chapter, the 3D Green’s function for the Chen-Holm equation is derived
91
in terms of symmetric stable distributions [49]. Symmetric stable distributions have
previously been used in polymer physics modeling [104], space-fractional diffusion
equations [45, 82], and relaxation processes [105]. The 1D and 2D Green’s func-
tions are also calculated explicitly in terms of the stable PDF and GDP, respectively.
Finally, the 3D Green’s function is given physical justiﬁcation by incorporating micro-
heterogeneity into a viscous model governed by the Stokes wave equation.
5.2 Chen-Holm Equation
5.2. 1 Formulation
In 3D, the Chen-Holm wave equation is given by [20]
V2 _i&_ﬂ(_ 2)y/2§E=
6t
0 (5.1)
where the power law exponent y ranges from 0 to 2. Eq. (5.1) is second-order in time
and order y in space. The loss operator, contained in the third term of Eq. (5.1),
contains a fractional Laplacian, or Riesz fractional differentiation, which is deﬁned
by spatial Fourier transforms via Eq. (A5) in Appendix A. In the viscous case
(y = 2), the fractional Laplacian in Eq. (5.1) reduces to the classical Laplacian and
Eq. (5.1) reduces to the Stokes wave equation given by Eq. (2.1) with 7 = Zoo/coo
and co = coo. In the frequency-independent case (y = 0), the fractional Laplacian in
Eq. (5.1) reduces to the identity operator and Eq. (5.1) reduces to the telegrapher’s
equation given by Eq. (4.6). Hence, the Chen-Holm equation interpolates between
the telegrapher’s equation and the Stokes’ equation.
Before solving the Chen-Holm equation, the dispersive and causal properties of
Eq. (5.1) are established.
92
5.2.2 Dispersion Relationship
Applying a four-fold Fourier transform to Eq. (5.1) yields the dispersion relationship
1:2 — Luz/c3 - 2iaow|k|y/c(l)-y = 0 (5.2)
For y 79 1, 2, Eq. (5.2) is transcendental in k and therefore cannot be solved in
closed form. Therefore, Eq. (5.2) is considered in the limit of low frequencies w
and small attenuation slopes do. In order to ﬁnd an approximate expression for
the attenuation coefﬁcient oz(w) and phase velocity C(w), the wavenumber is written
16(0)) = w/c(w) + ia(w) and inserted into Eq. (5.2)
= 0 (5.3)
w 2iwa 2 w 21'an w . y
a ( + to)
0
Assuming 01(0)) << 1 within the bandwidth of interest allows the a2(w) terms to be
neglected. Applying the binomial theorem and grouping the real and imaginary terms
yields
2 2
[”22— — “—2] + 2i [% - “WM? = 0 (5.4)
CO 63—34431
Solving for the real part gives C(w) z c0, indicating a lack of dispersion. Solving for
the imaginary part yields
at») = “gffch-y (5.5)
which yields a power law function of frequency given by Eq. (1.1). Based upon this
analysis of Eq. (5.1), the Chen-Holm equation has power law dissipation in the limit
of small (10 yet does not have a phase velocity predicted by the Kramers-Kronig given
by Eq. (4.15). Hence, the Chen-Holm model does not properly model the theoretical
and experimental phase velocity observed in power law media.
5.2.3 Causality
Since the Chen—Holm equation interpolates between the telegrapher’s equation and
the Stokes equation, both which are causal, the Chen-Holm equation is expected to
93
be causal for all 0 S y _<_ 2. To establish causality, Eq. (5.2) is solved for w as a
function of k using the analytical technique developed in [17], yielding
wiUC) = —z'aolklyc3+1 i Colk|Xy (5.6)
where
_ 2
Xy = \/1— ag|k|2y 2003’. (5.7)
Noting that Im [wi] < 0 establishes that all poles of 6' (k, w) lie in the lower half-plane.
Thus for t < O, the inverse temporal Fourier transform
G(k, t) = r1 L ’63 I (5.8)
w — w+)(w - 00—)
_ _1_ exp(—z’wt)
‘ 2w c . (5.18)
(s + (marl/(kw)? + cars? 3 + aocg‘wlkly)? + car?
Expanding the second factor as a geometric series yields
A w A
G(k, s) = Z anus, s), (5.19)
n=0
where 2 2 2
agncon+ + nyl 16'2”"
A
Gn(k, s) =
(5.20)
[(3 + ()z()c(l,+y|lc|y)2 + 63k2
The inverse Laplace transform of Eq. (5.20) is computed via Eq. (2.10), and the
]n+l°
constants /\ = cok and a = aoc6+y|k|y are deﬁned, yielding
n+1
t 2n n+2+2ny k 2yn
u( )CYO CO I l exp (_aoc(1)+ylk|yt) tkn jn(COtk) (5.21)
27171!
Gn(k:t) =
As with the Stokes wave equation, the 3D Green’s function is recovered via the three-
fold inverse Fourier transform in Eq. (2.13), yielding the inﬁnite series given by Eq.
(2.14), where each term is speciﬁed by
u(t)agn n+2+2ny
0 CO n+1 00 _ 1+1] y - 2yn+1—n -
2n+1nirr2R t /0 exp( 01000 k t) Jn(cotk)k 8111 (1:12) dk.
(5.22)
911(R1t) =
Eq. (5.22) is recognized as an inverse sine transform; unfortunately, this expression
cannot be evaluated in closed form for general 11. However, the ﬁrst term go may be
calculated explicitly by utilizing the symmetric stable distributions introduced in the
next subsection.
5.4.2 Symmetric Stable Distributions
The symmetric stable distribution wy(t) [88, 89, 106] is a special case of the stable dis-
tributions deﬁned by Eq. (8.3) with skewness ﬂ = 0. These distributions, which are
deﬁned in Appendix B, have been studied extensively [49] and possess the following
analytical properties
101
1. wy(t) is inﬁnitely continuously differentiable (or smooth).
2. wy(t) has support on (—00, oo).
3. wy(t) is symmetric (even) with respect to t and unimodal.
4. As |t| -—> oo, wy(t) ~ [fl—3"1 for 0 < y < 2.
5. fff‘;O wy(t)dt= 1.
6. wy(t) is bell shaped (the k—th derivative possesses exactly k simple zeros)
The symmetric stable distribution wy(t) has several analytical representations. The
special cases 3; = l, 3/ 2, and 2 correspond to Cauchy, Holtsmark, and Gaussian
distributions. In the special cases of y = 1/2, 2/3, 1, and 2, wy(t) is expressed in
closed-form [95, 105]
1 1 1 1
“ll/2U) = \/2—7r|tl3/2 [ <2 - 0( 27r|t|) cost-1'70)
+
2 4
102/30 )= 7—317r—I2tl exp (7) W—1/2,1/6 (W) (5231))
151(1) = «(1:12) (5.23c)
2
w2(t) = % exp (—€l) ‘ (5.23d)
where C(z) and S (2) are the Fresnel cosine and sine integrals, respectively and
Wk,m(z) is the Whittaker function [49]. For y = 3/ 2, rational approximations may
be utilized [107]. For other values of y, wy(t) is represented via the Fox H -function
[81,95] via
_ 1 1 (— 1,—1,)( 1/2,1/2) .
M(t)—EH” (I l(— 1/y.1/y>(—1/2.1/2>) ”<1 (“43’
102
_ 1 1,1 (1_1/y11/y)3(1/211/2) -
wy(t) — 9H2’2 (tl (0,1),(1/211/2) ) if y > 1 (5.24b)
Figure 5.3 displays plots of the symmetric stable PDFs wy(t) for y = 1 / 2, 1, 3/ 2,
and 2. These PDFs were generated via the STABLE toolbox [3] and veriﬁed against
the closed form expressions given by Eq. (5.23). As shown in Fig. 5.3, the densities
display heavy tails for y < 2, with the tails decaying slower as y descends from two to
zero. Also, the densities become more peaked as y decreases from two to zero, thus
yielding a larger fourth moment.
5.4.3 Leading Term (11 = 0)
As in the viscous case (y = 2), the dominant behavior is analyzed by examining the
leading order term in Eq. (2.14). For n = 0, the expresion for the spherical Bessel
function j0(z) in Eq. (2.16) is inserted into Eq. (5.22), yielding
oo
90(R,t) = 02015;) sin (kR) exp (—a0cg+1kyt) sin (00kt) dk. (5.25)
7’ 0
Eq. (5.25) is evaluated using Eq. (B.7). Evaluating the inverse transform via Eq.
(B.8) by taking a = 00%“ and b = cot yields
90(R, t) ~ u(t) )l/y [my ($5161) _ my (MN . (5.26)
~ 47rR_ 1) Terms
Evaluation of Eq. (5.22) for general 11 2 1 and y < 2 is not possible in closed form.
However, an estimate may be obtained which yields the decay properties for t > 1.
For large .2, the asymptotic form of the spherical Bessel function is employed
sin(z — n7r/2)
jn(z) -—> z (5.28)
Inserting Eq. (5.28) and utilizing [sin 2| S 1 yields the rough estimate
2n n+l+2ny 00
0‘0 C0 1+1: 2 —1
971(R, t) S 2n+l7r2an tn./0 exp (_GOCO kyt) k( y )n dk (5.29)
The integral in Eq. (5.29) is evaluated in terms of the Gamma function, yielding
2n n+1+2ny
a c —l+n—2ny y 1 _ 2
games §n+i,,2,,.R t"y‘1(aocé+yt)( V 1“(——"—f—"’-) (5.30)
Rearranging terms yields the ﬁnal estimate
Cn n/y —n(1—1/y) —n+(n—1)/y
gn(R, t) S 47rR(atc0)1/yao t co , (5.31)
where
0.. = r (2n — (n — 1m». (532)
2"“1nl7r
104
Eq. (5.31) shows that gn(R, t) --3 0 as t —» 00. Moreover, gn decays faster with
respect to t than go by a factor of t‘"(1’1/y). In addition, Eq. (5.31) scales with
n/y-th power of 010, thus becoming negligible for em small and n/ y >> 1.
5.4.5 Linear Attenuation Media (y = 1)
In the special case of linear attenuation media (y = 1), the temporal Fourier transform
is evaluated by letting s = —z'w in Eq. (4.23), yielding
(k2 — (122/0(2) — 2350mm) 6105,51) = 1. (5.33)
Factoring the dispersion polynomial in terms of w and solving for G(k, w) yields
(5.34)
where
wiUC) = -iao|k|63 i c0|k|X1 (5-35)
and X1 = ‘/1 — agcg. Performing an inverse temporal Fourier transform over w via
Eq. (5.9) and evaluating via contour integration yields
G(k, t) = 613:) exp (—aoc3|k|t) sin (col/chat) (5.36)
Inserting Eq. (5.36) into Eq. (2.13) and performing the same operations as in the
general case yields
g(R’t) = 47TR(:~
0.3 ‘
0.15
95
Figure 5.3. Plots of symmetric stable PDFs wy(t) for four cases: y = 0.5 , 1 (Cauchy) ,
1.5 (Holtsmark), and 2 (Gaussian) evaluated with the STABLE toolbox [3].
The linear attenuation media (y = 1) case is very similar to the non-causal solution
derived by Kak and Dines [108], suggesting the Kak and Dines solution is a non-causal
approximation to the Chen-Holm equation.
5.5 3D Green’s Function: Spatially Dispersive Wave Equa-
tion
The same solution technique is applied to Eq. (5.9) in this section. The spatially
dispersive wave equation, subject to an impulse point-source at the origin, is given
by
1 829 2C'0 v Q = —6(t)6(R). (5-38)
V29 — — — V
03 8t2 cos(1ry/2)c(l)—y S at
106
Fourier-Laplace transforming Eq. (5.38) from the space-time domain (R, t) to the
spectral-frequency domain (k, 3), yields
((153)? + s2/c3 + 255090.92) G(k,, s) = 1 (5.39)
Cé—y cos(7ry/ 2
Solving for G (k, s) and expanding in an inﬁnite series following the same steps used
in Section 5.4.1 yields Eq. (5.19), where
seCZn(7ry/2)(1gnc2n+2+2ny(_ ik3)2yn
[*~~ 1) given by Eq. (4.38). Since fy(t) > 0 for all t for y > 1
[89], Eq. (5.49) is non-zero for t S 0, the power law Green’s function is interpreted
as a non-causal approximation to the spatially dispersive Green’s function.
5.6 1D and 2D Green’s Functions
Asymptotic Green’s functions valid in one and two dimensions are calculated for
both the Chen-Holm wave equation and the spatially dispersive wave equation in this
section.
5.6.1 1D Green’s Functions
The 3D and 1D Green’s functions for any linear operator are related via [18]
1 691D
9 D R t — __ _ . 5.50
109
Inverting this relation yields
livl
91502.0 = — / 24R935(R,t)dR. (5.51)
—00
The mapping given by Eq. (5.51) is now applied to the asymptotic 3D Green’s
functions for the Chen-Holm equation in Eq. (5.27) and the spatially dispersive
equation in Eq. (5.47).
CHEN-HOLM EQUATION
To compute the 1D Chen-Holm Green’s function, Eq. (5.27) is inserted into Eq.
(5.51), and applying a change of variables yields
can t °°
91D(;r,t) = ——2—(—)- 1/ 7.119(1)) dv. (5.52)
(t—|a:|/c0)/(a0cot) 3"
Introducing the symmetric, stable cumulative distribution function (cdf) Wy(t) de-
ﬁned via Eq. (3.11) in Appendix B allows Eq. (5.52) to be evaluated as
c0W?) [ (t— [ail/Co )]
33,15 2 —— l — W —— 5.53
Closed-form expressions exist for Wy(t) for the special cases of y = 1 and y = 2 [49]:
1 1 _1
W1(t) = — + — tan (t) (5.54a)
2 7r
1 t
W2(t) = 5 + erf (5) (5.545)
SPATIALLY DISPERSIVE WAVE EQUATION
To compute the 1D spatially dispersive Green’s function, Eq. (5.47) is inserted into
Eq. (5.51), and applying a change of variables yields
CM(t) /°°
g D(a:, t) = —— w (v)dv. (5.55)
1 2 (t—lwl/60)/(ﬁ060t)1/y y
110
Identifying the skewed, stable cumulative distribution function (cdf) Fy(t) deﬁned
via Eq. (B.10) in Appendix B allows Eq. (5.55) to be evaluated as
915040 = @391 [1— F, ($95)] (556)
5.6.2 2D Green’s Functions
The 2D Green’s function is calculated by integrating the 3D result along the z—axis.
Putting R = \/ p2 + 22 yields
00
9250). t) = 2 f0 9350, 2. 0 42. (5.57)
which is applied to both the Chen-Holm and spatially dispersive equations.
CHEN-HOLM EQUATION
Inserting Eq. (5.27) into Eq. (5.57) and applying a change of variable yields
_ u(t) t-p/co v
925(p,t)— f wy(E————),—/,;)q(v)dv (5.58)
2740106001” —oo aocot
where the “wake function” q(v) is deﬁned via
9(2)) = 1/ \/(t - v)2 - (p/CO)2 (559)
Since wy(t) has support on (—00, 00), Eq. (5.58) is non-zero for t < p/co. Unlike the
lossless Green’s function given by
the lossy 2D Green’s function is ﬁnite at the arrival time t = p/co due to the smoothing
effect of wy(t). In addition, Eq. (5.58) is nonzero fort < p/co since wy(t) has support
on (-OO,OO).
111
SPATIALLY DISPERSIVE WAVE EQUATION
Inserting Eq. (5.47) into Eq. (5.57) yields
_ u(t) t’p/CO 1}
9250.0 — m 50401 h, f... fy(———(ﬁocot)1/y)q(v)dv (5.61)
As with Eq. (5.58), Eq. (5.61) is ﬁnite at the arrival time t = p/co due to the
smoothing effect Of fy(t). In fact, since both wy(t) and fy(t) are 00°, it follows
that the 2D Green’s functions for both the Chen-Holm and spatially dispersive wave
equations are also C°°.
5.7 Numerical Results
5.7.1 Chen-Holm Wave Equation
Figure 5.4 displays snapshots of the 3D approximate Chen-Holm Green’s function in
Eq. (5.27) for y = 1.0, 1.5, and 2.0 with em = 0.05 mm-IMHz-y at times t = 20, 30,
40, and 50 as. Like the power law Green’s functions shown in Fig. 4.4, the Green’s
functions are broader and are attenuated as time progresses. The approximate Chen-
Holm Green’s function and power law Green’s function are approximately the same
for y = 2, and both models have slowly decaying tails for y < 2. However, the Chen-
Holm Green’s functions are symmetric about the radial coordinate R = cot, whereas
the power law Green’s functions shown in Fig. 4.4 are skewed to the left. Thus,
the two models yield different behavior for 0 < y < 2, with the Chen-Holm Green’s
function lacking the characteristic skewness of a dispersive medium. Moreover, the
Chen-Holm model predicts a signiﬁcant amount of acoustic energy traveling faster
than the reference speed of sound c0, which is not observed in the power law model
developed in Chapter 3.
To explore broadband pulse propagation under the assumptions of the Chen-Holm
model, the velocity potential is computed by convolving a broadband pulse given by
112
Eq. (3.11) with the approximate Chen-Holm Green’s function. A liver-like and fat-
like medium are utilized with the same parameters used. in Fig. (4.5). The resulting
velocity potentials are displayed in Fig. 5.5 at distances R = 10 mm, 25 mm, 50
mm, and 100 mm. For small distances (R = 10 mm and 25 mm), both the power
law results in Fig. 4.5 and Chen-Holm results in Fig. 5.5 are very similar. For larger
distance (R = 50 mm and 100 mm), there are noticeable differences between the
two models. In particular, the Chen-Holm velocity potentials have a longer duration
than the power law model, which is caused by symmetric heavy tails present in the
asymptotic Green’s function given by Eq. (5.27).
5.7.2 Spatially Dispersive Wave Equation
Figures 5.6 and 5.7 compare the spatially dispersive and power law Green’s functions
in fat-like and liver-like media, respectively. Fig. 5.6, which evaluates Eqns. (5.47)
and (5.49) at observation points a)R = 1 mm and b) R = 10 mm, displays the simi-
larity between the spatially dispersive and power law Green’s functions. Both Green’s
functions quickly rise to a peak at time t as 0.69 us, after which they both decay very
slowly as t —> 00. At R = 1 mm, the power law Green’s function underestimates the
peak at 0.69 us. Also, the power law Green’s function is slightly delayed relative to
the spatially dispersive Green’s function to the left of the peak. However, for times
occurring after the peak, the spatially dispersive and power law Green’s functions are
indistinguishable. At R = 10 mm, the two Green’s function display similar behavior,
although there is closer correspondence between the two models for times preceding
the peak at t z 6.2 ps.
Fig. 5.7 displays a greater disparity between the power law and spatially dispersive
models, especially for observation points near the source R = 1 mm and small times.
Although both Green’s function display the same qualitative behavior, the power law
model is noticeably delayed relative to the spatially dispersive model. This delay
113
4x10”? 4x10?”
-y=2 —y=2
_---y=1.5 _---y=1.5
3 """" Y= 3 """" Y=
—L
d
Green's Function g(R,t)
N
O
Green's Function g(R,t)
N
O
\\
‘0
26_ 28 80 32 34 36 40 .42 44.46 48 50 52
radial coordinate R (mm) radial coordinate Fl (mm)
(a) t = 20 118. (b) t = 30 113.
x10:3
—Y=2
---y=1.5
....... y._..1
A
A
‘9
9°
—L
Green's Function g(Fl,t)
N
Green's Function g(R,t)
-* N
-(
W
\I
o
s‘
..L A .
54 58 , 62 66 68 2 76 "A 80 84
radial coordinate H (mm) radial coordinate R (mm)
O
O
(c)t=40ps. (d)t=50ps.
Figure 5.4. Snapshots of the 3D approximate Chen-Holm Green’s function in Eq. (5.27)
for y = 1.0, 1.5, and 2.0 for (10 = 0.05 Np/mm/MHz”. Snapshots of the Green’s function
are shown for t = 20, 30, 40, and 50 118.
114
U‘l
velocity potential (mm/us)
8 o
0.‘
'velocity potential (mm/11s)
5. o
x1Q .
1 . —fat .
5': ---liver
. i . .
. 6 7 8
t1me(us)
(a) R = 10 mm
x103 .
—fat
5 ---liver.
iii?"
is:-
53
i
30 82 34 36
time(p.s)
(c) R = 50 mm
velocity potential (mm/us)
velocity potential (mm/us)
x 10"] +
2’ —fat
f; ---liver
1' EE 1
c .5: 5'! ”TV A
55 =5
-1- ES
14 16 18
time ((16)
(b) R = 25 mm
x 10"4 '
2. —fat
; ---Iiver
1 . E .
0——~’.§’:.’ ‘
iii \l
-1- as ‘
-2. l
62’ 64 , 66 68 70
time (us)
(d) R = 100 mm
Figure 5.5. Velocity potential produced by a point source excited by a broadband pulse
deﬁned by Eq. (3.11) in two media governed by the Chen-Holm equation: 1) a fat-like
medium with y = 1.5, co = 1.432 mm/ps, and a0 = 0.086 Np/MHzl's/cm and 2) a liver-
like medium with y = 1.139, co = 1.569 mm/ps, and a0 = 0.0459 Np/MHzl'l39/cm. The
velocity potential is displayed at radial distances 10 mm, 25 mm, 50 mm, and 100 mm.
115
arises from the scaling parameters of the Green’s functions. The spatially dispersive
model has a scale parameter given by (ﬂocot)1/y, whereas the power law model has a
scale parameter given by (16013)” 3’. Since the spatially dispersive Green’s function has
a time-dependent scale parameter, Eq. (5.47 ) possesses greater temporal localization
for t < R/co. Thus, for t < R/co, the spatially dispersive Green’s function decays
faster for times preceding R/co relative to the power law Green’s function, which has
a time-independent scale.
As discussed in Chapter 4, the power law Green’s functions display heavy-tailed
asymptotic behavior due to the slow decay of fy(t), which decays with 0(t’y‘1).
This slow decay is shared by the spatially dispersive Green’s function, which has al-
most identical long—term behavior as the power law Green’s function. As Figs. 5.6
and 5.7 show, the spatially dispersive Green’s function also displays this heavy-tailed
behavior, which yields a slowly decaying wake that trails the primary wavefront. Since
large times correspond to the low-frequency information, the similar behavior of the
power law and spatially dispersive models indicates that both models correctly model
low-frequency power law behavior. On the other hand, the early—time behavior of
the power law and spatially dispersive models differs slightly, due to differing high-
frequency behavior. The power law model has a power law attenuation coefﬁcient
dependence my over all frequency bands, whereas the spatially dispersive model has
power law behavior only for low frequencies. At high frequencies, the spatially disper-
sive attenuation coefﬁcient grows more slowly than 02. As displayed in Figs. 5.6 and
5.7, this differing high-frequency behavior is manifested in different early-time behav-
ior of the power law and spatially dispersive Green’s functions. In fact, this slower
growth of the spatially dispersive attenuation coefficient is a necessary condition for
causality due to the Paley-Wiener condition.
116
5.8 Discussion
5.8.1 Stochastic Model for the Chen-Holm Wave Equation
The Chen-Holm approximate Green’s function given by Eq. (5.27) is derived by con-
sidering tissue heterogeneity on the macromolecular level. One possible explanation
of tissue absorption is Class 0 scattering on the length scale of 1 pm. At this length
scale, classical absorption is modeled by the Stokes wave equation with a variable
viscous relaxation time
1 829 6
v29 — 7? + 7(r)—V29 = —6(t)6(R), (5.52)
Co
at
where the viscous relaxation time 7(r) is a function of position to account for tissue
heterogeneity. However, the speed of sound is taken to be a deterministic constant.
That is, the fundamental loss mechanism is assumed to be a variable viscosity [109]
with a constant background c0. However, at larger spatial scales, the tissue is assumed
to be homogeneous. That is, there exists a volume such that the small scale variations
of 7(r) average out to a constant. The size of this volume is taken to be on the order
of Rayleigh scattering (10 — 100 pm). Equivalently, the relaxation time consists of two
components: 1) a constant '7 valid on the scale of Rayleigh, diffractive, and specular
scattering and 2) a high frequency component 7’ (r) valid on the micro-scale. Since
the oscillations of ’y’ (r) are on a small scale, the spatial Fourier transform 7’ (k) is
approximately zero for |k| less than some constant.
Assume that 7(r) is a random variable with distribution p(y) satisfying the fol-
lowing conditions: 1) 7(r) is strictly positive (e.g. 19(7) = 0 if 'y S 0), and 2) 7(r)
is an inﬁnitely-divisible distribution [89], resulting from the summation of many sim-
ilar random variables. Condition 1 is required since the coefﬁcient of viscosity in a
dissipative medium is strictly positive. Condition 2, although not required, is desired
since Class 0 scattering is modeled as the summation of many independent scatter-
ing events. In other words, 7(r) is the sum of a large number of independent and
117
identically distributed (i.i.d.) random variables. The only distributions that satisfy
both conditions are the maximally skewed stable distributions with PDFs given by
Eq. (4.41). A spread parameter A > 0 is incorporated into Eq. (4.41), yielding the
desired distribution of '7
p( )= —/2sz11 ( ,\1 /z) . (5.63)
Recalling the asymptotic Green’s function to the Stokes wave equation, the averaged
Green’s function g(R, t) is computed by taking the expected value:
86.0 a Elg 1 R)...
‘ 47m. 000(tT—/2)1/2f2((17/2)1/2)Al/Zfz(A1/z)d
The integration in the last line is computed in [89], yielding
_ u(t) t—R/co
90” t): 47rR(At/2)1/(2z) 2w" ((At/2)1/(22)) (5'65)
where wy(t) is the symmetric stable distribution deﬁned by Eq. (B.7). Eq. (5.65)
has the same form as the Chen-Holm Green’s function given by Eq. (5.27) with a
spread parameter A = 2a0c0 and a power-law exponent of y = 22:. The evaluation
of the integral in Eq. (5.65) is particularly straightforward when 2 = 1 / 2, yielding a
Cauchy distribution in Eq. (5.65) with a linear frequency dependence.
Thus, the assumption of a spatially varying viscous relaxation time yields the
Chen-Holm Green’s function. Note, however, that the sound speed Co was assumed
constant on all scales, which neglects the effect of dispersion. This analysis shows
that the Chen-Holm model correctly accounts for absorption due to viscous micro-
heterogeneity but neglects the dispersive effect of micro—heterogeneity. This model
agrees with Eq. (5.4), which predicts a frequency-independent speed of sound for the
Chen-Holm model. 2
118
5.8.2 Stochastic Interpretation of the Spatially Dispersive Wave Equation
Like the power law and Chen—Holm Green’s functions, Eqns. (5.47) and (5.48) are
probability distributions. Following the analysis in Section 4.5, the convolution of
a pulse u(t) with Eq. (5.47) may be expressed in the form of Eq. (4.66), where
the random variable T in Eq. (4.66) is distributed by the shifted and scaled stable
distribution
_ u(t) ~ 7' — R/co
gL(t,T — R/CO) _ (OOCOt)1/yfy ((aocot)1/y) (5.66)
where P(T(t) S 1') = f—Too gL(t,T' — R/c0)d7". Thus, T = T(t) is a continuous
time random process which is interpreted as randomized delays. Like the power law
model, the process T(t) consists of two components: 1) a bulk delay R/co and 2)
additional, random delays due to the micro-heterogeneity of the intervening medium.
As time evolves, the spread parameter in T(t) grows, indicating a greater probability
of encountering delays far from R/co. Since ﬂ, is skewed to the right, there is a
much larger probability that the outwardly propagating spherical wave encounters
delay larger than R/cO. Note that the same stochastic interpretation applies to the
1D and 2D spatially dispersive Green’s functions given by Eqns. (5.56) and (5.61),
respectively.
119
3.5
3_
23:25
C»
C
8 2-
O
C
3
u-1.5~
sn
5
9 1'
o
0.5
—spatially dispersive
- - -power law
0.08
0.07-
3;? 0.06-
‘5’.
g 0.05-
‘3
g 0.04»
5% 0.03~
5 0.02-
0.01-
8.
90
a;
a:
0.68 017 0.12 0.74
time (us)
(a) R: 1 mm.
-spatially dispersive
- - -power law
6.9 % 7:1 7.2 7.3
time (us)
(b) R = 10 mm.
Figure 5.6. Comparison between the spatially dispersive Green’s function given by Eq.
(5.47) and the power law Green’s function given by Eq. (5.49) in a fat-like medium (y = 1.5)
with attenuation coefﬁcient 0.086 Np/cm/MHzl's. The two Green’s functions are evaluated
as a function of time t at observation points a) R = 1 mm and b) R = 10 mm.
120
—spatially dispersive
15 ' - - -power law 4
.1314.
E
0124
8
1310-
S
11. 85
SD
5 6- '
9
(D 4-
2’ J
oo
.62 0.625 0.63 0.635 0.64 0.645 0.65
time (Its)
(a) R = 1 mm
0.25 . .
—spatially dispersive
- - -power law
I; 0.2”
If
6
,5 0.155
‘6
C
3
u.
g) 0.1 -
C
8
‘5 0.05-
6925 63 6.85 6.4 6.45 6.5
time ((15)
(b) R = 10 mm
Figure 5.7. Comparison between the spatially dispersive Green’s function given by Eq.
(5.47) and the power law Green’s function given by Eq. (5.49) in a liver-like medium (y =
1.139) with attenuation coefﬁcient 0.0459 Np /cm/ MHz1'139. The two Green’s functions are
evaluated as a function of time t at observation points a) R = 1 mm and b) R = 10 mm.
121
CHAPTER 6
Fractal Ladder Models and the
Caputo-Wismer Equation
6. 1 Introduction
The power law models discussed in Chapters 4 and 5, such the Szabo [19] and Chen—
Holm [20] wave equations are phenomenological insofar as these equations are not de-
rived from an underlying physical principle and / or constitutive equation (as opposed
to the Stokes wave equation studied in Chapter 2). Rather, the order of the fractional
derivative and coefficients are ﬁtted to measured attenuation measurements, but are
not derived from equations of state, continuity, and momentum. Although the FPDE
model may accurately predict wave propagation in biological media, the underlying
FPDE does not provide a physical picture or mechanism of the underlying absorption
processes. Therefore, this chapter focuses on deriving power law attenuation models,
and the corresponding FPDE, from a fractal description of tissue rooted in observed
biophysics. I
In order to derive a biophysically-based constitutive equation, the viscoelastic
propertias of both cells [110—112] and bulk tissue [113, 114] should be considered.
In particular, such a constitutive equation should account for 1) the tissue micro-
heterogeneity and 2) the observed viscoelastic response of tissue. Within the vis-
122
coelastic and biomechanics communities, lumped parameter networks, such as the
Maxwell and Voigt models [115], are commonly employed to generate tractable mod-
els that capture the salient properties of a material [112]. These lumped parameter
networks have been extended to include inﬁnite ladder networks consisting of alternat-
ing elastic and viscous elements [115—118], which generate time-fractional rheological
constitutive equations [119] for polymers. As discussed in the above references, the
time-fractional derivative in the constitutive equation captures the l) elastic, 2) vis-
cous, and 3) self-similar pr0perties described by these inﬁnite networks. To date,
however, these ladder networks have not been applied to the propagation of pulsed
ultrasound through soft tissue.
In this chapter, a time-fractional constitutive equation, based on the lumped pa—
rameter methodology, is proposed to model the dissipative properties of soft tissue.
This constitutive equation is used to derive FPDE models (the Caputo-Wismer and
Szabo wave equations) found in the literature [19, 21, 36, 120] . In Section 6.2, a con-
stitutive equation is formulated by considering fractal networks of springs and dash-
pots [115, 117, 118]. Using this constitutive equation, the Caputo-Wismer equation is
derived from basic conservation laws in Section 6.3 for linear macro-homogeneous me-
dia and inhomogeneous media. In Section 6.5, the ladder model is analyzed in terms
of previous biomechanical and fractal models, and extensions to nonlinear media are
discussed.
6.2 Fractal Ladder Model and Fractional Constitutive Equa-
tion
This section introduces a lumped parameter, fractal ladder network to model the
stress-strain relationship in biological media. From this fractal ladder network, a
time-fractional derivative constitutive equation is derived, thus providing a physical
123
basis for time fractional FPDE such as the Caputo—Wismer equation [21, 36].
A linear constitutive equation postulates a functional relationship between the
time-dependent stress tensor TU-(t) and strain tensor 5,-3- (t) (or, equivalently, the ve-
locity gradient Bui/ 8231-) via a differential, integral, or integro-differential relationship
that satisﬁes the principle of superposition. Familiar examples of constitutive equa-
tions, such as as Hooke’s Law for an elastic solid and Newton’s Law for a viscous
fluid, fail to predict the behavior of many viscoelastic solids; therefore, generalized
viscoelastic models, involving fractional derivatives and integrals, have been proposed
[119, 121]. This section proposes an constitutive equation for biological media using a
fractal ladder network as a lumped—parameter model. A qualitative biological model
is ﬁrst postulated, followed by the basic theory of Viscoelasticity.
6.2.1 Biological Motivation
A mechanical model for the loss mechanism in mammalian biological tissue is consid-
ered in this section. As will be shown, this model satisﬁes a power law attenuation
coefﬁcient 01(0)) given by Eq. (1.1) over an appropriate bandwidth of frequencies for
power law exponents 1 < y S 2. Speciﬁc tissue, such as breast, consist of hierarchical
arrangements of elastic and fluid-like components. For instance, breast consists of
interlobular stroma, which has ﬂuid like properties. Embedded in the stroma are
mammary lobules, which have elastic membranes. The mammary lobules contain
intra-lobular stroma (viscous) and smaller elastic structures such as basal lamina.
This heterogeneity continues down to the cellular and sub—cellular levels. Individual
tissue is highly heterogeneous and composed of over a hundred distinct cell types
These tissues consists of aggregates of cells suspended by an ﬂuid-like extra-cellular
matrix (ECM). The ECM is often modeled as an aqueous solution of viscoelastic poly-
mers, which possess both solid and ﬂuid-like properties. Individual cells are modeled
as elastic membranes containing ﬂuid-like cytoplasm [112]. Within the cytoplasm are
124
distributed organelles, such as the nucleus, endoplasmic reticulum, and lysosomes,
which in turn have an elastic membrane containing with a ﬂuid-like interior [4].
This hierarchical arrangement is displayed at several scales in Figure 6.1. For
reference, consider a medium with average sound speed c0 = 1.5 mm/ps excited at
7 .5 MHz, yielding an acoustic wavelength of A = 200 um. Panel a), which is on the
scale of an acoustic wavelength /\ ~ 200 pm, contains an ensemble of mammalian cells,
each bounded by an elastic membrane, suspended in a viscoelastic ECM. Both the
ECM and cytoplasm consist of complex polymers (e.g. collagen) dissolved in a viscous
ﬂuid, thus yielding a viscoelastic material. In panel b), which is on the scale of A/ 10 ~
20 pm, an individual cell is shown at a higher level of magniﬁcation. Inside the elastic
membrane is the cytoplasm, which has similar viscoelastic properties to the ECM.
Panel c) displays the cell nucleus on the scale of A/4O ~ 5 pm, consisting of a double
membrane, a ﬂuid-like nucleoplasm, and an elastic nucleolus in the interior, which
contains chromatin [4]. From this brief description, both the viscoelastic properties
and self-similar, or fractal properties, of tissue are reasonable simplifying assumptions
for a mathematical model. Applying this biological picture, tissue may be visualized
as a recursive arrangement of ﬂuid substrates containing elastic membranes. Similar
models, known as liquid drop models, have been proposed within the biomechanics
community to describe the deformation of eukaryotic cells [110, 112]. These liquid
drop models typically model the cell membrane as a cortical layer with a characteristic
surface tension, whereas the cytoplasm in the cell interior is modeled as a viscous,
incompressible ﬂuid with a characteristic coefﬁcient of viscosity. The cell nucleus may
also be included as an additional elastic component embedded within the viscous ﬂuid.
In the following fractal model for the viscoelastic properties of tissue, the liquid
drop and compound liquid drop models discussed in Ref. [112] are extended to
include an inﬁnite number of nested elastic membranes, each containing a viscous,
125
“g. u
v_——--—_--—-
\ : : \ : :/ .
~200 pm Cell / ~ 20 pm Nuclear ~ 5 pm Nucleolus
Membrane Membrane
Figure 6.1. Schematic showing tissue structure at three different spatial scales (tissue,
cellular, and sub-cellular). The ﬁrst panel (A ~ 200 pm) displays an ensemble of mammalian
cells, each bounded by an elastic membrane, shown suspended in viscoelastic ECM. The
second panel (A/ 10 ~ 20 pm) displays an individual cell at a. higher level of magniﬁcation.
The third panel (A / 40 ~ 5 pm) displays the cell nucleus, consisting of a double membrane, a
ﬂuid-like nucleoplasm, and an elastic nucleolus in the interior, which contains chromatin [4].
Although the speciﬁc biological structures vary at each successive spatial scale, the essential
features are the same: ﬂuid substrates containing elastic compartments. This self-similar
pattern forms the basis for the fractal structure shown in Fig. 6.2.
126
compressible, ﬂuid. By allowing an inﬁnite number of layers, larger structures (e.g.
ensembles of cells) and smaller structures (e.g. cell nuclei) may be included within
the lumped parameter framework. By extending the number of structural compo-
nents indeﬁnitely, the self-similarity of biological media is revealed. This topology
is depicted in Figure 6.2, where the alternating elastic and viscous components are
visualized as a self-similar hexagonal packing of spheres within spheres. Each of the
three panels in Fig. 6.2 corresponds to the three panels in Fig. 6.1. That is, the left
panel of Fig. 6.2 models the tissue level, the center panel models the cellular level,
and the right panel models the sub-cellular level. Comparing the three panels of Ref.
6.2, the self-similar nature of this fractal structure is immediately evident. Hence,
this fractal model captures the three essential features of the biological picture shown
in Fig. 6.1: 1) elastic membranes 2) ﬂuid compartments, and 3) self-similarity over
a range of spatial scales.
In order to capture these three salient properties of biological media, a fractal
network of springs and dashpots is proposed in Figure 6.3. The elastic membranes
displayed in Fig. 6.2 are represented by springs with Young’s moduli E, while the
viscous compartments are represented by dashpots with coefﬁcients of viscosity 1].
Each level in Fig. 6.3 corresponds to each pair of elastic/ viscous layers shown in Fig.
6.2 and each level of magniﬁcation depicted in Fig. 6.1. Similar fractal networks have
previously been used in lumped parameter models of viscoelastic systems [116, 122]
such as cross-linked polymers and gels [118].
6.2.2 Stress-Strain Relationships
In this section, the well—known constitutive equation for a viscous Newtonian ﬂuid
relating the strain 6,-3- and stress Ti]- tensors is generalized. The strain tensor is deﬁned
via eij = wi/asj, where w, denotes the i-th component of displacement. The stress
tensor Tij denotes the i-th component of stress per unit area along a surface normal
127
to the j-th direction, where l S 2', j _<_ 3 denote the ac, y, and 2 directions. The viscous
stress tensor Tij for a compressible, viscous ﬂuid with pressure p and velocity u is
given by [14] . 3
Ti]- = —p6,-j — 311V - udij + p (3%; + $3) (6.1)
where p is the coefﬁcient of shear viscosity, 6,5 is the Kronecker delta operator, and
i, j =1,2,3. For homogeneous gases, a may be computed via kinetic theory. For more
complicated ﬂuids, )4 must be measured experimentally. Eq. (6.1) contains three
terms: 1) an elastic term involving the thermodynamic pressure p, 2) an isotropic
frictional term and 3) a shearing term. In Eq. (6.1), the coefﬁcient of bulk viscosity
is assumed to be zero. Eq. (6.1) describes the viscous behavior of monoatomic gases
quite well, yet does not properly describe the dissipative behavior of more complex
polyatomic molecules.
Although Eq. (6.1) may describe the stress-strain relationship on a sufficiently
small micro-scale with a variable viscosity [1, Eq. (6.1) fails to predict dissipative be-
havior on the observable macro-scale in most biological media. To obtain a constitu-
tive equation on a scale commensurate with an acoustic wavelength (macro-scale), Eq.
(6.1) is averaged over a suﬂiciently large volume to achieve an constitutive equation
with constant coefficients. This averaging, or up-scaling procedure, should account
for the signature micro-heterogeneity and hierarchical micro-structure of biological
media. One simple up-scaling procedure employs a lumped parameter model dis-
cussed earlier, wherein the individual components of the medium (cells, membranes,
organelles, etc.) are represented via hierarchical arrangements of springs and dash-
pots. By assuming tissue is isotropic on the macro-scale, the 3D constitutive equation
can easily be generalized from the 1D expression.
On the macro-scale, the normal stress is decomposed according to
Tij(ry t) = -P(r, 19% + 0ij(rat) (65-?)
where 0(r, t) is the component of stress responsible for dissipation. The dissipative
128
component of Eq. (6.1) is generalized to include memory effects by relating each com-
ponent of stress and strain to a causal, stationary, hereditary integral (or Boltzmann
superposition integral) [115]:
t 2
aij(r,t) =£mg( (t—t’ )[-— Egg-[(r,t )6z-j+e,'j(r, t)+ej,-(r, t) dt' (6.3)
where g(t) is a relaxance, or memory, function [115] which relates the present state
of the material to its previous history. Since Eq. (6.3) is a convolution integral, the
stress—strain relationship becomes multiplication in the Laplace domain. In 1D, Eq.
(6.3) is represented in the Laplace domain as
5(1‘, 8) = 9(8)€(r, 8) (5-4)
6.2.3 Fractal Ladder Model
From a mechanical point of view, the combined viscous and elastic components are
modeled as springs and dashpots, respectively. Springs, which model energy storage,
represent the nested elastic membranes shown is Fig. 6.2, while dashpots, which
model dissipation, represent the viscous components, such as cytoplasm. The self-
similar structure is realized as a fractal ladder in Fig. 6.3, which provides a lumped
parameter description of the geometric model shown in Fig. 6.2. For simplicity, all
the springs have the same spring constant, or Young’s modulus, E and the dashpots
have the same coefficient of viscosity 17. Utilizing basic mechanics, the stress-strain
relationship given by Eq. (6.4) is evaluated an an inﬁnite, periodic continued fraction:
1
—l 1
E + —s+—j——
’7 E—I+...
= [713: E-linsiE—l: ' ° ]
—n/Es + \/n/Es(77/Es + 4)
2/E ’
6(a) ——- as + (6.5)
129
Elastic Membrane (Young's Modulus E)
Q-o
(D.
Q.
I
~200 um
Viscous Fluid (coefficient of viscosity 7])
Figure 6.2. Layered fractal model for biological tissue based on the schematic shown in
Fig. 6.1. The ﬁrst panel displays an inﬁnite number of thin elastic membranes with Young’s
modulus E alternating with viscous compartments that have coefﬁcients of viscosity 1]. The
second panel zooms in on the ﬁrst panel, thus showing the self-similar layered structure.
130
“Ll
1"I
Figure 6.3. Fractal ladder model for tissue micro-structure. The continuum model de-
picted in Fig. 6.2 is described with springs with Young’s modulus E and coefficients of
viscosity 7).
131
where the continued fraction has been evaluated in closed form [118]. For 317/ E < 1,
the binomial approximation is applied, yielding the low-frequency approximation
6 z nEs. (6.6)
Inserting Eq. (6.6) into Eq. (6.4) and performing an inverse Laplace transform by
applying Eq. (A.3) from Appendix A yields
81/26
where the fractional derivative operator is deﬁned by Eq. (A.2) in Appendix A.
Generalizing to 3D (under the macro-isotropic assumption) yields
E812/ 31/2
0ij = -— -\/’f] Ear/23.21%]. (lij + TIE—ﬂ atl/z (623’ + Eji) (6.8)
Constitutive equations similar to Eq. (6.7) have been previously proposed for vis-
coelastic materials within the geology community [36]. To cast Eq. 6.8 as a function
of velocity u, the following relationship between the i-th component of velocity u,-
and strain EU is utilized
Bu,- _ 86,-]-
axj — at
(6.9)
yielding
2 Eva—U2 6-1/2 62.,- an,
TU: 1.3-(19+ VTIE ——1—/2Vat_ ) 6ij + "Eat—1” (61“) + 8221'). (6.10)
Thus, a time-fractional stress-strain relationship follows from the fractal ladder model
originally proposed in Ref. [118].
6.2.4 Recursive Ractal Ladder Models
The ladder model can also be considered as a fundamental mechanical component (a
“springpot” [123]), allowing more complicated fractal networks, or recursive ladders,
to be constructed. For instance, consider a ladder model constructed by replacing the
132
viscous damper in Fig. 6.3 with a fractal ladder, producing the arrangement shown
in Figure 6.4. That is, a simple ladder with modulus K 37, where K = \/E_ﬁ and 7 =
1/2 is embedded within a larger ladder along with springs with elastic coefﬁcients
E. Computation of g(s) for this model using the low-frequency approximation given
by Eq. (6.6) yields g(s) as E3/ 4111/ 431/ 4. This construction may be generalized by
embedding N — 1 level ladders alternating with springs to create an N—level ladder-
N l N l N 1
z E1—1/2 + 171/2 + 31/2 + . By performing an
spring network, yielding g(s)
inverse Laplace transform, fractional derivative stress-strain relationships of order
1 / 2, 1/4, 1/8, are generated. As the depth of the ladder increases (N —> 00),
{1(3) -+ E, yielding a purely elastic response.
A similar recursive mechanical network is constructed by replacing the springs
with dashpots in the recursive ladder-spring network. Replacing the springs in
Fig. 6.3 with ladders yields 9(3) z E1/4n3/4s3/4. This model may also be gener-
alized to dampers alternating with N — 1 level ladders (or springpots), producing
9(3) z E” 2N+1n1'1/2N+131‘1/2N+1. By performing an inverse Laplace transform,
fractional derivative stress-strain relationships of order 1 / 2, 3/ 4, 7 / 8, are gener-
ated. As the depth of the ladder increases (N —+ 00), 6(3) —> 17.9, yielding a purely
viscous response.
In order to generate fractional derivatives of all orders within the unit interval,
fractal ladders containing alternating damper-M ladders and N ladder-spring net-
works are considered. Computing the frequency-domain modulus for this network
yields
1+1_111_1+111_1+1
27W+l 2N+l n2 2M+l 2N+l 3’2 2M+l 27V+l (
t0>
E
22
t1]
Lav—I
A
6.11)
Putting B = % (l — 5171:1- + 21713) and performing an inverse Laplace transform
(6.12)
(KY)
(KY)
(Kill
(K?!)
ml ..
Figure 6.4. Recursive fractal ladder model for tissue micro-structure which generalizes
the ladder shown in Fig. 6.3. The dashpots in the simple ladder model are replaced with
springpots with a modulus K .97, where each of the springpots are built from recursive
arrangements of ladders, springs, and dashpots.
134
which is a generalization of Eq. (6.7) for all 0 < ,8 g 1. For instance, putting
M = O and N = 1 yields ,6 = 3/ 8. Hence, arbitrary rational-order stress-strain
relationships of the form Eq. (6.12) are modeled by alternating ladder-spring/damper—
ladder networks. As 6 —+ 0, the response becomes more elastic, and as 6 —-> O, the
response becomes more viscous. Thus, Eq. (6.12) interpolates between Hooke’s law
for an elastic solid and Newton’s law for a viscous ﬂuid.
6.2.5 Constitutive Equation
The generalization of Eq. (6.12) to 3D yields
3
01'3“ = - 5E0 1756—”- ECU 515+ E0 "€37 (EU + éji) (6.13)
i=1
Eq. (6.13) consists of a two viscoelastic terms involving the fractional derivative of
strain with respect to time, where the fractional derivative term is responsible for the
coupled processes of attenuation and dispersion. For [3 < 1, Eq. (6.13) displays a
temporal non-locality commonly utilized in phenomenological Viscoelasticity [36, 119].
Theoretical justiﬁcation for constitutive equations similar to Eq. (6.13) were also
established [121] on the basis of dilute solutions of polymers [124] within a Newtonian
ﬂuid medium. In the Bagley-Torvik—Rouse theory [121, 124], a polymer solute in a
homogeneous, Newtonian solvent was considered. The polymer solute is modeled as
a chain of sub-molecules, each with a characteristic relaxation time. These relaxation
times are related to physical and chemical pr0perties of the macromolecules, such as
molecular weight and concentration.
By identifying a generalized coefﬁcient of viscosity
1—
p = E0 37763 (6.14)
and repeating the steps performed in the ﬁ = 1/2 case, the following averaged con-
135
stitutive equation involving velocity gradients is computed:
2 63-1 65-1 Bu,- an,-
Tij : —P6ij - Bylaw-TV - U5z'j + #atﬂ_1 (axj + 55;) (6.15)
Eq. (6.15) may also be written in dyadic notation as
—l aﬁ—l
I = —pl — War—IV - u; + [lam-:1- (Vu + uV) (6.16)
where underlines denote rank-2 tensors and l is the rank-2 identity tensor. The
dyad Vu has components Bu,- / 6123', whereas the dyad uV = (Vu)T has components
(911]- /8:1:,-
For 6 = 1, Eqns. (6.15) and (6.16) reduce to the viscous stress tensor for a com-
pressible, Newtonian ﬂuid given by Eq. (6.1). Eq. (6.15) contains three terms: 1)
an elastic term involving the thermodynamic pressure p, 2) an isotropic frictional
term and 3) a shearing term. The frictional terms tend to diffuse momentum through
the ﬂow. In a viscous ﬂuid (,6 = 1), the frictional term involves only a spatial
derivative and is purely local. For a homogeneous ﬂuid with simple molecular struc-
ture, this relation properly accounts for momentum diffusion. For biological tissue,
however, viscous loss does not properly account for observed dissipation. Tissue is
both heterogeneous and has a complex molecular structure, which is best modeled
as a viscoelastic media. Physically, momentum may diffuse faster and/or slower in
some directions due to the heterogeneity of tissue. To describe these effects, the
local constitutive relationship is generalized by a global relation that incorporates
memory into the ﬂow. The temporal operator 3(3—1/(9135—1 is the Riemann-Liouville
fractional derivative for O < ,6 < 1. Since the order of differentiation is negative,
the Riemann-Liouville fractional derivative is equivalent to a fractional integration.
Thus, the viscous shearing term in a standard Newtonian ﬂuid is replaced with a
memory term which relates the stress at time t to the entire history of the velocity
gradient. Similar hereditary constitutive equations, which utilize time-fractional and
time-convolutional operators, are widely used in theoretical Viscoelasticity [125, 126].
136
6.3 Derivation of the Caputo-Wismer Equation
This section demonstrates how the fractal ladder model and fractional constitutive
equation pr0posed in Section 6.2 lead to previously pr0posed F PDE for dispersive
wave propagation in biological media. In particular, the 3D Caputo-Wismer equation
[21], which models power law attenuation via a time-fractional derivative, is derived
for in a linear, compressible, isotropic, and adiabatic medium governed by Eq. (6.15).
Both homogeneous and inhomogeneous media are considered.
6.3.1 Homogeneous Media
Longitudinal wave motion is considered in a homogeneous medium with density p0,
sound speed c0, generalized coefﬁcient of viscosity )1, and exponent )6. To simplify the
analysis, the adiabatic hypothesis, whereby entropy is assumed constant, is adopted,
which neglects the dissipative effects of thermal conduction. In addition, nonlinearity
is ignored by assuming pressure and density are linearly related.
The constitutive equation given by Eq. (6.15) is complemented by 1) the linearized
Cauchy’s equation, 2) the linearized, adiabatic equation of state, and 3) the linearized
equation of continuity. The linearized Cauchy equation, which neglects the inertial
term u - Vu, restates Newton’s second law of motion for an arbitrary continuum. In
dyadic notation, the Cauchy equation is written
Bu
— = V - T. .17
PO at _ (6 )
The linearized, adiabatic equation of state is given by
_ 2
p — COP) (6.18)
where p denotes excess density and Co is the adiabatic speed of sound. Finally, the
linearized equation of continuity is given by
5p
a + pOV ° 11 — 0, (619)
137
which accounts for local mass conservation. First, the divergence of Eq. (6.15) is
taken and inserted into Eq. (6.17) using the identities [127] V - (Vu) = V (V - u) and
V-(uV) = V(V-u)—V x V xu, yielding
an 86" 4
p09? — —Vp+ MW (§V(V ' 11)- V X V X 11) . (6.20)
Any vector ﬁeld may be decomposed into a longitudinal component, for which V x u =
0, and a transverse component, for which V - u = 0. Since the transverse component
attenuates rapidly away from any boundaries, only the longitudinal component is
considered, thus causing the V x V x u term in Eq. (6.20) to vanish, yielding
Eq. (6.21) is interpreted as a force balance equation, where the left hand side is
mass per unit volume multiplied by acceleration. The right hand side consists of two
terms: 1) a pressure gradient responsible for longitudinal vibrations and 2) a frictional
term responsible for attenuation and dispersion. The frictional term consists of 1) a
diffusive operator which tends to average momentum variations in the ﬂow and 2) a
temporal integration operator which incorporates the past history of the ﬂow into the
ﬂow’s present state.
The divergence operator is applied to both sides of Eq. (6.21), yielding
3V ' u — _V2 +6. ﬂ
v . (V(V . u)). (6.22)
In order to eliminate particle velocity in favor of pressure, Eq. (6.19) is combined
with Eq. (6.18) to produce
1 0p
0
Inserting Eq. (6.23) into Eq. (6.22) and interchanging spatial and temporal differen-
tiation, yields
1 32p 2 4,11. 6B 2
—— = V ——.—V 6.24
6(2) 6t2 19+ 3C3p0 atla p ( )
138
Identifying the relaxation time
4
75 = 2“ (6.25)
3CoP0
yields the Caputo—Wismer equation [21]
162
V21) — ——I—) + Ty—
633112
where y = 5 +1 is the power law exponent. Frequency-dependent loss is incorporated
1 3161
2 _
via the Riemann-Liouville fractional derivative deﬁned by Eq. (A.1) in Appendix B.
For y = 2, Eq. (6.26) reduces to Stokes wave equation given by Eq. (2.1), which
models wave propagation in a homogeneous, viscous medium. As shown in Ref. [21],
for y > 1, Eq. (6.26) admits an attenuation coefﬁcient with a power-law dependence
in the low-frequency limit. For y = 1, however, the loss operator reduces to a spatial
Laplacian, thereby failing to model power law attenuation. For this reason, the
exponent is restricted to l < y s 2. The 1D version of Eq. (6.26) has previously
been studied within the context of Viscoelasticity [36]. Hence, the Caputo—Wismer
equation arises naturally as the wave equation modeling small-amplitude, longitudinal
disturbances in media governed by the constitutive equation given by Eq. (6.15). In
addition, the parameter T appearing in Eq. (6.26) is given physical meaning by Eq.
(6.25), which depends both on the micro-structural properties of the medium E0 and
770, as well as the macroscopic properties on and p0.
6.3.2 Inhomogeneous Medium
In general, biological media is inhomogeneous on both the microscopic scale (~ 1 11m)
and macroscopic scale (~ 1 mm). The fractional derivative operator in Eq. (6.26)
accounts for the effect of micro-heterogeneity on the macroscopic scale, yet does
not account for macro-heterogeneity responsible for coherent scattering (e.g. tissue
boundaries). To incorporate macro-heterogeneity, the material properties of density
p0(r), adiabatic compressibility n0(r) and generalized viscosity p(r) are assumed to
139
be functions of space. The resulting force balance equation, equation of state, and
equation of continuity are
an 4 85—1
P0095 - ‘VP + ‘3‘#(F)5t7_TV(V ' ‘1) (5373)
p = 60(r)p0(r)P (6271))
gt?- + p0(r)V - u = 0. (6.270)
Dividing Eq. (6.27a) and taking the divergence yields
6V - u _ Vp 4 aﬁ-l u(r) )
—at—— — V (p0(r)) + 3atﬁ-1V (p0(r)V(V u) . (6.28)
To eliminate particle velocity u, Eqns. (6.27b) and (6.270) are combined to yield
5’2
at (6.29)
V-u= —r.0(r)
Inserting Eq. (6.29) into Eq. (6.28) and moving the time derivative outside the spatial
derivatives yields
2
.El’. 54.93.1412 Hr _,,r§_P=
v (m(r,)+3at),v (p0(r)V(0()P)) 0()8t2 0. (6.30)
Finally, in analogy with Eq. (6.25), deﬁne a variable relaxation time via
75(r) = M2
3 , (6.31)
yielding
.& 23._T"(1_,,,. _,,,.§32_
V (p0(,,)+a,gv (po(r)nov( 0‘ )1”) 0‘ ’22 ‘0' (6'32)
Recognizing ﬂ = y—l, Eq. (6.32) matches Eq. (11) in [21] for constant compressibility
media. From Eq. (6.32), three sources of macro-scattering are recognized: compress—
ibility variations, density variations, and variations in the relaxation time T(r). The
power-law exponent y, on the other hand, is responsible for micro-scattering and the
associated attenuation and dispersion. For 6 = 1 (y = 2), the medium is homo-
geneous on the micro—scale. For 6 < 1, micro-heterogeneity manifests itself via a
temporal integration and the parameter deﬁned by Eq. (6.31).
140
6.3.3 Simpliﬁed Equations
Eq. (6.32) can be simpliﬁed considerably in several important special cases. In most
biological tissue, the density is nearly constant p0(r) z p0. Deﬁning the adiabatic
speed of sound via
600‘) = ——i— (6.33)
M(t)/20
and applying to Eq. (6.32) yields
2 9.6.. . 2 2 ,. _P_ _ __1_22_P _
v p+ atﬁv (600% ( )V(C(2)(r))) 02(1') 6:2 _ 0. (6.34)
If the local sound speed varies smoothly with position, then the second term in Eq.
(6.34) may be approximated:
1 (92p
ﬂ
6.35
Finally, if the relaxation time T(l‘) z r is nearly constant (that is, the only mechanism
for macro-scattering is sound speed inhomogeneity), then
1 0217
2 2.53 2 ___
Vp+r Vp c2(r)6t2
atﬁ = 0. (6.36)
6.4 Qualitative Properties of the Caputo-Wismer Equation
Although the Caputo-Wismer equation has been analyzed in Ref. [21], several facts
are brieﬂy reviewed in the following subsections. First, the dissipative and dispersive
properties of Eq. (6.26) and shown to satisfy both a power law attenuaton coefﬁ-
cient and a dispersive phase velocity. Secondly, causality is established by examining
the high-frequency limit. Finally, the Szabo wave equation, given by Eq. (4.4), is
interpreted as a plane-wave approximation to Eq. (6.26).
6.4.1 Low Frequency Limit
To derive a power law attenuation coefﬁcient for Eq. (6.26), the dispersion relation-
ship between angular frequency w and spatial wavenumber k is calculated. Applying
141
a space—time Fourier transform to Eq. (6.26) yields
2
—k2 + w? — k2Tﬁ(—iw)5 = 0. (6.37)
0
Solving for the wavenumber k(LU) yields
(6.38)
2(6) = w .
c0(/ 1 + (—z'm)ﬁ .
In the low frequency limit, the binomial approximation is applied, yielding
~ w TB ‘ ﬂrr ﬂ 3rﬂ , ,61r ﬂ
k(w) ~ g (1— —2— co (7) w + 2—2— sm (7) w . (6.39)
The attenuation coefﬁcient 0(a)) and phase velocity C(w) are calculated by inserting
ﬂ = y — 1 such that
y—l
r cos(7ry/2)wy.
2co
a(w) = Im h(w) = (6.40)
Deﬁning
20 = T’Hl C0~°»(7ry/2)|/(260) (6-41)
yields a power-law coefﬁcient given by Eq. (1.1). Likewise, the phase velocity C(w)
is computed by taking the real part of Eq. (6.39) yielding Eq. (4.15), where co is
taken as low frequency limit of phase velocity. Thus, the Caputo-Wismer equation
possesses a phase velocity predicted by the Kramers-Kronig relations [5—7]. In the
low frequency limit, the Caputo-Wismer equation shares the same dissipative and
dispersive properties of the both the power law model (Chapter 4) and the spatially
dispersive model (Chapter 5)
6.4.2 High Frequency Limit and Causality
For high-frequencies, the attenuation coefﬁcient assoCiated with Eq. (6.26) decays
slower than the low-frequency approximation given by Eq. (6.40). Letting w —-> 00,
Eq. (6.38) is approximated via
(6.42)
Thus, the attenuation coefﬁcient approaches
__) Sin("ﬂ/4) l—ﬁ/2
0(6)) ———COTﬁ/2 a) . (6.43)
To establish causality, an impulse point source 6(t)6(r) is applied to the Caputo-
Wismer equation, yielding an FPDE in terms of the Green’s function g(R, t)
1 62 _, 62-1
v29 — @5522 + Ty 16t2-1V29 = —6(t)6(r). (6.44)
Applying a temporal Fourier transform to Eq. (6.44) with the aid of Eq. (A.4) yields
the Helmholtz equation
v29 + k(w)§1 = — (1+ ((S—(ZW‘U' (6.45)
where 16(6)) is given by Eq. (6.38). Solving Eq. (6.45) subject to zero initial conditions
yields
eik(w)R
‘ R = . 6.46
9( ’w) 41rR(l + (—z‘w)y-1) ( )
Taking the magnitude of Eq. (6.46) in the limit as 0) —> 00 yields
, ex -a a) R
19024) —» pl ( ’ 1 (6.47)
47rRTy‘1w3/‘1 ’
where (1(0)) is speciﬁed by Eq. (6.43). The Paley-Wiener theorem, which relates
causality to the decay of the Fourier transform, states that g(m,t) is causal if and
only if the Fourier transform given by Eq. (6.47) decays slower than exp(—|w|) as
w -—> 00. Since 0 < B S 1, the high-frequency bound is deﬁned as
a(w) < C(r)wr, (6.48)
where C is independent of 0) (but possibly dependent on r) and r < 1. Application
of the Paley-Wiener theorem [128] thereby guarantees causality for solutions of Eq.
(6.44) for 1 < y S 2.
143
6.4.3 Relationship to the Szabo Wave Equation
In this subsection, the Szabo wave equation is shown to be a low-frequency approxi-
mation to the Caputo-Wismer Equation, just as the Blackstock equation in Eq. (4.10)
is the low-frequency approximation to the Stokes wave equation given by Eq. (2.1).
For observation points acoustically far from the radiator, the wavefront is an approx-
imately planar wave that propagates and propagating in the direction k = R/ R.
Under this approximation, which is valid in the limit of low frequency,
3%
CO at'
Taking the divergence of Eq. (6.49) and inserting Eq. (6.50) into the third term in
Eq. (6.26) yields
Vp z (6.49)
all-1‘72 ~ 163’“);
Bty‘l “23%:
where the composition rule for fractional derivatives has been employed [129]. Thus,
Ty—l
(6.50)
Eq. (6.26) is approximated as
_ 1% + 11222
Cgat2 6(2) aty+1
[which corresponds to the Szabo equation given in Eq. (4.4) by choosing r =
We = 0, (6.51)
(2a0c0/ cos(7ry/2))1/(y‘1). Due to the plane wave approximation, Eq. (6.51) is invalid
for 1) observation points acoustically close to radiating sources, 2) high frequencies,
and/ or 3) lossy media where 7' is large.
6.5 Discussion
6.5.1 Bio-Mechanical Interpretation
The physical signiﬁcance of the power law exponent is explored in this section. The
fractal ladder models covers four special cases, each of which is discussed below.
Case I: Micro-homogeneous Media (6 = 1): Let E = 0 and 17 > 0. The viscous
theory is recovered using a single dashpot, yielding a power law coefﬁcient of y = 2.
144
Due to the lack of springs, the media is homogeneous at all scales much smaller than
a wavelength, thereby indicating a lack of micro-heterogeneity.
Case 11: Simple Ladder Model (6 = 1 / 2): Let M = N = 0. In this case, there are
the same number of spring and dashpots at all levels of the mechanical arrangement,
yielding a power law exponent of y = 3/ 2.
Case 111: 1/2 < 6 < 1: Put M < N. In this case, there are more springs in the
network than dampers, indicating that the medium has a greater elastic component
than viscous component. The power law exponent y ranges form 1 to 1.5 in this case,
which is typical for most soft tissue [13].
Case IV: 0 < 6 < 1 / 2: Put M > N. In this case, there are more dampers than
springs, indicating a greater viscous component than elastic component. The power
law y ranges from 1.5 to 2 in this case, which is typical for complex ﬂuids like castor
oil and silicone ﬂuid [7].
Since springs correspond to elastic structures such as cellular and nuclear mem-
branes, while dampers correspond to inter-cellular and intra-cellular ﬂuids such as
cytoplasm, the exponent 6 measures the relative mechanical contributions of elastic
versus viscous structures. For instance, anatomical media such as liver have ,8 close
to zero due to the relatively complex tissue structure. Simpler media such as fat
have larger 6 near 0.5 due to a more ﬂuid-like structure. Due to the complexity of
biological tissue, most media have 6 between zero and 0.5, corresponding to power-
law exponents between 1 and 1.5. The compound ladder model also sheds light on
the dependence of the power law exponent on the pathological state of tissue. For
example, in [8], the power law exponent y in normal liver exhibited y z 1, whereas
y ranged from 1.25—1.4 in fatty liver. The increase in power law exponent was ex-
plained in terms of an increase in Rayleigh scattering in fatty liver relative to healthy
liver. Within the context of the present compound ladder model, the increase in y is
explained as an increase in viscous micro-structure of the tissue, relative to healthy
145
liver, which has a greater elastic component.
The recursive ladder model provides a simple explanation for the combined ef-
fects of absorption and incoherent scattering in biological media. First, note that the
(local) speed of sound is proportional to the square root of the spring coefﬁcient E.
Thus, the fractal arrangement of springs qualitatively accounts for sound speed inho-
mogeneity at multiple spatial scales, resulting in incoherent scattering of an incident
sound ﬁeld. However, sound speed inhomogeneity is not solely responsible for the
observed power-law dependence of the attenuation coefﬁcient. In addition, a viscous
mechanism is required to dissipate both the incident and incoherently scattered sound
ﬁelds. This viscous mechanism, like the sound-speed inhomogeneity, is represented at
multiple spatial scales by this model. The interactions between these two mechanisms
is mediated by the hierarchical arrangement of springs and dashpots.
6.5.2 Fractal Geometry Interpretation
Since the pioneering work of Mandelbrot, fractal geometry has been applied to ex-
plained the self-similar structure (e.g. alveolar surfaces, cell membranes, etc.), found
in biological systems [47]. More recently, fractal geometry has been applied to un-
derstanding the pathological architecture of tumors [130]. Since the vasculature of
tumors is more tortuous than healthy tissue, the measured fractal dimension of tu-
mor vessels is signiﬁcantly larger than normal veins and arteries [130]. To interpret
the power law exponent within the context of fractal geometry, which was discussed
qualitatively in Ref. [20], a quantitative relationship is provided by applying the
analysis presented in Ref. [118]. The exponent ﬂ appearing in the constitutive equa-
tion given by Eq. (6.15) is expressed in terms of the spectral dimension d3 of the
spring-dashpot network 6 = 1 — d3 / 2, where l S d3 3 2. The spectral dimension (or
fracton dimension) d, is related to the vibrational properties of the underlying fractal
network, such as the density of normal modes in the low-frequency limit [131, 132].
146
Intuitively, d3 measures the connectivity of a fractal network. Applying this relation
to the N -level recursive ladder-spring network yields d3 = 2 — 1/2N, which reveals
the fractal nature of the recursive ladder models presented in Sec. II.
The power law exponent y may also be related to the spectral dimension by noting
y = B +1, yielding
' y = 2 — d3/2. (6.52)
Furthermore, the dynamical pr0perties are related to the underlying fractal geometry
of the network, which is characterized by fractal dimension df. As a speciﬁc example,
consider a spring-dashpot network in the form of the Sierpinskigasket [46], which is
displayed in Figure 6.5, where each node corresponds to a dashpot, and each edge
corresponds to a spring. To construct such a fractal, the initial structure is an equilat-
eral triangle, which is then subdivided into four identical triangles, then the interior
triangle is removed, and this process is applied to the three remaining triangles. The
Sierpinski gasket displayed in Fig. 6.5 has a spectral dimension of d f z 1.37, which
yields a power law exponent of y a: 1.31 by Eq. (6.52). Higher dimensional analogs,
such as the 3D Sierpinski gasket may also be constructed [118], which has a fractal
dimension of df = 2 and a spectral dimension of d3 = 1.5474. The Sierpinski gasket
is an analogous model for both the self-similar and porous nature of tissue. Applying
Eq. (6.52) yields a power law exponent of y = 1.23, which is similar to the measured
power law exponents of soft tissue such as brain (y = 1.3) and liver (y = 1.14).
These calculations suggest a connection between the power law exponent and
the underlying fractal geometry of the tissue. As the spectral dimension increases
from 1 for the ladder network to 2 for more higher-connectivity network such as
the Sierpinski gasket or carpet, the power law exponents descends from 1.5 to 1.
Therefore, the power law exponent may describe the underlying connectivity of tissue
structures, with simple tissue such as breast fat has an exponent near 1.5, while more
complex tissue, such as liver or brain, having an exponent closer to 1. Moreover, this
147
Figure 6.5. Spring-dashpot network constructed in a Sierpinski gasket network. This
fractal exhibits self-similarity at ﬁve levels of magniﬁcation. Each node corresponds to a
dashpot, while each edge corresponds to a spring.
analysis suggests that Eq. (6.26) provides the governing equation for dissipative wave
propagation on fractal structures, much as the fractional diffusion equation describes
diffusion on fractals [97].
6.5.3 Stochastic Interpretation
Eq. (6.15) may also be viewed as a macroscopically averaged expression. Although
Eq. (6.3) assumes isotropy and spatial homogeneity on the macro-scale, spatial in-
homogeneity may be manifested by the temporal memory given by Eq. (6.3). At a
sufﬁciently ﬁne spatial scale, biological tissue may be regarded as a heterogeneous,
Newtonian ﬂuid with spatially varying material properties such as the coefﬁcients of
shear viscosity and adiabatic compressibility. That is, the viscous stress tensor given
by Eq. (6.1) with a variable ,u holds approximately on the micro—scale. Assuming
148
the medium to be random with zero mean, the stress tensor is regarded as a random
process Tij(r, t; f) unfolding on some microscopic time-scale 5. Assuming this random
process to be ergodic allows the statistical expectation to be interchanged with a tem-
poral ensemble average. Thus, Eq. (6.15) is regarded as a weighted, temporal average
of Eq. (6.1), where the temporal weight l/tf’ results from the statistical distribution
of a and (possibly) the compressibility K. and density p. Physically, the dissipative
stress at time t is proportional to the weighted averaged of the entire time history of
velocity gradients (or, equivalently, stresses). In short, the spatial inhomogeneity on
the micro-scale is mediated to the observable macro-scale by including this temporal
memory.
6.5.4 Nonlinear Media
The Caputo-Wismer equation, as derived in Section 3, assumes small amplitude os-
cillations and negligible heat conduction by utilizing a linear, adiabatic equation of
state. Although the adiabatic hypothesis is justiﬁed in most biological media due
to the negligible thermal conductivity, the linear assumption is not justiﬁed in many
biomedical applications where large amplitude effects must be modeled [133]. How-
ever, most nonlinear models, such as Burger’s equation and Westervelt’s equation,
assume a thermoviscous dissipation mechanism, resulting in an attenuation coeffi-
cient with frequency-squared dependence. In order to combine the effects of power
law attenuation with nonlinearity, several authors have formulated nonlinear FPDEs
[19, 20, 134]. Finite amplitude effects may be incorporated into the Caputo—Wismer
model by augmenting the equation of state with a quadratic term. Utilizing the
stress tensor given by Eq. (6.15), a nonlinear generalization of the Caputo-Wismer
equation may be rigorously derived for both homogeneous and inhomogeneous media.
The competing effects of nonlinearity and power law dissipation may then be studied
within the presented framework.
149
CHAPTER 7
Lossy Impulse Response in Power Law
Dispersive Media
7 .1 Introduction
This chapter returns to the problem of ﬁnite aperture radiation in power law media
discussed in Section 3.6.3. We are now equipped to solve this problem using the
Green’s function theory developed in the previous Chapters. Although many papers
have been published dealing with transient propagation in power law media [33, 34,
40], these studies considered plane-wave propagation for simplicity. Hence, the study
of time-domain radiation by ﬁnite apertures in power law media is limited. Time-
domain modeling in attenuating media was considered in [102] using the dispersive
tissue model developed in [30], which only considers linear attenuation media (y = 1)
within the Field II program. However, the impulse response was not presented in
closed form, and this expression required 1) numerical inverse Fourier transforms
and 2) approximated frequency-dependent attenuation as a constant for all points on
the surface of the radiator. A general 3D numerical model incorporating power-law
loss and ﬁnite apertures was developed in Ref. [35], while the combined effects of
diffraction, dissipation, and nonlinearity were studied numerically in Ref. [55] under
the parabolic approximation.
150
In this Chapter, the analytical tools developed in Chapters 3, 4, and 5 are syn-
thesized to provide an analytical description of pulsed ﬁnite apertures in power law
media. In particular, the lossy impulse response developed in Chapter 3 is general-
ized to power law media by deriving the loss functions mentioned in Section 3.6.3. To
derive these loss functions, the Green’s function for the spatially dispersive wave equa-
tion discussed in Chapter 5 are utilized. This Green’s function satisﬁes the following
desirable criteria:
1. In the low-frequency limit, the attenuation coefﬁcient satisﬁes a power law with
respect to frequency.
2. The frequency—dependent phase speed is commensurate with the Kramers—
Kronig relations.
3. The Green’s function is causal for all power law exponents between 0 and 2.
4. The scale paramter (ﬂocot)1/y does not depend on distance, allowing a closed-
form evaluation of the power law impulse response.
Following the analysis applied to the Stokes wave equation in Chapter 2, the causal
Green’s function given by Eq. (5.47) is decomposed into loss and diffraction compo-
nents. As in the special case of viscous media (y = 2), a mapping between the lossless
impulse response and the power law impulse response is developed. This mapping is
applied to our four canonical geometries: 1) the uniform circular piston the nearﬁeld,
2) the uniform circular piston the farﬁeld, 3) the rectangular piston in the nearﬁeld,
and 4) the focused, spherical shell. Since all of the essential features of the disper-
sive impulse response are exhibited by the circular piston, emphasis is placed on case
1. For biomedical applications, however, cases 3 and 4 are of greater interest. The
qualitative properties of ﬁelds produced by each geometry are discussed, followed by
numerical evaluation of each expression over both space and time.
151
7.2 Green’s Function Decomposition: Power Law Dispersive
Media
Similar to the viscous case discussed in Chapter 2, the asymptotic 3D Green’s function
to the spatially dispersive wave equation given by Eq. (5.47) admits the following
decomposition for y 74 1:
g(R, t) = L: [W — :ﬂ'éR/Col] [(60:31/96 (—(ﬁ0c:t)1/y)] dt'. (7.1)
Thus, the spatially dispersive Green’s function is given by the non-stationary convo—
lution
g(R, t) = 9P6, 13’) ® 90(R, t') (7-2)
where the power-law loss function is given by
u(t) t’
gp 2, the
diffusion equation is recovered, while as y —> 1, Eq. (7.4) approaches the one-way
wave equation (or the advection equation).
152
For y = 1, an alternate parameterization of the stable distribution f~1 (t) is utilized:
gpat'): “’l f.( " ) (7.5)
aoclt aotcl
7 .2.1 Power Law Impulse Response
As in the viscous case, the decomposed Green’s function given by Eq. (7.2) is inte-
grated over the radiating aperture S and multiplied by a factor of 2 to account for
the inﬁnite bafﬂe in the z = 0 plane, yielding
hp(r, t) = 2/Sg(r — r',t) d2r’. (7.6)
Inserting the decomposed Green’s function given by Eq. (7.2) into Eq. (7.6) and rec-
ognizing the deﬁnition of the lossless impulse response (Eq. (3.3) with the apodization
function set to unity) yields the mapping
hp(r. t) = gp(t. t’) e h(r. t’) (7.7)
where the convolution is taken over t’. Instead of a Gaussian, the loss function gp(t, t’)
for dispersive media is given in terms of the maximally skewed stable distribution
fy(t). Eq. (7.7), which maps boundary-value solutions of the lossless wave equation
to corresponding boundary value solutions of the spatially dispersive wave equation,
is a generalization of the earlier viscous result given by Eq. (3.6).
7.3 Circular Piston: Nearﬁeld
The baffled circular piston shown in Fig. 3.1 is considered in a homogeneous, power
law media with attenuation slope 0:0 and power law exponent y. Inserting Eq. (3.7)
into Eq. (7.7) and applying a substitution yields
czau 7' ’10“) I I A; (11)) I I
hp(a:,z,t)= (”/0 Ma(a:,¢) [[1 fy(t)dt —/ 2 fy(t)dt] d.) (7.8)
71' —oo —oo
153
where fm(7/)) = (t — rm)/([30tc0)1/y for m = 1 or 2. Noting that fy(t) is the PDF of a
maximally skewed y-stable distribution, the cumulative distribution function (CDF)
Fy(t) is identiﬁed as a deﬁnite integral of the Wright function
1 t 1 1 1 1
F2“) 712471-17) ’2 "13"("'1}’1”)’ (’9)
which follows from Eq. (C5) in the Appendix, yielding
hp(:r,z,t)= COW“) 7r1r/0 Mama) [17,, (W)— Fy (WM dz).
(7.10)
In the special case y = 2, Eq. (C.6b) is utilized, yielding
hp(x,z,t)=&:u?(QAﬂMa(x,w)[erf(2m)— er f(2—t—t\/E_)] 611/) (7.11)
Using the relation '7 = 200C; yields Eq. (3.9). Likewise, letting do —> 0 yields the
lossless case given by Eq. (3.7). Eq. (7.10) therefore generalizes the viscous case to
more general power law media. Finally, letting a0 —1 0 causes Fy(t) —+ u(t), thus
reducing Eq. (7.10) to the lossless expression given by Eq. (3.7). On axis (:1: = 0),
the integration in Eq. (7.10) is trivial, yielding
hp(0,z,t) = u(t)c0 [17,, (I’LL/Cl) — Fy (’ ‘ V z + “ /C°)] (7.12)
ﬂow-AW (ﬁotcowy
For y = 1.5, which corresponds to breast fat, the Wright function may be expressed
in terms of the generalized hypergeometric function in Eq. (C.6c), allowing the eval-
uation of Eqns. (7.10) and (7.12). In the special case of y = 1, the loss function given
by Eq. (7.5) is utilized, yielding
hp(a:,z,t) = 9%“2 [on Ma(a:,1/)) [61662) — 660161)] at (7.13)
where E1 (t) is the CDF for the maximally skewed stable distribution of index 1 using
the parameterization described in Ref. [87]. In Eq. (7.13) corresponds to the phase
velocity at a) = 1. On-axis, Eq. (7.13) reduces to
F1(t—z/cl)—F1(t_m/Cl)] (7.14)
h 0,z,t = to
P( ) U()1 aotq aOtCl
Both Eq. (7.13) and Eq. (7.14) are evaluated using the STABLE toolbox [3].
154
7 .3.1 Numerical Results
To verify Eqns. (7.13) and (7.14), the dispersive impulse response is compared
to numerical results generated by Field II [66, 135]. The Field II program’, which
is widely used within the medical ultrasonics community, generates apertures using
combinations of rectangular, triangular, and/ or bounding-line sub-elements. Disper-
sion is incorporated using the Hilbert dispersive model outlined in Refs. [30] and [102]
for linear attenuation media (y = 1). In addition, Field II decomposes the attenu-
ation coefﬁcient into a frequency independent component, evaluated at a reference
frequency f0, and a frequency dependent component. General power law attenuation
media is not implemented within Field II.
A circular piston of radius 10 mm was generated using 100 triangular sub-elements,
and ﬁrst tested in lossless media using an artiﬁcially large sampling rate of 800 MHz
to avoid aliasing. The computed impulse response was in good agreement with known
analytical results [48]. Linear with frequency attenuation was then activated using 0.5
dB / cm / MHz, and a sound speed of c1 = 1540 m/s. Impulse responses were calculated
on-axis at (0, O, 50) mm and (0, 0 100) mm and off-axis at (10, O, 50) mm and (10,
0,100)rnnr
The Field II dispersion model requires the following input parameters: a frequency
independent attenuation, a reference frequency, a frequency dependent attenuation
slope, and a minimum phase delay. In the following simulation, attenuation is initial-
ized with the following commands
set_fie1d (’att’, 50); % frequency independent attenuation
set_field (’att_f0’, 166); 2 reference frequency
set_fie1d (’freq_att’, 5e-5); % frequency dependent attenuation slope
'set_fie1d(’tau_m’,20.0); % minimum phase delay
lField II version 3.00 for MATLAB, http:Hm.es.oersted.dtu.dk/staff/jaj/field/
index.html
155
set_field (’use_att’, 1); % activate attenuation
The dispersive impulse response was then computed using Eqns. (7.13) and (7.14)
and compared with the Field 11 results. To evaluate the integral in Eq. (7.13), Gauss
quadrature was utilized. Since Field 11 incorporates additional delays in the dispersion
computation, the sound speed was adjusted to c1 = 1536 m/s in Eqns. (7.13) and
(7.14). Figure 7.1 compares these two methods at four observation points. In all
cases, the two models are in good agreement.
Several features of the impulse responses shown in Fig. 7.1 are notable. In all
cases, the impulse response is smooth and does not possess the breaks and corners
seen in the lossless impulse response shown in Fig. 3.3. Thus, like the viscous impulse
response discussed in Chapter 3, the dispersive impulse response is bandlimited and
not subject to the severe aliasing problems encountered in the lossless case. Unlike
the viscous impulse response, the dispersive response has a much longer temporal
duration due to the slowly decaying tail, which extends far beyond the duration of
the viscous impulse response. Hence, the dispersive impulse response is characterized
by a wake which trails the initial wavefront.
To compare the previously derived lossy impulse response developed in Chapter
3 with the power law impulse response derived in the preceding section, Eq. (7.10) is
evaluated in Figure 7.2 for a piston with radius a = 10 mm in media with attenuation
slopes a0 = 0.086 Np /mm/Msz, zero-frequency sound speed c0 = 1.5 mm/us.
Three power law media are considered with exponents y = 1.139, 1.5, and 2 (viscous).
Figure 7 .2(a) shows the on-axis impulse response for z = 50 mm, while Figure 7.2(b)
displays the impulse response for z = 100 mm. While all three media display spreading
of the impulse response due to dissipation, there are signiﬁcant differences between
the three media. In viscous media, the impulse response decays exponentially as
t —-> 00, while for y = 1.5 and y = 1.139, the decay is (an inverse power of t. Hence,
for y < 2, the impulse response decays much slower than the impulse response in
156
-—-Field ll ' ' 7 ' —Field ll w"
- - -power law _ - - -power Ia
d
v
.0
a:
.0
co
11L
.0
#1
.°
.5
impulse response (mm/us)
o
9
impulse response (mm/us)
.0 .0
'9 9
p
N
9
_a
O
C)
65 65.5 66 66.5 67 67.5 68 33
34 35 36 37 38
time (us) time 018)
(a) (:1), y, z) = (0, 0, 50) mm. (b) (3:,y, z) = (10, 0,50) m.
—Fleld ll ' ' ' ' ' —Field n w]
- - -power law ' ' 'POWOV '3
d
I
P
a:
477
impulse response (mm/us)
p o
a a:
impulse response (mm/us)
O
to
.0
N
L
65 65.5 66 66.5 67 67.5 68 65 66 67 66 69 70 71 72
time (us) time (us)
O
(c) (2:, y, z) = (0,0, 100) mm. (d) (a), y, z) = (10,0, 100) mm.
Figure 7.1. Comparison of the power law impulse response with the Field II model.
A baffled circular piston with radius a = 10 mm radiates in linear attenuation media
(y = 1) with an attenuation slope 0.5 dB/MHz/cm. Impulse responses are shown at a)
(a2,y,z) = (0, 0, 50) mm, b) (m,y,z) = (10, 0, 50) mm, c) (:r,y,z) = (0,0, 100) mm, and
(:r,y,z) = (10, O, 100) mm.
157
viscous media (y = 2). The slow decay, or long-term memory, is characteristic of
dispersive media. As y decreases from 1.5 to 1.139 with (10 ﬁxed, the decay of the
impulse response slows, and the duration of the wake increases. In Fig. 7.2(a), which
lies within the nearﬁeld of the radiator (z = 50 mm), the impulse responses in all
three media display both the “stretching” due to the ﬁnite extent of the aperture,
and the “spreading” due to dissipation. However, in Fig. 7.2(b), which is in the
farﬁeld of the radiator (24 = 100 mm), the impulse response is characterized mainly
by the “spreading” due to dissipation and dispersion, while the “stretching” due to
diffraction in Fig. 7.2(a) is absent. The amplitude decrease in all three media is
similar.
The on-axis velocity potential generated by a circular piston of radius a = 10
mm is evaluated by convolving Eq. (7.12) with the broadband pulse deﬁned by Eq.
(4.68) using the same parameters utilized in Fig. 4.5. As with the power law model
analyzed in Chapter 4, two power law media are evaluated: 1) a fat-like medium
with y = 1.5, CO = 1.432 mm/us, and a0 = 0.086 Np/MHz1'5/cm and 2) a liver-like
medium with y = 1.139, c0 = 1.569 mm/ps, and do = 0.0459 Np/MHzl'l39/cm. The
velocity potential in Figure 7.3 is displayed on-axis at a) z = 50 mm and b) 2 = 100
mm. In panel a) of Fig. 7.3, which lies within the nearﬁeld of the piston, the velocity
potential has a signiﬁcantly longer duration than the potential generated by a point
source displayed in panel a) of Fig. 4.5. In the case of liver, the interference between
edge and direct waves dominates Fig. 7.3(a), while the dissipative and dispersive
effects of the medium are minimal. The fat-like medium, which has a much larger
attenuation slope than liver, exhibits the combined effects of diffraction and dispersion
at depth 2 = 50 mm. In panel b) of Fig. 7.3, the temporal stretching of the pulse
caused by the ﬁnite extent of the aperture, is not pronounced. However, the frequency
down-shifting and amplitude attenuation are evident. Hence, for observation points
in the farﬁeld of the radiator, the effects of diffraction are much smaller, while the
158
---y=1.5
3 ""1.139
E
E
Q)
m
C
3.
w
9
0)
g
3
0.
.§
35 35.5 36
1.5 T 1 2
_y=
---y=1.5
a ""1.139
E
e 1-
Q
01
C
3
m
9
g 0.5-
3
0.
.§
065 69 70
(b) 2 = 100 mm.
Figure 7.2. Evaluation of the power law impulse response generated by a circular piston
with radius a = 10 mm with attenuation slope a0 = 0.086 Np/mm/MHz”, zero-frequency
sound speed c0 = 1.5 mm/us, and power law exponents y = 1.139, 1.5, and 2. Panel a)
shows the on-axis impulse response for z = 50 mm, while panel b) displays the result for
z z 100 mm.
159
dispersive properties of both liver and fat are more pronounced. Comparing the
velocity potential generated by a point source in Fig. 4.5(d) with Fig. 7.3(b), the
shapes of the two potentials are very similar in both media, although the magnitude
of traces shown in Fig. 7.3(b) are larger by a factor of 1ra2, which accounts for the
size of the radiating aperture.
7 .4 Uniform Circular Piston: Farﬁeld
Utilizing the dispersive loss function given by Eq. (7.5) and following the same steps
that produced the expression for the viscous case in Eq. (7.7) yields
ﬂ - o
hp(r,6,t) = 600746)] Fy (t r/CO + asmdcos ¢/CO) cosrpdtp. (7.15)
771' 31nd 0 ([306001/3/
The on—axis case is computed by convolving Eq. (3.18) with Eq. (7.3), yielding
2
a u(t) ( t— z/co )
h ,0,t = —— . 7.16
P“ ’ 2T(ﬂ0¢0t)1/yf ’ (aerial/2 ( ’
In the special case of linear attenuation media (y = 1), the parameterization described
in Ref. [87] is utilized with ,80 replaced by C10.
7 .4. 1 Numerical Results
The nearﬁeld and farﬁeld impulse response for a circular piston are compared in
Figure 7.4. A circular piston with radius a = 10 mm in a power law medium with
exponent y = 1.5, attenuation slope 0.086 Np/cm/MHz1'5, and phase velocity c0 =
1.5 mm / as is considered. In panel a), the nearﬁeld impulse response and the farﬁeld
impulse response were evaluated at r = 50 mm on—axis (6 = 0) and off-axis (9 = 77/ 12
and 05 = 0). Panel b) shows the nearﬁeld and farﬁeld impulse responses evaluated at
7‘ = 200 mm on—axis (6 = 0) and off-axis (6 = 77/12 and 05 = 0). For 7' = 50 mm, there
is signiﬁcant disparity between the nearﬁeld and farﬁeld responses, especially on-axis,
thus indicating that the farﬁeld approximation in Eqns. (7 .15) and (7.16) is invalid
160
0.08~ j j ' ' ' -‘_.fat .
---ﬁver
’23: 0.06- 7
g 0.04- :E
E :. i:
g 0.025 :: :: q
E O :: 52"‘v' Alia-P
$4.02 E .1 5 a 1
8 —0.04- :3: W 1
9 t: l:
—0.06 H I: ‘
l: -
-0.08- "' ' . ‘ .
31 32 33 34 35 36 37 38
- time (us)
(a) z = 50 mm.
—fat
0.06~ ---liver
"5 0.04-
B
E 0.02 , 5
s -': -=
8 l :
'5 -0.04: i:
-0.06* :g
I.
I. . 1 .
64 66 68 7O 72
time (us)
(b) z = 100 mm.
Figure 7.3. Velocity potential generated by a circular piston of radius a = 10 mm in
power law media: 1) a fat-like medium with y = 1.5, co = 1.432 mm/us, and a0 = 0.086
Np/MHzl'S/cm and 2) a liver-like medium with y = 1.139, co = 1.569 mm/us, and a0 =
0.0459 Np/MHzl'139/cm. The input pulse is deﬁned by Eq. (4.68). The velocity potential
is displayed on—axis at a) z = 50 mm and b) 2 = 100 mm.
161
this close to the aperture. For 7' = 200 mm, however, the nearﬁeld and farﬁeld results
are in closer agreement. Unlike the nearﬁeld and farﬁeld viscous impulse responses
depicted in Fig. 3.6, the responses in power law media are skewed toward later arrival
times and have much longer temporal durations.
Although the farﬁeld impulse response given by Eq. (7 .15) is limited to observation
points far from the aperture, Eq. (7 .15) requires only half as many CDF evaluations of
the nearﬁeld formula given by Eq. (7.11). Since the evaluation of stable CDFs requires
the evaluation of numerical integrals and series expansions, this reduction in CDF
evaluations will provide a signiﬁcant computational speed-up for large simulations.
7 .5 Rectangular Aperture
The baffled rectangular piston displayed in Fig. 3.7 radiating in a power law medium
is considered in this section. As in the lossless and viscous cases, the aperture is
sub—divided into four sub-apertures. To compute the dispersive impulse response,
Eq. (7.3) is convolved with Eq. (3.28), which yields
_7- _Z/
111—31 118—3 .)
P...1(z,t) =u(t)s/ (30600 " (Bocot) 2
0
c72+s2
d0 (7.17)
The contribution from each sub—rectangle to the lossy impulse response is then spec-
iﬁed by
h (t=ﬂP tP t 718)
P,s,l z) ) 27r( 8,1(z1 )+ l,s(za )) ° ( °
The dispersive impulse response is formed by Eqns. (7.18) and (3.27 ) Letting y = 2
in Eq. (7.17) yields the viscous impulse response given by Eq. (3.29). As with the cir-
cular piston, linear attenuation media (y = 1) in handled using the parameterization
described in Ref. [87] and replacing ,60 with 070.
162
3P ; I —On-a1XiS (near) ﬂ
:.'. ...on—axis (far)
A : : a... off—axis (near) -
€2.57 : : off-axis (far)
E l l
E, 2. : :
m I I
c I I
I I
§ 1.5- : l
9 I
E 1 - E
3 I
Q I
E I
-- 0.5. '
941 32 33 34 35 36
time (115)
(a) r = 50 mm.
0.35 p I I I on—aXiS (far)
i "'-'oﬂ-axis (near)
E 0.3 . .mno Oﬂ—aXis (far)
E
3 0.25~ 4
c
g 0.2-
E
g 0.15“
3
o. .1 -
.§ 0
0.05 '
930 . ' 138
time (us)
(b) r = 200 mm-
Figure 7.4. Lossy impulse response for a circular piston (a = 10 mm) in a power law
medium with exponent y = 1.5, attenuation slope 0.086 Np/cm/MHzl'5, and phase velocity
co = 1.5 mm/ps. In panel a), the nearﬁeld impulse response and the farﬁeld impulse
response are evaluated at r = 50 mm both on-axis (0 = O) and off-axis (0 = 1r/12 and
43 = 0). Panel b) shows the nearﬁeld and farﬁeld impulse responses evaluated at r = 200
mm on-axis (0 = 0) and off-axis (0 = 1r/12 and 45 = 0).
163
7.5.1 Numerical Results
Eq. (7.17) is evaluated for a rectangular piston with half-width a = 2.5 mm and
half-height b = 10 mm with attenuation slope do = 0.086 Np/mm/Msz and zero-
frequency sound speed c0 = 1.5 mm/ps. Three power law media are considered with
exponents y = 1.139, 1.5, and 2 (viscous). Figure 7.5(a) shows the on-axis impulse
response for z = 50 mm, while Figure 7.5(b) displays the impulse response for z =
100 mm. In panel a), which lies within the nearfield of the radiator, characteristics
of both the lossless impulse response and the loss function gp(t, t’) are evident. For
example, the corners in the lossless impulse response, which correspond to the corners
of the rectangular aperture, are smoothed by the dissipation in the media. In panel
b), which is in the farﬁeld of the radiator, the impulse response is dominated by the
dissipative properties of the loss function, such as the slowly decaying tail for media
with y = 1.5 and 1.139.
7 .6 Spherical Shell
The lossy impulse response for a spherical shell depicted in Fig. 3.11 evaluated in a
power-law media is considered in this section. The lossless impulse response, given
by Eq. (3.31), is convolved via Eq. (7.7), yielding a lossy impulse response
a1:Co f0“ Na,R(’3z’¢) [Fy (W) _ Fy (ﬁgﬂadE)
hp(r, z, t) = u(t)
where Fy(t) is the stable CDF given by Eq. (7.9). Eq. (7.19) is valid for all points
excluding the geometric focus, where
(7.20)
_ 2u(t)(R — VR2 — a?) t - R/co
121903 2, t) — (ﬁocoﬂVy y (_(ﬁocot)1/y) .
As in the viscous case, Eq. (7.20) is non—singular, thus indicating that the high fre-
quency components of the impulsive excitation have been low-pass ﬁltered by the
164
A 0.5 = 2
% IIIy =1.5
E 04 ""'1.139
7;
(I)
C
3 o 3
(I)
‘2
g 0.2-
3
‘s‘
-- 0.1
(92.5 L. . 35 35.5 35
—y = 2
0.25’ ...y=1.5
'-'-'1.139
.0
N
impulse response (mm/us)
.0 ‘3.
d 01
.0
o
01
(95.5 66 55.5 57 67.5 68 68.5
time (us)
(b) 2 = 100 mm.
Figure 7 .5. Lossy nearﬁeld impulse response for a rectangular piston with half-width a =
2.5 mm and half-height b = 10 mm evaluated in power law media with attenuation slope
a0 = 0.086 Np/mm/MHz”, zero-frequency sound speed co = 1.5 mm/ps, and power law
exponents y = 1.139, 1.5, and 2. Panel a) shows the on-axis impulse response for z = 50
mm, while panel b) displays the results for z = 100 mm.
165
medium. On-axis, the integration in Eq. (7.19) simpliﬁes to the closed—form expres-
sion
we. at) = u(t)—Rz—C°
F, (t — \/R2 + z2 +2zx/R2 _ (Ir/CO) _ F, (t — (R+z)/c0)]
(ﬂocot)1/y (ﬂocot)1/y
(7.21)
which reduces to Eq. (3.37) for y = 2. Linear attenuation media (y = 1) is handled
by using the parameterization described in Ref. [87] and replacing ,60 with 070.
7.6.1 Numerical Results
In the following computation, a Spherical shell with radius a = 10 mm and radius
of curvature R = 70 mm centered at (0,0, —70) mm radiates in power law media
with attenuation slope a0 = 0.086 Np/mm/Msz, zero-frequency sound speed c0 =
1.5 mm/ps, and power law exponents y = 1.139, 1.5, and 2. Figure 7.6 displays the
on-axis impulse response in the a) pre-focal region 2 = —20 mm and the focal region
2 = —1 mm. In a lossless medium, the impulse response increases without bound
to a scaled and delayed Dirac delta function as the observation point approaches the
focus, and then decreases in magnitude as the observation point moves past the focus.
In contrast, the magnitude of the impulse response is greater in the pre—focal region
relative to the focal region in all three lossy media considered, indicating a severe
degradation in the focusing ability of the radiator due to power law dissipation. In
addition, in media with exponents y < 2, the impulse response has an extended
duration due to the dispersion in the medium, which further degrades the focus.
7 .7 Discussion
7 .7 . 1 Stationary Approximation
For each point on the aperture S, a time-domain spherical wave given by Eq. (5.47)
is generated. At a given observation point, these spherical waves arrive with varying
166
_ 3. —y=2
6 :‘i . I-Iy=1.5
A 5'! :' “1.139
g~5‘ g E l
E s '5'
E s 2:
m4” _5 "
m I
c i
§3~ :
9 i
g .5
'5 2’ 5
O- i
.E i
1' E
. .5
0 32.5
6" —y=2 J
---y=1.5
35 ----1.139
E
E.
m 4'
(D
C
§3-
93
33
5 2'
Q.
.E
1.
0 45 45.5 46 46.5 47 47.5
time (us)
(b)z=-1 mm.
Figure 7.6. Lossy impulse response generated by a baﬂ‘led spherical shell with radius a
= 10 mm and radius of curvature R = 70 mm in power law media with attenuation slope
a0 = 0-086 Np/mm/Msz, zero-frequency sound speed c0 = 1.5 mm/ps, and power law
exponents y = 1.139, 1.5, and 2. Panel a) shows the on-axis impulse response for z = -20
mm, while panel b) displays the result for z = -1 mm.
167
phases and amplitudes. The amplitudes of the transient spherical waves differ on
account of both the diffraction of the wave and the attenuation of the wave. Note
that Eq. (7.7) allows for varying attenuated path lengths due to the ﬁnite extent of
the aperture. Hence, the power law loss and diffraction mechanisms are coupled in
Eq. (7.7).
Loss and diffraction may be decoupled by ﬁxing a particular value of time
t = R*/c0 in Eq. (3.6), where R* denotes a representative distance between source
and observer. For instance, R* is chosen to be the mean distance between the aper-
ture and the observer. Physically, this approximation assumes that the attenuated
path between the observer and all points on the source is approximately identical.
Approximating Eq. (7.7) yields a stationary convolution over t:
hL(I‘,t) % 9P(R*/Co,t')®h(l‘,t')
z gp(R*/co, t) * h(r, t) (7.22)
where the convolution * is now performed over the t variable. For y 2 l, gp(R*/c(), t)
is noncausal since u(R*/c0) = 1 and fy(t) > 0 for all t. However, as discussed in
Section 4.7.1, the loss of causality has little practical impact for observation points
more than one wavelength from the aperture. As shown below, Eq. (7.22) is an
excellent approximation for many problems. Also, since standard convolutions are
involved, Eq. (7.22) is easy to implement using FFTs.
Eq. (7.22) decouples the attenuation and diffraction problems, approximating
the lossy impulse response by cascading the lossless impulse response and a 1D loss
function. As discussed in Ref. [31] and Ref. [102], this approximation captures the
essential ﬁeld properties “far” from the radiating aperture. Then, using Eq. (7.22),
the velocity potential (I) is expressed via
@(r, t) z u(t) * gp(R*/c0, t) * h(r, t). (7.23)
168
7 .7.2 Physical Interpretation
The physical interpretation of the impulse response expression for the uniformly ex-
cited circular piston discussed in Ref. [74], which decomposes the impulse response
into direct and edge waves, is also applicable to transient propagation in power law
media. To explain the source of the maximally skewed CDF within the integrand of
Eq. (7.10), a random process T = T (t) is deﬁned, corresponding to randomized delay
times with a PDF given by Eq. (5.66). The nearﬁeld impulse response given by Eq.
(7.10) is then expressed as
hpcr. 2. t) = @ﬁf—(t) f0” Mae. «7) [P(T s n) -. P(T s 72)] cm. (7.24)
Physically, the ﬁrst term in Eq. (3.9) corresponds to the edge wave generated by
the discontinuity at the of the piston, whereas the second term corresponds to the
direct wave due to the bulk motion of the piston. The delay 71 represents the time
for sound on the rim of the piston at angle 2;» to travel to the observation point
(r, 0, z) in a homogeneous medium with sound speed c0, while the delay 72 represents
the shortest travel time from the piston to an observation point within the paraxial
region of the piston. Since time has been randomized in Eq. (7.24), the edge wave
in the lossless impulse response given by Eq. (3.7) is replaced with the probability
that the delay T is less than the ﬁxed delay 71. Likewise, the direct wave is replaced
with the probability that the delay T is less than 72. In the lossless case, where time
ﬂows deterministically, P(T S 73-) = u(t — Ti) for the edge ( = 1) wave or the direct
(2' = 2) wave. For power law media, however T may assume a continuum of values
according to the distribution described by Eq. (5.66). In the special case of viscous
(y = 2) media, T is symmetrically distributed about R/co, while for y < 2, T is
skewed toward larger delays due to the dispersion in the medium.
Numerical evaluations of the dispersive impulse response show distinct behavior
within the nearﬁeld, the transition region, and the farﬁeld. For the values of 010 and
169
y evaluated in Fig. 7.2, the diffraction component dominates in the nearﬁeld, with
only a small amount of smoothing and skewing in the vicinity of the head and tail
of the response. For observation points in the transition region between the nearﬁeld
and farﬁeld, the effects of diffraction and loss are both apparent. In this transition
region, the impulse response is both smooth and asymmetric. Finally, in the farﬁeld,
the effects of viscous loss predominate, yielding increasingly broad responses with a
reduced directivity.
7.7 .3 Complex Apertures and B-Mode Imaging
Ultrasound transducers often employ more complicated, non-canonical geometries in
dispersive media, such as arrays of curved strips [136], apodized circular pistons [69],
or concave annular arrays [137]. These more realistic, albeit more complicated, radi-
ator geometries may be handled by either 1) applying Eq. (7 .2) to a single-integral
representation of the lossless impulse response, or 2) constructing the desired geom-
etry using circular or rectangular sub-elements and applying the principle of super-
position, similar to the approach utilized by Field 11 [66]. Since complex geometric
shapes are easily constructed from triangular sub-elements, the fast nearﬁeld expres—
sion derived for uniformly excited triangular aperture in [138] can be utilized within
the present dispersive impulsive response framework to analyze complicated aperture
geometries in lossy media. Ti'ansient ﬁelds due to electronically and / or geometrically
focused 2D and 3D phased arrays may then be synthesized by superposition, allow—
ing the point-spread functions generated by non-circular geometries in viscous media
to be computed. Finally, the analytical results in this chapter may be utilized for
the modeling of diffraction and scattering of phased array imaging systems. Thus,
the analytical tools developed below may be used to analyze the effect of power law
attenuation on B-mode imaging systems.
170
CHAPTER 8
Conclusion
8.1 Summary
Chapter 2 develops the transient Green’s function theory for the Stokes wave equation.
A previously derived approximate Green’s function, given by Eq. (2.20), is derived
as the leading-order term in an inﬁnite series solution. The error incurred via the
approximate Green’s function given by is analyzed as function of the relaxation time
7 and the distance from the source R, yielding the simple relation between 7 and
R for a 1% error threshold. The asymptotic Green’s function is decomposed into
diffraction and loss operators, where the diffraction component satisﬁes the lossless
wave equation and the loss component satisﬁes the diffusion equation.
Chapter 3 develops the lossy impulse response for circular, rectangular, and spher-
ical pistons using the Green’s function decomposition developed in Chapter 2. This
decomposed Green’s function relates the solutions to the Stokes wave solution and
the lossless wave equation via Eq. (3.6), thereby facilitating the derivation of the
nearﬁeld lossy impulse response for a circular piston in Eq. (3.9). The resulting ex-
pressions are strongly causal and straightforward to implement numerically. Eq. (3.9)
also eliminates the aliasing problems associated with the lossless impulse response.
Expressions are also developed for the farﬁeld of a circular piston in Eqns. (3.21) and
(3.22), the nearﬁeld of a rectangular piston in Eq. (3.30), and the bafﬂed spherical
171
shell in Eq. (3.35). Velocity potential computations display a signiﬁcant attenuation
of the direct wave within the paraxial region of the piston due to viscous diffusion.
Both the physical and numerical properties of the resulting lossy impulse responses
are analyzed, revealing the critical role of dissipation in the farﬁeld. Applications to
phased array ultrasonic imaging show a reduction in both axial and lateral resolution.
In Chapter 4, 3D Green’s functions in power law media are analytically computed
in the time-domain using the tools of fractional calculus. Chapter 4 shows that the
Green’s functions for a medium with an attenuation coefﬁcient given by Eq. (1.1) are
maximally-skewed stable distributions, shifted by a delay R/co, scaled by (aOR)1/y,
and multiplied by a spherical diffraction factor of 1/(477R). Mathematically, the power
law Green’s functions are represented by a transient spherical wave, which embodies
the diffractive process, convolved with a loss function gL(R, t), which embodies both
the attenuation and dispersion of the medium. Physically, therefore, the solutions to
the power law wave equation are represented as a linear coupling between diffraction
and dispersive loss. The time-domain power law Green’s functions presented here
are exact solutions to the power law wave equation in Eq. (4.19) and approximate
solutions to the Szabo wave equation given by Eq. (4.4). 2D Green’s functions in
power law media are also calculated in terms of stable distributions.
Efficient numerical evaluation of stable distributions is available via the STABLE
toolbox [3]. For all y < l, the Green’s functions are expressed in terms of Fox H -
functions, while for y > 1, the Green’s functions are expressed in terms of the H
function and Wright functions. For y = O, 1/3, 1/2, 2/3, 3/2, and 2, the Green’s
function solutions are expressed in terms of widely-known special functions. These
analytical Green’s functions were numerically veriﬁed against the MIRF result [31]
and Hilbert dispersive model [30], and example ﬁelds were computed. The results
show that the power law Green’s function is causal for y < 1 and noncausal for
'1 S y S 2. Despite being noncausal for 1 S y S 2, the power law Green’s function
172
provides an excellent causal approximation even for observation points close to the
radiating source. For 0 < y < 2, the Green’s function decays as t—y'l'l, leaving a
slowly decaying wake behind the primary wavefront.
The power law model discussed above suffers from one major drawback: solutions
are noncausal for power law exponents greater than or equal to one. To address this
deﬁciency, alternate power law FPDE based on fractional spatial Operators are con-
sidered in Chapter 5: the Chen-Holm wave equation and a spatially dispersive wave
equation. Green’s functions are derived for both equations, which yield causal solu-
tions for all applicable power law exponents. The Chen-Holm equation is shown to be
non-dispersive, while the spatially dispersive wave equation supports a phase velocity
predicted by the Kramers—Kronig relations. 3D Green’s functions are computed for
both the Chen-Holm and spatially dispersive wave equations. For the Chen-Holm
model, the leading order term is a shifted and scaled symmetric stable distribution
multiplied by a spherical spreading factor. For the spatially dispersive model, the
3D Green’s function is a shifted and scaled maximally skewed stable distribution. In
both cases, the Green’s function is causal, yet the Chen-Holm model does not cor-
rectly account for dispersion. 1D and 2D Green’s functions for both models are also
calculated, followed by a stochastic interpretation.
To model the relationship between viscoelastic parameters and dispersion, Chapter
6 proposes a fractal ladder network of springs and dashpots as a model for power
law attenuation media. In order to derive the constitutive equation, both a simple
and a recursive ladder model are considered in order to capture the viscoelastic,
self-similar, and hierarchical properties of biological tissue. These fractal ladders
capture the hierarchical arrangement of elastic and viscous components present in
biological media. The fractal ladder network produces a stress-strain relationship
with a fractional derivative of order 1 / 2, while the recursive ladder produces fractional
derivatives of all orders between zero and one. Hence, the resulting constitutive
173
equation interpolates between a Hookean solid and Newtonian ﬂuid via the non-
local fractional derivative operator. The attenuation coefficient computed from this
constitutive equation follows a power law in the low frequency limit, thus agreeing
with previous phenomenological theories and experimental data.
From the fractional constitutive equation in Eq. (6.15) and a linear equation of
state, the Caputo-Wismer equation, which models longitudinal wave propagation in
power law media via a time-fractional derivative, is derived. By applying a plane
wave approximation to the loss operator in the Caputo—Wismer equation, the Szabo
wave equation is recovered [19] as a noncausal approximation to the Caputo-Wismer
equation. Hence, the Szabo and Caputo—Wismer equations, which were originally
proposed as phenomenological models for power law attenuation and dispersion, are
derived from a fractal ladder model. The power law exponent y is related to the
spectral dimension d3 of the underlying fractal geometric structure, thereby supplying
the power law exponent with physical meaning.
Chapter 7 extends the lossy impulse response methodology developed in Chapter
3 to power law media. The Green’s function to the spatially dispersive wave equation,
which is obtained in Chapter 5, is utilized. The Green’s function is decomposed into
diffraction and loss components, and impulse response expressions for the nearﬁeld
circular piston (Eq. (7.10)), farﬁeld circular piston (Eqns. (7.15) and (7.16)), rectan-
gular piston (Eq. (7.17)), and spherical shell (Eq. (7.19)) are derived in power law
media. In the special case of linear attenuation media (y = 1), the nearﬁeld circular
piston is numerically veriﬁed against the Field 11 program with good agreement. The
coupled effects of power—law attenuation and diffraction are evaluated in the time-
domain for liver and fat, thus forming the basis for medical ultrasound simulations
in biological media.
174
Spatially Caputo-
Dispersive Wismer Chen-Holm
I y=2 I y=2 >
1 Stokes
y=2
Plane Plane
Wave Wave Ham Wm
Approx Approx Approx
V V
Szabo ] W > Blackstock
Figure 8.1. Diagram showing the relationships between the spatially dispersive, Caputo—
Wismer, Chen-Holm, Szabo, Stokes, and Blackstock wave equations studied in this thesis.
8.2 Relationships between FPDE Models
Figure 8.1 summarizes the relationship between the Stokes (Chapter 2), Szabo (Chap-
ter 4), Blackstock (Chapter 4), Chen-Holm (Chapter 5), spatially dispersive (Chapter
5), and Caputo—Wismer (Chapter 6) wave equations. The tOp row (spatially disper-
sive, Caputo—Wismer, Chen-Holm, and Stokes equations) are causal for all power law
exponents between 0 and 2, while the bottom row (Szabo and Blackstock) are non-
causal in general due to the plane wave approximation. Finally, in the special case
of viscous media (y = 2), the spatially dispersive, Caputo-Wismer, and Chen-Holm
wave equations all reduce to the Stokes wave equation, while the Szabo equation
reduces to the Blackstock equation.
175
APPENDIX A
Fractional Derivatives
The notion of fractional derivative was originally postulated by L’Hospital in 1695,
although a rigorous formulation was not published until Liouville and Reimann’s work
in the 1830’s and 1840’s. Until the 1970’s, however, fractional calculus was limited to
the arena of pure mathematics [129]. In the late 1970’s and early 1980’s, physicists
and engineers, spurred largely by the visionary work of B. Mandelbrot in the theory
of fractals [46], began to apply fractional calculus to problems in anomalous diffusion
and relaxation.
A.1 Riemann-Liouville Fractional Derivatives
The Riemann-Liouville fractional derivative is formally deﬁned via a hyper-singular
integral [129]
dyf__ 1 ‘ f(t’) ,
3f; — 1“(-y) ./—00 (t - t’)1+y dt' (Al)
The non-integrable singularity is removed from Eq. (A.1) by recognizing a differenti-
ation operator, yielding
1 d2 t 99') I -
1 d3 t g(t’ .
I‘Zz—ngtﬁf—oom‘ylrrdt' 1f1 < y < 2.
In the literature, the 0 < y < 1 case is known as the Reimann-Liouville form, whereas
(A.2)
the 1 < y < 2 case is the Liouville form. The following Laplace transform relationship
176
for fractional derivatives [129] is essential for the derivations presented herein:
dyg y
— = . A03
£(dty) sag) < )
Letting s = —iw yields the Fourier transform relationship
dyg _ _ y
.7: (rt?) — (-—zw) T(g). (A.4)
In the 0 < y < 1 case, Eq. (A.4) assumes g(t) vanishes for t < 0. However, for
l S y S 2, g(t) may be nonzero for t < 0.
A.2 Fractional Laplacian
The fractional Laplacian, or Riesz fractional derivative, is deﬁned via spatial Fourier
transforms [129]
(472)”2 we) a 5“ mama]. (45)
where \I'(k) is the spatial Fourier transform of 2,1)(r) deﬁned via Eq. (1.4b). Eq. (A.5)
may also be stated as a symmetric convolution over 723 (see Eq. (12) in [20] with d =
3):
(—V2)y/22,b(r) = Ay / ———V—29—(rI—)—d3r', (A.6)
R3 |r - r’|1+y
where
1" (1/2 + y/2)
(2x/7rlz“yP (1 - y/Z)’
4,, = — (A.7)
A.3 Skew Fractional Laplacian
More general fractional Laplacians may be constructed by introducing skewness. In
the analysis of the spatially dispersive wave equation, we deﬁne the skew fractional
Laplacian via
$46) = F1 [wax—7km] (A.8)
177
where k3 = |k|sgn (cos kg) is the signed wavenumber which ranges over R. That
is, k3 > 0, if k lies in upper half-space of k-space, while k3 is negative if k lies in
the lower half-space of k-space. In 1D, the skew Laplacian reduces to the Reimann-
Louiville fractional derivative by the relation Eq. (A.4). Thus, Eq. (A.8) generalizes
the Reimann-Louville derivative to higher-dimensions.
178
APPENDIX B
Stable Distributions
Two classes of stable distributions, namely the maximally skewed stable distribution
fy(t) and the symmetric stable distribution wy(t), are introduced in Chapters 4 and
5, respectively. As noted in these chapters, these stable distributions are actually
special cases of more general stable distributions, which are discussed in detail in Ref.
[49, 87, 89]. Although the mathematical theory of stable distributions was developed
independently of fractional calculus, the deep link between these two subjects has
recently been noted [45]. In brief, the solutions to fractional diffusion equations are
stable distributions, while stable distributions provide a stochastic description of the
underlying physics that give rise to fractional diffusion.
B.1 Deﬁnitions
Stable random variables are deﬁned via the characteristic function w(k) [87]
‘I’(k) = exp {—|k|aoa [1 —iﬁsgn(k)tan(1ra/2)] +ikp}, (B.l)
where a # 1 is the index, lﬁl S 1 is the skewness, o is the spread, and a is the location.
In the special case of a = 1 [87],
\i/(k) = exp (-—|k|o [1 + 2iﬁ/7rsgn(k) ln lkl] + wk). (13.2)
179
The density ib(:z:) is given by the 1D inverse Fourier transform (see Eq. (1.3b):
¢3(t)=f-l [xi/(4)] — 1 [00 \‘I’J(—k)e-iktdk (3.3)
_ 5; _00
Two special cases are utilized extensively: 1) the maximally skewed stable distribution
(6 = 1), denoted by fy(t), and the symmetric stable distribution (6 = 0), denoted
by iDy(t).
Like Fourier transforms, alternate conventions are utilized to characterize stable
distributions. The convention deﬁned by [89] is widely used in the literature; in this
alternate convention, the symmetric stable is denoted by wy(t) while the maximally
skewed stable is denoted by fy(t). Fortunately, the symmetric stable distributions
coincide for both conventions (iDy(t) = wy(t)). For y # 1,
as) = Isec(7ry/2)|’/yfy (I sec(7ry/2)|’/yt) . (8.4)
As noted in Chapters 4 and 5, the stable distributions fy(t) and fy(t) are represented
as Fox H -functions and Wright functions for y 75 1. These higher transcendental
functions are deﬁned and analyzed in Appendix C.
B.1.1 Maximally Skewed Stable Distributions
The maximally skewed stable distribution fy(t) utilized in Chapters 4 and 5 may be
represented via by the one—sided Fourier integral [89]
fy(t) = %Re‘/:0 exp ((—ik)y) exp(—ikt) dk. (B5)
For y = 0, f0(t) = 6(t — 1). For y =1,f1(t)= 6(t — 1). For all other 0 < y < land
1 < y S 2, fy(t) is a smooth (Coo) function of t. By application of Euler’s formula,
fy(t) = l AGO exp (cos (I?) Icy) cos (kt — sin (g) k”) dk. (B6)
77
180
B.1.2 Symmetric Stable Distributions
The PDF 103/(t) may be expressed as a one-sided inverse Fourier transform [89]
wy(t) = -71;/000 exp(—wy)cos(wt) dw. (B.7)
By exploiting the cosine addition formula and the scaling properties of Fourier trans-
forms, the following identity is obtained
;/0 exp (—aky) sm(Rk) Sln(bk) dk = —2a1/y [my ( al/y ) _ wy ( al/y )] '
B.2 Cumulative Distribution Functions (CDFs)
The computation of 1D Green’s functions and lossy impulse responses requires the
evaluation of cumulative distribution functions (CDFs). For an arbitrary probability
distribution function p(x), the CDF is deﬁned via
2
P(x) =/ p(rc’) dx'. (B9)
-00
The maximally skewed stable CDF is deﬁned as
t
Fy(t) =/ fy(t’)dt', (B.10)
—00
while the symmetric stable CDF is deﬁned via
t
mm = / wy(t’) dt’. (8.11)
—00
Special cases of Eqns. (8.10) and (BM) are given in the main text as needed.
181
APPENDIX C
Special Functions
The Fox H -function and Wright function, while commonly utilized in the fractional
calculus and stable random variable communities, are not included within the common
parlance of electrical engineers and applied physicists. Deﬁnitions and properties of
these transcendental functions are therefore tabulated below.
C.1 The Fox-H Function
The Fox H -function [81, 129] is deﬁned via a complex contour integral
mn )A 3
HM} (z 8:23:13) ) E 2—17r_i/LX(S)Z ds, (C.1)
where
(010: AP) = (a1) A1): (a21A2)9 ' ' ' i (apt/1P)
(bqqu) = (bltBl)!(b2iBz)i°”t(bQiBQ)i
and the integral density x(s) is deﬁned as
= H.911 P(bj - Bjs) 113-1:1 P(l -- aj + Ajs)
H§=m+1F(l — bi + Bjs) Hg=n+1 P(aj — Ajs)
Since the Gamma function is singular when the argument is a negative integer,
x(8)
(0.2)
Eq. (C.2) contains an inﬁnite number of poles in general. Therefore, the contour
L is chosen to separate the poles of Eq. (C.2). Note that Eq. (C.1) is the inverse
182
Mellin transform of Eq. (C.2). The Fox H -function, which generalizes most of the
special functions of mathematical physics, is discussed in great detail in Refs. [81]
and [83].
C.2 The Wright Function
The Wright function ¢(oz, 6,2) is deﬁned by an inﬁnite series (see Eq. (1.11.1) in
Ref. [129]):
Z]:
¢(aaIB;Z) E Igm-ﬁ. (0.3)
The Wright function and the more general H -function are related (see Eq. (6.2.66)
Hi’? (2] E303 ) = ¢(—a,b; —z). (C.4)
Eq. (0.4) is computed by evaluating the contour integral deﬁnition of the H -function,
in Ref. [129]) via
which possesses an inﬁnite number of poles, with Cauchy’s residue theorem. The
resulting inﬁnite series is identiﬁed as Eq. (C3).
The Wright function is differentiated and integrated via the following relation
[129]
d n
(a) ¢(a,ﬂ; z) = ¢(a, a + nﬁ; z), ((3.5)
where n is a positive integer. Three special cases, which were calculated using Math-
ematica, are
1 1 1 z2
4) —§,§,Z) —-\7_;exp (—-:1-) (0.68.)
1 z
(b (—§,1,z) =1+erf(§) (C.6b)
1 15244.22 22 274542
_—)1;2 =1+ F -)_;—:_;— —— _3—1—)_1
¢( 3 ) r(1/3)22(3 633 27)+I‘(—l/3)2F2(3633 272)
(C.6c)
In general, the Wright function may be reduced to sums of generalized hypergeometric
functions when y is rational.
183
BIBLIOGRAPHY
[1] S. A. Goss, L. A. Frizzell, and F. Dunn, “Ultrasonic absorption and attenuation
in mammalian tissues”, Ultrasound Med. Biol. 5, 181—186 (1979).
[2] R. J. McGough, “Rapid calculations of time-harmonic nearﬁeld pressures pro-
duced by rectangular pistons”, J. Acoust. Soc. Am. 115, 1934—1941 (2004).
[3] J. P. Nolan, “Numerical calculation of stable densities and distribution func-
tions”, Commun. Statist. Stochastic Models 13, 759—774 (1997).
[4] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, Molecular
Biology of the Cell ( Fourth Edition) (Garland Science, New York) (2002).
[5] M. O’Donnell, E. T. Jaynes, and J. G. Miller, “Kramers-Kronig relationship
between ultrasonic attenuation and phase velocity”, J. Acoust. Soc. Am. 69,
696—701 (1981).
[6] T. L. Szabo, “Causal theories and data for acoustic attenuation obeying a fre-
quency power-law”, J. Acoust. Soc. Am. 97, 14—24 (1995).
[7] K. R. Waters, M. S. Hughes, J. Mobley, G. H. Brandenburger, and J. G. Miller,
“On the applicability of Kramers—Kréinig relations for ultrasonic attenuation
obeying a frequency power law”, J. Acoust. Soc._Am. 108, 556—563 (2000).
[8] P. A. N arayana and J. Ophir, “On the frequency dependence of attenuation in
normal and fatty liver”, IEEE Trans. Sonics Ultrason. SU-306, 379—383 (1983).
[9] T. Lin, J. Ophir, and G. Potter, “Frequency-dependent ultrasonic differentiation
of normal and diffusely diseased liver”, J. Acoust. Soc. Am. 82, 1131—1138
(1987).
[10] F. T. D’Astous and F. S. Foster, “Frequency dependence of ultrasound atten-
uation and backscatter in breast tissue”, Ultrasound Med. Biol. 12, 795-808
(1986).
[11] K. J. Parker and W. C. Waag, “Measurement of ultrasonic attenuation within
regions selected from B-scan images”, IEEE Trans. Biomed. Eng. 30, 431—437
(1983).
184
[12] P. He, “Experimental veriﬁcation of models for determining dispersion from at-
tenuation”, IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 46, 706—714 (1999).
[13] F. A. Duck, Physical Properties of Tissue { First Edition), 99—124 (Academic
Press, London) (1990).
[14] P. M. Morse and K. U. Ingard, Theoretical Acoustics, 270—300 (Princeton Uni-
versity Press, Princeton, New Jersey) (1968).
[15] M. A. Biot, “Theory of propagation of elastic waves in ﬂuid-saturated porous
solid. 1. Low-frequency range”, J. Acoust. Soc. Am. 28, 168—178 (1956).
[16] A. I. Nachman, J. F. Smith, and R. C. Waag, “An equation for acoustic prop—
agation in inhomogeneous media with relaxation losses”, J. Acoust. Soc. Am.
88, 1584—1595 (1990).
[17] M. J. Buckingham, “Causality, Stokes’ wave equation, and acoustic pulse pr0p-
agation in a viscous ﬁuid”, Phys. Rev. E 72, 026610 (2005).
[18] P. M. Morse and H. Feshbach, Methods of Theoretical Physics (McGraw-Hill,
New York) (1953).
[19] T. L. Szabo, “Time-domain wave-equations for lossy media obeying a frequency
power-law”, J. Acoust. Soc. Am. 96, 491—500 (1994).
[20] W. Chen and S. Holm, “Fractional Laplacian time-space models for linear and
nonlinear lossy media exhibiting arbitrary frequency power-law dependency”,
J. Acoust. Soc. Am. 115, 1424—1430 (2004).
[21] M. G. Wismer, “Finite element analysis of broadband acoustic pulses through
inhomogeneous media with power law attenuatiOn”, J. Acoust. Soc. Am. 120,
3493—3502 (2006).
[22] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. S. F. Edition), Fundementals
of Acoustics (Fourth Edition ), 210—213 (John Wiley and Sons, Inc., New York)
(2000).
[23] M. F. Hamilton and C. L. Morfrey, “Model equations”, in Nonlinear Acoustics,
edited by M. F. Hamilton and D. T. Blackstock, 41-63 (Academic Press) (1998).
[24] R. 0. Cleveland, M. F. Hamilton, and D. T. Blackstock, “Time-domain mod-
eling of ﬁnite-amplitude sound in relaxing fluids”, J. Acoust. Soc. Am. 99,
3312—3318 (1996).
[25] T. L. Szabo, Diagnostic Ultrasound Imaging: Inside Out (Elsevier Academic
Press, Burlington, MA) (2004).
185
[26] P. A. J. Bascom and R. S. C. Cobbold, “On a fractal packing approach for
understanding ultrasonic backscattering from blood”, J. Acoust. Soc. Am. 98,
3040—3049 (1995).
[27] A. Hanyga and M. Seredyr'iska, “Power-law attenuation in acoustic and isotropic
anelastic media”, Geophys. J. Int. 155, 830—838 (2003).
[28] E. J. Rothwell and M. J. Cloud, Electromagnetics, 216—231 (CRC Press, Boca
Raton, FL) (2001).
[29] J. D. Holmes, W. M. Carey, S. M. Dediu, and W. L. Siegmann, “Nonlinear
frequencydependent attenuation in sandy sediments” , J. Acoust. Soc. Am. 121,
EL218—EL222 (2007).
[30] K. V. Gurumurthy and R. M. Arthur, “A dispersive model for the propagation
of ultrasound in soft tissue”, Ultrason. Imaging 4, 355—377 (1982).
[31] T. L. Szabo, “The material impulse response for broadband pulses in lossy
media”, Proceedings of the IEEE Ultrasonics Symposium, Honolulu, H1 748—
751 (2003).
[32] P. He, “Simulation of ultrasound pulse propagation in lossy media obeying a
frequency power law”, IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 45, 114-
125(1998)
[33] N. V. Sushilov and R. S. C. Cobbold, “Frequency-domain wave equation and
its time-domain solutions in attenuating media”, J. Acoust. Soc. Am. 115,
1431—1436(2004)
[34] R. S. C. Cobbold, N. V. Sushilov, and A. C. Weathermon, “Transient prOp—
agation in media with classical or power-law loss”, J. Acoust. Soc. Am. 116,
3294—3303 (2004).
[35] M. G. Wismer and R. Ludwig, “An explicit numerical time domain formulation
to simulate pulsed pressure waves in viscous fluids exhibiting arbitrary fre-
quency power law attenuation”, IEEE Trans. Ultrason. Ferroelect. Freq. Contr.
42, 1040—1049 (1995).
[36] M. Caputo, “Linear models of dissipation whose Q is almost frequency
independent—II”, Geophys. J. R. Astr. Soc. 13, 529—539 (1967).
[37] R. Schumer, D. A. Benson, M. M. Meerschaert, and S. W. Wheatcraft, “Eu-
lerian derivation of the fractional advection-dispersion equation”, Journal of
Contaminant Hydrology 48, 69—88 (2001).
186
[38] S. Leeman, “Ultrasound pulse propagation in dispersive media”, Ultrasound
Med. Biol. 25, 481— 488 (1980).
I [39] J. M. Blackledge and S. Leeman, “Green’s functions for acoustic ﬁelds in dis-
persive media”, J. Phys. D: Appl. Phys. 16, L247— L250 (1983).
[40] G. V. Norton and J. C. N ovarini, “Including dispersion and attenuation directly
in the time domain for wave propagation in isotropic media”, J. Acoust. Soc.
Am. 113, 3024—3031 (2003).
[41] G. Pinton and G. Trahey, “Full-wave simulation of ﬁnite-amplitude ultrasound
in heterogeneous media”, Proceedings of the IEEE Ultrasonics Symposium,
New York, NY 130—133 (2007).
[42] D. T. Blackstock, “Transient solution for sound radiated into a viscous ﬂuid”,
J. Acoust. Soc. Am. 41, 1312 - 1319 (1967).
[43] M. Liebler, S. Ginter, T. Dreyer, and R. E. Riedlinger, “Full wave modeling of
therapeutic ultrasound: Efﬁcient time-domain implementation of the frequency
power-law attenuation”, J. Acoust. Soc. Am. 116, 2742—2750 (2004).
[44] M. Giona and H. E. Roman, “Fractional diffusion equation on fractals: One
dimensional case and asymptotic behavior”, J. Phys. A: Math. Gen. 25, 2093—
2105 (1992).
[45] R. Metzler and J. Klafter, “The random walk’s guide to anomalous diffusion:
A fractional dynamics approach”, Physics Reports 339, 1—77 (2000).
[46] B. B. Mandelbrot, The Fractal Geometry of Nature (Freeman, New York)
(1983).
[47] E. R. Weibel, “Fractal geometry: A design principle for living organisms”, Am.
J. Physiol. Lung Cell Mol. Physiol. 261, L361—L369 (1991).
[48] J. C. Lockwood and J. G. Willette, “High-speed method for computing the
exact solution for the pressure variations in the nearﬁeld of a baffled piston”,
J. Acoust. Soc. Am. 53, 735—741 (1973).
[49] V. M. Zolotarev, One-Dimensional Stable Distributions (American Mathemat-
ical Society, Providence, RI) (1986).
[50] P. M. Jordan, M. R. Meyer, and A: Puri, “Causal implications of viscous damp-
ing in compressible fluid ﬂows”, Phys. Rev. E 62, 7918 (2000).
187
[51] R. Ludwig and P. L. Levin, “Analytical and numerical treatment of pulsed
wave propagation into a viscous ﬂuid”, IEEE Trans. Ultrason. Ferroelect. Freq.
Contr. 42, 789—792 (1995).
[52] G. C. Gaunaurd and G. C. Everstine, “Viscosity effects on the propagation of
acoustic transients”, J. Vib. Acoust. 124, 19—25 (2002).
[53] G.-P. J. Too, “New phenomena on King integral with dissipation”, J. Acoust.
Soc. Am. 101, 119—124 (1997).
[54] P. T. Christopher and K. J. Parker, “New approaches to linear propagation of
acoustic ﬁelds”, J. Acoust. Soc. Am. 90, 507-521 (1991).
[55] Y.-S. Lee and M. F. Hamilton, “Time-domain modeling of pulsed ﬁnite-
amplitude sound beams”, J. Acoust. Soc. Am. 97, 906—917 (1995).
[56] D. G. Crighton, Modern Methods in Analytical Acoustics: Lecture Notes
(Springer-Verlag, London) (1996).
[57] L. E. Kinsler, A. R. Frey, A. B. Coppens, and J. V. Sanders, Fundamentals of
Acoustics (Fourth Edition), 210—245 (John Wiley and Sons, New York) (2000).
[58] I. Podlubny, “The Laplace rIfansform Method for Linear Differential Equations
of the Fractional Order”, in eprint arXiv:funct-an/9710005, 10005-+ (1997).
[59] J. F. Kelly and R. J. McGough, “The causal impulse response for cicular sources
in viscous media”, J. Acoust. Soc. Am. 123, 2107—2116 (2008).
[60] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, with
Formulas, Graphs, and Mathematical Tables, 887—889 and 916—919 (Dover Pub-
lications, Inc., New York) (1972).
[61] I. Stakgold, Green’s Functions and Boundary Value Problems ( Second Edition),
198 (John Wiley and Sons, Inc., New York) (1998).
[62] F. Oberhettinger, “On transient solutions of the “bafﬂed piston” problem”,
Journal of Research of the National Bureau of Standards-B. Mathematics and
Mathematical Physics 653, 1-6 (1961).
[63] G. R. Harris, “Review of transient ﬁeld-theory for a baffled planar piston”, J.
Acoust. Soc. Am. 70, 10—20 (1981).
[64] J. Djelouah, N. Bouaoua, A. Alia, H. Khelladi, and D. Belgrounde, “Theoretical
and experimental study of the diffraction of focused ultrasonic waves in viscous
ﬂuids”, 2003 World Congress on Ultrasonics, Paris 1347—1350 (2003).
188
[65] P. Chadwick and G. E. Tupholme, “Generation of an acoustic pulse by a bafﬂed
circular piston”, Proc. Edinburgh Math. Soc. 15, 263—277 (1967).
[66] J. A. Jensen and N. B. Svendsen, “Calculation of pressure ﬁelds from arbitrarily
shaped, apodized, and excited ultrasound transducers”, IEEE Trans. Ultrason.
Ferroelect. Freq. Contr. 39, 262—267 (1992).
[67] G. R. Harris, “Transient ﬁeld of a bafﬂed planar piston having an arbitrary
vibration amplitude distribution”, J. Acoust. Soc. Am. 70, 186—204 (1981).
[68] R. Burridge, “An apodized spherical transducer: Analytical solution”, Wave
Motion 26, 13—24 (1997).
[69] J. F. Kelly and R. J. McGough, “An annular superposition integral for axisym-
metric radiators”, J. Acoust. Soc. Am. 121, 759—765 (2007).
[70] K. H. Lee and G. Xie, “A new approach to imaging with low frequency electro—
magnetic ﬁelds”, Geophysics 58, 780-796 (1993).
[71] J. F. Kelly and R. J. McGough, “A time-space decomposition method for cal-
culating the nearﬁeld pressure generated by a pulsed circular piston”, IEEE
Trans. Ultrason. Ferroelect. Freq. Contr. 53, 1150—1159 (2006).
[72] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, with
Formulas, Graphs, and Mathematical Tables, 295—309 (Dover Publications, Inc.,
New York) (1972).
[73] P. J. Davis and P. Rabinowitz, Numerical Integration, 139—140 (Academic, New
York) (1975).
[74] G. E. Tupholme, “Generation of acoustic pulses by bafﬂed plane pistons”, Math-
ematika 16, 209—224 (1969).
[75] P. M. Morse and K. U. Ingard, Theoretical Acoustics, 388 (McGraw—Hill, New
York) (1968).
[76] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, with
Formulas, Graphs, and Mathematical Tables, 360 (Dover Publications, Inc.,
New York) (1972).
[77] M. R. Jennings and R. J. McGough, “A fast nearﬁeld method for calculating
steady-state and transient presssures generated by focused radiators”, Preprint
(2008).
189
[78] D. P. Oroﬁno and P. C. Pedersen, “Multirate digital signal-processing algorithm
to calculate complex acoustic pressure ﬁelds”, J. Acoust. Soc. Am. 92, 563—582
(1992).
[79] R. N. Bracewell, The Fourier Transform and Its Applications, 144 (McGraw-
Hill, New York) (1978).
[80] R. Gorenﬂo, Y. Luchko, and F. Mainardi, “Analytical properties and applica-
tions of the Wright function”, Frac. Calc. Appl. Anal. 2, 383—414 (1999).
[81] A. M. Mathai and R. K. Saxena, The H-function with Applications in Statistics
and Other Disciplines (John Wiley and Sons, New York) (1978).
[82] R. Gorenﬂo and F. Mainardi, “Fractional calculus and stable probability dis-
tributions”, Arch. Mech. 50, 377—388 (1998).
[83] W. G. Glbckle and T. F. Nonnenmacher, “Fox function representation of non-
Debye relaxation processess”, J. Stat. Phys. 71, 741—756 (1993).
[84] W. C. Chew, Waves and Fields in Inhomogeneous Media, 232—235 (IEEE Press,
New York) (1995).
[85] T. Kobayashi and S. Imai, “Spectral analysis using generalized cepstrum”, IEEE
Trans. Acoustics, Speech, and Signal Processing 32, 1087—1089 (1984).
[86] M. J. Lighthill, Introduction to Fourier analysis and the Theory of Distributions
(Cambridge U. Press, Cambridge) (1959).
[87] G. Samorodnitsky and M. Taqqu, Stable non-Gaussian Random Processes
(Chapman and Hall, New York) (1994).
[88] J. P. Nolan, Stable Distributions - Models for Heavy Tailed Data
(Birkhauser, Boston) (2009), in progress, Chapter 1 online at
academic2 . american.edu / ~ j pnolan.
[89] W. Feller, An Introduction to Probability Theory and Its Applications, 165—173
and 548—549 (John Wiley and Sons, New York) (1966).
[90] M. Kanter, “Stable densities under change of scale and total variation inequal-
ities”, Ann. Probab. 3, 697—707 (1975).
[91] H. Pollard, “The representation of e_z’\ as a Laplace integral”, Bull. Amer.
Math. Soc. 52, 908—910 (1946).
[92] M. Dishon, J. T. Bendler, and G. H. Weiss, “Tables of the inverse Laplace
transform of the function e’sﬂ”, Journal of Research of the National Bureau
of Standards-B. Mathematics and Mathematical Physics 95, 433—467 (1990).
190
[93] A. Hanyga and V. E. Rok, “Wave propagation in micro-heterogeneous porous
media: A model based on an integral-differential wave equation”, J. Acoust.
SOC. Am. 107, 2965-2972 (2000).
[94] A. Hanyga, “Propagation of pulses in viscoelastic media”, Pure Appl. Geophys.
159, 1749 — 1769 (2002).
[95] W. R. Schneider, “Stable distributions: Fox function representation and gen-
eralization”, in Stochastic Processes in Classical and Quantum Systems, edited
by S. Albeverio, G. Casati, and D. Merlini (Spring-Verlag, Berlin) (1985).
[96] F. Mainardi and G. Pagnini, “The Wright functions as solutions of the time-
fractional diffusion equation”, Appl. Math. Comput. 51—62 (2003).
[97] R. Metzler, W. G. Gléckle, and T. F. Nonnenmacher, “Fractional model equa-
tion for anomalous diffusion”, Physica A. 211, 13—24 (1994).
[98] M. Galassi, GNU Scientiﬁc Library Reference Manual (2nd Ed.) (Network The-
ory Ltd., Bristol, UK) (2006).
[99] G. Robinson, “Rapid computations concerning log maximally-skewed stable
distributions, with potential use for pricing Options and evaluating portfolio
risk”, Private Report (2005).
[100] N. Denisenko, G. Scarano, M. Matteucci, and M. Pappalardo, “An approximate
solution of the transient acoustic ﬁeld”, IEEE Trans. Ultrason. Ferroelect. Freq.
Contr. 46, 821—827 (1985).
[101] J. L. San Emeterio and L. G. Ullate, “Diffraction impulse response of rectan-
gular transducers”, J. Acoust. Soc. Am. 92, 651—662 (1992).
[102] J. A. Jensen, D. Gandhi, and W. D. O’Brien, “Ultrasound ﬁelds in an attenu-
ating medium”, Proceedings of the IEEE Ultrasonics Symposium, Baltimore,
MD 943—946 (1993).
[103] J. F. Kelly and R. J. McGough, “Tfansient acoustic ﬁelds produced by rectan-
gular apertures and linear arrays in viscous media”, Proceedings of the IEEE
Ultrasonics Symposium, New York, NY 1553—1556 (2007).
[104] J. F. Douglas, “Polymer science applications of path-integration, integral equa-
tions, and fractional calculus”, in Applications of Fractional Calculus in Physics,
edited by R. Hilfer, 241—330 (World Scientiﬁc) (2000).
[105] M. Dishon, G. H. Weiss, and J. T. Bendler, “Stable law densities and linear re-
laxation phenomena”, Journal of Research of the National Bureau of Standards-
B. Mathematics and Mathematical Physics 90, 27—39 (1985).
191
[106] W. Gawronski, “On the bell-shape of stable densities”, Annals of Probability
12, 230—242 (1984).
[107] D. G. Hummer, “Rational approximations for the Holtsmark distribution, its
cumulative and derivative”, Journal of Quantitative Spectroscopy and Radiative
Transfer 36, 1—5 (1986).
[108] A. C. Kak and K. A. Dines, “Signal-processing of broad-band pulsed ultrasound
- Measurement of attenuation of soft biological tissues”, IEEE Trans. Biomed.
Eng. 25, 321—344 (1978).
[109] W. J. Fry, “Mechanism of acoustic absorption in tissue”, J. Acoust. Soc. Am.
24, 412—415 (1952).
[110] A. Yeung and E. Evans, “Cortical shell-liquid core model for passive ﬂow of
liquid-like spherical cells into micropipets”, Biophys. J. 56, 139-149 (1989).
[111] E. M. Darling and M. Topel and S. Zauscher and T. P. Vail and F. Guilak,
“Viscoelastic properties of human mesenchymally—derived stem cells and pri-
mary osteoblasts, chondrocytes, and adipocytes”, Journal of Biomechanics 41,
454—464 (2008).
[112] C. T. Lim, E. H. Zhou, and S. T. Quek, “Mechanical models for living cells—A
review”, Journal of Biomechanics 39, 195—216 (2006).
[113] J. D. Humphrey, “Continuum biomechanics of soft biological tissues”, Proc. R.
Soc. Lond. A 459, 3—46 (2003).
[114] C. Verdier, “Rheological properties of living materials. From cells to tissues”,
J. Theoretical Medicine 5, 67—91 (2003).
[115] N. W. Tschoegl, The Phenomenological Theory of Linear Viscoelastic Behavior:
An Introduction (Springer-Verlag, Berlin) (1989).
[116] B. Gross and R. M Fuoss, “Ladder structures for representation of viscoelastic
systems”, J. Polymer Sci. 19, 39—50 (1956).
[117] H. Schiessel and A. Blumen, “Hierarchical analogues to fractional relaxation
equations”, J. Phys. A.: Math. Gen. 26, 5057-5069 ( 1993).
[118] H. Schiessel and A. Blumen, “Mesoscopic pictures of the sol-gel transition:
Ladder models and fractal networks”, Macromolecules 28, 4013—4019 (1995).
[119] H. Schiessel, R. Metzler, A. Blumen, and T. F. Nonnenmacher, “Generalized
viscoelastic models: Their fractional equations with solutions”, J. Phys. A:
Math. Gen. 28, 6567—6584 (1995).
102
[120] J. F. Kelly, M. M. Meerschaert, and R. J. McGough, “Analytical time-domain
Green’s functions for power-law media”, J. Acoust. Soc. Am. 124, 2861—2872
(2008).
[121] R. L. Bagley and P. J. Torvik, “A theoretical basis for the application of frac-
tional calculus to Viscoelasticity”, J. Rheol. 27 , 201—210 ( 1983).
[122] B. Gross, “Ladder structures for representation of viscoelastic systems. 11”, J.
Polymer Sci. 20, 123—131 (1956).
[123] N. Heymans and J .-C. Bauwens, “Fractal rheological models and fractional
differential equations for viscoelastic behavior”, Rheol. Acta 33, 1994 (210-
219).
[124] P. E. Rouse, “A theory of linear viscoelastic properties of dilute solutions of
coiling polymers”, J. Chem. Phys. 21, 1272—1280 (1953).
[125] A. Kreis and A. C. Pipkin, “ViscOelastic pulse propagation and stable proba-
bility distributions”, Q. Appl. Math. 44, 353—360 (1986).
[126] T. L. Szabo and J. Wu., “A model for longitudinal and shear wave propagation
in viscoelastic media”, J. Acoust. Soc. Am. 107, 2437—2446 (2000).
[127] P. M. Morse and K. U. Ingard, Theoretical Acoustics, 278 (McGraw-Hill, New
York) (1968).
[128] A. Papoulis, The Fourier Integral and its Applications, 213-217 (McGraw-Hill,
New York) (1987).
[129] A. A. Kilbas, H. M. Srivastava, and J. J. 'Ifuhillo, Theory and Applications of
Fractional Differential Equations (Elsevier, Amsterdam) (2006).
[130] J. W. Baish and R. K. Jain, “Fractals and cancer”, Perspectives in Cancer
Research 60, 3683—3688 (2000).
[131] R. Orbach, “Dyanamics of fractal networks”, Science 231, 814—819 (1986).
[132] R. Rammal and G. Toulouse, “Random walks on fractal structures and perco-
lation clusters”, J. Physique 44, L13—L22 (1983).
[133] G. Wojcik, J. Mould, F. Lizzi, N. Abboud, M. Ostromogilsky, and D. Vaughn,
“Nonlinear modeling of therapeutic ultrasound”, Proceedings of the IEEE Ul-
trasonics Symposium, Cannes, France 1617— 1622 (1995).
[134] M. Ochmann and S. Makarov, “Representation of the absorption of nonlinear
waves by fractional derivatives”, J. Acoust. Soc. Am. 94, 3392—3399 (1993).
193
[135] J. A. Jensen, “Ultrasound ﬁelds from triangular apertures”, J. Acoust. Soc.
Am. 100, 2049-2056 (1996).
[136] P. Faure, D. Cathignol, J. Y. Chapelon, and V. L. Newhouse, “On the pressure
ﬁeld of a transducer in the form of a curved strip”, J. Acoust. Soc. Am. 95,
628—637 (1994).
[137] M. Arditi, F. S. Foster, and J. W. Hunt, “Transient ﬁelds of concave annular
arrays”, Ultrason. Imaging 3, 37—61 (1981).
[138] D. Chen, J. F. Kelly, and R. J. McGough, “A fast near-ﬁeld method for calcula-
tions of time-harmonic and transient pressures produced by triangular pistons”,
J. Acoust. Soc. Am. 120, 2450—2459 (2006).
194
~~