é... ; :21! L. 4.; 25w e. C. 3.1,.353 s . . a . m»? h. .r. .a , “w. i : x. .. . .. 3.3. $- nugn s unflnvuxs ".‘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 .firfifafiafl ..(‘...oxf.~o.1¢ :. 1.1:! 0.7.0, )I II Ilfir i,‘4¥:r§!”ll 2.! i. 3'! 4.? .3415! l. ...‘>I$tc,’z 9t 1‘ .. .58.:2."avfl.uv (a 01?. 11’! all} 5.. {$1. (2... ..v,« I... I! .1112 17v. Iii. «v.51! kl. In l 173!!!» aa‘mfiizé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 fulfillment of the requirements for the PhD. degree in Electrical Engineen‘gq Mo V. I 3r 2008 Date MSU is an Affirmative 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 fulfillment 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 coefficient 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 fields 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 verified using the material impulse response function (MIRF) approach. A lossy impulse response for the Stokes wave equation is derived for calculations of transient fields generated by finite aperture radiators. Expressions for the uniform circular aper- ture (in both the nearfield and the farfield), the uniform rectangular aperture in the nearfield, and the spherical shell in the nearfield 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 infinite- 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 finite 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: Nearfield ................... 3.2.1 Lossless Media .......................... 3.2.2 Viscous Media ........................... 3.2.3 Numerical Results ........................ 3.3 Uniform Circular Piston: Farfield .................... 3.3.1 Lossless Media .......................... 3.3.2 Viscous Media ........................... 3.3.3 Numerical Results ........................ 3.4 Rectangular Sources: Nearfield ..................... 3.4.1 Lossless Media .......................... 3.4.2 Viscous Media ........................... CDCDQOTCfiuD-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 Simplified 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: Nearfield ........................ 153 7.3.1 Numerical Results ........................ 155 7.4 Uniform Circular Piston: Farfield .................... 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 Definitions ................................. 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 coefficients for kidney, liver, and tendon taken from [1]. The data was fit 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 infinite 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 verification. ................................ 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 field 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 nearfield impulse response and the farfield 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 nearfield 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 infinite rigid baffle in the z = 0 plane. Reprinted from [2]. ..................................... Lossy nearfield 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 infinite 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 baffled 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 defined 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 coefficient 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 defined 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 coefficient 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 coefficient 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 first 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- nification. 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 specific biological structures vary at each successive spatial scale, the essential features are the same: fluid 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 first panel displays an infinite 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 first 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 coefficients 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 five levels of magnification. 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 defined 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 nearfield impulse response and the farfield 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 nearfield and farfield impulse responses evaluated at r = 200 mm on-axis (t9 = O) and off-axis (6 = 1r/ 12 and 43:0) .................................... Lossy nearfield 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 ramifications 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 field 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 coefficient 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 coefficient 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 confirm that the attenuation coefficient 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 coefficient 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 finite 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 fields 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 fluids (except superfluids) possess an internal friction called viscosity, which results from the diffusion of momentum between fluid particles with different velocities [14]. Since momentum diffusion is a localized near velocity gradi- ents, a viscous fluid is memoryless and only weakly dispersive. As shown in standard textbooks [14, 22], viscous dissipation of sonic energy results in an attenuation coef- ficient 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 rarefied 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 fluid, 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 fitted to experimental data. Although relaxation models may be fitted 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 conflicts with observed attenuation and backseatter data. To resolve this conflict, 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 coefficient. 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 coefficient. In addition, tissue possesses both solid-like and fluid-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 fluid-like extra-cellular matrix (ECM) [4]. The ECM is often modeled as an aqueous solution of viscoelastic polymers, which possess both solid and fluid-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 coefficient, 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 coefficient, on is the attenuation slope, and y is the power law exponent. Measured attenuation coeflicients 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. Stratified 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 coefficients as a function of frequency [1]. Attenuation coefficients, 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 fitted 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 coeflicient’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 filters, 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 coefficients for kidney, liver, and tendon taken from [1]. The data was fit 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 coefficient 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 flexible, 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 defined 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 finite-difference time-domain (FDTD) [40, 41] or finite 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 first two initial conditions. To explore the ramifications 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 defined 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 defined 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 fields produced by finite 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 fields. 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 fluids 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], infinite 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 sufficiently far from the source and 2) moderate values of the attenuation coefficient. The Green’s function is then decomposed into diffraction and loss components, where the diffraction component satisfies the lossless wave equation and the loss component satisfies 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 fluid 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 coefficient 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 first 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 infinite series. 1 The first 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 coefficients [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) satisfies 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 first 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 defined 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 dflk 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 infinite propagation speed, yet is strictly causal. As the following section shows, Eq. (2.20) is an excellent, causal approximation except in the extreme nearfield 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 first 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 definition 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 specified 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, defined 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 significant 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 fixing R and identifying the range of 7 satisfying the 1% error threshold, yielding the relation 7 S 2 x 10-4R. By fixing 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 coeflicient a at a fixed frequency f via the relation [57] 7 = COO/(2W2f2). For instance, soft tissue at f = 3 MHz with an attenuation coefficient 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) flimexp 0%)] dt’, (2.30) where the sifting property of the Dirac delta function applied to Eq. (2.30) yields Eq. (2.20). Identifying the first 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 finite 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 - fl 1.2 ' ' ' M —relerence _'°f°'°"_°° 0.35 a - uappmxlmafion " ' '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 difl'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 finite 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 field radiated by finite 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 field points off-axis were not considered. In particular, the effect of viscous loss on the impulse response of baffled 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 fields produced by finite apertures [66]. The (lossless) impulse response theory assumes a homogeneous, linear acoustic medium with a finite radiating aperture S surrounded by an infinite, rigid baffle. 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 defined by integrating Eq. (2.20) over the surface of the radiating aperture and multiplying by two to account for the infinite, 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 finite 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 nearfield solution given by Eq. (3.7) is an exact solution to the transient wave equation assuming a uniform piston in an infinite rigid baffle. 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 first 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 infinitely 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 nearfield 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 field 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 nearfield, 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 nearfield 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 field becomes progressively closer to that generated by an omni-directional source as the directivity of the aperture is weakened far from the source. The nearfield velocity potential field 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 significant 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 first 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: Farfield Although Eq. (3.9) is valid for all field points, a simpler expression exists in the farfield, where the effects of attenuation are more pronounced. Assuming the defi- nition of the farfield utilized by Morse and Ingard [75], the observation distance is assumed to be much larger than the radius of the aperture (7' >> a) . This definition of the farfield is distinguished from the classical farfield distance r > a2/A, which is only valid over a finite frequency band. 3.3. 1 Lossless Media A single-integral representation of the farfield lossless impulse response function is derived from the well-known frequency-domain expression. In the farfield 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-field pattern and k = w/c is the (lossless) wavenumber. Eq. (3.16) is computed by inverse Fourier transformng the classical frequency-domain result. The farfield 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 first 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 final result h(r, 6, t) = f“ costb u(t — r/co + asin dcos rb/co) dip, (3.16) firmnd 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 nearfield 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 nearfield and the farfield, 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 farfield 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_/fl 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/yfl), (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 infinite 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-infinite region of support implies wave components traveling with phase speeds varying from zero to infinity, but the infinite phase speeds are infinitely attenuated [17], so the result is causal. 3.3.3 Numerical Results Figure 3.6 displays both the nearfield and farfield 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 farfield approximation differs significantly from the nearfield 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 farfield approximation agrees with the nearfield solution, indicating that the simplified farfield 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 nearfield, the transition region, and the farfield. For the values of 7 evaluated in Fig. 3.6, the diffraction component dominates in the nearfield, 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 farfields, the effects of diffraction and loss are both apparent. In this transition region, the impulse response is both smooth and asymmetric. Finally, in the farfield, the effects of viscous loss predominate, yielding increasingly symmetric and broad responses with a reduced directivity. 3.4 Rectangular Sources: Nearfield 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 fields 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 nearfield 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 eilcvrzgflrifisQ _ 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 nearfield 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 specified 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 nearfield solution is provided for all field points. The solution is derived for a spherical shell centered at (0,0, —R) and surrounded by an infinite, 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 infinite 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 + zzJRfl—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 first 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 significantly 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 five-cycle Hamming-weighted toneburst in Eq. (3.11) with center frequency f0 = 5 MHz. For liver, the attenuation coefficient 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 finite 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 fields 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 efficient 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 benefit 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 field 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 artifi- 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 magnified on-axis and in the farfield, 43 where the lossless impulse response is represented by a short-duration rect function in the nearfield and a scaled delta function in the farfield. However, these numerical difficulties 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 filter and removing the high-frequency components of the impulse response. This filtering ef- fect is reflected in the lossy impulse response expressions, which are infinitely 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 difficulties 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 first 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 fields 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 benefits 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 finite 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 verification. 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 fl 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 field 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 nearfield impulse response and the farfield 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 nearfield and farfield 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 infinite rigid baffle 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 nearfield 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 infinite 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. ... ”fill/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 coefficient 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 coefficient 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 identified 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 reflection 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 first and third temporal derivative, respectively. For fractional y, the third term is defined by the Riemann-Lionville fractional derivative defined 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 definition 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 infinitely high frequency, and do is the attenuation coefficient 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 modified 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 coefficient 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, defined 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 wavefield is large near the source, Eq. (4.10) fails to capture the high-frequency content in the 62 extreme near-field 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 field radiates away from the source. However, Eq. (4.13) is noncausal, which means that Eq. (4.13) specifies non-zero field values fort < 0. However, for (10 small and R sufficiently 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 defined 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 infinite 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 coefficient and 2) a frequency—dependent phase speed. The power law dispersion relationship that satisfies 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 coefficient 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 fi = ;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 coefficient 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————(fi—-——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. Defining 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 —» --+-2fl’-(lw|‘—1) c1 we where the limit tanz —> z is used in the last line. To finish 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 satisfies 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-lflflw) - 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 fi(R,w) = [exP(:j’,fi/C°’] [exp («04:0wa — Haney/244142540] (4.36) where the first 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) filters, Eq. (4.38) combines the effects of dispersive loss and diffraction as a cascade of two system functions. The first 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 defined via gL(R, t) = 3:4 [exp (—aoR(|w|y — itan(7ry/2)w|w|y_l))] . (4.39) Interestingly, the inverse Fourier transform defined 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 identified 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 defined 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 infinitely 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 first 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 defined 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 defined 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 confluent 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 coefficients 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 infinitely 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 first 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/fisl/y) ), (4.48) where 1 < y S 2. Eq. (4.48) may be simplified 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 confluent 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 infinitely 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 infinite 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 first 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 (fiop008h¢)’/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 (flopcosh ¢)1/yfy ( (flopcosh ¢)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 first 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 finite range of integration u (t —- p/coo) /°°Sh_1(C°Ot/p) 1 (t — pcosh (D/coo) dtb. 27r 0 gm” = (flopcosW/y ” (fiopcosh 4)”?! (4.63) Although the range of integration is infinite 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 finite. 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 definition 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 fl, 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 defined 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 Scientific 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 verified 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 fields 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 figures, the impulsive excitation has been stretched and smoothed by the filtering 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 field 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 (—fit) 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 significant 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 significant 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 modified 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 justification 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&_fl(_ 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 defined 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 find an approximate expression for the attenuation coefficient 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 defined, 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 infinite series given by Eq. (2.14), where each term is specified 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 first 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 defined by Eq. (8.3) with skewness fl = 0. These distributions, which are defined in Appendix B, have been studied extensively [49] and possess the following analytical properties 101 1. wy(t) is infinitely 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 verified 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 final 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 infinite 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- fined 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)/(fi060t)1/y y 110 Identifying the skewed, stable cumulative distribution function (cdf) Fy(t) defined 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 defined 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 finite 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(———(fiocot)1/y)q(v)dv (5.61) As with Eq. (5.58), Eq. (5.61) is finite 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 significant 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 defined 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 (flocot)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 coefficient 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 coefficient 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 infinitely-divisible distribution [89], resulting from the summation of many sim- ilar random variables. Condition 1 is required since the coefficient 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 defined 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 fl, 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 coefficient 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 coefficient 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 fitted 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 infinite 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 infinite 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 satisfies 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 first 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 satisfies a power law attenuation coefficient 01(0)) given by Eq. (1.1) over an appropriate bandwidth of frequencies for power law exponents 1 < y S 2. Specific tissue, such as breast, consist of hierarchical arrangements of elastic and fluid-like components. For instance, breast consists of interlobular stroma, which has fluid 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 fluid-like extra-cellular matrix (ECM). The ECM is often modeled as an aqueous solution of viscoelastic poly- mers, which possess both solid and fluid-like properties. Individual cells are modeled as elastic membranes containing fluid-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 fluid-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 fluid, 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 magnification. 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 fluid-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 fluid 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 fluid with a characteristic coefficient of viscosity. The cell nucleus may also be included as an additional elastic component embedded within the viscous fluid. 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 infinite 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 first 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 magnification. The third panel (A / 40 ~ 5 pm) displays the cell nucleus, consisting of a double membrane, a fluid-like nucleoplasm, and an elastic nucleolus in the interior, which contains chromatin [4]. Although the specific biological structures vary at each successive spatial scale, the essential features are the same: fluid substrates containing elastic compartments. This self-similar pattern forms the basis for the fractal structure shown in Fig. 6.2. 126 compressible, fluid. By allowing an infinite 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 indefinitely, 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) fluid 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 coefficients 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 magnification 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 fluid relating the strain 6,-3- and stress Ti]- tensors is generalized. The strain tensor is defined 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 fluid 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 coefficient 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 fluids, )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 coefficient 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 sufliciently 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 infinite, 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 first panel displays an infinite number of thin elastic membranes with Young’s modulus E alternating with viscous compartments that have coefficients of viscosity 1]. The second panel zooms in on the first 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 defined 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—fl 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_fi and 7 = 1/2 is embedded within a larger ladder along with springs with elastic coefficients 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 fluid. 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 justification for constitutive equations similar to Eq. (6.13) were also established [121] on the basis of dilute solutions of polymers [124] within a Newtonian fluid 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 coefficient of viscosity 1— p = E0 37763 (6.14) and repeating the steps performed in the fi = 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 + #atfl_1 (axj + 55;) (6.15) Eq. (6.15) may also be written in dyadic notation as —l afi—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 fluid 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 flow. In a viscous fluid (,6 = 1), the frictional term involves only a spatial derivative and is purely local. For a homogeneous fluid 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 flow. 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 fluid 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 coefficient 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 field 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 flow and 2) a temporal integration operator which incorporates the past history of the flow into the flow’s present state. The divergence operator is applied to both sides of Eq. (6.21), yielding 3V ' u — _V2 +6. fl 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 defined 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 coefficient 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 afi-l u(r) ) —at—— — V (p0(r)) + 3atfi-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), define 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 fl = 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 defined by Eq. (6.31). 140 6.3.3 Simplified Equations Eq. (6.32) can be simplified considerably in several important special cases. In most biological tissue, the density is nearly constant p0(r) z p0. Defining 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+ atfiv (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 fl 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 atfi = 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 briefly reviewed in the following subsections. First, the dissipative and dispersive properties of Eq. (6.26) and shown to satisfy both a power law attenuaton coeffi- 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 coefficient 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? — k2Tfi(—iw)5 = 0. (6.37) 0 Solving for the wavenumber k(LU) yields (6.38) 2(6) = w . c0(/ 1 + (—z'm)fi . In the low frequency limit, the binomial approximation is applied, yielding ~ w TB ‘ flrr fl 3rfl , ,61r fl k(w) ~ g (1— —2— co (7) w + 2—2— sm (7) w . (6.39) The attenuation coefficient 0(a)) and phase velocity C(w) are calculated by inserting fl = y — 1 such that y—l r cos(7ry/2)wy. 2co a(w) = Im h(w) = (6.40) Defining 20 = T’Hl C0~°»(7ry/2)|/(260) (6-41) yields a power-law coefficient 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 coefficient 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 coefficient approaches __) Sin("fl/4) l—fi/2 0(6)) ———COTfi/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 specified 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 defined 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 significance 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 coefficient 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 fluids like castor oil and silicone fluid [7]. Since springs correspond to elastic structures such as cellular and nuclear mem- branes, while dampers correspond to inter-cellular and intra-cellular fluids 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 fluid-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 coefficient 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 field. However, sound speed inhomogeneity is not solely responsible for the observed power-law dependence of the attenuation coefficient. In addition, a viscous mechanism is required to dissipate both the incident and incoherently scattered sound fields. 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 significantly 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 fl 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 specific 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 five levels of magnification. 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 sufficiently fine spatial scale, biological tissue may be regarded as a heterogeneous, Newtonian fluid with spatially varying material properties such as the coefficients 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 justified in most biological media due to the negligible thermal conductivity, the linear assumption is not justified 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 finite 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 finite 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 finite 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 finite 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 satisfies the following desirable criteria: 1. In the low-frequency limit, the attenuation coefficient satisfies 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 (flocot)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 nearfield, 2) the uniform circular piston the farfield, 3) the rectangular piston in the nearfield, 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 fields 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 — :fl'éR/Col] [(60:31/96 (—(fi0c: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 infinite baffle 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 definition 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: Nearfield 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 identified as a definite 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?(QAflMa(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) flow-AW (fiotcowy 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 coefficient 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 first tested in lossless media using an artificially 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 significant 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 fixed, the decay of the impulse response slows, and the duration of the wake increases. In Fig. 7.2(a), which lies within the nearfield of the radiator (z = 50 mm), the impulse responses in all three media display both the “stretching” due to the finite extent of the aperture, and the “spreading” due to dissipation. However, in Fig. 7.2(b), which is in the farfield 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 defined 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 nearfield of the piston, the velocity potential has a significantly 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 finite extent of the aperture, is not pronounced. However, the frequency down-shifting and amplitude attenuation are evident. Hence, for observation points in the farfield 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: Farfield 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 fl - 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(fl0¢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 nearfield and farfield 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 nearfield impulse response and the farfield 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 nearfield and farfield 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 significant disparity between the nearfield and farfield responses, especially on-axis, thus indicating that the farfield approximation in Eqns. (7 .15) and (7.16) is invalid 160 0.08~ j j ' ' ' -‘_.fat . ---fiver ’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 defined 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 nearfield and farfield results are in closer agreement. Unlike the nearfield and farfield 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 farfield 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 nearfield 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 significant 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- ified by h (t=flP 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) fl :.'. ...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 "'-'ofl-axis (near) E 0.3 . .mno Ofl—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 nearfield impulse response and the farfield 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 nearfield and farfield 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 farfield 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 (figfladE) 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) — (fiocoflVy y (_(fiocot)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 filtered 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 nearfield 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) simplifies 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)] (flocot)1/y (flocot)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 bafl‘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 finite 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 fixing 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 field 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 defined, corresponding to randomized delay times with a PDF given by Eq. (5.66). The nearfield impulse response given by Eq. (7.10) is then expressed as hpcr. 2. t) = @fif—(t) f0” Mae. «7) [P(T s n) -. P(T s 72)] cm. (7.24) Physically, the first 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 fixed 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 flows 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 nearfield, the transition region, and the farfield. For the values of 010 and 169 y evaluated in Fig. 7.2, the diffraction component dominates in the nearfield, 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 nearfield and farfield, the effects of diffraction and loss are both apparent. In this transition region, the impulse response is both smooth and asymmetric. Finally, in the farfield, 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 nearfield 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 fields 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 infinite 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 satisfies the lossless wave equation and the loss component satisfies 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 nearfield 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 farfield of a circular piston in Eqns. (3.21) and (3.22), the nearfield of a rectangular piston in Eq. (3.30), and the baffled spherical 171 shell in Eq. (3.35). Velocity potential computations display a significant 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 farfield. 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 coefficient 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 verified against the MIRF result [31] and Hilbert dispersive model [30], and example fields 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 deficiency, 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 fluid 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 nearfield circular piston (Eq. (7.10)), farfield 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 nearfield circular piston is numerically verified 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 defined 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—ngtfif—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 defined via spatial Fourier transforms [129] (472)”2 we) a 5“ mama]. (45) where \I'(k) is the spatial Fourier transform of 2,1)(r) defined 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 define 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 Definitions Stable random variables are defined via the characteristic function w(k) [87] ‘I’(k) = exp {—|k|aoa [1 —ifisgn(k)tan(1ra/2)] +ikp}, (B.l) where a # 1 is the index, lfil 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 + 2ifi/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 defined 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 defined 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 defined via 2 P(x) =/ p(rc’) dx'. (B9) -00 The maximally skewed stable CDF is defined as t Fy(t) =/ fy(t’)dt', (B.10) —00 while the symmetric stable CDF is defined 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. Definitions and properties of these transcendental functions are therefore tabulated below. C.1 The Fox-H Function The Fox H -function [81, 129] is defined 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 defined 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 infinite 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 defined by an infinite series (see Eq. (1.11.1) in Ref. [129]): Z]: ¢(aaIB;Z) E Igm-fi. (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 definition of the H -function, in Ref. [129]) via which possesses an infinite number of poles, with Cauchy’s residue theorem. The resulting infinite series is identified as Eq. (C3). The Wright function is differentiated and integrated via the following relation [129] d n (a) ¢(a,fl; z) = ¢(a, a + nfi; 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 nearfield 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 verification 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 fluid-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 fiuid”, 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 finite-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 fields 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 finite-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 fluid”, 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: Efficient 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 nearfield 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 flows”, 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 fluid”, 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 fields”, J. Acoust. Soc. Am. 90, 507-521 (1991). [55] Y.-S. Lee and M. F. Hamilton, “Time-domain modeling of pulsed finite- 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 “baffled 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 field-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 fluids”, 2003 World Congress on Ultrasonics, Paris 1347—1350 (2003). 188 [65] P. Chadwick and G. E. Tupholme, “Generation of an acoustic pulse by a baffled circular piston”, Proc. Edinburgh Math. Soc. 15, 263—277 (1967). [66] J. A. Jensen and N. B. Svendsen, “Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers”, IEEE Trans. Ultrason. Ferroelect. Freq. Contr. 39, 262—267 (1992). [67] G. R. Harris, “Transient field of a baffled 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 fields”, Geophysics 58, 780-796 (1993). [71] J. F. Kelly and R. J. McGough, “A time-space decomposition method for cal- culating the nearfield 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 baffled 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 nearfield method for calculating steady-state and transient presssures generated by focused radiators”, Preprint (2008). 189 [78] D. P. Orofino and P. C. Pedersen, “Multirate digital signal-processing algorithm to calculate complex acoustic pressure fields”, 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. Gorenflo, 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. Gorenflo 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’sfl”, 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 Scientific 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 field”, 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 fields 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 fields 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 Scientific) (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 flow 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 fields 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 field 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 fields of concave annular arrays”, Ultrason. Imaging 3, 37—61 (1981). [138] D. Chen, J. F. Kelly, and R. J. McGough, “A fast near-field method for calcula- tions of time-harmonic and transient pressures produced by triangular pistons”, J. Acoust. Soc. Am. 120, 2450—2459 (2006). 194