TIME-DOMAIN ANALYSIS OF FRACTIONAL WAVE EQUATIONS AND IMPLEMENTATIONS OF PERFECTLY MATCHED LAYERS IN NONLINEAR ULTRASOUND SIMULATIONS By Xiaofeng Zhao A DISSERTATION Submitted to Mi higan State University in partial fulllment of the requirements for the degree of Ele tri al Engineering - Do tor of Philosophy 2018 ABSTRACT TIME-DOMAIN ANALYSIS OF FRACTIONAL WAVE EQUATIONS AND IMPLEMENTATIONS OF PERFECTLY MATCHED LAYERS IN NONLINEAR ULTRASOUND SIMULATIONS By Xiaofeng Zhao The attenuation of ultrasound propagating in human tissue follows a power law with respe t to frequen y that is modeled by several dierent fra tional partial dierential equations. These models for the power law attenuation of medi al ultrasound have been developed using fra tional al ulus, where ea h ontains one or more time-fra tional or spa e-fra tional derivatives. To demonstrate the similarities and dieren es in the solutions to ausal and non ausal fra tional partial dierential equations, time-domain Green's fun tions are al ulated numeri ally for the fra tional wave equations. For three time-fra tional wave equations, namely the power law wave equation, the Szabo wave equation, and the Caputo wave equation, these Green's fun tions are evaluated for water with a power law exponent of y = 2, of liver with a power law exponent of y = 1.5. y = 1.139, and breast with a power law exponent Simulation results show that the non ausal features of the numeri ally al ulated time-domain response are only evident in the extreme neareld region and that the ausal and the non ausal Green's fun tions onverge to the same time-domain waveform in the fareld. When non ausal time-domain Green's fun tions are onvolved with nite-bandwidth signals, the non ausal behavior in the time-domain is eliminated, whi h suggests that non ausal time-domain behavior only appears in a very limited set of ir umstan es and that these time-fra tional models are equally ee tive for most numeri al al ulations. For the al ulation of spa e-fra tional wave equations, time-domain Green's fun tions are numeri ally al ulated for two spa e-fra tional models, namely the Chen-Holm and TreebyCox wave equations. Numeri al results are omputed for these in breast and liver. results show that these two spa e-fra tional wave equations are ausal everywhere. The Away from the origin, the time-domain Green's fun tion for the dispersive Treeby-Cox spa efra tional wave equation is very similar to the time-domain Green's fun tions al ulated for the orresponding time-fra tional wave equations, but the time-domain Green's fun tion for the nondispersive Chen-Holm spa e-fra tional wave equation is quite dierent. To highlight the similarities and dieren es between these, time-domain Green's fun tions are ompared and evaluated at dierent distan es for breast and liver parameters. When time-domain Green's fun tions are onvolved with nite-bandwidth signals, the phase velo ity dieren e in these two spa e-fra tional wave equations is responsible for a time delay that is espe ially evident in the fareld. The power law wave equation is also utilized to implement a perfe tly mat hed layer (PML) for numeri al al ulations with the Khokhlov - Zabolotskaya - Kuznetsov (KZK) equation. KZK simulations previously required a omputational grid with a large radial distan e relative to the aperture radius to delay the ree tions from the boundary. To de rease the size of the omputational grid, an absorbing boundary layer derived from the power law wave equation. Simulations of linear pressure elds generated by a spheri ally fo used transdu er are evaluated for a short pulse. Numeri al results for linear KZK simulations with and without the absorbing boundary layer are ompared to the numeri al results with a su iently large radial distan e. Simulation results with and without the PML are also evaluated, where these show that the absorbing layer ee tively attenuates the wavefronts that rea h the boundary of the omputational grid. Copyright by XIAOFENG ZHAO 2018 ACKNOWLEDGEMENTS The last six years have been a very rewarding experien e. I an only make su h a statement due to the diligen e of others who have made a profound impa t in my life. First and foremost, I would like to express my sin ere gratitude to Dr. Robert M Gough for his gra ious support throughout my entire Ph.D. program. Dr. M Gough's reativity, enthusiasm, and problem-solving a umen always inspired me in pursuing every resear h goal. He has been a great mentor and has always made sure I understand the importan e of good s ienti work. He has made a great deal of eort to ensure my su ess as a Ph.D. student. I am very grateful for his never-ending help, patien e, understanding, and en ouragement in past few years. I would also like to thank my ommittee members: Dr. Edward Rothwell, Dr. Selin Aviyente, and Dr. Brian Feeny for their mentorship and invaluable advi e. I really appre iate their guidan e in my Ph.D. program and instru tion in my ourse study of Ele tromagneti Fields and Time-Frequen y Wavelet Analysis. Furthermore, I would like to thank my olleagues from the Biomedi al Ultrasoni s and Ele tromagneti s Laboratory: James Kelly, Drew Murray, Kirk Sales, Peter Beard and Leslie P. Thomas. It is my honor to have the opportunity to work with all of them. In addition, I would like to give spe ial thanks to Bin Fan, Yiqun Yang, Pedro Nariyoshi and Jie Li. Thank you for your support, ompany, and help over the past few years. Finally, I would like to thank my family for all of their en ouragement and support throughout my life and studies. To my mom and dad, the things I have done in my life are only possible as result of your sa ri e. I moved to an unknown ountry and ould not a ompany you in the last few years. I have not always been the best son or herished this sa ri e. I want to let you know that I have learned from my mistakes that I have made in life and that I love you from the bottom of my heart. v TABLE OF CONTENTS LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 1 1.1 Ultrasound attenuation in soft tissue . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Nonlinear ultrasound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Fra tional derivative operators . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Dispersion relations for fra tional wave equations 1.5 . . . . . . . . . . . . . . . 5 . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4.1 First order approximation 1.4.2 Se ond order approximation . . . . . . . . . . . . . . . . . . . . . . . 7 Thesis stru ture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Chapter 2 Time fra tional wave equations . . . . . . . . . . . . . . . . . . . . 9 2.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Power law attenuation and dispersion . . . . . . . . . . . . . . . . . . . . . . 10 2.3 The Szabo wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.4 The power law wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.5 The Caputo wave equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.6 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.7 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.7.1 Time-domain Green's fun tions for a ousti propagation in water . . 18 2.7.2 Time-domain Green's fun tions for a ousti propagation in breast . . 20 2.7.3 Time-domain Green's fun tions for a ousti propagation in liver . . . 21 2.7.4 Verti al axis s aling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 2.7.5 Comparisons between time-domain Green's fun tions for a ousti prop- 2.7.6 Convergen e of the Green's fun tions for a ousti propagation in the agation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . time domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7.7 . 2.7.8 Time-domain Green's fun tions onvolved with a three y le Hanning- 2.7.9 Chara terizing the dispersion of weighted pulse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v(t) ∗ 4πrg(r, t) 26 28 with the full width at half maximum (FWHM) of the envelope . . . . . . . . . . . . . . . 32 Dis ussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.8.1 Causal and non ausal time-domain Green's fun tions for a ousti prop- 2.8.2 Convolving time-domain Green's fun tions with 3 y le . . . . . . . . . . . . . . . . . . . . . . . . 35 2.8.3 Causality in a ousti wave propagation . . . . . . . . . . . . . . . . . 35 Con lusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 agation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Hanning-weighted pulses 2.9 24 Chara terizing the non ausal omponent of the time-domain Green's fun tions for the Bla ksto k, Szabo, and power law wave equations 2.8 23 vi 33 Chapter 3 3.1 Spa e fra tional wave equations . . . . . . . . . . . . . . . . . . . . 38 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 3.2 The Chen-Holm spa e-fra tional wave equation . . . . . . . . . . . . . . . . 40 3.3 The Treeby-Cox spa e-fra tional wave equation . . . . . . . . . . . . . . . . 42 3.4 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.5 3.4.1 Dispersion Relations . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4.2 The Pantis Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4.3 Time windows for omputed time-domain Green's fun tions . . . . . 46 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.5.1 47 Phase velo ity and attenuation in breast and liver . . . . . . . . . . 3.5.2 Time-domain Green's fun tions al ulated for breast 3.5.3 Time-domain Green's fun tions al ulated for liver . . . . . . . . . 50 . . . . . . . . . . 52 3.5.4 Amplitude and full width at half maximum (FWHM) values in breast and liver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.5 Time-domain Green's fun tions onvolved with a three y le Hanningweighted pulse 3.6 3.7 54 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Dis ussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.6.1 59 Numeri al evaluations of the inverse 3D Fourier transform . . . . . . 3.6.2 Improved approximations for the attenuation and phase velo ity . . 61 3.6.3 Time-domain Green's fun tions . . . . . . . . . . . . . . . . . . . . . 63 3.6.4 Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.6.5 Convolving time-domain Green's fun tions with short pulses 65 3.6.6 Comparisons with time-domain Green's fun tions al ulated with 3D . . . . . FFTs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 Con lusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Chapter 4 Perfe tly mat hed layers for nonlinear ultrasound simulations with the KZK equation . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.1 73 4.3 4.4 4.5 4.6 The 2D wave equation with power law attenuation . . . . . . . . . . 4.2.2 Coordinate stret hing . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.3 The KZK equation 75 4.2.4 PML derivation for the KZK equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.3.1 80 Error al ulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Finite dieren e al ulations with the KZK equation 4.3.3 Muir's method Results . . . . . . . . . 80 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.4.1 KZK simulations for a linear lossless medium . . . . . . . . . . . . . 84 4.4.2 KZK simulations for a nonlinear lossy medium . . . . . . . . . . . . . 88 Dis ussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.5.1 Computation time . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.5.2 Continuous wave (CW) KZK al ulations . . . . . . . . . . . . . . . . 92 Con lusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 vii Chapter 5 Con lusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 APPENDIX A Derivation of the nonlinear wave equations . . . . . . . . . . . . 103 APPENDIX B Simulations of ultrasound wave propagation . . . . . . . . . . . 107 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 viii LIST OF FIGURES Figure 2.1: (a-d) Simulated time-domain Green's fun tions of the Caputo, Szabo, and power law wave equations al ulated for water with y = 2, α0 = 2.5328 × 10−4 Np/ m/MHz2 , and c0 = 1500 m/s and s aled by 4πr at (a) r=1 nm, (b) r = 100 nm, ( ) r=1 m, and (d) r = 10 m. (e-h) Simulated time-domain Green's fun tions of the Caputo, Szabo, and power law wave equations al ulated for breast with y = α0 = 0.086 Np/ m/MHz1.5 , and c0 = 1450 m/s and s aled by at (e) r = 10 nm, (f ) r = 1 µm, (g) r =1 m, and (h) r = 10 1.5, 4πr m. (i-l) Simulated time-domain Green's fun tions of the Caputo, Szabo, and power law wave equations al ulated for liver with y = 1.139, α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s and s aled by 4πr at (i) r = 100 zm, (j) r = 100 am, (k) r = 1 m, and (l) r = 10 m. Figure 2.2: 17 a) The per ent dieren e between the time-domain Green's fun tions for the Bla ksto k and Stokes wave equations and the timedomain Green's fun tion for the power law wave equation as a fun tion of distan e al ulated for water with y = 2, α0 = 2.5328 × 10−4 Np/ m/MHz2 , with c0 = 1500m/s. b- ) The per ent dieren e between the time-domain Green's fun tions for the Szabo and Caputo wave equation and the time-domain Green's fun tions for the power law wave equation as a fun tion of distan e al ulated for breast 1.5 with y = 1.5, α0 = 0.086 Np/ m/MHz , and c0 = 1450 m/s and 1.139 al ulated for liver with y = 1.139, α0 = 0.0459 Np/ m/MHz , and Figure 2.3: c0 = 1540 m/s. . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Relative non ausal ontributions of the Bla ksto k, Szabo, and power law wave equations hara terized with 20 log10 {g(r, t = 0)/ max [g(r, t)]} r . This quantity is al ulated for a) water y = 1.139, α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s 1.5 from r = 1 nm to r = 1 m, b) breast with y = 1.5, α0 = 0.086 Np/ m/MHz , and c0 = 1450 m/s from r = 1 nm to r = 1 m, and ) liver with y = 1.139, α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s from r = 1 zm to r = 1 pm. . . . . . . . . . . . . . . . . . . . . . . . . . 27 as a fun tion of distan e with Figure 2.4: Simulated three- y le Hanning-weighted pulse with a enter frequen y f0 = 7.5 MHz onvolved with time-domain Green's fun tions multiplied by 4πr al ulated for water at (a) r = 100 µm, (b) r = 1 mm, ( ) r = 1 m, and (d) r = 10 m, al ulated for breast at (e) r = 100 µm, (f ) r = 1 mm, (g) r = 1 m, and (h) r = 10 m, and al ulated for liver at (i) r= 100 µm, (j) r = 1 mm, (k) r = 1 m, and (l) r = 10 m. of ix 29 Figure 2.5: The FWHM of the envelope of b) breast, and ) liver, where pulse and g(r, t) v(t)∗4πrg(r, t) al ulated for a) water, v(t) is a three y le Hanning-weighted is the time-domain Green's fun tion for the Stokes, Bla ksto k, Caputo, Szabo, or power law wave equation. . . . . . . . Figure 3.1: Phase velo ity and attenuation in breast and liver obtained from the dispersion relation in Eq. 3.6 for the Chen-Holm wave equation and two dierent approximations (• and ◦) (∗) to the dispersion relation for the Chen-Holm wave equation. . . . . . . . . . . . . . . . . . . . Figure 3.2: 32 47 Phase velo ity and attenuation in breast and liver obtained from the dispersion relation in Eq. 3.12 for the Treeby-Cox wave equation (∗), the approximations to the dispersion relation for the Treeby-Cox wave equation given in Eqs. 3.14-3.15 (•), and the attenuation and phase velo ity for the power law wave equation given in Eqs. 2.1-2.2 Figure 3.3: (◦). . 48 Time-domain Green's fun tions s aled by 4πr al ulated for breast 1.5 with y = 1.5, α0 = 0.086 Np/ m/MHz , and c0 = 1450 m/s at (a) (e) r = 1 nm, r = 100 µm, r = 10 nm, ( ) r = 100 nm, (d) r = 1 µm, r = 1 mm, (g) r = 1 m, and (h) r = 10 m with (b) (f ) the power law (solid line), Chen-Holm (dashed line), and Treeby-Cox (dot-dashed line) wave equations. At all distan es, the time-domain Green's fun tions for the Chen-Holm and Treeby-Cox wave equations evaluated for breast are ausal while the time-domain Green's fun tion for the power law wave equation is learly non- ausal for and r = 10 nm. Beyond about r = 100 µm, r = 1 nm the time-domain Green's fun tions for the power law wave equation and the Treeby-Cox wave equation are nearly indistinguishable, while the time-domain Green's fun tion for the Chen-Holm wave equation is distin t from the timedomain Green's fun tions for the time-fra tional power law wave equation and the spa e-fra tional Treeby-Cox wave equation at all distan es. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 50 Figure 3.4: Time-domain Green's fun tions s aled by 4πr al ulated for liver with y = 1.139, α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s at (a) (e) r = 100 zm, (b) r = 1 am, ( ) r = 10 am, (d) r = 100 am, r = 100 µm, (f ) r = 1 mm, (g) r = 1 m, and (h) r = 10 m with the power law (solid line), Chen-Holm (dashed line), and TreebyCox (dot-dashed line) wave equations. At all distan es, the time- domain Green's fun tions for the Chen-Holm and Treeby-Cox wave equations evaluated for liver are ausal while the time-domain Green's fun tion for the power law wave equation is learly non- ausal for r = 100 z m. Beyond about r = 100 µm, the time-domain Green's fun tions for the power law wave equation and the Treeby-Cox wave equation are nearly indistinguishable, while the time-domain Green's fun tion for the Chen-Holm wave equation is onsistently distin t from the time-domain Green's fun tions for the time-fra tional power law wave equation and the spa e-fra tional Treeby-Cox wave equation at all distan es. Figure 3.5: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 The (a, ) amplitudes and (b,d) FWHM values of the time-domain Green's fun tions al ulated for the power law, Chen-Holm, and TreebyCox wave equations in (a,b) breast and ( ,d) liver. The amplitudes of all three time-domain Green's fun tions de rease as the distan e in reases while the FWHM values of all three time-domain Green's fun tions in rease as the distan e in reases. The amplitudes of the time-domain Green's fun tions for all three of these fra tional wave equations are very similar at ea h distan e, and the FWHM values are all approximately the same at longer distan es, although there is a small dieren e in the FWHM values at shorter distan es that diminishes with in reasing distan e. Figure 3.6: . . . . . . . . . . . . . . . . . . 54 Simulated three- y le Hanning-weighted pulses with enter frequen ies (a) f0 = 7.5 MHz and (b) f0 = 29 MHz onvolved with time- domain Green's fun tions for the power law, Chen-Holm, and TreebyCox wave equations multiplied by 1 4πr evaluated for breast at r = m. The onvolution results for the power law wave equation and the Treeby-Cox wave equation are very similar while the onvolution result for the Chen-Holm wave equation learly shows a time delay. Signi ant attenuation and waveform spreading are observed for all three signals in (b) produ ed by inputs with f0 = 29 MHz, whereas a moderate amount of attenuation and waveform spreading is observed for all three signals in (a) produ ed by inputs with xi f0 = 7.5 MHz. . 56 Figure 3.7: Simulated three- y le Hanning-weighted pulses with enter frequen ies (a) f0 = 7.5 MHz and (b) f0 = 29 MHz onvolved with time- domain Green's fun tions for the power law, Chen-Holm, and TreebyCox wave equations multiplied by 4πr evaluated for liver at r = 1 m. The onvolution results for the power law wave equation and the Treeby-Cox wave equation are very similar while the onvolution result for the Chen-Holm wave equation learly shows a time delay. A moderate amount of attenuation and waveform spreading are observed for all three signals in (b) produ ed by inputs with f0 = 29 MHz, whereas a smaller amount of attenuation and waveform spreading is observed for all three signals in (a) produ ed by inputs with Figure 3.8: f0 = 7.5 MHz. . . . . . . . . . . . . . . . . . . . . . . . . . . . Computed time-domain Green's fun tion s aled by 4πr 58 evaluated for r = 1 m omputed with the Pantis method using (a) 2,000 Filon abs issas and m = 401, (b) 2,000 Filon abs issas and m = 101, and ( ) 500 Filon abs issas and m = 401. . . . . . . . . . . . . . . . breast at Figure 3.9: 4πr evaluated for m = 5 (a) without Computed time-domain Green's fun tions s aled by breast at r=1 nm using 500 Filon abs issas and and (b) with the Pantis term. Figure 3.10: 60 . . . . . . . . . . . . . . . . . . . . . 60 4πr al ulated for breast y = 1.5, α0 = 0.086 Np/ m/MHz1.5 , and c0 = 1450 m/s evaluated at (a) r = 10 nm, (b) r = 100 nm, and ( ) r = 1 µm with Time-domain Green's fun tions s aled by with the Pantis method (solid line) and the 3D FFT approa h (dot-dashed line) using Figure 3.11: dx = dy = dz = 0.5 nm in (a), (b), and ( ). . . . . . . . 67 4πr al ulated for liver 1.139 Np/ m/MHz , and c0 = 1540 m/s Time-domain Green's fun tions s aled by with y = 1.139, α0 = 0.0459 r = 100 zm, (b) r = 1 evaluated at (a) am, and ( ) r = 10 am with the Pantis method (solid line) and 3D FFT approa h (dot-dashed line) using Figure 3.12: dx = dy = dz = 50 zm in (a), (b), and ( ). . . . . . . . . 67 4πr al ulated for liver y = 1.139, α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s evaluated at r = 10 m with the Pantis method (solid line) and with 3D FFTs (dot-dashed line) using (a) dx = r/100, (b) dx = r/200, ( ) dx = r/300, (d) dx = r/400, and (e) dx = r/500. In all simulations with 3D FFTs evaluated here, dx = dy = dz . . . . . . . 69 Time-domain Green's fun tions s aled by with Figure 4.1: Comparison of simulated on-axis waveforms obtained from nite differen e KZK al ulations (solid line) and Muir's method (dashed line) evaluated in a linear lossless medium at xii z=6 m. . . . . . . . . . . 85 Figure 4.2: Comparison between on-axis waveforms generated by nite dieren e KZK al ulations in a linear lossless medium at z=6 m with and without a PML using dierent radial boundaries. (a) KZK simulation without a PML that denes a radial boundary at solid line) and at rmax = 6a rmax = 2a (bla k (red dashed line). (b) KZK simulation rmax = 6a (red y = 0 single term PML that denes a radial boundary at rmax = 2a (blue solid line), and with a y = 2 single term PML that denes a radial boundary at rmax = 2a (green solid line). without a PML that denes a radial boundary at dashed line), with a Figure 4.3: 86 Simulated 2D pressure eld and dieren es between KZK al ulations without and with PMLs, where the radial boundaries are lo ated at rmax = 9 m and at rmax = 3 m in a linear lossless medium. (a) The peak pressure distribution for the KZK simulation without a PML that denes a radial boundary at rmax = 9 m. (b) The dieren e between the KZK simulation without a PML that denes a radial boundary at rmax = 9 m and the KZK simulation without a PML that denes a radial boundary at rmax = 3 m. ( ) The dieren e between the KZK simulation without a PML that denes a radial boundary at rmax = 9 y = 0 = 3 m. m and the KZK simulation with a single term PML that denes a radial boundary at rmax (d) The dieren e between the KZK simulation without a PML that rmax = 9 denes a radial boundary at y = 2 = 3 m. with a rmax Figure 4.4: m and the KZK simulation single term PML that denes a radial boundary at . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Comparison between on-axis waveforms generated by nite dieren e KZK al ulations in a nonlinear medium at z = 6 m with and without a PML for dierent radial boundaries. The attenuation 2 −3 parameter is α = 2.2 × 10 dB/ m/MHz , and the nonlinearity parameter is β = 3.5. (a) KZK simulation without a PML that denes a radial boundary at rmax = 2a rmax = 6a (b) KZK simulation without a PML (red dashed line). that denes a radial boundary at y=0 (bla k solid line) and at rmax = 6a (red dashed line), with a single term PML that denes a radial boundary at (blue solid line), and with a radial boundary at rmax = 2a y =2 single term PML that denes a (green solid line). xiii rmax = 2a . . . . . . . . . . . 89 Figure 4.5: Simulated 2D pressure eld and dieren es between KZK al ulations without and with PMLs, where the radial boundaries are lo ated at rmax = 9 rmax = 3 m in a nonlinear medium. The α = 2.2 × 10−3 dB/ m/MHz2 and the nonlinearity parameter is β = 3.5. (a) The peak pressure distribution m and at attenuation parameter is for the KZK simulation without a PML that denes a radial boundary at rmax = 9 m. (b) The dieren e between the KZK simulation rmax = 9 without a PML that denes a radial boundary at m and the KZK simulation without a PML that denes a radial boundary at rmax = 3 m. ( ) The dieren e between the KZK simulation rmax = 9 without a PML that denes a radial boundary at y = 0 = 3 m. the KZK simulation with a radial boundary at rmax m and single term PML that denes a (d) The dieren e between the KZK simulation without a PML that denes a radial boundary at rmax = 9 m and the KZK simulation with a that denes a radial boundary at Figure 4.6: rmax = 3 y=2 single term PML m. . . . . . . . . . . . . 91 On-axis omparisons between ontinuous wave nite dieren e KZK simulations with the results al ulated with Muir's method evaluated in a linear lossless medium. (a) The result obtained with Muir's method (red solid line) and the nite dieren e KZK simulation results without a PML that denes a radial boundary at 10.5 rmax = m (blue dashed line). (b) The result with Muir's method (red solid line) and the nite dieren e KZK simulation results with a y = 0 PML that denes a radial boundary at rmax = 3 m (blue dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 4.7: The ontinuous-wave 2D pressure distribution for a spheri ally-fo used transdu er with a = 1.5 m, R=6 m, and f =1 a linear lossless medium with Muir's method. Figure 4.8: 93 MHz al ulated in . . . . . . . . . . . . 95 The 2D pressure dieren e between linear lossless nite dieren e KZK numeri al results for a spheri ally-fo used transdu er with 1.5 m, R=6 m, and f =1 a= MHz and the results for the same on- guration evaluated with Muir's method. (a) The dieren e between the nite dieren e KZK simulation without a PML that denes a radial boundary at rmax = 10.5 m and the results obtained with Muir's method. (b) The dieren e between the nite dieren e KZK simulation with a rmax = 3 y = 0 PML that denes a radial boundary at m and the results obtained with Muir's method. xiv . . . . . 96 Figure 4.9: The rst four harmoni s generated by a spheri ally-fo used transdu er with a = 1.5 m, R=6 m, and f =1 MHz for on-axis nite dieren e simulations of the ontinuous wave KZK equation in water. 2 −3 The attenuation parameter is α = 2.2 × 10 dB/ m/MHz , and the nonlinearity parameter is β = 3.5. (a) The nite dieren e KZK simulation results without a PML that denes a radial boundary at rmax = 10.5 with a y = 0 Figure 4.10: m. (b) The nite dieren e KZK simulation results PML that denes a radial boundary at rmax = 3 m. . 97 The 2D pressure dieren e between the nite dieren e KZK numeri al results with and without a y = 0 PML in water with α0 = 2.2 × 10−3 dB/ m/MHz2 and β = 3.5 evaluated for rst four harmoni s produ ed by a spheri ally-fo used transdu er with 6 m, and f =1 MHz. a = 1.5 m, R= . . . . . . . . . . . . . . . . . . . . . . . . . xv 98 Chapter 1 Introdu tion 1.1 Ultrasound attenuation in soft tissue As sound waves propagate, the medium is temporarily displa ed in a dire tion parallel (longitudinal wave) or perpendi ular (transverse wave) to the dire tion of energy transport and then the medium returns to the equilibrium state. When ultrasound travels through a medium, the intensity diminishes with distan e. In a lossless medium, the amplitude is only redu ed by the spreading of the wave. However, when ultrasound propagates through soft tissue, the amplitude is redu ed as a fun tion of propagation distan e, and the enter frequen y of the signal is also downshifted by attenuation. As indi ated by Laugier and Haïat [1℄, Goss et al. [2℄, and Parker [3℄, the two main me hanisms that ontribute to ultrasound attenuation are absorption and s attering. Absorption is the onversion of the sound energy to other forms of energy [4, 5℄, espe ially heat as a result of fri tion between the vibrating parti les that transmit the a ousti wave within soft tissue. In homogeneous vis ous media, the vis ous for es between neighboring parti les moving with dierent velo ities are the major sour e of a ousti wave absorption. A ousti wave attenuation is also aused by s attering, whi h des ribes the redire tion of the in ident wave in multiple dire tions [6, 7℄. In heterogeneous media, where the physi al properties su h as density or sound speed are dierent from those of the surrounding medium, s attering also redire ts the a ousti energy. 1 In three-dimensional (3D) spa e, the amplitude de ay and attenuation of ultrasound are mathemati ally des ribed by where p is the pressure, α0 p(r) = p0 e−α(f )r for the attenuation given by r is the attenuation oe ient, r α(f ) = α0 |f |y , is the distan e in 3D, f frequen y, and y y However, in most biologi al tissues, the measured power law exponents is equal to 2. is the power law exponent. For instan e, in water, the power law exponent within the range of 0.7 ≤ y ≤ 1.5 y = 1.5 y are for the range of frequen ies utilized in medi al ultrasound [8℄. For example, measured values for the power law exponent are [9℄ and is the y = 1.139 in human liver in human breast [10℄. The power law exponents and attenuation oe ients vary for dierent tissues, and measurements of these parameters have been widely evaluated for medi al ultrasound in human tissue [11, 12, 13℄. 1.2 Nonlinear ultrasound The fundamental equations of nonlinear ultrasound are derived from the three onstitutive relations, namely the equation of motion, the ontinuity equation, and the equation of state [14℄. For small pressure amplitudes, the linearized versions of these three fundamental equations are ombined to produ e a linear wave equation. However, when the pressure amplitudes are su iently large, the se ond order terms in these fundamental equations must be retained, and the ombination of the three onstitutive relations yields a nonlinear wave equation. The amount of nonlinearity in a material through whi h a nite-amplitude ultrasoni wave propagates is expressed by the nonlinearity parameter and B B/A. The values of A are the oe ients of the rst and se ond order terms of the Taylor series expansion for the equation of state, whi h relates the pressure to the density. Some values of B/A in biologi al tissues are given by Wells [15℄. Some ommon models that des ribe nonlinear ultrasound propagation in lude the Westervelt equation, the Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation, and Burgers equation. A general wave equation that a ounts for nonlinearity up to se ond-order is given by the Westervelt equation [16℄. After a paraboli approximation is applied, the Westervelt equation 2 redu es to the KZK equation, whi h a ounts for the ombined ee ts of nonlinearity, dira tion, and absorption in dire tional sound beams. ommonly used to model problems in nonlinear a ousti s. Solutions to this equation are Several numeri al approa hes that solve the KZK equation have been proposed by Lee and Hamilton [17℄, Cleveland [18℄, and Berntsen [19℄. When the dira tion term is dis arded, the KZK equation redu es to Burgers equation, whi h des ribes the ombined ee ts of nonlinearity and attenuation on the propagation of progressive plane waves. Solutions to Burgers equation an be obtained with several dierent methods [20, 21℄. Nonlinear wave propagation has been widely analyzed in the medi al ultrasound eld. Two ommon appli ations in lude high intensity fo used ultrasound (HIFU) in therapeuti ultrasound [22, 23, 24℄ and harmoni imaging in diagnosti ultrasound [25, 26, 27℄. HIFU generates high intensity pressure elds in the fo al zone to heat tumors or break up kidney stones [28℄. Compared to diagnosti ultrasound, HIFU uses higher energies and lower frequen ies. In harmoni imaging, sin e s attering and rst ree tions are redu ed in the se ond harmoni , the resulting images provide better ontrast, better resolution, and diminished ee ts of undesirable sidelobes. 1.3 Fra tional derivative operators Fra tional derivative operators are applied widely in the eld of fra tional al ulus. A variety of fra tional derivatives are dened to repla e the integer order derivative, in luding the RiemannLiouville fra tional derivative [29℄, the Caputo fra tional derivative [30℄, the AtanganaBaleanu derivative [31℄, the Katugampola fra tional derivative [32℄, and so on. Two of these fra tional derivative operators are utilized here, namely the RiemannLiouville fra tional derivative and the Caputo fra tional derivative [33℄. 3 The RiemannLiouville fra tional derivative is dened by DLy f where (t) = y > 0, t > a,     dn 1 Γ(n−y) dtn f (τ ) dτ, a (t−τ )y+1−n dn f dtn    and ´t y, a, t ⊂ R. n − 1 < y < n ⊂ N, (t) , (1.1) y = n ⊂ N, The gamma fun tion, whi h frequently appears in fra tional al ulus, is dened as Γ (z) = ˆ ∞ tz−1 e−t dt. (1.2) 0 The Caputo fra tional derivative takes the following form: DCy f (t) =        1 Γ(n−y) f (n) (τ ) dτ, a (t−τ )y+1−n ´t dn f dtn (t) , n − 1 < y < n ⊂ N, (1.3) y = n ⊂ N. These two fra tional derivatives are similar, where the main dieren e is the order in whi h the dierentiation and integration operations are performed. For a fra tional derivative operator, the derivative of a fun tion evaluated at a point is no longer a lo al property, so additional knowledge of previous states is required in either time or spa e. The properties of integer derivatives for Fourier and Lapla e transforms are readily extended to fra tional derivatives [33℄. The 1D Lapla e transform for the integer derivative is dened as L {f (t)} = ∞ ˆ e−st f (t) dt. (1.4) 0 This is extended to the fra tional Riemann Liouville derivative as [29℄ L {DLy f y (t)} = s F (s) − 4 n−1 X k=0 h i sk DLy−k−1 f (t) t=0 (1.5) and to the Caputo derivative as [30℄ L {DCy f y (t)} = s F (s) − n−1 X sy−k−1 f (k) (0) . (1.6) k=0 For the Fourier transform in 1D F {f (t)} = ˆ ∞ f (t) e−jωt dt, (1.7) −∞ both RiemannLiouville and Caputo derivatives are the same [34, 35℄: F {DLy f (t)} = (jω)y F (ω) , (1.8) F {DCy f (t)} = (jω)y F (ω) . (1.9) 1.4 Dispersion relations for fra tional wave equations When the ee ts of attenuation are in luded through the term εLr,t (p), the wave equation be omes ∇2 p − where the fra tional operator εLr,t (p) 1 ∂2p − εLr,t (p) = 0, c20 ∂t2 (1.10) may be either time-fra tional or spa e-fra tional or both. After Fourier transforms are evaluated in both time and spa e domains, the dispersion relation for Eq. 1.10 is given by ω2 −k + 2 − εL̃ (k, ω) = 0. c0 2 5 (1.11) When Lr,t (p) ontains only time-fra tional derivatives, εL̃ (k, ω) = εL̃ (ω), and the analyti al expression for the wavenumber is then represented by k (ω) = s ω2 − εL̃ (ω). c20 (1.12) The relationship of phase velo ity and attenuation is then obtained from the real and imaginary parts of Eq. 1.12. For the general ase when the loss term εL̃ (k, ω) ontains at least one spa e-fra tional operator, the wavenumber in Eq. 1.11 is al ulated with the binomial approximation   ω c20 c40 2 2 c60 3 3 k (ω) ≈ 1 − 2 εL̃ (k, ω) − 4 ε L̃ (k, ω) − ε L̃ (k, ω) + . . . . c0 2ω 8ω 16ω 6 The right-hand side of Eq. 1.13 in ludes several terms that ontains fun tion of k. To obtain an expression for the wavenumber L̃ (k, ω), (1.13) whi h is a k (ω) that is independent of k on the right-hand side, further approximations are required. 1.4.1 First order approximation If O (ε2 ) and higher order terms in Eq. 1.13 are dis arded, and the rst-order approxi- mation for Eq. 1.13 is then given by k (ω) ≈ where k+ is obtained by setting  ω c0 − εL̃ k + , ω , c0 2ω L̃ (k, ω) = 0, substituted ba k to Eq. 1.14 and the O (ε) whi h yields (1.14) k+ ≈ ω . c0 This expression is terms are ignored, whi h yields c0 ω − εL̃ k (ω) ≈ c0 2ω 6   ω ,ω . c0 (1.15) 1.4.2 Se ond order approximation To obtain a more a urate approximation for the wavenumber, third order and higher terms are dis arded from Eq. 1.13. Then, the se ond-order approximation is given by  c0 c30 2 2 +  ω + − ε L̃ k , ω − ε L̃ k , ω . k (ω) ≈ c0 2ω 8ω 3 Similarly, k+ is approximated by substituting ω/c0 ω c0 k ≈ − εL̃ c0 2ω + and then terms that are third order or higher in into Eq. 1.16  ε (1.16)  ω ,ω , c0 (1.17) are dis arded. 1.5 Thesis stru ture For ertain time-fra tional and spa e-fra tional models, exa t and approximate timedomain Green's fun tions have been derived and evaluated numeri ally. More a urate expressions for the phase velo ity and attenuation are also derived for several fra tional al ulus models. To demonstrate some of the similarities and dieren es in these fra tional partial dierential equations, ausality is analyzed for ea h of these, the time-domain Green's fun tions are ompared, and full width at half maximum (FWHM) values for ea h timedomain Green's fun tion are evaluated for breast and liver models. Chapter 2 numeri ally evaluates time-domain Green's fun tions for three time-fra tional models, namely the power law wave equation, the Szabo wave equation, and the Caputo wave equation. These Green's fun tions are evaluated for water with a power law exponent of y = 2, breast with a power law exponent of y = 1.5, and liver with a power law exponent of y = 1.139. The ausality of ea h fra tional wave equation is analyzed, and the time-domain Green's fun tions for these three time-fra tional models are ompared at dierent distan es. To demonstrate the ee ts of power law attenuation and dispersion on transient ex itations, a 7 three- y le Hanning-weighted pulse is also onvolved with the time-domain Green's fun tions for these three time domain Green's fun tions. Chapter 3 evaluates improved approximations for the frequen y-dependent phase velo ity and attenuation that were derived from two spa e-fra tional models, namely the Chen-Holm and Treeby-Cox spa e-fra tional wave equations, and these are evaluated using parameters for breast and liver. After the ausality of the two spa e-fra tional models is established, the amplitudes and FWHM values of the time-domain Green's fun tions are evaluated at short distan es from the origin. In addition, a three- y le Hanning weighted pulse is onvolved with ea h time-domain Green's fun tion to show how dieren es in these Green's fun tions inuen e the results for a nite bandwidth ex itation. Chapter 4 introdu es new expressions that des ribe perfe tly mat hed layers (PML) for numeri al simulations with the transient KZK equation. Arti ial attenuation in these new PMLs is implemented through terms derived from the power law wave equation with and y = 2. y=0 These expressions are further simplied by retaining only one term, whi h is su ient to redu e ree tions from the radial boundary. For a spheri ally fo used transdu er with aperture radius a = 1.5 m and radius of urvature R = 6 m, simulations in both linear lossless and nonlinear media validate the ee tiveness of these new PMLs. Similar simulations are then evaluated for the ontinuous wave KZK equation in both linear lossless and nonlinear media. 8 Chapter 2 Time fra tional wave equations 1 2.1 Introdu tion The attenuation of ompressional ultrasound waves in soft tissue is des ribed by a power law of the form α(f ) = α0 |f |y , onstant in Np/m/Hz y where f is the frequen y in MHz, y or dB/m/Hz , and y y = 1.5 is the attenuation is the power law exponent. measured values for the power law exponent are liver [9℄, and α0 y = 2 in water, Examples of y = 1.139 in human in human breast [10℄. Additional values for mammalian tissues with various power law exponents are tabulated in the book by Du k [8℄, and other attenuation values are ompiled in papers by Goss et al. [36, 37℄. The orresponding wave equations that des ribe power law attenuation in soft tissue utilize fra tional derivatives, whi h are non-integer order derivatives. These fra tional derivatives are often time-fra tional [38, 39, 40℄, although spa e-fra tional derivatives are also used [41, 42℄. Examples of time-fra tional wave equations that model the attenuation and dispersion of ultrasound in soft tissue in lude the Szabo wave equation [38℄, the Caputo wave equation [30℄, and the power law wave equation [39℄. The Szabo and power law wave equations were developed for medi al ultrasound appli ations, and the Caputo wave equation [30℄ was originally dened for appli ations in geophysi s and then independently onsidered by Wismer as a model for attenuation and dispersion in soft tissue [40℄. 1 Reprodu ed from X. Zhao and R. J. M Gough, Time-domain omparisons of power law attenuation in ausal and non ausal time-fra tional wave equations, The Journal of the A ousti al So iety of Ameri a, 139(5):30213031, 2016, with the permission of the A ousti al So iety of Ameri a. 9 These time-fra tional wave equations are parti ularly amenable to analyti al methods for analyzing ausality, in luding the Paley-Wiener riterion [43℄, Kramers-Kronig analysis [44℄, and a time ausal theory [38℄; however, in onsistent on lusions are often rea hed with dierent methods, espe ially for power law exponents y ≥ 1. In an eort to resolve some of these apparent in onsisten ies, time-domain Green's fun tions are al ulated numeri ally for the Bla ksto k wave equation, the Stokes wave equation, the Szabo wave equation, the Caputo wave equation, and the power law wave equation. In addition, a three- y le Hanningweighted pulse is onvolved with ea h of these to show the ee ts of ausal and non ausal Green's fun tions on the al ulated signals. The results show that non ausal behavior is only evident very lose to the sour e in time-domain Green's fun tion al ulations, that this non ausal behavior is no longer evident after onvolution with a short pulse, and that time-domain al ulations with these ausal and non ausal time-fra tional models onverge a short distan e from the sour e. 2.2 Power law attenuation and dispersion The frequen y-dependent attenuation α(ω) of ultrasound in soft tissue is des ribed by the power law [45℄ α (ω) = α0 |ω|y , where y is the power-law exponent, α0 (2.1) is the attenuation onstant, and ω is the angular fre- quen y in radians/se ond. The orresponding frequen y-dependent sound speed (dispersion) c(ω) satises [45℄  πy  1 1 + α0 tan |ω|y−1 . = c (ω) c0 2 In Eq. 2.2, c0 ω = ∞ 0 ≤ y < 1. for is the sound speed at When ω = 0 y = 2, for 1 < y ≤ 2, and (2.2) c0 is the sound speed at Eq. 2.2 is nondispersive be ause the ω dependen e in Eq. 2.2 disappears. For the numeri al al ulations that follow, the attenuation onstant α0 with units Np/ m/MHz y is multiplied by 100 and divided by 10 106y and (2π)y to onvert m into m, MHz into Hz, and frequen y in Hz into angular frequen y in radians/se ond, respe tively. 2.3 The Szabo wave equation α0 For an attenuation onstant y with units Np/m/Hz , the Szabo wave equation [38℄ is given by ∇2 p − where p represents 2α0 ∂ y+1 p 1 ∂2p − = 0, c20 ∂t2 c0 cos (πy/2) ∂ty+1 (2.3) t is the time in se onds. The Szabo wave equation the pressure in Pa and is a time-fra tional extension of the Bla ksto k wave equation [46℄, where the third term approximately des ribes the ee ts of power law attenuation and dispersion in Eqs. 2.1 and 2.2 over the range of frequen ies where the smallness approximation [38℄ holds. There is no known exa t time-domain Green's fun tion for the Szabo wave equation, but the 3D frequen y-domain Green's fun tion for the Szabo wave equation is − cr G(r, ω) = for frequen ies ω ≥ 0, r= where origin to an observation point at p e 0 q 2α c 0 0 (jω)y+1 −ω 2 + cos(πy/2) 4πr x2 + y 2 + z 2 (2.4) is the distan e from a point sour e at the (x, y, z). The phase velo ity and attenuation are derived by solving the dispersion relation 2α0 ω2 (−jω)y+1 k = 2 − c0 c0 cos (πy/2) 2 for k. (2.5) By taking the square root of Eq. 2.5 and utilizing the binomial approximation, an approximate expression for the wave number is obtained k≈ + α20 c0 2 ω c0 {1 + jα0 c0 [1 − j tan (πy/2)] ω y−1 } [1 − j2 tan (πy/2) − tan2 (πy/2)] ω 2y−1 . 11 (2.6) The approximate phase velo ity is then extra ted from the real part of the wavenumber divided by ω  1 1 1 1 − tan2 (πy/2) α02 c0 ω 2y−2 ≈ + tan (πy/2) α0 ω y−1 + c (ω) c0 2 (2.7) and the approximate attenuation is the imaginary part of the wavenumber α (ω) ≈ α0 ω y − tan (πy/2) α02 c0 ω 2y−1 . When the power law exponent y (2.8) is equal to 2, the Szabo wave equation redu es to the Bla ksto k equation ∇2 p − 1 ∂ 2 p 2α0 ∂ 3 p + = 0, c20 ∂t2 c0 ∂t3 (2.9) where the 3D frequen y-domain Green's fun tion for the Bla ksto k equation is given by − cr G(r, ω) = e 0 √ −ω 2 −2α0 c0 (jω)3 . 4πr (2.10) 2.4 The power law wave equation The power-law wave equation [39℄, whi h is losely related to the Szabo wave equation, is given by ∇2 p − 2α0 α02 ∂ y+1 p ∂ 2y p 1 ∂2p − − = 0. c20 ∂t2 c0 cos (πy/2) ∂ty+1 cos2 (πy/2) ∂t2y (2.11) The rst three terms in the power law wave equation also appear in the Szabo wave equation, where the fourth time-fra tional term yields a omplex wavenumber that exa tly satises Eqs. 2.1 and 2.2 for all frequen ies ω. The 3D frequen y-domain Green's fun tion for the power law wave equation is y e−jωr/c0 e−α0 (jω) G(r, ω) = 4πr 12 r/ cos(πy/2) . (2.12) By expanding the argument of the se ond exponential fun tion in Eq. 2.12 after applying Euler's formula with j y = ejπy/2 and olle ting real and imaginary terms, the attenuation and dispersion relations in Eqs. 2.1 and 2.2 are exa tly re overed. Furthermore, unlike the other time-fra tional wave equations evaluated here, the power law wave equation has an exa t losed form 3D time-domain Green's fun tion, whi h is # " t − cr0 1 1 . f˜y g (r, t) = 4πr (α0 r)1/y (α0 r)1/y In Eq. 2.13, f˜y (2.13) is the probability density fun tion (pdf ) for a maximally skewed stable distribution [47℄ with parameter y. Sin e the power law wave equation exa tly satises Eq. 2.1 and Eq. 2.2, demonstrating whether the time-domain Green's fun tion in Eq. 2.13 is ausal or non ausal is equivalent to demonstrating whether the ombination of Eqs. 2.1 and 2.2 is ausal or non ausal. For the power law exponent y = 2, the 3D frequen y-domain Green's fun tion in Eq. 2.12 redu es to 2 e−jωr/c0 e−α0 ω r , G(r, ω) = 4πr whi h is a Gaussian fun tion multiplied by the 1/ (4πr) (2.14) geometri spreading fa tor and a omplex exponential delay term. The inverse Fourier transform of Eq. 2.14 is exa tly equal to the time-shifted Gaussian fun tion g (r, t) = 2 1 1 √ e−(t−r/c0 ) /(4α0 r) , 4πr 4πα0 r (2.15) whi h is equivalent to the time-domain Green's fun tion in Eq. 2.13 with the power law exponent y = 2. Although the expressions in Eqs. 2.14 and 2.15 are only appli able to a few materials with frequen y-squared attenuation su h as water and air, these expressions are nevertheless onvenient for preliminary evaluations and omparisons. 13 2.5 The Caputo wave equation The Caputo wave equation [30℄, whi h is a time-fra tional extension of the Stokes wave equation [48℄, is given by ∇2 p − where τ y−1 1 ∂2p y−1 ∂ + τ ∇2 p = 0 2 2 y−1 c0 ∂t ∂t (2.16) is the fra tional relaxation time. The Caputo wave equation approximately satises the attenuation and dispersion relations in Eqs. 2.1 and 2.2, respe tively, and no exa t losed form time-domain Green's fun tion is available for the Caputo wave equation. However, there is an exa t 3D frequen y-domain Green's fun tion for the Caputo wave equation, whi h is − jωr c G(r, ω) = for frequen ies e 1 y−1 1 + (jωτ ) 0 √ 1 1+(jωτ )y−1 4πr . (2.17) ω ≥ 0. To obtain an expression that relates the value of the power law attenuation onstant the fra tional relaxation time α0 to τ , the power law wave equation and the Caputo wave equation are Fourier-transformed in time and spa e. After solving for the square of the wavenumber and taking the square root of both sides, the smallness approximation [38℄ is applied to the expression obtained from the Caputo wave equation. The resulting onversion fa tor [41℄ is τ y−1 = −2α0 c0 / cos (πy/2) . The expression in Eq. 2.18 is singular at are also singular at y = 1, y = 1, (2.18) the Szabo and power law wave equations and the Caputo wave equation is non-attenuating at only values of the power law exponent that satisfy numeri al evaluations. 14 1 τ c0 k = (2/τ c0 )1/(y−1) , ĝ(k, t) = tc20 e−τ k and t. y c2 t/2 0 u (t), In all three of these expressions for so ĝ(k, t) ĝ(k, t), u (t) guarantees that the time-domain response is equal to 0 (3.4) is ontinuous with respe t to is the unit step fun tion, whi h for all times t<0 and that ĝ(k, t) is ausal. The time-domain Green's fun tion is then obtained by evaluating 4π g(r, t) = (2π)3 r where r= point at p x2 + y 2 + z 2 (x, y, z). ˆ ∞ ĝ(k, t) sin (kr) kdk, (3.5) 0 is the distan e from a point sour e at the origin to an observation The integral in Eq. 3.5 also ontains the unit step fun tion 41 u (t) through the expressions for ĝ(k, t), so the time-domain Green's fun tion spa e-fra tional wave equation is also ausal. expressions, the time-domain Green's fun tion g(r, t) for the Chen-Holm Although Eqs. 3.3 and 3.4 are analyti al g(r, t) in Eq. 3.5 is evaluated numeri ally. Approximate expressions for the phase velo ity and attenuation are derived from the dispersion relation k2 = with τ = 2α0 cy−1 0 . ω2 + jωτ k y c20 (3.6) After taking the square root of both sides of Eq. 3.6, fa toring out ω/c0 on the right hand side, and keeping the rst three terms in the binomial approximation, the wavenumber is rewritten as side, ky is repla ed by the approximation the approximation to α0 k ≈ ω/c0 + jα0 cy0 k y + α02 c2y+1 ω −1 k 2y /2. 0 (ω/c0 )2y , (ω/c0)y (1 + jyα0 c0 ω y−1 ) and On the right hand k 2y is repla ed by respe tively, and all third order and higher terms with respe t are dis arded. The approximate phase velo ity is then obtained from the real part of the wavenumber divided by ω,   1 1 1 α02 c0 ω 2y−2 . ≈ − y− c (ω) c0 2 The attenuation α (ω) = α0 ω y , approximation for k, (3.7) whi h is obtained from the imaginary part of the resulting is the same expression given in Eq. 2.1. In Treeby and Cox [42℄, the attenuation for the Chen-Holm wave equation is also α (ω) = α0 ω y , but the approximate phase velo ity for the Chen-Holm wave equation is instead expressed as c (ω) ≈ c0 . 3.3 The Treeby-Cox spa e-fra tional wave equation Treeby and Cox [42℄ re ognized that the Chen-Holm spa e-fra tional wave equation orre tly models the power law attenuation des ribed by Eq. 2.1 but that the dispersion for the Chen-Holm wave equation diers from the desired expression given in Eq. 2.2. To address this problem, Treeby and Cox inserted a se ond spa e-fra tional term that a ounts 42 for the dispersion, and the resulting expression is given by [42℄     1 ∂2p ∂ 2 y/2 2 (y+1)/2 ∇ p − 2 2 + −τ p = 0, −∇ + η −∇ c0 ∂t ∂t 2 where τ = 2α0 cy−1 0 and η = 2α0 cy0 tan (πy/2). (3.8) In Eq. 3.8, the third and fourth terms a ount for the ee ts of attenuation and dispersion, respe tively. Similar to the Chen- Holm wave equation, there is no known exa t analyti al losed-form time-domain Green's fun tion for the Treeby-Cox wave equation, so numeri al evaluations are required. To derive an approximate expression for these numeri al al ulations, the transfer fun tion for an impulsive for ing fun tion applied to the Treeby-Cox wave equation is expressed in terms of the spatial Fourier transform variable k and the Lapla e transform variable c20 G(k, s) = s, whi h is . (3.9)   p sin kc0 t 1 − τ 2 k 2y−2 c20 /4 − ηk y−1 (3.10) 2 (s + τ k y c20 /2) + k 2 c20 − τ 2 k 2y c40 /4 − ηc20 k y+1 After the inverse Lapla e transform of Eq. 3.9 is evaluated, the result is ĝ(k, t) = e−τ k k √ y c2 t/2 0 c0 1−τ 2 k 2y−2 c20 /4−ηk u (t) y−1 for k<κ and ĝ(k, t) = e−τ k k where κ = h expression for η+ y c2 t/2 0   p sinh kc0 t τ 2 k 2y−2 c20 /4 + ηk y−1 − 1 c0 u (t) τ 2 k 2y−2 c20 /4+ηk y−1 −1 √  i−1/(y−1) p η 2 + τ 2 c20 /2 . ĝ(k, t) At ĝ(k, t) k and indi ates that (3.11) k > κ, k = κ, ĝ(k, t) = tc20 e−τ κ is ontinuous with respe t to in all three of these expressions for for t. y c2 t/2 0 u (t), so the The unit step fun tion ĝ(k, t) u (t) is ausal. The time-domain Green's fun tion for the Treeby-Cox spa e-fra tional wave equation is al ulated with Eq. 3.5 ombined with the expressions for ĝ(k, t), and due to the presen e of 43 u (t) in ea h expression for ĝ(k, t), the time-domain Green's fun tion g(r, t) for the Treeby-Cox spa e-fra tional wave equation is also ausal. The phase velo ity and attenuation for the Treeby-Cox spa e-fra tional wave equation are derived from the dispersion relation k2 = ω2 y y y+1 + j2α0 cy−1 , 0 ωk + 2α0 c0 tan (πy/2) k c20 (3.12) starting with the assumption that the attenuation and dispersion are both small relative to ω/c0 . The binomial approximation is then applied to Eq. 3.12, and the rst three terms are retained, yielding k≈ ω c0 + α0 cy0 (jk y + tan (πy/2) c0 ω −1 k y+1 ) (3.13) 2 ω −1 (jk y + tan (πy/2) c0 ω −1 k y+1 ) . − 21 α02 c2y+1 0 On the right hand side, is repla ed with ky is repla ed with (ω/c0 )y [1 + yα0c0 ω y−1 (j + tan (πy/2))] and k y+1 (ω/c0 )y+1 [1 + (y + 1) α0 c0 ω y−1 (j + tan (πy/2))], order and higher terms with respe t to α0 respe tively, and all third are dis arded. The approximate phase velo ity is obtained from the real part of the wavenumber divided by ω,     1 1 1 1 2 y−1 tan (πy/2) α02 c0 ω 2y−2 , ≈ + tan (πy/2) α0 ω + −y + + y + c (ω) c0 2 2 (3.14) and the approximate attenuation is the imaginary part of the wavenumber, α (ω) ≈ α0 ω y + 2y tan (πy/2) α02 c0 ω 2y−1 . (3.15) Treeby and Cox [42℄ instead express the approximate phase velo ity and attenuation for the Treeby-Cox wave equation as 1/c (ω) ≈ 1/c0 + tan (πy/2) α0 ω y−1 respe tively. 44 and α (ω) ≈ α0 ω y , 3.4 Methods 3.4.1 Dispersion Relations The attenuation and phase velo ity of the power law wave equation are al ulated with Eqs. 2.1 and 2.2, respe tively. For the Chen-Holm and the Treeby-Cox spa e-fra tional wave equations, sin e there is no analyti al solution, the phase velo ity and attenuation are both al ulated using the fsolve routine in Matlab. This Matlab fun tion applies the LevenbergMarquardt algorithm to numeri ally determine the omplex value of the dispersion relations in Eqs. 3.6 and 3.12 using the initial value k that solves k = ω/c0 . 3.4.2 The Pantis Method Numeri al al ulations of the time-domain Green's fun tions for the Chen-Holm and Treeby-Cox wave equations are hallenging be ause Eq. 3.5 is a highly os illatory improper integral. When applied to this problem, most standard numeri al integration te hniques perform poorly, so an alternative approa h is required. The Pantis method [56, 57℄ provides an ideal solution to this problem, where the improper integral in Eq. 3.5 is rewritten as the sum of a proper integral and an improper integral a ording to 4π g(r, t) = (2π)3 r ˆ ∞ kĝ(k, t) sin (kr) dk = 0 4π [I1 (r, t) + I2 (r, t)] , (2π)3 r (3.16) with I1 (r, t) = ˆ mπ/r kĝ(k, t) sin (kr) dk, (3.17) 0 and I2 (r, t) = The integral I1 ˆ ∞ kĝ(k, t) sin (kr) dk. (3.18) mπ/r (r, t) is evaluated with Filon's method [58, 59, 60℄, whi h approximates kĝ(k, t) with a pie ewise se ond order polynomial. Filon's method is implemented within a Matlab 45 routine [59℄, where the input parameters in lude the Matlab fun tion handle, the lower and upper limits of integration, and the number of points at whi h the integrand is evaluated. Ea h numeri al evaluation of I1 (r, t) with Filon's method is al ulated with 70,000 abs issas, whi h is su ient for onvergen e at all temporal and spatial points evaluated here for both breast and liver. integral I2 (r, t) After integrating by parts and re ognizing that − For these al ulations, m is Also, when m odd, m m cos (kr) r ˆ ∞ mπ/r I1 (r, t) k=mπ/r − sin (kr) ∂ [kĝ(k, t)] ∂k r2 value the k=mπ/r sin (kr) ∂ [kĝ(k, t)] dk. ∂k r2 (3.19) whi h is an ee tive approximation for these al ula- su iently large and when the ontribution from I2 (r, t) is relatively small. is an odd integer, the se ond term in Eq. 3.19 disappears, and when is onsistently positive for the expressions onsidered here. is advantageous for these al ulations be ause an even value of value for = 0, is an odd integer, and only the rst term in Eq. 3.19 is retained. I2 (r, t) ≈ −mπĝ(mπ/r, t)/r 2 , tions when k=∞ over the upper subinterval is rewritten as I2 (r, t) = kĝ(k, t) Thus, [kĝ(k, t)] I1 (r, t), m = 1601 whi h is undesirable be ause yields ee tive results for all g(r, t) ≥ 0 (r, t) an yield a negative pairs evaluated here. sinh (·) is An odd value of for all values of problems with numeri al overow when the argument of the exp (−a) sinh (b) m m r and t. The Also, to avoid fun tion grows large,  eb−a − e−b−a /2. is instead al ulated using 3.4.3 Time windows for omputed time-domain Green's fun tions The start and end times for all time-domain Green's fun tion al ulations are adjusted manually for ea h lo ation su h that ea h waveform lls a signi ant portion of the time window. For both the Chen-Holm wave equation and the Treeby-Cox wave equation, at t=0 be ause ĝ(k, t) = 0 for all values of obtained through dire t substitution of t=0 k when t = 0, g (r, t) = 0 whi h is the analyti al result into Eqs. 3.3-3.4 and Eqs. 3.10-3.11. The unit 46 2 Attenuation (Np/m) Phase velocity (m/s) liver 1540.005 c(ω) ≈ c 0 1450.2 1450.1 1450 0.1 10 20 30 40 (a) Frequency (MHz) 2500 1/c(ω) ≈ α(ω) ≈ α 0ωy 2000 1500 1000 500 0 0.1 (b) 10 20 30 1540.004 c(ω) ≈ c 0 1540.002 1540.001 1540 0.1 10 20 30 40 (c) Frequency (MHz) Frequency (MHz) Chen-Holm (Eq. 8) 2y-2 2 1/c0-(y-1/2)α 0c0ω 1540.003 40 liver 400 Chen-Holm (Eq. 8) Chen-Holm (Eq. 8) 1/c(ω) ≈ 1/c0-(y-1/2)α 0c0ω2y-2 1450.3 breast 3000 Attenuation (Np/m) breast Chen-Holm (Eq. 8) Phase velocity (m/s) 1450.4 α(ω) ≈ α 0ωy 300 200 100 0 0.1 (d) 10 20 30 40 Frequency (MHz) Figure 3.1: Phase velo ity and attenuation in breast and liver obtained from the dispersion relation in Eq. 3.6 for the Chen-Holm wave equation (• and ◦) (∗) and two dierent approximations to the dispersion relation for the Chen-Holm wave equation. step fun tion u (t) also guarantees that ĝ(k, t) = 0 for all t < 0, so g (r, t) = 0 for all t<0 in the Chen-Holm and Treeby-Cox wave equations. For the power law wave equation, the time-domain Green's fun tion is al ulated in Matlab with the STABLE toolbox [47, 53℄, whi h numeri ally evaluates the expression for the stable pdf in Eq. 2.13. The time window dened for time domain Green's fun tion al ulations very lose to the sour e in ludes negative time values to apture the non ausal omponents of the response in these lo ations. In all plots, the arrival time r/c0 in a lossless medium is also indi ated to provide a time referen e. 3.5 Results 3.5.1 Phase velo ity and attenuation in breast and liver Fig. 3.1 shows the exa t and approximate frequen y-dependent phase velo ity and attenuation for the Chen-Holm spa e-fra tional wave equation using values for human breast with y = 1.5, α0 = 0.086 with Np/ m/MHz y = 1.139, α0 = 0.0459 1.5 c0 = 1450 , and 1.139 Np/ m/MHz , and m/s and using values for human liver c0 = 1540 m/s. The frequen y range is 0.1 to 40 MHz, and the step size on the horizontal axis is ∆f = 0.3 MHz for ea h of the phase velo ity and attenuation plots. For the Chen-Holm wave equation, the result obtained with the Matlab fsolve routine applied to Eq. 3.6 is indi ated by dashed lines with star markers 47 1460 1455 1450 0.1 10 20 30 40 (a) Frequency (MHz) 2500 1560 Treeby-Cox (Eq. 14) Treeby-Cox approx. (Eq. 17) power law wave eqn (Eq. 1) 2000 1500 1000 500 0 0.1 (b) 10 20 30 Frequency (MHz) 40 liver Treeby-Cox (Eq. 14) Treeby-Cox approx. (Eq. 16) power law wave eqn (Eq. 2) 1555 1550 liver 400 Attenuation (Np/m) 1465 Attenuation (Np/m) Phase velocity (m/s) 1470 breast 3000 Treeby-Cox (Eq. 14) Treeby-Cox approx. (Eq. 16) power law wave eqn (Eq. 2) Phase velocity (m/s) breast 1475 1545 0.1 10 20 30 40 (c) Frequency (MHz) Treeby-Cox (Eq. 14) Treeby-Cox approx. (Eq. 17) power law wave eqn (Eq. 1) 300 200 100 0 0.1 (d) 10 20 30 40 Frequency (MHz) Figure 3.2: Phase velo ity and attenuation in breast and liver obtained from the dispersion relation in Eq. 3.12 for the Treeby-Cox wave equation (∗), the approximations to the dispersion relation for the Treeby-Cox wave equation given in Eqs. 3.14-3.15 (•), and the attenuation and phase velo ity for the power law wave equation given in Eqs. 2.1-2.2 (∗), (◦). and the approximate solution from Eq. 3.7 is indi ated by dot-dashed lines with solid dot markers (•). The zero-order approximation to the phase velo ity by a solid line with ir le markers that the approximation to (◦) in Figs. 3.1(a) and 3.1( ). c (ω) ≈ c0 is indi ated Figs. 3.1(a) and 3.1( ) show c (ω) introdu ed in Eq. 3.7 demonstrates ex ellent agreement with the numeri al solution of Eq. 3.6. Figs. 3.1(a) and 3.1( ) also indi ate that the phase velo ity of the Chen-Holm wave equation is not quite equal to c0 be ause c (ω) in reases slightly as the frequen y in reases. In addition, Fig. 3.1(b) shows that the attenuation predi ted by the Chen-Holm wave equation follows the power law indi ated by Eq. 2.1 very losely in breast. Furthermore, Figs. 3.1(a) and 3.1( ) demonstrate that the hange in the phase velo ity in liver over this frequen y range is mu h less than the hange in the phase velo ity in breast, and Figs. 3.1(b) and 3.1(d) indi ate that the attenuation of liver is mu h smaller than that of breast. Fig. 3.2 shows the exa t and approximate values for the phase velo ity and attenuation of the Treeby-Cox spa e-fra tional wave equation and the exa t phase velo ity and attenuation of the power law wave equation using the same parameters evaluated in Fig. 3.1. For the Treeby-Cox wave equation, the numeri al results obtained with the Matlab fsolve routine applied to Eq. 3.12 are indi ated by dashed lines with star markers (∗), and the approximate solution from Eq. 3.14 is indi ated by dot-dashed lines with solid dot markers 48 (•). For the power law wave equation, whi h exa tly satises Eqs. 2.1 and 2.2, the results are indi ated by solid lines with ir le markers (◦). Figs. 3.2(a-d) show that the approximations introdu ed in Eqs. 3.14-3.15 losely mat h the numeri al values obtained from Eq. 3.12. Fig. 3.2 also shows that the phase velo ity and attenuation urves are similar for the power law wave equation and the Treeby-Cox wave equation; however, there is a small but noti eable dieren e between these urves that is asso iated with the higher order terms on the right hand side of Eqs. 3.14-3.15. Fig. 3.2(b) indi ates that the dieren e between the attenuation for the Treeby-Cox wave equation and the power law wave equation in reases slightly as the frequen y grows larger. For example, at f = 40 MHz, the phase velo ities obtained from Eq. 3.12 (the dispersion relation for the Treeby-Cox wave equation) with the fsolve routine in Matlab, Eqs. 3.14-3.15 (approximations to the phase velo ity and the attenuation for the Treeby-Cox wave equation), and Eqs. 2.1-2.2 (the exa t attenuation and phase velo ity for the power law wave equation) in breast are and the orresponding attenuations are respe tively. At f = 40 1553.1 2096.6 1468.2 Np/m, m/s, 1468.2 2093.7 m/s, and Np/m, and 2175.6 MHz, the orresponding phase velo ities in liver are 1553.0 m/s, and 306.59 Np/m, respe tively. m/s, and the attenuations are Thus, at f = 40 300.79 Np/m, 1468.4 300.68 m/s, Np/m, 1553.0 m/s, Np/m, and MHz and at all other frequen ies evaluated in Fig. 3.2, Eqs. 3.14-3.15 are better approximations for the phase velo ity and attenuation of the Treeby-Cox wave equation than Eqs. 2.2 and 2.1, where the largest dieren es are observed at f = 40 MHz. Also, the phase velo ity and attenuation are slightly smaller for the Treeby-Cox wave equation than for the power law wave equation. Figs. 3.2(a) and 3.2( ) also demonstrate that the hange in the phase velo ity in liver over this frequen y range is smaller than the hange in the phase velo ity in breast, and Figs. 3.2(b) and 3.2(d) show that the attenuation for liver is mu h less than that for breast over the range of frequen ies evaluated here. In Figs. 3.2(a) and 3.2( ), the dieren es in the phase velo ities for the power law wave equation and the Treeby-Cox wave equation are due to the ontribution from the third term in Eq. 3.14, and in Figs. 3.2(b) and 3.2(d), the dieren es in the attenuation for 49 the power law wave equation and the Treeby-Cox wave equation are due to the ontribution from the se ond term in Eq. 3.15. 3.5.2 Time-domain Green's fun tions al ulated for breast 15 ×10 10 2.5 0 r = 10nm 7.5 10 8 r = 100µ m 7.5 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 1 0.5 30 40 7 60 70 80 90 100 time (ns) 2.5 t=r/c -2.5 0 r = 1mm 100 ×10 3 1.5 300 400 6 r = 1cm 5 675 700 725 750 time (ns) -5 6.7 (g) 6.9 ×10 3 1 1.5 2 6 r = 10cm power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 2 1 0 t=r/c t=r/c 0.5 time (ns) 4 0 -1.5 650 (f) 0 (d) power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 10 4.5 200 time (ps) 15 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn t=r/c r = 1µ m power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 0 (c) 0 0 -0.5 50 (e) 20 time (ps) ×10 6 1.5 10 9 5 -1 0 ×10 t=r/c 4πr g(r,t) 5 7.5 0 -5 -10 (b) 4πr g(r,t) 2 2.5 r = 100nm power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn t=r/c time (ps) ×10 2.5 4πr g(r,t) 0 10 1 0 -2.5 -2.5 (a) ×10 2 5 t=r/c 3 3 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 10 4πr g(r,t) 5 4πr g(r,t) r = 1nm power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 4πr g(r,t) 11 4πr g(r,t) ×10 4πr g(r,t) 7.5 t=r/c 7.1 7.3 time (µ s) -1 68 68.5 69 69.5 70 70.5 (h) time (µ s) Figure 3.3: Time-domain Green's fun tions s aled by 4πr al ulated for breast with α0 = 0.086 Np/ m/MHz1.5 , and c0 = 1450 m/s at (a) r = 1 nm, (b) r = 10 nm, ( ) nm, (d) r = 1 µm, (e) r = 100 µm, (f ) r = 1 mm, (g) r = 1 m, and (h) y = 1.5, r = 100 r = 10 m with the power law (solid line), Chen-Holm (dashed line), and Treeby-Cox (dot-dashed line) wave equations. At all distan es, the time-domain Green's fun tions for the Chen-Holm and Treeby-Cox wave equations evaluated for breast are ausal while the time-domain Green's fun tion for the power law wave equation is learly non- ausal for Beyond about r = 100 µm, r = 1 nm and r = 10 nm. the time-domain Green's fun tions for the power law wave equation and the Treeby-Cox wave equation are nearly indistinguishable, while the timedomain Green's fun tion for the Chen-Holm wave equation is distin t from the time-domain Green's fun tions for the time-fra tional power law wave equation and the spa e-fra tional Treeby-Cox wave equation at all distan es. To demonstrate some of the similarities and dieren es between these three fra tional wave equations evaluated at various distan es, Fig. 3.3 shows the time-domain Green's fun tions multiplied by 4πr for the power law wave equation (solid line), the Chen-Holm wave equation (dashed line), and the Treeby-Cox wave equation (dot-dashed line) al ulated for human breast with are omputed at y = 1.5, c0 = 1450 r =1 nm, r = 10 nm, m/s, and r = 100 50 α0 = 0.086 nm, 1.5 Np/ m/MHz . The results r = 1 µm, r = 100 µm, r = 1 mm, r = 1 m, and r = 10 m. The units dened for the horizontal axis in Fig. 3.3 are pi ose onds, nanose onds, or mi rose onds. The verti al dashed line des ribes the time t = r/c0 in ea h subgure. Figs. 3.3(a) and 3.3(b) show that the Green's fun tion for the power law wave equation is non ausal at r=1 nm and r = 10 nm sin e the Green's fun tion for the power law wave equation is learly nonzero before time t = 0 in these two subgures. However, the Green's fun tions for the Chen-Holm wave equation and the Treeby-Cox wave equation are both ausal as these two Green's fun tions are always equal to zero for all times t ≤ 0, as demonstrated in se tions 3.2 and 3.3, respe tively. The amplitudes of the Green's fun tions for these three wave equations are similar in ea h individual subgure, and the amplitudes hange by a signi ant amount in ea h su essive subgure. the peaks in Fig. 3.3(a) are all approximately Fig. 3.3(b) are all around 9 × 1010 −1 s . 4 × 1011 −1 s For example, the values of , and the values of the peaks in Three distin t peaks are observed in Figs. 3.3(a) and 3.3(b), and then the peaks of the time-domain Green's fun tions for the power law wave equation and the Treeby-Cox wave equation start to move mu h loser together in Figs. 3.3( ) and 3.3(d) while the peak for the Chen-Holm wave equation remains distin t from the peaks for the other two fra tional wave equations in ea h subgure. In Figs. 3.3(e)-3.3(h), the waveforms for the Chen-Holm wave equation are onsistently quite dierent from those obtained with the power law and Treeby-Cox wave equations. This is expe ted be ause the phase velo ity for the Chen-Holm wave equation is approximately equal to c0 for all frequen y omponents; however, the phase velo ity demonstrates mu h greater variation as a fun tion of frequen y in the power law and Treeby-Cox wave equations. Furthermore, in Figs. 3.3(e)-3.3(h), the time-domain Green's fun tions for the power law wave equation and the Treeby-Cox wave equation are nearly indistinguishable, whi h suggests that these two time-domain Green's fun tions are onverging to the same result. Fig. 3.3 shows that the shape of the time-domain Green's fun tion is the same in all eight subgures for the power law wave equation, where the amplitude and time-s aling is dierent in ea h subgure. However, the shapes of the time-domain Green's fun tions for the other two fra tional wave equations 51 hange as a fun tion of distan e. For instan e, Fig. 3.3(a) shows that the waveforms for the Chen-Holm and Treeby-Cox wave equations both experien e abrupt hanges at t = 0, whereas the orresponding waveforms for these two fra tional wave equations are smooth at all time points shown in Fig. 3.3(e). 3.5.3 Time-domain Green's fun tions al ulated for liver 4.5 2 1 2 0 t=r/c (a) 200 300 ×10 r = 100µ m ×10 2 3 8 0 64 65 66 time (ns) ×10 4 1 68 -1 630 (f) 650 15 20 7 0 r = 1cm 6 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 2 670 690 6.3 (g) 100 150 ×10 6 r = 10cm power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 4 2 0 t=r/c 6.5 6.7 6.9 time (µ s) -2 64 64.5 65 65.5 66 66.5 (h) time (µ s) Figure 3.4: Time-domain Green's fun tions s aled by 4πr al ulated for liver with y = α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s at (a) r = 100 zm, (b) r = ( ) (h) r = 10 am, (d) r = 100 am, r = 10 m with the power law (e) r = 100 µm, 200 time (zs) t=r/c time (ns) 50 (d) 0 t=r/c 67 5 t=r/c 10 time (zs) r = 1mm 2 t=r/c 5 (c) 0 0 63 1 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 3 1 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 0 t=r/c time (zs) 4 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 4πr g(r,t) 2 9 0 (b) time (ys) ×1019 r = 100am -2 4πr g(r,t) 100 10 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 4 0 -1.5 0 ×1020 r = 10am 6 3 t=r/c 4πr g(r,t) 8 1.5 0 -1 (e) r = 1am power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 4πr g(r,t) 3 ×1021 4πr g(r,t) 4 4πr g(r,t) 6 power law wave eqn Chen and Holm wave eqn Treeby and Cox wave eqn 4πr g(r,t) ×1022 r = 100zm 4πr g(r,t) 5 (f ) r = 1 mm, (g) r = 1 1.139, 1 am, m, and (solid line), Chen-Holm (dashed line), and Treeby-Cox (dot-dashed line) wave equations. At all distan es, the time-domain Green's fun tions for the Chen-Holm and Treeby-Cox wave equations evaluated for liver are ausal while the time-domain Green's fun tion for the power law wave equation is learly non- ausal for r = 100 z m. Beyond about r = 100 µm, the time-domain Green's fun tions for the power law wave equation and the Treeby-Cox wave equation are nearly indistinguishable, while the time-domain Green's fun tion for the Chen-Holm wave equation is onsistently distin t from the time-domain Green's fun tions for the time-fra tional power law wave equation and the spa e-fra tional Treeby-Cox wave equation at all distan es. Fig. 3.4 shows the time-domain Green's fun tions multiplied by 4πr for the power law wave equation (solid line), the Chen-Holm wave equation (dashed line), and the Treeby-Cox 52 wave equation (dot-dashed line) al ulated for human liver with and α0 = 0.0459 are evaluated at mm, r = 1 Np/ m/MHz r = 100 m, and yo tose onds (ys or zm, r = 10 10−24 1.139 . y = 1.139, c0 = 1540 m/s, In Fig. 3.4, the three time-domain Green's fun tions r = 1 am, r = 10 am, r = 100 am, r = 100 µm, r = 1 m. The units dened for the horizontal axis in Fig. 3.4 are se onds), zeptose onds (zs or 10−21 se onds), nanose onds (ns), or mi rose onds (µs). The verti al dashed line des ribes the time t = r/c0 in ea h subgure. In Fig. 3.4, the rst row of subgures shows that at relatively short distan es, the time-domain Green's fun tions for the power law wave equation, the Chen-Holm wave equation, and the Treeby-Cox wave equation are again all noti eably dierent. Fig. 3.4(a) shows that the r = 100 time-domain Green's fun tion for the power law wave equation is non ausal at zm sin e this time-domain Green's fun tion is learly nonzero before time t = 0. However, the time-domain Green's fun tions for the Chen-Holm wave equation and the Treeby-Cox wave equation are both ausal as these are onsistently equal to zero for all times t ≤ 0. Figs. 3.4(a)-3.4(d) show that, ompared to the results al ulated for breast, a mu h shorter distan e is required to show the non ausal behavior of the power law wave equation in liver. The amplitudes of the Green's fun tions for these three wave equations are again similar in ea h individual subgure, and the amplitudes again de rease by a signi ant amount in ea h su essive subgure. For instan e, the values of the peaks in Fig. 3.4(a) are all approximately 2×1022 s−1 , and the values of the peaks in Fig. 3.4(b) are all near 3×1021 s−1 . Three distin t peaks are observed in Figs. 3.4(a)-3.4(d), while the peaks for the power law wave equation and the Treeby-Cox wave equation are very lose together in Figs. 3.4(e)-3.4(h). In ea h subgure, the lo ation of the peak for the Chen-Holm wave equation is distin t from the peak lo ations for the other two wave equations. Furthermore, in Figs. 3.4(e)-3.4(h), the time-domain Green's fun tions for the power law wave equation and the Treeby-Cox wave equation are nearly indistinguishable, and with in reasing distan e, these two time-domain Green's fun tions are onverging to the same result. 53 3.5.4 Amplitude and full width at half maximum (FWHM) values in breast and liver Figure 3.5: The (a, ) amplitudes and (b,d) FWHM values of the time-domain Green's fun tions al ulated for the power law, Chen-Holm, and Treeby-Cox wave equations in (a,b) breast and ( ,d) liver. The amplitudes of all three time-domain Green's fun tions de rease as the distan e in reases while the FWHM values of all three time-domain Green's fun tions in rease as the distan e in reases. The amplitudes of the time-domain Green's fun tions for all three of these fra tional wave equations are very similar at ea h distan e, and the FWHM values are all approximately the same at longer distan es, although there is a small dieren e in the FWHM values at shorter distan es that diminishes with in reasing distan e. Fig. 3.5 shows the amplitudes and full width at half maximum (FWHM) values for the time-domain Green's fun tions multiplied by 4πr in breast and liver for the power law wave equation (solid line), the Chen-Holm wave equation (dashed line), and the Treeby-Cox wave equation (dot-dashed line), where these four subgures summarize the ee ts of attenuation and dispersion in the time-domain. The initial distan es are r = 100 am for liver, and ea h plot ends at r = 10 m. r = 100 nm for breast and Figs. 3.5(a) and 3.5( ) indi ate that the amplitudes of the time-domain Green's fun tions for these three fra tional wave equations, whi h are plotted on a log-log s ale, de rease with frequen y as the distan e in reases due to the ee t of attenuation. Also, the amplitudes of the time-domain Green's fun tions for three wave equations are all very similar at all distan es for both breast and liver. The slopes of the amplitudes in Fig. 3.5(a) and 3.5( ) are −0.6677 and −0.8790, respe tively, where the absolute values of these two quantities are approximately equal to the re ipro als of the power law exponents in breast and liver, respe tively. Figs. 3.5(b) and 3.5(d), whi h are plotted 54 on a log-log s ale, show that the FWHM values of the time-domain Green's fun tions for these three fra tional wave equations in rease with distan e due to the ee t of dispersion. In Figs. 3.5(b) and 3.5(d), the FWHM values of the time-domain Green's fun tions for these three fra tional wave equations dier slightly lose to the sour e. For example, the FWHM values of the time-domain Green's fun tions for the power law, Chen-Holm, and Treeby-Cox wave equations are respe tively, and 4.4 × 10−11 2.0 × 10−21 s, s, 4.1 × 10−11 1.7 × 10−21 s, and 3.3 × 10−11 s at s, and 9.4 × 10−22 s at r = 100 nm in breast, r = 100 am in liver, respe tively. At larger distan es, the FWHM values of the time-domain Green's fun tions for all three wave equations yield nearly the same result even though the shape of the time-domain Green's fun tion for the Chen-Holm wave equation is distin t from the shapes of the time-domain Green's fun tions for the power law and Treeby-Cox wave equations. Figs. 3.5(b) and 3.5(d) also show that, at r = 10 m, the FWHM values of the time-domain Green's fun tions for the power law, Chen-Holm, and Treeby-Cox wave equations are equal to 4.4 × 10−7 1.9 × 10−7 s, s, and 4.1 × 10−7 2.2 × 10−7 the FWHM values are s, and 4.4 × 10−7 s in breast, respe tively, and 2.2 × 10−7 s, s in liver, respe tively. In Fig. 3.5(b) and 3.5(d), the slopes of 0.6711 and 0.8825, whi h are approximately equal to the re ipro als of the power law exponents in breast and liver, respe tively. 3.5.5 Time-domain Green's fun tions onvolved with a three y le Hanning-weighted pulse Fig. 3.6 shows the waveforms obtained when a three y le Hanning-weighted pulse [50, 51℄ is onvolved with the time-domain Green's fun tions for the power law wave equation (solid line), the Chen-Holm wave equation (dashed line), and the Treeby-Cox wave equation (dot-dashed line) multiplied by α0 = 0.086 Np/ m/MHz 1.5 , and 4πr c0 = 1450 evaluated at r = 1 m in breast with y = 1.5, m/s. This gure shows how power law attenua- tion and phase velo ity inuen e the shape of a short pulse. In Fig. 3.6, input pulses with two 55 Figure 3.6: Simulated three- y le Hanning-weighted pulses with enter frequen ies (a) 7.5 MHz and (b) f0 = 29 MHz onvolved with time-domain Green's fun tions for the power law, Chen-Holm, and Treeby-Cox wave equations multiplied by at r=1 f0 = 4πr evaluated for breast m. The onvolution results for the power law wave equation and the Treeby-Cox wave equation are very similar while the onvolution result for the Chen-Holm wave equation learly shows a time delay. Signi ant attenuation and waveform spreading are observed for all three signals in (b) produ ed by inputs with f0 = 29 MHz, whereas a moderate amount of attenuation and waveform spreading is observed for all three signals in (a) produ ed by inputs with f0 = 7.5 MHz. dierent enter frequen ies, f0 = 7.5 MHz and f0 = 29 MHz, are evaluated to highlight the dieren es in the resulting waveforms. Fig. 3.6 indi ates that the onvolution results for the power law wave equation and the Treeby-Cox wave equation are nearly indistinguishable for pulses with f0 = 7.5 MHz and f0 = 29 MHz while the onvolution results for the Chen-Holm wave equation evaluated at these two frequen ies ontain a signi ant time delay be ause of the substantial dieren e in the phase velo ity. In parti ular, when the time delay is dened as the time between the peaks, the time delay between the onvolution results for the power law wave equation and the Chen-Holm wave equation is 35 ns in Fig. 3.6(a), and the time delay between the onvolution results for the power law wave equation and the Chen-Holm wave equation is 41 ns in Fig. 3.6(b), whereas the time delay between the onvolution results for the power law wave equation and the Treeby-Cox wave equation is less than 56 1 ns in Figs. 3.6(a) and 3.6(b). In Fig. 3.6(a), the amplitudes of all three waveforms are omparable and the waveform shapes are also similar. This indi ates that the attenuation is nearly the same in all three waveforms and that the amount of waveform spreading (dispersion) for the propagating waveforms is also approximately the same. In Fig. 3.6(b), the amplitudes are again similar and the amount of spreading is similar for all three waveforms. The onvolution results for the power law wave equation and the Treeby-Cox wave equation are again nearly indistinguishable, but the waveform shape for the Chen-Holm wave equation onvolution result is somewhat dierent from the other two, where the sour e of this dieren e is again the phase velo ity. Also, relative to Fig. 3.6(a), the signal amplitude drops o onsiderably in Fig. 3.6(b), and the ltering and spreading of the signal in the time domain is mu h more signi ant in Fig. 3.6(b), where the initial duration of the three y le is 0.103 µs. these at 0.035 µs, Sin e the FWHM values are approximately equal to r=1 m, whi h is larger than one period of the 29 0.09 µs 29 MHz pulse for all three of MHz enter frequen y, namely the attenuation and waveform spreading are more signi ant in Fig. 3.6(b) than in Fig. 3.6(a). In ontrast, one period of the larger than the FWHM values 7.5 MHz enter frequen y is 0.133 µs, whi h is 0.09 µs of all three time-domain Green's fun tions at r = 1 m, so there is mu h less attenuation and waveform spreading in Fig. 3.6(a). Fig. 3.7 shows the results obtained when a three y le Hanning-weighted pulse is onvolved with the time-domain Green's fun tions for the power law wave equation (solid line), the Chen-Holm wave equation (dashed line), and the Treeby-Cox wave equation (dot-dashed line) multiplied by α0 = 0.0459 to f0 = 7.5 Np/ m/MHz MHz and 1.139 f0 = 29 4πr , and evaluated at c0 = 1540 r = 1 m for liver with y = 1.139, m/s. The enter frequen ies are again equal MHz. As indi ated in Fig. 3.7(a), the onvolution results for the power law wave equation and the Treeby-Cox wave equation also tra k very losely for both pulses while the onvolution result for the Chen-Holm wave equation ontains a noti eable time delay. For example, the time delay between the onvolution results for the power law wave equation and the Chen-Holm wave equation is 57 44 ns in Fig. 3.7(a), Figure 3.7: Simulated three- y le Hanning-weighted pulses with enter frequen ies (a) 7.5 MHz and (b) f0 = 29 MHz onvolved with time-domain Green's fun tions for the power law, Chen-Holm, and Treeby-Cox wave equations multiplied by r = 1 m. f0 = 4πr evaluated for liver at The onvolution results for the power law wave equation and the Treeby-Cox wave equation are very similar while the onvolution result for the Chen-Holm wave equation learly shows a time delay. A moderate amount of attenuation and waveform spreading are observed for all three signals in (b) produ ed by inputs with f0 = 29 MHz, whereas a smaller amount of attenuation and waveform spreading is observed for all three signals in (a) produ ed by inputs with f0 = 7.5 MHz. and the time delay between the onvolution results for the power law wave equation and the Chen-Holm wave equation is 52 ns in Fig. 3.7(b), whereas the time delay between the onvolution results for the power law wave equation and the Treeby-Cox wave equation is again less than 1 ns in Figs. 3.7(a) and 3.7(b). Fig. 3.7(a) also shows that the amplitudes of all three waveforms are almost the same and the shapes of these waveforms are again similar, whi h indi ates that the amount of attenuation and dispersion is nearly the same in all three waveforms. Fig. 3.7 shows that there is more attenuation and dispersion for the 29 MHz signal in Fig. 3.7(b) relative to the 7.5 MHz signal in Fig. 3.7(a) and that the shape of the three y le Hanning-weighted input pulse is still re ognizable in Fig. 3.7(b). 58 3.6 Dis ussion 3.6.1 Numeri al evaluations of the inverse 3D Fourier transform Numeri al evaluations of improper integrals with highly os illatory integrands are often hallenging. Certain methods, su h as the Matlab quadgk routine, whi h is based on adaptive Gauss-Kronrod quadrature, and the Matlab integral routine, whi h is based on global adaptive quadrature, are able to evaluate some improper integrals. Unfortunately, neither of these methods onsistently onverge to the orre t result when applied to Eq. 3.5. Although Filon's formula provides an ee tive method for evaluating integrals with highly os illatory integrands, Filon's formula yields in onsistent results when applied to Eq. 3.5 be ause the upper limit of integration is innite. However, the Pantis method [56, 57℄ solves this problem by applying Filon's formula to the integral I1 (r, t) with nite limits of integration, and then the rst term of an innite sum approximates the ontribution from the remaining improper integral I2 (r, t). To determine whether numeri al onvergen e is a hieved at a given distan e, the result is ompared to the result obtained with twi e as many Filon abs issas and to a value for m that is twi e as large. If the Eu lidean norm of 1 × 10−3 , the dieren e between the two results is within lose and onvergen e is a hieved. for liver, 900 Filon abs issas and For instan e, at m = 21 r = 100 µm points is a hieved with 35,000 Filon abs issas and r and t r = 1 nm for breast and r = 100 zm are su ient to a hieve onvergen e at all time points in Figs. 3.3(a) and 3.4(a). However, at values of then the two results are su iently for liver, onvergen e at all time m = 801. Convergen e is a hieved at all for al ulations in breast and liver with 70,000 abs issas and m = 1601, so instead of spe ifying these values at ea h distan e for ea h material, these two values are used throughout. Fig. 3.8 des ribes the time-domain Green's fun tions s aled by the Pantis method at r = 1 4πr al ulated with m with dierent numbers of Filon abs issas and values of 59 15 ×10 6 r = 1cm 15 ×10 6 Treeby and Cox wave eqn 15 5 0 6.8 6.9 7 7.1 7.2 -5 6.7 7.3 6.8 6.9 (a) 7 7.1 7.2 -5 6.7 7.3 4 7 time (µ s) (b) ( ) ×10 11 m = 101, 4πr r = 1nm 4 ×10 11 7.1 7.2 7.3 evaluated for breast at and ( ) 500 Filon abs issas and m = 401, m = 401. r = 1nm Treeby and Cox wave eqn Treeby and Cox wave eqn 3 4πr g(r,t) 3 4πr g(r,t) 6.9 m omputed with the Pantis method using (a) 2,000 Filon abs issas and (b) 2,000 Filon abs issas and 2 1 0 2 1 0 -1 -1 0 1 2 3 4 0 1 2 3 time (ps) time (ps) (a) (b) Figure 3.9: Computed time-domain Green's fun tions s aled by m. 6.8 time (µ s) Figure 3.8: Computed time-domain Green's fun tion s aled by r=1 5 0 time (µ s) r=1 r = 1cm 10 4πr g(r,t) 4πr g(r,t) 0 -5 6.7 6 Treeby and Cox wave eqn 10 5 ×10 Treeby and Cox wave eqn 10 4πr g(r,t) r = 1cm nm using 500 Filon abs issas and m=5 4πr 4 evaluated for breast at (a) without and (b) with the Pantis term. A omparison between Figs. 3.8(a) and 3.8(b) indi ates that when m is too small, arti ial os illations appear in the omputed Green's fun tion. This problem is addressed by in reasing the value of m. Also, a omparison between Figs. 3.8(a) and 3.8( ) shows that when the number of abs issas is insu ient, the amplitude of the waveform is diminished by a signi ant amount relative to the orre t result for the omputed time-domain Green's fun tion. Fig. 3.9 shows the time-domain Green's fun tions s aled by after (b) the Pantis method is applied at r = 1 60 4πr al ulated before (a) and nm with 500 Filon abs issas and m = 5. These two waveforms mat h losely for larger values of ontribution from lose to 0. For I2 (r, t) t ≤ 0.5 t. As indi ated in Fig. 3.9, the is parti ularly important at shorter distan es, espe ially when ps at r=1 nm in breast, the ontribution from 1% of the peak value of the time-domain Green's fun tion. only the ontribution from I1 (r, t) If the I2 (r, t) t is is larger than I2 (r, t) term is omitted and al ulated with the Filon's formula is in luded, arti ial os illations appear in the omputed Green's fun tion for small values of t. 3.6.2 Improved approximations for the attenuation and phase velo ity Figs. 3.1 and 3.2 show that the attenuation of the Chen-Holm wave equation is a urately represented by a power law and that the attenuation of the Treeby-Cox wave equation is very lose to a power law. Figs. 3.1 and 3.2 also show that a more a urate representation for the attenuation of the Treeby-Cox wave equation is a hieved when the se ond term in Eq. 3.15 is in luded. In addition, Figs. 3.1 and 3.2 demonstrate that more a urate representation for the phase velo ities of the Chen-Holm and Treeby-Cox wave equations are obtained when a se ond (Chen-Holm) or a third (Treeby-Cox) term is in luded in the expression for c (ω). The similarities and dieren es between these phase velo ities are further reinfor ed when binomial approximations are evaluated for ea h expression. For instan e, the phase velo ity obtained from the binomial approximation to Eq. 3.7 is c (ω) ≈ c0 + (y − 1/2) α02 c30 ω 2y−2 [(y − 1/2) α02 c20 ] 1012 −1/(2y−2) , radians/se ond (or whi h yields for the Chen-Holm wave equation when c (ω) ≈ 1450 + 9.0899 × 10−10 ω f ≪ 2.5388 × 1011 Hz) for breast. Thus, when c (ω) ω ≪ ω ≪ 1.5952 × in reases linearly as a fun tion of frequen y when the Chen-Holm wave equation is evaluated for breast, as indi ated in Fig. 3.1(a). In Fig. 3.1( ), the phase velo ity of the Chen-Holm wave equation is approximately ans/se ond (or c (ω) ≈ 1540 + 1.6049 × 10−5 ω 0.278 , f ≪ 8.2103 × 1027 whi h holds when ω ≪ 5.1587 × 1028 radi- Hz) for liver. When plotted on the same axes, the results obtained from these approximations are indistinguishable from Eq. 3.7 and are therefore not 61 shown. When the binomial approximation is applied to Eq. 3.14, this yields c (ω) ≈ c0 − tan (πy/2) α0 c20 ω y−1 + [tan2 (πy/2) − (−y + 1/2 + (y + 1/2) tan2 (πy/2))] α02 c30 ω 2y−2 Treeby-Cox wave equation when ω ≪ [|tan (πy/2)| α0 c0 ]−1/(y−1) , for the and similarly, the binomial approximation applied to the phase velo ity for the power law wave equation in Eq. 2.2 yields c (ω) ≈ c0 −tan (πy/2) α0 c20 ω y−1 + tan2 (πy/2) α02 c30 ω 2y−2 ω ≪ [|tan (πy/2)| α0 c0 ]−1/(y−1) . when the angular frequen y satises In Fig. 3.2(a), the frequen y-dependent phase velo ities of both the Treeby-Cox wave equation and the power law wave equation for breast are represented by the approximate expressions for the Treeby-Cox wave equation and c (ω) ≈ 1450+1.1481×10−3ω 0.5 −4.0367×10−25 ω c (ω) ≈ 1450 + 1.1481 × 10−3 ω 0.5 + 9.0899 × 10−10 ω for the power law wave equation, whi h hold when f ≪ 2.5388 × 1011 ω ≪ 1.5952 × 1012 radians/se ond (or Hz). In Fig. 3.2( ), the phase velo ities of the Treeby-Cox wave equation and the power law wave equation for liver are represented by the approximate expressions c (ω) ≈ 1540 + 0.8864ω 0.139 − 3.0994 × 10−4ω 0.278 for the Treeby-Cox wave equation and c (ω) ≈ 1540+0.8864ω 0.139 +5.1016×10−4ω 0.278 for the power law wave equation, respe tively, when ω ≪ 2.0356×1023 radians/se ond (or f ≪ 3.2398×1022 Hz). axes, the results obtained from these two approximations for from Eq. 3.14 and are therefore not shown. When plotted on the same c (ω) are also indistinguishable The rst two terms in these two binomial expressions for the phase velo ity of the Treeby-Cox and power law wave equations are the same, and the rst term in the expression for the attenuation of the Treeby-Cox wave equation in Eq. 3.15 is also the same as that for the power law wave equation. Thus, although Figs. 3.1 and 3.2 indi ate that there is a small but noti eable dieren e in α (ω) and c (ω) for the Treeby-Cox and power law wave equations over the frequen y range from mu h loser agreement in α (ω) is expe ted for smaller values of and ω. c (ω) 0.1 to 40 MHz, for the Treeby-Cox and power law wave equations Furthermore, the rst terms in the series expansions for both the attenuation and the phase velo ity are the same for all three of these fra tional wave equations. velo ity c (ω) However, the se ond term in the binomial approximation for the phase of the Chen-Holm wave equation is quite dierent from the se ond terms in 62 the binomial approximations for the other two fra tional wave equations, where these terms are responsible for the signi ant dieren es between the phase velo ities shown in Fig. 3.1 and Fig. 3.2. 3.6.3 Time-domain Green's fun tions Figs. 3.3 and 3.4 show several examples of time-domain Green's fun tions for the ChenHolm and Treeby-Cox fra tional wave equations, whi h are ausal for all values of the power law exponent 1 < y ≤ 2. Figs. 3.3 and 3.4 also show that the time-domain Green's fun tion for the power law wave equation is non ausal for time-domain at very short distan es. 1 < y ≤ 2, whi h is only evident in the At longer distan es, Figs. 3.3 and 3.4 suggest that the time-domain Green's fun tions for the Treeby-Cox wave equation and the power law wave equation onverge to the same result. Figs. 3.3 and 3.4 also show that amplitudes are similar for all three Green's fun tions in ea h subgure. The amplitudes for all three of these Green's fun tions are proportional to the s ale fa tor (α0 r)−1/y , whi h was also demonstrated for the three time-fra tional wave equations evaluated in Zhao and M Gough [53℄. The peak values are all about the same for ea h subplot in Figs. 3.3 and 3.4, whi h is onrmed by the results in Figs. 3.5(a) and 3.5( ) showing the amplitudes al ulated for the time-domain Green's fun tions of all three fra tional wave equations. Figs. 3.5(a) and 3.5( ) therefore demonstrate that the attenuations are all about the same for all three of these fra tional wave equations, whi h is expe ted sin e the attenuations of all three fra tional wave equations either exa tly or approximately follow the power law des ribed by Eq. 2.1. Figs. 3.3 and 3.4 indi ate that three distin t peaks are observed in the time-domain Green's fun tions at shorter distan es, and then the peak lo ations for the time-domain Green's fun tions of the power law wave equation and the Treeby-Cox wave equation onverge at larger distan es. In Figs. 3.3 and 3.4, temporal spreading is observed in all three waveforms. The amount of temporal spreading is also determined by (α0 r)1/y , whi h appears in the denominator of the argument for the stable distribution in Eq. 2.13. The temporal spreading is also about the same for all three of these fra tional wave equations, as shown in Figs. 3.5(b) and 3.5(d). 63 These plots, whi h des ribe the FWHMs of the three fra tional wave equations, all mat h very losely. Thus, the amplitude peaks, the temporal spreading, and the rst terms in the series expansions for α (ω) and c (ω) are nearly the same for all three of these fra tional wave equations, but the remaining terms in the expressions for the phase velo ity c (ω) of the Chen-Holm wave equation are distin t from the orresponding terms in the expressions for the phase velo ities of the other two fra tional wave equations. This suggests that the attenuation α (ω) is primarily responsible for both the de ay in the peak values and the temporal waveform spreading in these three fra tional wave equations. The phase velo ity c (ω), parti ularly the rst non- onstant term in the binomial approximation for c (ω), is responsible for the `skew' or symmetry/asymmetry of the time-domain Green's fun tions shown in Figs. 3.3 and 3.4. 3.6.4 Dispersion The Chen-Holm wave equation is `dispersionless' in the sense that the phase velo ity c (ω) is nearly equal to a onstant, whi h yields the symmetri time-domain Green's fun tions depi ted in Figs. 3.3 and 3.4, where the axis of symmetry is dened as t = r/c in ea h subplot. Two other examples of `dispersionless' wave equations are the Stokes and Bla ksto k wave equations. Both of these also have nearly onstant phase velo ity c (ω), and the time-domain Green's fun tions for both of these in the far eld are approximately represented by Gaussian fun tions that are entered at redu es to c (ω) = c0 , t = r/c. When y=2 in the power law wave equation, Eq. 2.2 whi h is also dispersionless. Also, the time-domain Green's fun tion for the power law wave equation with y=2 is a Gaussian entered at t = r/c. Thus, these four examples of `dispersionless' lossy wave equations all possess exa tly or approximately symmetri time-domain Green's fun tions. In a ousti s and medi al ultrasound, the term `dispersion' has two dierent meanings. In some ontexts, the dispersion refers to the temporal spreading, and in others, the dispersion des ribes the frequen y-dependent phase velo ity the attenuation α (ω), c (ω). However, the ontribution from not the phase velo ity, is responsible for the waveform spreading in 64 all three of these fra tional wave equations. Thus, these two meanings for the dispersion interestingly refer to two ompletely dierent aspe ts of ultrasound propagation in soft tissue. that In parti ular, the Chen-Holm wave equation, whi h is dispersionless in the sense c (ω) is nearly onstant, demonstrates a onsiderable amount of waveform spreading as the propagation distan e in reases in the time-domain Green's fun tions shown in Figs. 3.3 and 3.4 and in the FWHM values shown in Figs. 3.5(b) and 3.5(d). This is in ontrast to the Treeby-Cox and power law wave equations, whi h are also dispersive in the sense that the phase velo ity c (ω) varies with frequen y. However, all three of these fra tional wave equations have approximately the same amount of dispersion in terms of the FWHM values shown in Figs. 3.5(b) and 3.5(d). 3.6.5 Convolving time-domain Green's fun tions with short pulses Figs. 3.6 and 3.7 des ribe the waveforms obtained when a three y le Hanning-weighted pulse is onvolved with the time-domain Green's fun tions for the three fra tional wave equations evaluated in breast and liver. Although Figs. 3.1 and 3.2 indi ate that there is a slight dieren e in the phase velo ity for the Treeby-Cox and the power law wave equations at 7.5 MHz and 29 MHz, the time-domain Green's fun tions for the Treeby-Cox and the power law wave equations are nearly identi al in Figs. 3.3(g) and 3.4(g) at r=1 m, so the results obtained in Figs. 3.6 and 3.7 with the Treeby-Cox and power law wave equation are less sensitive to the dieren es in the attenuation, whi h are ltered out after propagating 1 m. Also, the rst two terms in the binomial expansions for the phase velo ities of the Treeby-Cox and the power law wave equations are identi al. These observations, along with the results shown in Figs. 3.3 and 3.4, suggest that the Treeby-Cox and the power law wave equations exhibit very similar behaviors at distan es r ≥ 100 µm, but at very short distan es and for very high frequen y ex itations, some dieren es are observed. In addition, the dieren es between the time-domain Green's fun tion for the Chen-Holm wave equation and the time-domain Green's fun tions for the other two wave equations shown in Figs. 3.3(g) 65 and 3.4(g) are also ree ted in Figs. 3.6 and 3.7, whi h is ree ted in the dieren es in the phase velo ities. 3.6.6 Comparisons with time-domain Green's fun tions al ulated with 3D FFTs Time-domain Green's fun tions are also al ulated with 3D fast Fourier transforms (FFTs) using the approa h des ribed in Treeby and Cox [61℄ and ompared to the Green's fun tions results omputed with the Pantis method shown in Figs. 3.3 and 3.4. Although the analyti al expressions evaluated with ea h approa h are similar, the numeri al performan e of these two methods is quite dierent. For example, a single al ulation with the Pantis method was ompleted in 35 se onds on a Mi rosoft Surfa e Pro with an Intel Core m3-6Y30 CPU  0.90 GHz, whereas the al ulation using the same parameters with 3D FFTs required ompute servers (Dual Intel Xeon E5-26xx  2.3 GHz, 2.4 GHz, and 2.7 GHz) with 48-88 ores and 384-768 GB RAM, whi h took anywhere between 15 minutes and a few hours depending on a variety of fa tors in luding the number of pro essors available and memory/CPU usage of other jobs running on the ompute servers. The al ulations with the Pantis method were performed serially, whereas the al ulations with 3D FFTs were performed in parallel using `parfor' al ulations in Matlab using up to 36 `workers.' Reasonably a urate results were obtained when al ulations with 3D FFTs were performed with 512 x 512 x 512 spatial points (i.e., 512 points in ea h dire tion). When 3D FFTs with more points were evaluated (only radix 2 FFTs were onsidered here), Matlab either rashed or ran out of memory, and when 3D FFTs with fewer points were evaluated, the omputed result usually demonstrated severe nonphysi al os illations or other undesirable numeri al artifa ts. We also found that time-domain Green's fun tion al ulations with 3D FFTs are very sensitive to the grid spa ing. For instan e, Gibbs os illations (whi h often in lude nonphysi al negative values for the time-domain Green's fun tion) were ommonly observed in the results 66 Figure 3.10: Time-domain Green's fun tions s aled by 4πr al ulated for breast with y = 1.5, α0 = 0.086 Np/ m/MHz1.5 , and c0 = 1450 m/s evaluated at (a) r = 10 nm, (b) r = 100 nm, and ( ) r = 1 µm with the Pantis method (solid line) and the 3D FFT approa h (dot-dashed dx = dy = dz = 0.5 nm in (a), (b), and ( ). line) using Figure 3.11: Time-domain Green's fun tions s aled by 4πr al ulated for liver with y = 1.139, α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s evaluated at (a) r = 100 zm, (b) r = 1 am, r = 10 am with the Pantis method (solid line) using dx = dy = dz = 50 zm in (a), (b), and ( ). and ( ) line) and 3D FFT approa h (dot-dashed obtained with 3D FFTs. These os illations were redu ed by de reasing the spa ing between adja ent points in the omputational grid (i.e., by de reasing dx, dy , and dz ). However, if the grid spa ing is too small in al ulations with 3D FFTs, then this auses problems with frequen y-domain aliasing, whi h introdu es errors in the `heavy tail' of the time-domain Green's fun tion. If the grid spa ing is too small by a small amount, then a small to moderate oset from the orre t value is observed in the `tail'. If the grid spa ing is too small by a larger amount, the `tail' grows with in reasing time instead of slowly de aying to zero as shown in Figs. 3.3 and 3.4. In extreme ases, when the value of the grid spa ing is mu h too 67 large, most if not all values in the omputed time-domain Green's fun tion are mu h larger than the orre t values. Usually, Gibbs os illations or errors in the `heavy tail' or both are observed in time-domain Green's fun tion al ulations with 3D FFTs, although we were often able to manually sele t appropriate parameters that a hieve a reasonable trade-o between these two types of numeri al artifa ts for the results shown in Figs. 3.3(a-h) and 3.4(a-d). Examples of typi al results obtained after some parameter tuning are shown in Fig. 3.10(b) and 3.11(b), whi h ontain Gibbs os illations at the very beginning and good agreement with the Pantis results elsewhere. Interestingly, when the grid spa ing is sele ted to redu e both types of numeri al artifa ts for one value of r, the time-domain Green's fun tion al ulated with 3D FFTs in other lo ations (e.g., r/10 and 10r ) onsistently yield signi ant errors. This ee t is shown in Fig. 3.10 with three dierent time-domain Green's fun tions al ulated with 3D FFTs evaluated for breast at value for the grid spa ing, namely r = 10 nm, r = 100 dx = dy = dz = 0.5 nm, and nm. r = 1 µm with the same Fig. 3.10(b) ( enter panel) shows some Gibbs os illations in the result omputed with 3D FFTs near t = 0 and also good agreement in the heavy tail. In Fig. 3.10(a) (left panel), the grid spa ing is too large. When the grid spa ing is redu ed to dx = dy = dz = 0.2 nm or dx = dy = dz = 0.1 nm in Fig. 3.10(a) (assuming that the number of spatial points is xed and equal to 512 x 512 x 512), the agreement is mu h better (not shown), although some small Gibbs os illations still remain in the result obtained with 3D FFTs near t = 0. In Fig. 3.10( ) (right panel), errors in the heavy tail are produ ed by frequen y-domain aliasing. The errors in Fig. 3.10( ) are alleviated by in reasing the grid spa ing to approximately small amount of ringing at time t = 0. al ulations evaluated for liver at 5 nm, whi h also introdu es a Similar results are also observed in Green's fun tion r = 100 zm, r=1 am, and r = 10 am, whi h are shown in Fig. 3.11. Although time-domain Green's fun tions al ulated with 3D FFTs produ e a 68 Figure 3.12: Time-domain Green's fun tions s aled by 4πr al ulated for liver with y = 1.139, α0 = 0.0459 Np/ m/MHz1.139 , and c0 = 1540 m/s evaluated at r = 10 m with the Pantis method (solid line) and with 3D FFTs (dot-dashed line) using (a) (b) 3D dx = r/200, ( ) dx = r/300, (d) dx = r/400, and (e) dx = r/500. FFTs evaluated here, dx = dy = dz . dx = r/100, In all simulations with large grid of values, the results are only a urate in ertain lo ations where the grid spa ing is sele ted to avoid problems with Gibbs os illations and with frequen y-domain aliasing. We also found that the time-domain Green's fun tion al ulations with 3D FFTs failed to onverge for the results shown in Figs. 3.4(e-h). This is demonstrated in Fig. 3.12 using a time-domain Green's fun tion al ulated with 3D FFTs evaluated for liver at r = 10 m with grid spa ings (dx = dy = dz) of 1 mm, 500 µm, 333 µm, 250 µm, and 200 µm. types of problems are also observed for r = 100 µm, r = 1 mm, and r=1 The same m. Thus, the Pantis method is also useful for determining whi h parameter ombinations yield reasonable results in time-domain Green's fun tion al ulations with 3D FFTs. 3.7 Con lusion Phase velo ities and attenuation values were evaluated over a range of ultrasound frequen ies, and time-domain Green's fun tions were al ulated for the Chen-Holm, Treeby-Cox, and power law wave equations in breast and liver. In addition, the amplitudes and FWHM values of the time-domain Green's fun tions for these three fra tional wave equations were al ulated, and three- y le Hanning-weighted pulses at two dierent enter frequen ies were onvolved with the time-domain Green's fun tions for three fra tional wave equations. An additional term was derived for the power series that des ribes 69 c (ω) for the Chen-Holm wave equation and additional terms were derived for the power series that des ribe α (ω) c (ω) and for the Treeby-Cox wave equation, where these new power series more losely mat h the results obtained by numeri ally evaluating the dispersion relation than the previous approximations. Causality was also demonstrated analyti ally in the time domain for both the Chen-Holm and Treeby-Cox wave equations. The Pantis method was introdu ed as an ee tive approa h for evaluating the highly os illatory improper integrals that arise in numeri al al ulations of the time-domain Green's fun tions for the Chen-Holm and Treeby-Cox spa e-fra tional wave equations. The time-domain Green's fun tions for all three time-domain Green's fun tions show a similar amount of temporal spreading and amplitude redu tion, but the time-domain Green's fun tions for the Treeby-Cox and power law wave equations are skewed whereas the time-domain Green's fun tion for the Chen-Holm wave equation is symmetri . The Chen-Holm spa e-fra tional wave equation is non-dispersive in the sense that the phase velo ity is nearly onstant, but the Chen-Holm wave equation is also dispersive in the sense that waveform spreading is learly evident in the omputed time-domain Green's fun tions. The Treeby-Cox wave equation is dispersive in the sense that the phase velo ity is non- onstant and also in the sense that waveform spreading is learly evident in the omputed time-domain Green's fun tions. The Chen-Holm and Treeby-Cox spa e-fra tional wave equations demonstrate approximately the same amount of attenuation and waveform spreading, but the phase velo ities, time-domain Green's fun tions, and onvolution results obtained with these time-domain Green's fun tions all dier signi antly. Also, although the attenuation and phase velo ity for the Treeby-Cox and power law wave equation dier slightly, the time-domain Green's fun tions and the onvolution results obtained with the Green's fun tions for these two fra tional wave equations onverge to the same result. 70 Chapter 4 Perfe tly mat hed layers for nonlinear ultrasound simulations with the KZK equation 4.1 Introdu tion High intensity pulses are often used in appli ations of therapeuti ultrasound. Linear assumptions do not always predi t the orre t results in these ases. The attenuation and dira tion are also important in models of nonlinear ultrasound propagation. Thus, a model ombining dira tion, attenuation, and nonlinearity is needed to analyze the propagation of ultrasoni waves. The Khokhlov-Zabolotskaya-Kuznetsov (KZK) equation [62, 63℄ is a paraboli approximation to the Westervelt equation. Lee and Hamilton [17℄ propose an operator-splitting based numeri al algorithm to solve the transient KZK equation numeri ally in the time domain, using a ir ular or spheri ally fo used tradu er. Cleveland [18℄ expands the appli ation of this model by in luding the ee t of relaxation. An approa h for simulating nonlinear ontinuous wave (CW) pressure elds with the KZK equation is presented by Berntsen [19℄. The validation of the KZK equation for axisymmetri ultrasound beams is evaluated by Soneson [64℄. Curra and Filonenko [65, 66℄ numeri ally al ulate solutions for the CW KZK equation to model heat deposition in biologi al tissue. Huijssen [67℄ ompares the Iterative 71 Nonlinear Contrast Sour e (INCS) method with the KZK model for the omputation of nonlinear, wide-angle, pulsed a ousti elds. Attempts to al ulate solutions to the KZK equation with nite dieren e simulations, however, are limited by the implementation of the boundary onditions. Sin e there is no absorbing boundary layer implemented in the KZK models of Berntsen [19℄ or Lee [17℄, large grids are needed for these al ulations. In ele tromagneti s, this problem is addressed with absorbing boundary onditions or perfe tly mat hed layers. For example, Mur [68℄ proposes an a urate absorbing boundary ondition for both 2D and 3D time-domain ele tromagneti eld equations. Berenger [69℄ introdu es a perfe tly mat hed layer (PML) for nite dieren e time domain (FDTD) simulations of ele tromagneti waves, whi h is validated by Katz [70℄ in the 2D ase and then extended to the 3D ase. The PML is applied to a ousti models with dierent oordinate systems by Yuan and Liu [71, 72, 73℄. Abarbanel and Hu [74, 75℄ have developed PML equations for 2D linearized Euler equations. Sheu [76℄ applies a PML to a photoa ousti s model in axisymmetri ylindri al oordinates. Based on rst-order auxiliary dierential equations, Ma [77℄ proposes an unsplit PML for a se ond-order a ousti wave equation in 3D Cartesian oordinates. Ehrli h [78℄ ombines the a ousti wave propagator and a PML to model a ousti propagation in the time domain. Pinton [79℄ implements a PML for a 3D nonlinear attenuating full-wave model in the time domain. Resear h on 3D a ousti s attering problems has also been performed in both the time and frequen y domains by Kaltenba her [80℄, Chen [81℄, Alles [82℄, and Katsibas [83℄. Duru [84℄ also uses a PML in a 2D s alar wave equation for heterogeneous and layered media. In the following se tion, a perfe tly mat hed layer implemented through terms from the power law wave equation with y=0 and y=2 is derived for the transient and CW KZK equations. Arti ial attenuation is introdu ed to redu e ree tions from the radial boundary. Numeri al validations are demonstrated for both linear and nonlinear media. With these new PML implementations, the size of the omputational grid is redu ed, whi h a elerates numeri al solutions of the KZK equation. 72 4.2 Theory 4.2.1 The 2D wave equation with power law attenuation The linear lossless wave equation in a two-dimensional ylindri al oordinate system is expressed as ∂2p = c20 ∂t2  ∂2p 1 ∂ + ∂z 2 r ∂r  ∂p r ∂r  . (4.1) When the ee t of attenuation is introdu ed through the power law wave equation, the wave equation in 2D ylindri al oordinates be omes  ∂2p 1 ∂ + ∂z 2 r ∂r   ∂p 1 ∂2p 2α0 α02 ∂ y+1 p ∂ 2y p r − 2 2 − − = 0. ∂r c0 ∂t c0 cos (πy/2) ∂ty+1 cos2 (πy/2) ∂t2y For power law exponents y = 0 y = 2, and (4.2) the power law wave equation in Eq. 4.2 is expressed as ∂2p = c20 ∂t2  ∂2p 1 ∂ + ∂z 2 r ∂r   ∂p ∂p r − 2α0 c0 − α02 c20 p y = 0, ∂r ∂t (4.3) and ∂2p = c20 ∂t2  ∂2p 1 ∂ + ∂z 2 r ∂r  ∂p r ∂r  + 2α0 c0 4 ∂3p 2 2∂ p − α c 0 0 ∂t3 ∂t4 y = 2. (4.4) 4.2.2 Coordinate stret hing The analyti al expressions that des ribe perfe tly mat hed layers are often derived from an augmented model for the gradient that stret hes ea h oordinate a ording to ∇ = x̂ 1 ∂ 1 ∂ 1 ∂ + ŷ + ẑ , sx ∂x sy ∂y sz ∂z 73 (4.5) where sx , sy , and sz des ribe the ee ts of oordinate stret hing in all three dire tions. In the x- dire tion, the oordinate stret hing variable is dened as  sx = σ (x) 1+ jω  , (4.6) and similar expressions are dened in the y- and z- dire tions. The lossless Helmholtz equation in axisymmetri ylindri al oordinates is expressed as 1 ∂ r ∂r k= where the wavenumber is   ∂P ∂2P r + + k 2 P = 0, ∂r ∂z 2 ω and c0 P (4.7) is the pressure eld. Applying oordinate stret hing to Eq. 4.7 in both the radial and axial dire tions, ∂ 1 ∂  , → ∂r 1 + σrjω(r) ∂r 1 1 1  , → r 1 + σrjω(r) r ∂ 1 ∂  , → ∂z 1 + σzjω(z) ∂z (4.8) yields 1 whi h be omes  1+ σr (r) jω 2  ∂2P 1 ∂P + r ∂r ∂r 2 + 1+ ∂2P 2 2 2 + k P = 0, ∂z σz (z) 1 σz (z) jω after ea h term is multiplied by  1+ 1 ∂P r ∂r σz (z) jω 2  1+ (4.9) jω 2   2 + ∂∂rP2 + 1 + σrjω(r) 2  2  1 + σzjω(z) k 2 P = 0 + 1 + σrjω(r)  1+ 2   σr (r) jω 2 ∂2P ∂z 2 . Multiplying Eq. 4.10 by (4.10) (jω)2 and inverse Fourier transforming in time, the result is the lossy time-domain wave equation ∂ ∂t + σz (z) − ∂ ∂t 2  1 ∂p r ∂r + 2 + σr (r) ∂2p ∂r 2 ∂ ∂t  + ∂ ∂t + σz (z) + σr (r) 2 1 p c20 2 ∂2p ∂z 2 = 0, whi h attenuates wave propagation in the radial and axial dire tions. 74 (4.11) 4.2.3 The KZK equation The Westervelt equation is a full wave nonlinear model that ombines the ee ts of dira tion, attenuation, and nonlinearity β ∂ 2 p2 1 ∂ 2 p 2α0 ∂ 3 p + + = 0. ∇ p− 2 2 c0 ∂ t c0 ∂t3 ρ0 c40 ∂t2 2 (4.12) After applying a paraboli approximation and a hange of variables, Eq. 4.12 redu es to the KZK equation, whi h is dened in ylindri al oordinates (r, z) as ∂2p c0 ∂3p β ∂ 2 p2 = ∇2⊥ p + α0 3 + , ∂z∂τ 2 ∂τ 2ρ0 c30 ∂τ 2 where p is the pressure, onstant, and β c0 is the sound speed, ρ0 is the density, (4.13) α0 is the attenuation is the nonlinearity parameter. The propagation dire tion of the sound beam is aligned with the z-axis, and ∇2⊥ = ∂2 ∂r 2 + 1 ∂ is the radial omponent of the Lapla ian r ∂r operator in ylindri al oordinates. A detailed derivation of the KZK wave equation is given in Appendix A. Due to axial symmetry, the spatial variables for the KZK equation are (r, z) in the ylindri al oordinate system. On the right-hand side of Eq. 4.13, the three terms from left to right represent the ee ts of dira tion, attenuation, and nonlinearity. Similar to the Westervelt equation, there is no analyti al solution for the KZK equation. Thus, the KZK model is evaluated numeri ally. Although the KZK equation is only valid in the far eld of the paraxial region, this model is still ommonly applied to simulations of medi al ultrasound due to the omputational advantages of the paraboli approximation. As developed by Lee and Hamilton [17℄ and implemented in the KZK Texas ode [17℄, the nite dieren e method numeri ally solves the transient KZK equation with operator splitting to separately a ount for these three ee ts at ea h step. 75 4.2.4 PML derivation for the KZK equation Stret hed oordinates for the nonlinear lossy KZK equation Sin e waves propagate in the positive z dire tion but not in the negative z dire tion with the KZK model, there is no need for a PML in the axial dire tion. Eq. 4.13 in the frequen y domain is given by jω c0 β ∂P (jω)2 P 2 = ∇2⊥ P + α0 (jω)3 P + ∂z 2 2ρ0 c30 (4.14) and when oordinate stret hing is applied to the KZK equation for a PML in the radial dire tion, this yields  2  2  σ(r) c0 ∂ 1 ∂ 1 + P+ + = jω ∂P 2 ∂z jω 2 ∂r r ∂r 2 2   β α0 (jω)3 P + 1 + σ(r) 1 + σ(r) (jω)2 P 2 . jω jω 2ρ0 c3 (4.15) 0 After expanding and inverse Fourier transforming with respe t to time, the result is ∂2p ∂z∂τ 3 ∂p + 2σ (r) ∂z + σ (r)2 2 ´τ ∂p dτ ′ −∞ ∂z 2 ∂ p ∂ p +α0 ∂τ 3 + 2σ (r) α0 ∂τ 2 + α0 σ (r) ∂p ∂τ + = c0 2 β ∂ 2 p2 2ρ0 c30 ∂τ 2  ∂2 ∂r 2 + + 1 ∂ r ∂r σ(r)β ∂p2 ρ0 c30 ∂τ  + p σ(r)2 β 2 p , 2ρ0 c30 (4.16) whi h is a hanllenging to evaluate with nite dieren e al ulations, in part due to the integral on the left-hand side. Stret hed oordinates for the lossless linear KZK equation When only the ontributions from the dira tion terms are onsidered, the lossless linear KZK equation is expressed in the frequen y domain as jω ∂P c0 = ∇2⊥ P. ∂z 2 76 (4.17) If oordinate stret hing for a PML is applied in the radial dire tion, Eq. 4.17 then be omes ∂P jω ∂z  σ (r) 1+ jω 2 c0 = 2  ∂2 1 ∂ + ∂r 2 r ∂r  P. (4.18) After expanding and inverse Fourier transforming, ∂2p ∂p + 2σ (r) + σ (r)2 ∂z∂τ ∂z ˆ τ −∞ ∂p ′ c0 dτ = ∂z 2  ∂2 1 ∂ + 2 ∂r r ∂r  p, (4.19) whi h also ontains an integral on the left-hand side. Combining the power law equation with the KZK equation When Eq. 4.6 is applied to the 1D wave equation, ∂2p 1 ∂ 2 p 2σ (x) ∂p σ (x)2 p, = + + ∂x2 c20 ∂t2 c20 ∂t c20 whi h, after substituting with y = 0. αP M L = c0 σ (x), (4.20) is re ognized as the 1D power law wave equation This suggests that there is another possible approa h for implementing a PML that utilizes the power law equation. The 3D power law wave equation for the y=0 ase is expressed as ∇2 p = 1 ∂ 2 p 2αP M L ∂p + + αP2 M L p. c20 ∂t2 c0 ∂t (4.21) After applying a paraboli approximation and a hange of variables in an axisymmetri ylindri al oordinate system, Eq. 4.21 be omes ∂2p c0 ∂p c0 2 = ∇2⊥ p − αP M L − αP M L p. ∂z∂τ 2 ∂τ 2 Similarly, the 3D power law wave equation for the ∇2 p = y=2 ase is expressed as 1 ∂ 2 p 2αP M L ∂ 3 p ∂4p 2 − + α , P ML c20 ∂t2 c0 ∂t3 ∂t4 77 (4.22) (4.23) and the orresponding expression that is obtained from a paraboli approximation and a hange of variables in an axisymmetri ylindri al oordinate system is c0 ∂ 3 p c0 ∂2p ∂4p = ∇2⊥ p + αP M L 3 − αP2 M L 4 . ∂z∂τ 2 ∂τ 2 ∂τ (4.24) Eqs. 4.22 and 4.24 provide two related yet dierent approa hes for dening a perfe tly mat hed layer for the transient KZK equation, where the distribution of αP M L values an vary spatially. Using the power law wave equation with PML y=0 (two terms) for the Sin e the KZK equation is a one-way wave equation in the axial dire tion, a PML is only needed in the radial dire tion. This suggests that stret hed oordinates are not required for the derivation of the PML and that the power law wave equation with y =0 or y =2 in the radial dire tion provide ee tive expressions for PMLs. This indi ates that the stret hed oordinate system is not entral to the derivation of an ee tive PML. The power law wave equation is all that is required, ombined with the on ept introdu ed by Berenger [69℄ that the PML should be several ells thi k with a slowly varying lossy impedan e that attenuates the in ident wave with minimal ree tions. By applying 1 ∂ c0 ∂t → 1 ∂ c0 ∂t + αP M L to the wave equation in the time domain and utilizing the paraboli approximation in retarded time, one su h PML for the transient KZK equation is obtained from ∂2p c0 ∂p c0 2 = ∇2⊥ p − αP M L − αP M L p ∂z∂τ 2 ∂τ 2 for y = 0. (4.25) For numeri al al ulations, the expression in Eq. 4.25 is solved with the Thomas algorithm. 78 Using the Telegrapher's equation y = 0 (one term) for the PML For the y = 0 ase, the se ond term c0 2 α p an be negle ted when the value of c20 αP2 M L 2 P ML is small. Thus, another expression that des ribes a PML for the transient KZK is given by ∂2p c0 ∂p = ∇2⊥ p − αP M L , ∂z∂τ 2 ∂τ (4.26) whi h is the retarded time paraboli approximation to the 3D Telegrapher's equation. Using the power law wave equation with PML y=2 (two terms) for the The key to an ee tive absorbing boundary layer appears to be independent of the stret hed oordinate system and is mu h more strongly inuen ed by other fa tors su h as the slowly in reasing attenuation that minimizes ree tions at the interfa e between any two adja ent ells in the omputational grid. This motivates the onstru tion of a third PML based on the retarded time paraboli approximation to the power law wave equation with y = 2, whi h is given by ∂2p c0 ∂ 3 p c0 ∂4p = ∇2⊥ p + αP M L 3 − αP2 M L 4 . ∂z∂τ 2 ∂τ 2 ∂τ (4.27) Numeri ally, the rst term on the right hand side ombined with the original attenuation term is solved with the Thomas algorithm. The addition of the se ond attenuation term requires a penta-diagonal matrix algorithm. Using the Bla ksto k wave equation with PML Similar to the y = 2, y=0 y =2 (one term) for the ase, the se ond attenuation term is negligible for small αP M L when whi h enables a fourth PML that is very losely related to the KZK equation. This 79 PML evaluates a straightforward modi ation of the KZK equation, ∂2p c0 ∂3p = ∇2⊥ p + αP M L 3 , ∂z∂τ 2 ∂τ (4.28) whi h is the retarded time paraboli approximation to the Bla ksto k wave equation. 4.3 Methods 4.3.1 Error al ulations To validate these PMLs for the KZK equation, omparisons between simulations with and without PMLs that in rease the radial boundary limit to avoid ree tions are evaluated. The formula for al ulating the dieren e between these two results is D(r, z) = max|p(r, z) − pref (r, z)| , max|pref | (4.29) where the denominator is the maximum value of the referen e pressure evaluated at all spatial and temporal points, and the numerator is the maximum value of the dieren e between the referen e and simulation results evaluated at one spatial point for all time points. dieren e D(r, z) The is dependent on both radial and axial oordinates. 4.3.2 Finite dieren e al ulations with the KZK equation Numeri al al ulations with the KZK equation are often evaluated with the nite dieren e method. The pressure eld is rst dis retized in both the radial and axial dire tions, after whi h nite dieren e al ulations are evaluated layer-by-layer in the axial dire tion. Within ea h layer, the ee ts of dira tion, attenuation, and nonlinearity are al ulated separately through operator splitting. For these al ulations, two types of nite dieren e al ulations are evaluated. For the rst several iterations, the impli it ba kward Euler nite dieren e method (IBFD) is used, after whi h the Crank-Ni olson nite dieren e method (CNFD) is applied. Both of these methods are numeri ally stable, and the CNFD method 80 results in a smaller lo al trun ation error than the IBFD method with the same step size. However, the CNFD method is sensitive to os illations, whi h means that the omputed result an ontain spurious os illations, espe ially in the region lose to the edge of the transdu er where there is a dis ontinuity in pressure amplitude. Thus, the IBFD method is applied as a low pass lter in this region. For these reasons, these two nite dieren e methods are ombined for the numeri al evaluations of the KZK equation. For the transient KZK equation, three ee ts are al ulated separately within one spatial step using operator splitting. Dira tion ee ts are modeled by ∂p = ∂z ˆ τ −∞ c0 2  ∂2 1 ∂ + 2 ∂r r ∂r  pdτ ′ , (4.30) whi h is obtained after the nonlinear and loss terms are negle ted and the remaining terms are integrated on both sides. The indi ies of the nite dieren e al ulation in the temporal, radial, and axial dire tions are i, j , and k . Thus, nite dieren e approximations for Eq. 4.30 are dened as pij,k+1 − pij,k ∂p , → ∂z (△z)k pij+1,k+1 − pij−1,k+1 1 ∂p → , r ∂r 2j (△r)2 pij+1,k+1 − 2pij,k+1 + pij−1,k+1 ∂2p , → ∂r 2 (△r)2 # " i−1 ˆ τ X 1 f (τ ′ )dτ ′ → (△τ ) fm + (f0 + fi ) . 2 τ min m=1 (4.31) (4.32) (4.33) When only the ee t of attenuation is onsidered, Eq. 4.13 be omes ∂p ∂2p = α0 2 . ∂z ∂τ (4.34) Similarly, the nite dieren e approximations for Eq. 4.34 are dened as pij,k+1 − pij,k ∂p → , ∂z (△z)k i−1 i pi+1 ∂2p j,k+1 − 2pj,k+1 + pj,k+1 . → ∂τ 2 (△τ )2 81 (4.35) When only the nonlinear ee t is onsidered, Eq. 4.13 redu es to ∂p β ∂p p . = ∂z ρ0 c30 ∂τ (4.36) As indi ated by Lee and Hamilton [17℄, the solution to Eq. 4.36 is given by pij,k+1 = pij,k 1−   , 3 i pi+1 j,k − pj,k /△τ β (△z)k / (ρ0 c0 ) 1−   , pij,k − pj,k /△τ β (△z)k / (ρ0 c30 )  pij,k ≥ 0 (4.37) pij,k ≤ 0. (4.38) and pij,k+1 = pij,k  i−1 For ontinuous wave propagation, the a ousti pressure is often des ribed as series expansion of dierent harmoni s p (τ, r, z) = N max X Cn (r, z) e−jn2πf0 τ , (4.39) n=−Nmax where f0 is the fundamental frequen y, Nmax is the total number of harmoni s, and Cn (r, z) is the omplex amplitude of the n-th harmoni . When Eq. 4.39 is ombined with the transient KZK equation, the amplitudes are expressed as dCn (r,z) dz   ∂ 2 Cn (r,z) + 1r ∂Cn∂r(r,z) − α0 (2πnf0 )2 Cn ∂r 2 NP max 0β − jn2πf Cm (r, z) Cn−m (r, z) , 3 2ρ0 c0 m=−Nmax c0 = j 4πnf 0 (r, z) (4.40) n = ±1, ±2, ..., ±Nmax . In Eq. 4.40, all of the harmoni s intera t through nonlinear mixing. The derivatives in Eq. 4.40 are also dened in Eq. 4.31, as for the transient KZK equation. 4.3.3 Muir's method To establish the a ura y of the numeri al results al ulated with linear nite dieren e implementations of the KZK wave equation, Muir's method is utilized as a referen e. Other 82 methods for validation are given in Appendix B. Muir's method, whi h is valid for small aperture angle and large ka, is ee tive for al ulating the linear pressure eld generated by a spheri ally fo used ontinuous-wave sour e. Muir's formula for the pressure distribution in ylindri al oordinates predi ted by the linear lossless KZK equation is given by 2 r −jkp0 exp(jkz + jk 2z ) P (r, z) = z ˆa " jk (r ′ )2 exp 2 0  1 1 − R z #   rr ′ r ′ dr ′ J0 k z (z 6= R) , (4.41)   2 ka 2   1 1 (z 6= R) , − z R   r 2 2J1 (kar/R) ka2 p0 exp jkz + jk , P (r, R) = −j 2R 2R kar/R P (0, z) = p0 Rexp(jkz) (R − z) 1 − exp j P (0, R) = −j where R p0 is the pressure at the sour e, is the radius of urvature, and order 0 and 1, ka2 p0 exp (jkR) , 2R k = ω/c J0 (·) and is the wavenumber, J1 (·) (4.42) (4.43) (4.44) a is the aperture radius, are Bessel fun tions of the rst kind of respe tively. Muir's method evaluates a single integral when al ulating the o-axis pressure. For transient al ulations, a Fourier transform operation is required before applying Muir's method. When al ulating the transient pressure, the input pulse on the surfa e of the transdu er is expressed as p0 (t) = N P Pn ejnω0 t , where Pn is the pressure for frequen y sample n=1 n, ω0 is the fundamental frequen y, and frequen y sample n with wavenumber N kn , is the number of frequen y samples. Thus, for Muir's method in Eqs. 4.41-4.44 al ulates the pressure distribution in 2D spa e for ea h frequen y omponent with Pn (r, z) = P (r, z, kn ) 83 n = 1, 2, 3..N. (4.45) Then, after evaluating the inverse Fourier transform in time, the solution in the time domain is p(r, z, t) = N X Pn (r, z)ejnω0 t . (4.46) n=1 Compared to nite dieren e KZK al ulations, Muir's method is more time- onsuming for transient al ulations be ause of the large number of frequen ies that are ne essary to re onstru t the waveforms at ea h spatial point. However, Muir's method provides an a urate referen e for validating solutions to the linear lossless KZK equation, whi h is useful for debugging nite dieren e al ulations. 4.4 Results These simulations were performed on a ompute server (Dual Intel Xeon E5-2670  2.5 GHz) with 256 GB RAM. The KZK simulation evaluates a nite dieren e ode written in C++/Mex that runs on 64-bit MATLAB R2017a. Simulations for both linear and nonlinear media are evaluated. The transient input pressure, whi h is a one- y le Gaussian weighted sine wave, is generated by a spheri ally fo used transdu er. surfa e of the transdu er is urvature is P0 = 1.5 MPa, the aperture radius is a = 1.5 m, the radius of R = 6 m, the density is ρ = 1000 kg/m3 , and the sound speed is c0 = 1500 m/s. The enter frequen y of the input pressure is frequen y is The input pressure on the λ = 0.15 f = 1MHz, m. The sampling frequen y is and the wavelength at the enter fs = 200 MHz. The nite dieren e KZK al ulation that is employed as a referen e utilizes a radial boundary at and the spatial step size is at rmax = 3 m. ends at λ/40. m, The KZK simulation with a PML utilizes a radial boundary For the KZK simulation with the PML, the PML starts at rmax = 3 rmax = 9 m. The thi kness of the layer is therefore equal to 0.75 r = 2.25 m and m. 4.4.1 KZK simulations for a linear lossless medium KZK simulations are rst performed without attenuation or nonlinearity. Fig. 4.1 ompares two on-axis waveforms at z=6 m, where one is produ ed by KZK nite dieren e 84 Pressure (Pa) ×106 6 4 2 0 -2 -4 -6 -8 -10 -12 38 r = 0, z = 6cm KZK without PML (rmax = 6a) Muir's method 40 42 44 46 48 50 52 54 time (µs) Figure 4.1: Comparison of simulated on-axis waveforms obtained from nite dieren e KZK al ulations (solid line) and Muir's method (dashed line) evaluated in a linear lossless medium at z=6 m. al ulations (red solid line) and the other is generated by Muir's method (blue dashed line). Fig. 4.1 demonstrates that the results obtained with these two methods mat h losely at all temporal points in this lo ation. Also, the dire t wave and the edge wave have merged in this lo ation. Fig. 4.2(a) ontains the on-axis waveforms evaluated at KZK al ulations without a PML for a radial boundary at without a PML for a radial boundary at for t < 48 µs. Near t = 50 µs, rmax = 2a. z=6 rmax = 6a and the KZK simulation These two waveforms mat h losely the ree tion from the radial boundary at whi h indi ates that either the boundary at rmax = 2a is two dierent PMLs that start at (y =0 rmax = 6a. r = 2.25 rmax = 2a arrives, too lose or that a PML is needed. Fig. 4.2(b) ompares two other on-axis waveforms evaluated at radial boundary at m using nite dieren e z = 6 m al ulated with m to the result without a PML that denes the The blue solid line des ribes the result al ulated with a PML with one term) and the green solid line gives the result al ulated with another PML 85 6 ×10 6 r = 0, z = 6cm KZK without PML (r max = 6a) 4 KZK without PML (r max = 2a) Pressure (Pa) 2 0 -2 -4 -6 -8 -10 -12 38 40 42 44 46 48 50 52 54 time (µ s) (a) 6 ×10 6 r = 0, z = 6cm KZK without PML (r max = 6a) 4 KZK with PML y = 0 (r max = 2a) Pressure (Pa) 2 KZK with PML y = 2 (r max = 2a) 0 -2 -4 -6 -8 -10 -12 38 40 42 44 46 48 50 52 54 time (µ s) (b) Figure 4.2: Comparison between on-axis waveforms generated by nite dieren e KZK al ulations in a linear lossless medium at z = 6 m with and without a PML using dierent radial boundaries. (a) KZK simulation without a PML that denes a radial boundary at rmax = 2a (bla k solid line) and at rmax = 6a (red dashed line). (b) KZK simulation without a PML that denes a radial boundary at rmax = 6a (red dashed line), with a y = 0 single term PML that denes a radial boundary at rmax = 2a (blue solid line), and with a y = 2 single term PML that denes a radial boundary at rmax = 2a (green solid line). 86 (a) (b) ( ) (d) Figure 4.3: Simulated 2D pressure eld and dieren es between KZK al ulations without and with PMLs, where the radial boundaries are lo ated at in a linear lossless medium. rmax = 9 m and at rmax = 3 m (a) The peak pressure distribution for the KZK simulation without a PML that denes a radial boundary at rmax = 9 m. (b) The dieren e between rmax = 9 m rmax = 3 m. the KZK simulation without a PML that denes a radial boundary at and the KZK simulation without a PML that denes a radial boundary at ( ) The dieren e between the KZK simulation without a PML that denes a radial boundary at rmax = 9 m and the KZK simulation with a boundary at rmax = 3 y=0 single term PML that denes a radial m. (d) The dieren e between the KZK simulation without a PML that denes a radial boundary at rmax = 9 m and the KZK simulation with a term PML that denes a radial boundary at rmax = 3 87 m. y=2 single (y =2 with one term). Fig. 4.2(b) shows that the ree tion from the boundary is removed by ea h of these PMLs. Fig. 4.3(a) shows the 2D peak pressure distribution al ulated with nite dieren e KZK al ulations without a PML that denes a radial boundary at rmax = 9 m (rmax = 6a), whi h is su iently large so that radial ree tions are avoided in most lo ations. on-axis fo al peak is lo ated at about 12 z=6 The m. The maximum pressure value is equal to MPa. Fig. 4.3(b) shows the 2D dieren e between the nite dieren e KZK al ulation without a PML for a radial boundary dened at rmax = 9 m and the nite dieren e KZK al ulation without a PML for a radial boundary dened at of the dieren e, whi h is about 10%, rmax = 3 m. The peak value is lo ated on-axis in the far eld region. O-axis, the dieren e is mu h smaller. Fig. 4.3( ) shows the 2D dieren e between the nite dieren e KZK al ulation without a PML that denes a radial boundary at dieren e KZK al ulation with a rmax = 3 m. y=0 rmax = 9 m and the nite single term PML that denes a radial boundary at The peak value of the dieren e, whi h is now only about 0.3%, is again lo ated on-axis in the far eld region. There is also some dieren e o-axis in the far eld region in Fig. 4.3( ) due to ree tions that arrive later and/or a small amount of mismat h in the PML. Fig. 4.3(d) shows the 2D dieren e between the nite dieren e KZK al ulation without a PML for a radial boundary dened at with a single term y=2 rmax = 9 m and the nite dieren e PML that denes a radial boundary dened at peak value of the dieren e, whi h is about 0.5% KZK al ulation rmax = 3 m. The for this PML, is also lo ated on-axis in the far eld region. Also, some small dieren es appear on-axis and o-axis in the far eld. Comparisons between Fig. 4.3( ) and Fig. 4.3(d) indi ate slightly better performan e for the y=0 single term PML. 4.4.2 KZK simulations for a nonlinear lossy medium Simulations are also performed with the attenuation and nonlinearity values for water, whi h are α0 = 2.2 × 10−3 2 dB/ m/MHz and between the on-axis waveforms evaluated at β = 3.5. z = 6 88 Fig. 4.4(a) shows a omparison m obtained from the nite dieren e 12 ×10 6 r = 0, z = 6cm KZK without PML (r max = 6a) 10 KZK without PML (r max = 2a) Pressure (Pa) 8 6 4 2 0 -2 -4 -6 -8 38 40 42 44 46 48 50 52 54 time (µ s) 12 ×10 6 r = 0, z = 6cm KZK without PML (r max = 6a) 10 KZK with PML y = 0 (r max = 2a) Pressure (Pa) 8 KZK with PML y = 2 (r max = 2a) 6 4 2 0 -2 -4 -6 -8 38 40 42 44 46 48 50 52 54 time (µ s) Figure 4.4: Comparison between on-axis waveforms generated by nite dieren e KZK al ulations in a nonlinear medium at radial boundaries. m with and without a PML for dierent 2 −3 The attenuation parameter is α = 2.2 × 10 dB/ m/MHz , and the nonlinearity parameter is boundary at z = 6 rmax = 2a β = 3.5. (a) KZK simulation without a PML that denes a radial (bla k solid line) and at rmax = 6a (b) KZK rmax = 6a (red dashed line), y = 0 single term PML that denes a radial boundary at rmax = 2a (blue solid line), with a y = 2 single term PML that denes a radial boundary at rmax = 2a (green solid simulation without a PML that denes a radial boundary at with a and (red dashed line). line). 89 rmax = 6a KZK simulation without a PML that denes a radial boundary at (red dashed line) and the nite dieren e KZK simulations without a PML that denes a radial boundary at rmax = 2a. When nonlinearity is in luded, the waveforms are tilted where a sho kwave is formed in the fo al zone. The small dieren e near t = 50 µs in the nite dieren e KZK simulation without the PML is again aused by the ree tion at the radial boundary. Fig. 4.4(b) shows a omparison between the on-axis waveforms at z =6 m for the nite dieren e KZK simulation without a PML that denes a radial boundary at rmax = 6a and the nite dieren e KZK simulation with two dierent PMLs that dene a radial boundary at rmax = 2a. The blue solid line shows the result with a green solid line des ribes the result with a y=0 single term PML and the y = 2 single term PML. Fig. 4.4(b) indi ates that the ree tion from the boundary is removed by ea h of these PMLs. Fig. 4.5 shows the simulated 2D pressure eld and the dieren es between KZK al ulations without and with PMLs in a nonlinear medium. Fig. 4.5(a) shows the 2D peak pressure distribution al ulated with the nite dieren e implementation of the KZK equation without a PML that denes a radial boundary at z=6 rmax = 9 m. The fo al peak is lo ated at about m on axis, and the maximum pressure value is equal to 12 MPa. Fig. 4.5(b) shows the 2D dieren e between the nite dieren e KZK al ulation without a PML that denes a radial boundary at rmax = 9 m and the nite dieren e KZK al ulation without a PML that denes a radial boundary at about 10%, rmax = 3 m. The peak value of the dieren e, whi h is is lo ated on-axis in the far eld region. Fig. 4.5( ) shows the 2D dieren e between the nite dieren e KZK al ulation without a PML that denes a radial boundary at rmax = 9 m and the nite dieren e KZK al ulation with the that denes a radial boundary at is about 0.3%, rmax = 3 m. y=0 single term PML The peak value of the dieren e, whi h is again lo ated on-axis in the far region. There is also some dieren e in the far o-axis region in Fig. 4.5( ). Fig. 4.5(d) shows the 2D dieren e between the nite dieren e KZK al ulation without a PML that denes a radial boundary at and the nite dieren e KZK al ulation with the 90 y = 2 rmax = 9 m single term PML that denes a (a) (b) ( ) (d) Figure 4.5: Simulated 2D pressure eld and dieren es between KZK al ulations without rmax = 9 m and at rmax = 3 m α = 2.2 × 10−3 dB/ m/MHz2 and the and with PMLs, where the radial boundaries are lo ated at in a nonlinear medium. The attenuation parameter is nonlinearity parameter is β = 3.5. (a) The peak pressure distribution for the KZK simulation without a PML that denes a radial boundary at rmax = 9 m. (b) The dieren e between rmax = 9 m rmax = 3 m. the KZK simulation without a PML that denes a radial boundary at and the KZK simulation without a PML that denes a radial boundary at ( ) The dieren e between the KZK simulation without a PML that denes a radial boundary at rmax = 9 m and the KZK simulation with a boundary at rmax = 3 y=0 single term PML that denes a radial m. (d) The dieren e between the KZK simulation without a PML that denes a radial boundary at rmax = 9 m and the KZK simulation with a term PML that denes a radial boundary at rmax = 3 91 m. y=2 single radial boundary at rmax = 3 m. The peak value of the dieren e, whi h is about 0.5% for this PML, is also lo ated on-axis in the far eld region. Again, a smaller dieren e is observed in Fig. 4.5(d) than in Fig. 4.5( ), whi h indi ates slightly better performan e for y=0 the single term PML. 4.5 Dis ussion 4.5.1 Computation time Figs. 4.3 and 4.5 indi ate that, for both linear and nonlinear KZK simulations, the PMLs y = 0 with and y = 2 are ee tive in suppressing ree tions from the radial boundary, espe ially in the fo al zone. The hoi e of when αP M L αP M L balan es the ee t of impedan e mismat h is large versus the ee t of boundary ree tions when αP M L is too small. For ea h type of PML, our experien e is that there is no signi ant dieren e in the PML when only one or two terms are onsidered in the power law wave equation. Thus, to further a elerate these KZK simulations, only one term is in luded in ea h simulation for the and y=2 y=0 PMLs. The omputation time for the KZK nite dieren e al ulation without a PML that denes a radial boundary at simulations using a y=0 the omputation time is rmax = 9 m shown in Fig. 4.3 is 2807s. single term PML that denes a radial boundary at 913s. denes a radial boundary at For KZK simulations using a rmax = 3 y =2 m, the omputation time is For KZK rmax = 3 m, single term PML that 1028s. 4.5.2 Continuous wave (CW) KZK al ulations The same approa h for dening a PML is also appli able to ontinuous-wave KZK simulations, as des ribed in Appendix A. For the ontinuous wave KZK equation, a PML is only ne essary in the radial dire tion. The strength of the PML, as implemented here, in reases proportionally to the radial dire tion ubed. The ee tiveness of the 2 y = 0 and y = single term PMLs is also demonstrated for CW al ulations with the spheri ally-fo used transdu er evaluated in se tion 4.4. The input pressure on the surfa e of the transdu er 92 5 ×106 Muir's method KZK without a PML (rmax = 7a) Pressure (Pa) 4 3 2 1 0 0 5 10 15 20 axial distance (cm) (a) 5 ×106 Muir's method KZK with y = 0 PML (rmax = 2a) Pressure (Pa) 4 3 2 1 0 0 5 10 15 20 axial distance (cm) (b) Figure 4.6: On-axis omparisons between ontinuous wave nite dieren e KZK simulations with the results al ulated with Muir's method evaluated in a linear lossless medium. (a) The result obtained with Muir's method (red solid line) and the nite dieren e KZK simulation results without a PML that denes a radial boundary at rmax = 10.5 m (blue dashed line). (b) The result with Muir's method (red solid line) and the nite dieren e KZK simulation results with a y=0 PML that denes a radial boundary at 93 rmax = 3 m (blue dashed line). is P0 = 0.5 MPa. The density is The ex itation frequen y is f =1 ρ = 1000 3 kg/m , and the sound speed is c = 1500 MHz, whi h orresponds to a wavelength of m/s. λ = 0.15 m. For the nite dieren e KZK simulation without a PML, the radial boundary is lo ated at rmax = 10.5 m. PML starts at to 0.75 For the nite dieren e KZK simulation with a r = 2.25 m and ends at rmax = 3 y = 0 single term PML, the m, so the thi kness of the PML is equal m. The nite dieren e solution to the ontinuous-wave KZK equation is rst omputed in a linear lossless medium. Only the rst harmoni is omputed in the simulation, where the spatial step size is equal to λ/40 in both dire tions. Fig. 4.6 des ribes the on-axis results for the ontinuous wave KZK simulations with and without a PML, whi h are ompared to the results al ulated with Muir's method. As shown in Fig. 4.6(a), the on-axis pressure waveform obtained from the nite dieren e KZK simulation without a PML losely mat hes the waveform omputed with Muir's method in the fo al zone; however, there is some dieren e in the far eld region, even for the large radial boundary that is dened at rmax = 10.5 m. Fig. 4.6(b) des ribes the on-axis waveform obtained from the nite dieren e KZK simulation with a y = 0 PML that denes a radial boundary at rmax = 3 m, whi h losely mat hes the result obtained with Muir's method both in the fo al zone and in the far eld region. This indi ates that the ree tion from the boundary is removed by the PML. Fig. 4.7 shows the ontinuous-wave 2D pressure distribution for a spheri ally-fo used transdu er with a = 1.5 m, R = 6 m, and f = 1 MHz al ulated in a linear lossless medium with Muir's method. In Fig. 4.7, the peak pressure in the fo al zone is approximately 4.3 MPa. Fig. 4.8(a) des ribes the dieren e between the nite dieren e KZK numeri al al ulation without a PML that denes a radial boundary at obtained with Muir's method. results The ree tion from the radial boundary is learly evident in Fig. 4.8(a), whi h starts at the edge near 0.5 rmax = 10.5 m and the z = 16.5 m. The peak dieren e is about MPa in the on-axis far eld region. Fig. 4.8(b) shows that, after introdu ing a 94 y =0 Figure 4.7: The ontinuous-wave 2D pressure distribution for a spheri ally-fo used transdu er with a = 1.5 m, R=6 m, and f =1 MHz al ulated in a linear lossless medium with Muir's method. PML that denes a radial boundary at rmax = 3 m, the on-axis dieren e is mu h smaller than the on-axis dieren e without a PML in Fig. 4.8(a). The largest dieren e that o urs in Fig. 4.8(b) is observed where the PML is applied. Simulations are then evaluated for the same transdu er geometry using the attenuation and nonlinearity values of water, whi h are respe tively. α0 = 2.2 × 10−3 The number of harmoni s omputed in this simulation is the spatial step size is λ/N harm /40. 2 dB/ m/MHz and β = 3.5, N harm = 50, and Fig. 4.9(a) shows the rst four harmoni s of the on-axis nite dieren e simulation results for the ontinuous wave KZK equation evaluated in water. This simulation denes a radial boundary at rmax = 10.5 m without a PML. In Fig. 4.9(a), there is a very strong os illation in the fundamental due to ree tions from the boundary. However, no su h os illations appear in the higher harmoni s. Fig. 4.9(b) shows the rst four harmoni s evaluated on-axis for the ontinuous-wave KZK equation in water with a PML that denes a radial boundary at rmax = 3 95 y=0 m. In Fig. 4.9(b), no os illations appear (a) (b) Figure 4.8: The 2D pressure dieren e between linear lossless nite dieren e KZK numeri al results for a spheri ally-fo used transdu er with a = 1.5 m, R= 6 m, and f = 1 MHz and the results for the same onguration evaluated with Muir's method. (a) The dieren e between the nite dieren e KZK simulation without a PML that denes a radial boundary at rmax = 10.5 m and the results obtained with Muir's method. (b) The dieren e between the nite dieren e KZK simulation with a rmax = 3 y =0 PML that denes a radial boundary at m and the results obtained with Muir's method. 96 5 ×106 KZK without a PML N=1 (r KZK without a PML N=2 (r Pressure (Pa) 4 KZK without a PML N=3 (r KZK without a PML N=4 (r max max max max = 7a) = 7a) = 7a) = 7a) 3 2 1 0 0 5 10 15 20 axial distance (cm) (a) 5 ×106 KZK with y = 0 PML N=1 (r KZK with y = 0 PML N=2 (r Pressure (Pa) 4 KZK with y = 0 PML N=3 (r KZK with y = 0 PML N=4 (r max max max max = 2a) = 2a) = 2a) = 2a) 3 2 1 0 0 5 10 15 20 axial distance (cm) (b) Figure 4.9: a = 1.5 The rst four harmoni s generated by a spheri ally-fo used transdu er with m, R = 6 m, and f = 1 MHz for on-axis nite dieren e simulations of the −3 ontinuous wave KZK equation in water. The attenuation parameter is α = 2.2 × 10 2 dB/ m/MHz , and the nonlinearity parameter is β = 3.5. (a) The nite dieren e KZK rmax = 10.5 m. (b) y = 0 PML that denes a radial boundary simulation results without a PML that denes a radial boundary at The nite dieren e KZK simulation results with a at rmax = 3 m. 97 Figure 4.10: The 2D pressure dieren e between the nite dieren e KZK numeri al results 2 −3 with and without a y = 0 PML in water with α0 = 2.2 × 10 dB/ m/MHz and β = 3.5 evaluated for rst four harmoni s produ ed by a spheri ally-fo used transdu er with 1.5 m, R=6 m, and f =1 a = MHz. in the fundamental or in any of the higher harmoni s. For the higher harmoni s, the nite dieren e solution to the ontinuous-wave KZK equation with or without a PML produ es exa tly the same on-axis result. Fig. 4.10 shows the rst four harmoni s of the dieren e between the nite dieren e KZK numeri al al ulation without a PML that denes a radial boundary at and the KZK numeri al al ulation with a rmax = 3 y = 0 rmax = 10.5 m PML that denes a radial boundary at m. Fig. 4.10(a) indi ates the main dieren e in the rst harmoni o urs lose to the entral axis, where the sour e of this dieren e is the ree tion from the radial boundary at rmax = 10.5 m. Figs. 4.10(b-d) indi ate that, for higher harmoni s, the dieren e is negligible sin e the ree tion from the radial boundary is mu h smaller for the higher frequen y omponents. 98 4.6 Con lusion A new perfe tly mat hed layer was implemented for simulations of nonlinear wave propagation based on the Khokhlov-Zabolotskaya-Kuznetsov equation. Instead of deriving a PML with stret hed oordinates, the power law wave equation is introdu ed as an alternative model for the attenuation that o urs within the PML. For ea h value of the power law exponent onsidered here, the two terms that are responsible for the attenuation are redu ed to a single term, whi h yields the Telegrapher's equation within the PML when vis ous wave equation when y = 2. y = 0 and the Bla ksto k Numeri al simulations are evaluated in both linear lossless and nonlinear lossy media for inputs generated by a spheri ally fo used transdu er. The simulation results are ompared to Muir's formula and to nite dieren e KZK solutions with a very large radial boundary. Comparisons show that the PMLs ee tively eliminate the ree tions from the radial boundary. In addition, the formulas for implementing PMLs are readily integrated into existing KZK simulation programs. With these new PMLs, ree tions from the radial boundary are eliminated, whi h enables a onsiderable redu tion in the omputation time for nite dieren e simulations of the KZK equation. 99 Chapter 5 Con lusion Chapter 2 numeri ally evaluates time-domain Green's fun tions for three time-fra tional wave equations, and the results are ompared at various distan es for water, breast, and liver. At larger distan es, the time-domain Green's fun tions for all three fra tional wave equations onverge to the same result in the water, breast, and liver models. The results also demonstrate that the Szabo and power law wave equations are non ausal and that the Caputo wave equation is ausal, where the distin tion between these is learly evident at distan es very lose to the sour e. However, beyond a ertain distan e, the non ausal ontributions are negligible for both the Szabo and power law wave equations. When these time-domain Green's fun tions are onvolved with a three- y le Hanning-weighted pulse, no non ausal behavior is observed in the time-domain results, and the FWHMs of the envelopes of the onvolution results are all approximately the same. In Chapter 3, improved approximations for the attenuation and phase velo ity are derived for the Chen-Holm and Treeby-Cox wave equations. Numeri al al ulations of the attenuation and phase velo ity for the Chen-Holm, Treeby-Cox, and power law wave equations in breast and liver are evaluated over a range of ultrasound frequen ies. New expressions for power series mat h the results obtained by numeri ally evaluating the dispersion relation more losely than previous approximations. The time-domain Green's fun tions for these three fra tional wave equations are al ulated at various distan es, and the amplitudes and FWHM values of the time-domain Green's fun tions are also evaluated. The results show some similarities and dieren es between these three fra tional wave equations. For instan e, 100 the attenuation terms in all three fra tional wave equations are very similar, while the phase velo ity for the Chen-Holm wave equation is nearly onstant. Causality is demonstrated analyti ally and numeri ally in the time domain for both the Chen-Holm and Treeby-Cox wave equations. At mu h larger distan es, the time-domain Green's fun tions for the Treeby-Cox and power law wave equations onverge to the same result while the time-domain Green's fun tion for the Chen-Holm wave equation learly diers from the other two. The Pantis method is introdu ed as an ee tive approa h for evaluating the highly os illatory improper integrals that arise in numeri al al ulations of the time-domain Green's fun tions for the Chen-Holm and Treeby-Cox spa e-fra tional wave equations. The Pantis method provides an a urate result when the number of Filon abs issas and the value of m are su iently large. Three- y le Hanning-weighted pulses with two dierent enter frequen ies are onvolved with the time-domain Green's fun tions for three fra tional wave equations. The onvolution results for the power law wave equation and the Treeby-Cox wave equation are very similar while the onvolution result for the Chen-Holm wave equation learly shows a time delay. The onvolution results also indi ate that there is more attenuation and waveform spreading in signals with higher enter frequen ies. In Chapter 4, a new PML, whi h is based on the power law wave equation with y=0 or y = 2, is implemented to a elerate nonlinear ultrasound simulations with the KZK equation. For ea h power law exponent, a single attenuation term is su ient to avoid radial ree tions. In addition, the nite dieren e stru ture of the KZK equation is also des ribed. Numeri al simulations for the transient and ontinuous-wave KZK equations are then evaluated for both linear lossless and nonlinear media, where the inputs are generated by a spheri ally fo used transdu er. These results are ompared to Muir's formula and to nite dieren e KZK al ulations with a large radial boundary. The omparisons indi ate that the new PML ee tively eliminates the ree tions from the radial boundary, whi h subsequently redu es the omputation time. 101 APPENDICES 102 APPENDIX A Derivation of the nonlinear wave equations The Westervelt equation A wave equation that des ribes nonlinear wave propagation is derived from three fundamental equations, namely the equation of motion, the ontinuity equation, and the equation of state. The equation of motion is given by ρ → Du − + ∇P = 0, Dt (A.1) ρ = ρ0 + ρa is the total density, P = P0 + p is the total pressure, u is the velo ity, and → − ∂ = ∂t + u · ∇ . A ousti quantities are usually very small ompared to the stati values. where D Dt The ontinuity equation is expressed as Dρ + ρ∇ · u = 0. Dt (A.2) The equation of state is given by P = P0 +  ∂P ∂ρ  1 ρa + 2 ρ0 ,s 103  ∂2P ∂ρ2  ρ2a + . . . .. ρ0 ,s (A.3) For all three onstitutive equations, all terms up to se ond order a ura y are retained with respe t to ρa , p, and u, whi h are su ient for most appli ations of nonlinear a ousti s in uids. A ordingly, Eqs. A.1-A.3 are rewritten as ρ0  − → → ∂u ∂u − + ∇p = −ρa − ρ0 u · ∇ u, ∂t ∂t (A.4) → − → − → − ∂ρa + ρ0 ∇ · u = −ρa ∇ · u − u · ∇ρa , ∂t ρa − where A = ρ0  ∂P ∂ρ  ρ0 ,s B = ρ20 and  ∂2P ∂ρ2 (A.5) p B p2 ≈ − , c20 2A ρ0 c40  ρ0 ,s (A.6) . The right hand sides of Eqs. A.4-A.6 ontain the se ond order terms. Combining Eqs. A.4-A.6 yields the Westervelt equation ∇2 p − where δ and β = 1 + B/2A δ ∂3p β ∂ 2 p2 1 ∂2p = − − , c20 ∂t2 c40 ∂t3 ρ0 c40 ∂t2 (A.7) are the attenuation and nonlinearity parameters. The KZK equation To obtain an approximate nonlinear wave equation that des ribes one way wave motion in the axial dire tion, let z represent the dire tion of propagation, where (x, y) indi ates the oordinates perpendi ular to the z axis. Assuming that, for a sour e with radius and z > 0.5ka2 are both satised, the ee ts of dira tion are O(ε̃2) a, ka ≫ 1 in ea h dire tion are s aled by dierent amounts a ording to p = p(x1 , y1 , z1 , t′ ), (x1 , y1, z1 ) = (ε̃1/2 x, ε̃1/2 y, ε̃z), t′ = t − z/c0 . (A.8) The Lapla ian that appears in the Westervelt equation is then rewritten as ∂2 ∂2 + ∇ = ε̃ ∂x21 ∂y12 2   + ε̃2 ∂2 2 ∂2p 1 ∂2 − ε̃ + . ∂z12 c0 ∂z1 ∂t′ c20 ∂t′2 104 (A.9) If only O(ε̃) terms are retained, the left side of the Westervelt equation be omes 1 ∂2p ∇ p − 2 2 = ε̃ c0 ∂t 2 Let ∇⊥ = ∂2 ∂x2 + ∂2 and repla e ∂y 2  ∂2 ∂2 + ∂x21 ∂y12 (x1 , y1 , z1 ) with  p − ε̃ (x, y, z) 2 ∂2p . c0 ∂z1 ∂t′ (A.10) to obtain the KZK (Khokhlov- Zabolotskaya-Kuznetsov) equation ∂2p c0 δ ∂3p β ∂ 2 p2 = ∇ p + + . ⊥ ∂z∂t′ 2 2c30 ∂t′3 2ρ0 c30 ∂t′2 (A.11) This yields a simplied model that in ludes the ee ts of dira tion, attenuation, and nonlinearity. When the sour e is a ir ular transdu er, the pressure eld is symmetri in the radial dire tion. In axisymmetri ylindri al oordinates, the KZK equation is given by ∂2p c0 = ′ ∂z∂t 2  ∂ 2 p 1 ∂p + ∂r 2 r ∂r  + δ ∂3p β ∂ 2 p2 + . 2c30 ∂t′3 2ρ0 c30 ∂t′2 (A.12) Burgers' equation Burgers' equation, whi h is a one dimensional nonlinear equation, an be derived dire tly from the Westervelt equation. After the operators in Eq. A.7 are fa tored, this yields  1 ∂ δ ∂2 βp ∂ ∂ − + + 2 ∂z c0 ∂t 2c0 ∂t ρ0 c30 ∂t In Eq. A.13,  ∂ 1 ∂ δ ∂2 βp ∂ + − − 2 ∂z c0 ∂t 2c0 ∂t ρ0 c30 ∂t  p = 0. (A.13) βp ∂ δ ∂2 and are higher order terms, and the produ t of these are dis arded 2c0 ∂t2 ρ0 c30 ∂t when transforming Eq. A.13 ba k to Eq. A.7. Assuming one way approximation and retaining only the forward propagation terms gives 1 ∂p δ ∂2p βp ∂p ∂p − + + = 0. 2 ∂z c0 ∂t 2c0 ∂t ρ0 c30 ∂t 105 (A.14) applying the hange of variables z ′ = z − c0 t , Burgers' equation is obtained ∂p δ ∂2p βp ∂p + + = 0. ∂z ′ 2c0 ∂t2 ρ0 c30 ∂t 106 (A.15) APPENDIX B Simulations of ultrasound wave propagation The fast neareld method The fast neareld method (FNM) simulates the linear lossless pressure eld generated by transdu ers of various shapes. In time domain, the pressure eld produ ed by a ir ular piston is given by: ρ0 ca p(r, z; t) = π τ1 = where ˆ p π 0 r2 r cos ψ − a × [v(t − τ1 ) − v(t − τ2 )]dψ, + a2 − 2ar cos ψ z 2 + r 2 + a2 − 2ar cos ψ/c, a is the radius of the ir ular piston, τ1 and τ2 τ2 = z/c, are the delay times, and (B.1) (B.2) v is the normal velo ity for sour e points on the piston. The fast neareld method is an a urate method for omputing pressures in both the near eld and the far eld. The Cole-Hopf Model The Cole-Hopf model gives an exa t solution to Burgers' equation for given values of the nonlinearity and attenuation oe ients. For an arbitrary sour e pressure on the piston, p(0, t) = p0 F (t), 107 (B.3) the Cole-Hopf solution in luding both nonlinearity and attenuation is given by ′ p(z, t ) = p0 ´∞ F (t′′ )eEζ e−EG dt′′ ´∞ E , ζ e−EG dt′′ e −∞ −∞ βp0 Eζ (t ) = ρ0 δ ′′ EG = where t′ = t − z/c0 oe ient. When is retarded time, β equals zero, Eζ the linear ase. For ea h frequen y a distan e z. β ω, ˆ (B.4) t′′ F (t′′′ )dt′′′ , (B.5) −∞ c30 (t′ − t′′ )2 , 2zδ (B.6) is the nonlinearity oe ient, and δ is the attenuation equals zero, and the Cole-Hopf solution simplies to the amplitude is attenuated by exp(−ω 2 zδ/2c30 ) after Thus, for a single frequen y ex itation, the linear Cole-Hopf solution with attenuation only is given by p(z, t′ ) = p0 e−ω 2 zδ/2c3 0 sin (ωt′) . (B.7) Fay and Fubini Model When only nonlinearity is onsidered in the simulation, the Fay and Fubini solutions an be used for omparison. For periodi waves, the expression for the pressure an be expanded as a Fourier series, whi h learly shows how the harmoni s grow as the periodi waves propagate. In Fubini's model, the sour e pressure is given by p(0, t) = p0 sin (ωt) . (B.8) For this input, the Fubini solution for a single frequen y sour e is ∞ X 2 Jn (nσ) sin (nωt′ ) , p(σ, t ) = p0 nσ n=1 ′ 108 (B.9) where σ region σ ≤ 1. is a dimensionless distan e. In this region, when σ The Fubini solution is only valid in the pre-sho k in reases, the amplitude of the fundamental omponent de reases and the energy is transferred to higher harmoni omponents. After the sho k wave is fully developed, the Fay solution is hosen instead of the Fubini solution sin e the Fay solution is valid for larger values of σ. When both attenuation and nonlinearity are in luded, the Fay solution is expressed as ∞ p(σ, t′ ) = p0 where Γ is a parameter that des ribes 2 X sin (nωt′ ) , Γ n=1 sin[n(1 + σ)/Γ] (B.10) the attenuation. When only the ee ts of nonlinearity are in luded, the Fay solution simplies to a sawtooth wave ∞ 2 X sin (nωt′ ) , p(σ, t ) = p0 1 + σ n=1 n ′ whi h is valid for the region where σ > 3. Thus, by ombining the Fubini and Fay solutions, results obtained with Burgers' equation an be validated. needed in the region 1 < σ ≤ 3, (B.11) If additional omparisons are the solution in the transition region is needed [85℄. 109 BIBLIOGRAPHY 110 BIBLIOGRAPHY [1℄ P. Laugier and G. Haïat. Introdu tion to the physi s of ultrasound. In Ultrasound, pages 2945. Springer, 2011. [2℄ S. A. Goss, L. A. Frizzell, and F. Dunn. mammalian tissues. [3℄ K. J. Parker. Bone Quantitative Ultrasoni absorption and attenuation in Ultrasound in Medi ine & Biology, 5(2):181186, 1979. Ultrasoni attenuation and absorption in liver tissue. Medi ine & Biology, 9(4):363369, 1983. [4℄ A. S. Khimunin. Ultrasound in Numeri al al ulation of the dira tion orre tions for the pre ise measurement of ultrasound absorption. A ta A usti a United with A usti a, 27(4): 173181, 1972. [5℄ N. Purdie and M. M. Farrow. kineti s. The appli ation of ultrasound absorption to rea tion Coordination Chemistry Reviews, 11(3):189226, 1973. [6℄ K. K. Shung. On the ultrasound s attering from blood as a fun tion of hemato rit. IEEE Transa tions on Soni s and Ultrasoni s, 29(6):327330, 1982. [7℄ L. Y. L. Mo, I. Y. Kuo, K. K. Shung, L. Ceresne, and R. S. C. Cobbold. Ultrasound s attering from blood with hemato rits up to 100%. Engineering, 41(1):9195, 1994. [8℄ F. A. Du k. IEEE Transa tions on Biomedi al Physi al Properties of Tissue (First Edition), pages 99124. A ademi Press, London, 1996. [9℄ T. Lin, J. Ophir, and G. Potter. normal and diusely diseased liver. Frequen y-dependent ultrasoni dierentiation of J. A oust. So . Am., 82(4):11311138, 1987. [10℄ F. S. Foster and J. W. Hunt. Transmission of ultrasound beams through human tissue fo using and attenuation studies. Ultrasound Med. Biol., 5(3):257268, 1979. [11℄ C. Li, N. Duri , and L. Huang. Breast imaging using transmission ultrasound: Re onstru ting tissue parameters of sound speed and attenuation. In on BioMedi al Engineering and Informati s, 2008., International Conferen e volume 2, pages 708712. IEEE, 2008. [12℄ P. Laugier, P. Giat, and G. Berger. Broadband ultrasoni attenuation imaging: a new imaging te hnique of the os al is. Cal ied Tissue International, 54(2):8386, 1994. [13℄ M. Ribault, J.Y. Chapelon, D. Cathignol, and A. Gelet. Dierential attenuation imaging for the hara terization of high intensity fo used ultrasound lesions. Ultrasoni Imaging, 20(3):160177, 1998. [14℄ L. Bjørnø. Introdu tion to nonlinear a ousti s. 111 Physi s Pro edia, 3(1):516, 2010. [15℄ P. N. Wells. Ultrasoni imaging of the human body. Reports on Progress in Physi s, 62 (5):671, 1999. [16℄ P. J. Westervelt. Parametri a ousti array. Ameri a, 35(4):535537, 1963. The Journal of the A ousti al So iety of [17℄ Y. S. Lee and M. F. Hamilton. Time-domain modeling of pulsed nite-amplitude sound beams. The Journal of the A ousti al So iety of Ameri a, 97(2):906917, 1995. [18℄ R. O. Cleveland, M. F. Hamilton, and D. T. Bla ksto k. Time-domain modeling of niteamplitude sound in relaxing uids. The Journal of the A ousti al So iety of Ameri a, 99(6):33123318, 1996. [19℄ J. A. R. L. E. Berntsen. Numeri al al ulations of nite amplitude sound beams. Frontiers of Nonlinear A ousti s, pages 191196, 1990. [20℄ E. Hopf. The partial dierential equation ut+ uux= and Applied Mathemati s, 3(3):201230, 1950. µxx. Communi ations on Pure [21℄ S. Kutluay, A. R. Bahadir, and A. Özde³. Numeri al solution of one-dimensional Burgers equation: expli it and exa t-expli it nite dieren e methods. and Applied Mathemati s, 103(2):251261, 1999. Journal of Computational [22℄ R. O. Illing, J. E. Kennedy, F. Wu, G. R. Ter Haar, A. S. Protheroe, P. J. Friend, F. V. Gleeson, D. W. Cranston, R. R. Phillips, and M. R. Middleton. The safety and feasibility of extra orporeal high-intensity fo used ultrasound (HIFU) for the treatment of liver and kidney tumours in a Western population. British Journal of Can er, 93(8): 890895, 2005. [23℄ L. Poissonnier, J. Y. Chapelon, O. Rouviere, L. Curiel, R. Bouvier, X. Martin, J. M. Dubernard, and A. Gelet. patients. Control of prostate an er by transre tal HIFU in 227 European Urology, 51(2):381387, 2007. [24℄ C. C. Coussios, C. H. Farny, G. Ter Haar, and R. A. Roy. Role of a ousti avitation in the delivery and monitoring of an er treatment by high-intensity fo used ultrasound (HIFU). International Journal of Hyperthermia, 23(2):105120, 2007. [25℄ F. Tranquart, N. Grenier, V. Eder, and L. Pour elot. Clini al use of ultrasound tissue harmoni imaging. Ultrasound in Medi ine & Biology, 25(6):889894, 1999. [26℄ R. S. Shapiro, J. Wagrei h, R. B. Parsons, A. Stan ato-Pasik, H. C. Yeh, and R. Lao. Tissue harmoni imaging sonography: onventional sonography. evaluation of image quality ompared with AJR. Ameri an Journal of Roentgenology, 171(5):12031206, 1998. [27℄ E. L. Rosen and M. S. Soo. Tissue harmoni imaging sonography of breast lesions: improved margin analysis, onspi uity, and image quality ompared to onventional ultrasound. Clini al Imaging, 25(6):379384, 2001. 112 [28℄ A. Hirs hberg, J. Gilbert, R. Msallam, and A. P. J. Wijnands. Sho k waves in trombones. The Journal of the A ousti al So iety of Ameri a, 99(3):17541758, 1996. [29℄ A. Kilbas, H.M. Srivastava, and J.J. Trujillo. Theory and Appli ations of Fra tional Dierential Equations, Elsevier, North-Holland Mathemati s Studies, 204. Cal ulus and Applied Analysis, 9(1):71, 2006. Fra tional [30℄ M. Caputo. Linear models of dissipation whose Q is almost frequen y independent-II. Geophys. J. R. Astr. So ., 13:529539, 1967. [31℄ A. Atangana and D. Baleanu. New fra tional derivatives with nonlo al and non-singular kernel: theory and appli ation to heat transfer model. arXiv preprint arXiv:1602.03408, 2016. [32℄ U. N. Katugampola. A new approa h to generalized fra tional derivatives. Anal. Appl, 6(4):115, 2014. [33℄ M. Ishteva. Properties and appli ations of the Caputo fra tional operator. Bull. Math. MS . Thesis, 2005. Fra tional Dierential Equations: an Introdu tion to Fra tional Derivatives, Fra tional Dierential Equations, to Methods of Their Solution and Some of Their Appli ations, volume 198. A ademi Press, 1998. [34℄ I. Podlubny. [35℄ K. Singh, R. Saxena, and S. Kumar. Caputo-based fra tional derivative in fra tional Fourier transform domain. IEEE Journal on Emerging and Sele ted Topi s in Cir uits and Systems, 3(3):330337, 2013. [36℄ S. A. Goss, R. L. Johnston, and F. Dunn. ultrasoni properties of mammalian tissues. Comprehensive ompilation of empiri al J. A oust. So . Am., 64(2):423457, 1978. [37℄ S. A. Goss, R. L. Johnston, and F. Dunn. Compilation of empiri al ultrasoni properties of mammalian tissues. II. J. A oust. So . Am., 68(1):93108, 1980. [38℄ T. L. Szabo. Time-domain wave-equations for lossy media obeying a frequen y powerlaw. J. A oust. So . Am., 96(1):491500, 1994. [39℄ J. F. Kelly, R. J. M Gough, and M. M. Meers haert. Analyti al time-domain Green's fun tions for power-law media. [40℄ M. G. Wismer. J. A oust. So . Am., 124(5):28612872, 2008. Finite element analysis of broadband a ousti pulses through inho- mogenous media with power law attenuation. J. A oust. So . Am., 120(6):34933502, 2006. [41℄ W. Chen and S. Holm. Fra tional Lapla ian time-spa e models for linear and nonlinear lossy media exhibiting arbitrary frequen y power-law dependen y. 115(4):14241430, 2004. 113 J. A oust. So . Am., [42℄ B. E. Treeby and B. T. Cox. Modeling power law absorption and dispersion for a ousti propagation using the fra tional Lapla ian. J. A oust. So . Am., 127(5):27412748, 2010. [43℄ P. He. Simulation of ultrasound pulse propagation in lossy media obeying a frequen y power law. IEEE Trans. Ultrason. Ferroele t. Freq. Contr., 45(1):114125, 1998. [44℄ K. R. Waters, M. S. Hughes, J. Mobley, G. H. Brandenburger, and J. G. Miller. On the appli ability of Kramers-Kronig relations for ultrasoni attenuation obeying a frequen y power law. J. A oust. So . Am., 108(2):556563, 2000. [45℄ T. L. Szabo. power-law. Causal theories and data for a ousti attenuation obeying a frequen y J. A oust. So . Am., 97(1):1424, 1995. [46℄ D. T. Bla ksto k. Transient solution for sound radiated into a vis ous uid. So . Am., 41(5):13121319, 1967. [47℄ J. P. Nolan. J. A oust. Numeri al al ulation of stable densities and distribution fun tions. Commun. Statist. Sto h. Models, 13(4):759774, 1997. [48℄ M. J. Bu kingham. Causality, Stokes' wave equation, and a ousti pulse propagation in a vis ous uid. Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 72(2 Pt 2):026610, 2005. [49℄ J. P. Nolan. User Manual for STABLE 5.1. 2009. [50℄ D. Chen, J. F. Kelly, and R. J. M Gough. A fast neareld method for al ulations of time-harmoni and transient pressures produ ed by triangular pistons. Am., 120(5):24502459, 2006. J. A oust. So . [51℄ J. F. Kelly and R. J. M Gough. A time-spa e de omposition method for al ulating the neareld pressure generated by a pulsed ir ular piston. Ferroele t. Freq. Contr., 53(6):11501159, 2006. IEEE Trans. Ultrason. [52℄ T. L. Szabo. The material impulse response for broadband pulses in lossy media. Symposium on Ultrasoni s, 1:748751, 2003. IEEE [53℄ X. Zhao and R. J. M Gough. Time-domain omparisons of power law attenuation in ausal and non ausal time-fra tional wave equations. So iety of Ameri a, 139(5):30213031, 2016. The Journal of the A ousti al [54℄ S. Holm and R. Sinkus. A unifying fra tional wave equation for ompressional and shear waves. J. A oust. So . Am., 127(1):542548, 2010. [55℄ S. Holm and S. P. Näsholm. A ausal and fra tional all-frequen y wave equation for lossy media. [56℄ G. Pantis. J. A oust. So . Am., 130(4):21952202, 2011. The evaluation of integrals with os illatory integrands. Computational Physi s, 17(2):229233, 1975. 114 Journal of Handbook of Computational Methods for [57℄ P. K. Kythe and M. R. S häferkotter. Integration. CRC Press, 2004. [58℄ L. N. G. Filon. III.-On a Quadrature Formula for Trigonometri Integrals. of the Royal So iety of Edinburgh, 49:3847, 1930. [59℄ J. E. T. Penny and G. R. Lindeld. Numeri al Methods Using Matlab. Pro eedings Upper Saddle River, NJ : Prenti e Hall, 2nd edition, 2000. ISBN 0130126411 (paper). Previous ed.: 1995. [60℄ P. J. Davis and P. Rabinowitz. Methods of Numeri al Integration. Courier Corporation, 2007. [61℄ B. E. Treeby and B. T. Cox. A k-spa e Green's fun tion solution for a ousti initial value problems in homogeneous media with power law absorption. J. A oust. So . Am., 129(6):36523660, 2011. [62℄ E. A. Zabolotskaya. Quasiplane waves in the nonlinear a ousti s of onned beams. Sov. Phys. A oust., 15:3540, 1969. [63℄ V. P. Kuznetsov. Equation of nonlinear a ousti s. Sov. Phys. A oust., 16(4):467470, 1971. [64℄ J. E. Soneson. A parametri study of error in the paraboli approximation of fo used axisymmetri ultrasound beams. The Journal of the A ousti al So iety of Ameri a, 131 (6):EL481EL486, 2012. [65℄ F. P. Curra, P. D. Mourad, V. A. Khokhlova, R. O. Cleveland, and L. A. Crum. Numeri al simulations of heating patterns and tissue temperature response due to high-intensity fo used ultrasound. IEEE Transa tions on Ultrasoni s, Ferroele tri s, and Frequen y Control, 47(4):10771089, 2000. [66℄ E. A. Filonenko and V. A. Khokhlova. Ee t of a ousti nonlinearity on heating of biologi al tissue by high-intensity fo used ultrasound. A ousti al Physi s, 47(4):468475, 2001. [67℄ J. Huijssen and M. D. Verweij. An iterative method for the omputation of nonlinear, wide-angle, pulsed a ousti elds of medi al diagnosti transdu ers. A ousti al So iety of Ameri a, 127(1):3344, 2010. The Journal of the [68℄ G. Mur. Absorbing boundary onditions for the nite-dieren e approximation of the time-domain ele tromagneti -eld equations. Compatibility, (4):377382, 1981. IEEE Transa tions on Ele tromagneti [69℄ J. P. Berenger. A perfe tly mat hed layer for the absorption of ele tromagneti waves. Journal of Computational Physi s, 114(2):185200, 1994. [70℄ D. S. Katz, E. T. Thiele, and A. Taove. Validation and extension to three dimensions of the Berenger PML absorbing boundary ondition for FD-TD meshes. 115 IEEE Mi rowave and Guided Wave Letters, 4(8):268270, 1994. [71℄ X. Yuan, D. Borup, J. W. Wiskin, M. Berggren, R. Eidens, and S. A. Johnson. Formulation and validation of Berenger's PML absorbing boundary for the FDTD simulation of a ousti s attering. IEEE Transa tions on Ultrasoni s, Ferroele tri s, and Frequen y Control, 44(4):816822, 1997. [72℄ Q. H. Liu and J. Tao. media. The perfe tly mat hed layer for a ousti waves in absorptive The Journal of the A ousti al So iety of Ameri a, 102(4):20722082, 1997. [73℄ Q. H. Liu. oordinates. Perfe tly mat hed layers for elasti waves in ylindri al and spheri al The Journal of the A ousti al So iety of Ameri a, 105(4):20752084, 1999. [74℄ S. Abarbanel, D. Gottlieb, and J. S. Hesthaven. Well-posed perfe tly mat hed layers for adve tive a ousti s. Journal of Computational Physi s, 154(2):266283, 1999. [75℄ F. Q. Hu. On absorbing boundary onditions for linearized euler equations by a perfe tly mat hed layer. Journal of Computational Physi s, 129(1):201219, 1996. [76℄ Y. L. Sheu and P. C. Li. Simulations of photoa ousti wave propagation using a nitedieren e time-domain method with Berenger's perfe tly mat hed layers. of the A ousti al So iety of Ameri a, 124(6):34713480, 2008. The Journal [77℄ Y. Ma, J. Yu, and Y. Wang. A novel unsplit perfe tly mat hed layer for the se ond-order a ousti wave equation. [78℄ J. H. Ehrli h. Ultrasoni s, 54(6):15681574, 2014. Time domain modeling of a ousti propagation with a ousti wave propagator and absorbing boundary onditions. Ameri a, 123(5):3530, 2008. Journal of the A ousti al So iety of [79℄ G. F. Pinton, J. Dahl, S. Rosenzweig, and G. E. Trahey. attenuating full-wave model of ultrasound. A heterogeneous nonlinear IEEE Transa tions on Ultrasoni s, Ferro- ele tri s, and Frequen y Control, 56(3), 2009. [80℄ B. Kaltenba her, M. Kaltenba her, and I. Sim. A modied and stable version of a perfe tly mat hed layer te hnique for the 3-d se ond order wave equation in time domain with an appli ation to aeroa ousti s. Journal of Computational Physi s, 235:407422, 2013. [81℄ Z. Chen and X. Wu. Long-time stability and onvergen e of the uniaxial perfe tly mat hed layer method for time-domain a ousti s attering problems. Numeri al Analysis, 50(5):26322655, 2012. [82℄ E. J. Alles and K. W. van Dongen. SIAM Journal on Perfe tly mat hed layers for frequen y-domain integral equation a ousti s attering problems. IEEE Transa tions on Ultrasoni s, Ferroele tri s, and Frequen y Control, 58(5):10771086, 2011. [83℄ T. K. Katsibas and C. S. Antonopoulos. A general form of perfe tly mat hed layers for for three-dimensional problems of a ousti s attering in lossless and lossy uid media. 116 IEEE Transa tions on Ultrasoni s, Ferroele tri s, and Frequen y Control, 51(8):964 972, 2004. [84℄ K. Duru. A perfe tly mat hed layer for the time-dependent wave heterogeneous and layered media. equation in Journal of Computational Physi s, 257:757781, 2014. [85℄ D. T. Bla ksto k. Conne tion between the Fay and Fubini solutions for plane sound waves of nite amplitude. The Journal of the A ousti al So iety of Ameri a, 10191026, 1966. 117 39(6):