## 2 Theory and Experiments
### The Numerical Algorithm  
In the numerical analysis, both the KZK model and its modified form were solved with the finite difference time domain (FDTD) method for the acoustic field emitted from a single-element self-focusing transducer. For the traditional model, coordinate transformation \(Z=z/F\), \(R=r/a\), \(T=\omega t\), \(P=p/P_{0}\) were introduced, with \(F\), \(r\), \(a\), \(P_{0}\) being the geometrical focal length, the radial coordinate, the aperture radius of the transducer, and the surficial acoustic pressure, respectively. The following assumptions were then made,  
\[G=\frac{ka^{2}}{2F}, \tag{7}\]  
\[A=\frac{\omega^{2}\delta}{2\rho_{0}c_{0}{}^{3}}\,F=\alpha^{r}F, \tag{8}\]  
\[N=\frac{F}{\rho_{0}c_{0}{}^{3}/(P_{0}\beta\omega)}=\frac{F}{l_{d}}, \tag{9}\]  
in order that the KZK equation could be normalized to the following form,  
\[\frac{\partial^{2}P}{\partial T\partial Z}=\frac{1}{4G}\bigtriangleup_{\bot }P+A\frac{\partial^{3}P}{\partial T^{3}}+\frac{N}{2}\frac{\partial^{2}P^{2}}{ \partial T^{2}}. \tag{10}\]  
The values of the physical constants used for acoustic modeling were \(\rho_{0}\) = 1000 kg/m\({}^{3}\), \(c_{0}\) = 1486 m/s, \(\beta\) = 3.5, \(\alpha\) = 0.025 Np/m at 1 MHz, and \(\mu\) = 2 for fresh water [45]. The FDTD algorithm adopted here was generally the same as that used in a previous study [45; 46], in which Equation (10) was decomposed into three independent equations, accounting for the diffraction, attenuation and nonlinearity, respectively. Based on an orthogonal spatial grid, the discretized forms of these equations were expressed as,  
\[\frac{P_{i,j+1}^{n}-P_{i,j}^{n}}{dZ}=\frac{1}{4G}\int\limits_{T_{\min}}^{T}{ \left(\frac{P_{i+1,j+1}^{n}-2P_{i,j+1}^{n}+P_{i-1,j+1}^{n}}{dR^{2}}+\frac{P_{i+1,j+1}^{n}-P_{i-1,j+1}^{n}}{2 kdR^{2}}\right)dt}, \tag{11}\]  
\[\frac{P_{i,j+1}^{n}-P_{i,j}^{n}}{dZ}=A\frac{P_{i,j+1}^{n+1}-2P_{i,j+1}^{n}+P_{i,j+1}^{n-1}}{dT^{2}}, \tag{12}\]  
And  
\[P_{j+1}^{n}=\left\{\begin{array}{ll}P_{j}^{n}{\left(1-N\frac{P_{j}^{n+1}-P_{ j}^{n}}{dT}dZ\right)^{-1}}&P_{j}^{n}\geq 0\\ P_{j}^{n}{\left(1-N\frac{P_{j}^{n}-P_{j}^{n-1}}{dT}dZ\right)^{-1}}&P_{j}^{n}< 0\end{array}\right.. \tag{13}\]  
Here \(i\) and \(j\) were the spatial coordinate indexes in the radial and axial directions, respectively, and \(n\) was the time step, i.e., \(P_{i,j}^{n}=P(i\cdot dr,j\cdot dz;n\cdot dt)\).  
In the simulations for the mKZK equation, the major difference was the differential form of the fractional derivative in Equation (6),  
\[\frac{\partial^{y+1}p}{\partial\tau^{y+1}}=\frac{A}{\Delta\tau^{y}}\times{ \left[\Delta\tau^{-1}\sum\limits_{r=0}^{n}\omega_{r}^{(y)}\left(P_{i,j}^{n-r+ 1}+P_{i,j}^{n-r}\right)+\omega_{n+1}^{(y)}P^{\prime}(0)_{i,j}^{k}\right]}, \tag{14}\]  
where  
\[A=2\Gamma(-y)\Gamma(y+1)\cos[(y+1)\pi/2]/\pi, \tag{15}\]  
\[\omega_{r}^{(y)}=(-1)^{r}y(y-1)\Gamma(y-r+1)/r!. \tag{16}\]  
The FDTD simulations were then accomplished through a self-developed FORTRAN code package running on a \(\times\)64 PC platform. In the calculation, the boundary condition was symmetrical and had equal amplitude acoustic pressure driving conditions, which was also the boundary condition generally used to calculate HIFU [29; 41].