On nonlinear effects in holographic-modulated ultrasound

Focused ultrasound (FU) is emerging noninvasive and nonionizing therapeutic technology that has the potential to treat a variety of medical conditions. FU can heat up, destroy, or change target tissue while minimizing damage to tissue outside the target area by precisely focusing ultrasonic beams.^1,2 FU has been shown as an effective tool for tumor ablation,^3,4 lithotripsy,^5,6 histotripsy,^7–9 drug delivery,^10–12 and contactless ultrasound power transfer.¹³ The current FU technology focuses ultrasound energy using a curved single-element or phased-array transducers (PATs). Single-element transducers can only target a single location, whereas PATs have limited fidelity and require expensive and complex hardware. Holographic acoustic lenses (HALs), i.e., acoustic holograms, have been proposed as a simple, robust, and cost-effective method for creating highly localized ultrasound fields.^14–19 HALs, which are often 3D printed, act as phase plates and can be used to reconstruct ultrasound onto a desired spatial distribution deep inside living tissue. As a result, acoustic holograms are being rapidly investigated for various biomedical applications with the potential to enable low-cost therapeutic procedures such as transcranial neuromodulation,²⁰ blood–brain barrier opening,²¹ ultrasound-induced hyperthermia,²² and burst-wave lithotripsy.⁵ Furthermore, HALs have also been proposed for producing precise cavitation patterns²³ for various therapeutic applications. Many FU therapeutic techniques capitalize on the rise of nonlinear effects at high intensities to enhance the desired thermal and mechanical effects. At higher intensities, nonlinear propagation effects cause a sharp increase in ultrasound energy absorption resulting in accelerated and more localized thermal tissue ablation.²⁴ Short ultrasound pulses with shocks can also be used to generate purely mechanical liquefaction of tissue.^25,26 Acoustic holograms, on the other hand, are currently computed using linear propagation techniques and have only been demonstrated in the linear regime. If the acoustic amplitude is sufficiently large, the linearity assumptions are no longer valid and will result in significant errors in the predicted sound field. Studies that highlight the interaction of complex physics between the constructed holographic sound field and nonlinearities are absent. As a result, a significant potential role for acoustic holograms has yet to be explored: creating precise and complex high-intensity focal regions, particularly for clinical treatments. The current study numerically investigates the effects of increasing acoustic intensity on holographic modulated ultrasound. Modeling three-dimensional nonlinear ultrasound is a physically and computationally demanding problem. As such several difficulties coexist: the complex diffraction structure necessitates the use of accurate diffraction models and a fine spatial numerical grid. Furthermore, strong nonlinear effects necessitate the inclusion of a large number of harmonics in simulations with fine temporal and spatial grids.^27,28 In general, three-dimensional nonlinear simulations require large memory (RAM) and lengthy computation time. In this work, the k-space pseudospectral method was used to perform homogeneous three-dimensional nonlinear acoustic simulations with power law absorption.²⁸ The numerical model solves a set of coupled partial differential equations equivalent to a generalized Westervelt equation. The accuracy and validity of this numerical model have been previously verified with analytical solutions and experimental measurements.^29,30 The simulations were run in an optimized C++ implementation that makes use of the OpenMP interface for parallelized execution across multiple cores. The C++ implementation is an extension of the open-source K-wave MATLAB toolbox, which is available online.³¹ See the supplementary material for verification of the present numerical approach by previously published experimental data. The simulations were carried out on a high-performance computing cluster hosted at the Virginia Tech campus (Blacksburg, VA, USA). The used computational node has 128 cores (two AMD EPYC 7702 64-core processors) and 1024 GB of RAM. The largest simulation in this study took 56 h to run with a spatial discretization of 60 grid points per wavelength (supporting up to 30 harmonics) and a Courant–Friedrichs–Lewy (CFL) number of 0.3. The calculations were carried out in water using the appropriate physical parameters: $ρ0=1000$ kg$/$m³, $c0=1481$ m$/$s, $β=3.6, α0=0.0022$ dB$/$MHz, and γ = 2 for an excitation frequency of $f0=1$ MHz. Here, ρ₀, c₀, and β are the density, ambient speed of sound, and nonlinearity coefficient, respectively. While α₀ and γ, respectively, are the absorption coefficient and the power law exponent. The size of the spatial domain is $50×50×30$ mm³ in addition to the size of the perfectly matched layers (PMLs). To ensure that a steady state is reached throughout the computational domain, the temporal domain is extended to 12 000 time steps (50 time periods with respect to the excitation frequency). All simulations in this paper used the same set of parameters. To investigate the nonlinear effects in holographic-modulated ultrasound, the iterative angular spectrum algorithm (IASA)¹⁴ was first used to create a star-shaped amplitude pattern 25 mm from the hologram plane, as shown in Fig. 1. The algorithm was initialized with a uniform phase and amplitude input and then used to iteratively calculate the required phase map at the hologram plane (50 iterations). Figure 1 shows the computed phase map at the hologram plane and the simulated linear pressure field at the target plane. The computed phase was then used as a Dirichlet boundary condition in the nonlinear simulations. The linear sound field was first calculated with a uniform source pressure amplitude of $P0=0.1$ MPa. The source pressure was then increased to $P0=1.8$ MPa to investigate the emergence of nonlinearities and their effects on the sound field. The compressional phase in a propagating nonlinear wave travels at a speed of $c0+βu$ faster than the speed of the rarefaction phase $c0−βu$⁠, where u is the particle velocity.^32,33 As a result, the sinusoidal shape is distorted to a steeper waveform as the wave propagates further. Harmonics at integer multiples of the fundamental frequency are generated as a result of waveform distortion, resulting in enhanced heating and streaming effects. Furthermore, in the nonlinear problem, diffraction and focusing effects cause a relative phase shift in the harmonics with respect to the fundamental component, resulting in an asymmetric waveform and a discrepancy between the peak positive pressure (PPP) and the peak negative pressure (PNP). The effects mentioned above are depicted schematically in Fig. 1. The PPP is connected to the thermal effects in therapeutic ultrasound, whereas the PNP controls the cavitation effect. Understanding and characterizing such effects are essential to ensure that therapies are administered at the proper place and intensity, as required by the Food and Drug Administration (FDA) and International Electrotechnical Commission (IEC) rules.³⁴ The spatial distributions of important acoustic parameters at the target plane in the linear and nonlinear regimes are depicted in Fig. 2. For linear $P0=0.1$ MPa (a) and (b) and nonlinear $P0=1.8$ MPa (c)–(f) conditions, the PPP (⁠$P+$⁠), PNP (⁠$P−$⁠), intensity, and heat deposition are shown. The PPP and PNP are spatially identical in the linear case. Furthermore, the linear spatial structure of the wave intensity and heat deposition are identical. As a result, only the pressure amplitude is shown in (a), whereas the intensity and heat deposition are both shown in a single normalized plot in (b). The spatial structure of the acoustic parameters deviates from the initial linear prediction due to nonlinear distortion at higher intensities. When compared to the linear field, the PPP (c) and especially heat deposition (d) are more localized and less uniform in space in the nonlinear regime. The nonlinear spatial structure of the PNP (e) deviates from the linear prediction as well, though to a lesser extent. In contrast, the intensity distribution differs only slightly between the linear (b) and nonlinear (f) conditions. The spatial structures and magnitudes of the PPP and PNP are no longer identical in the nonlinear regime. It is also clear that the intensity cannot be used to estimate heat deposition in the nonlinear regime. Furthermore, due to the nonlinear distortion and diffraction, the nonlinear waveform shown in Fig. 2(g) exhibits top-bottom asymmetry, and the wave is no longer monochromatic as the corresponding spectrum shows, Fig. 2(h). As a consequence of the waveform's distorted asymmetric profile, the real gain of the PPP (⁠$P+/P0$⁠) increases while the gain of the PNP (⁠$P−/P0$⁠) decreases with respect to the linear value. Figure 3(a) compares the gain of the PPP and PNP to the linear gain as the amplitude of the source pressure increases. The values of the maximum peak pressures at the target plane were normalized by the amplitude of source pressure and, thus, represent the real gain. The gain of the PPP and PNP at low source pressure is that of linear propagation. However, under nonlinear conditions, the PPP gain grows larger than the linear gain to a maximum of 11. Meanwhile, as the source pressure increases, the gain for PNP decreases. As a result, where there are highly distorted waveforms, the maximum PPP grows faster while the maximum PNP grows more slowly. In the nonlinear regime, this is the underlying mechanism that results in a spatial structure that is more localized for the PPP and flattened for the PNP. Since higher harmonics increase wave energy absorption, it is critical to investigate the spectral distribution of acoustic energy at the target plane. Total harmonic distortion (THD) is defined as the ratio of all harmonic component powers to fundamental frequency power³⁵ 

Here, A is the amplitude, and n is an integer representing multiples of the fundamental frequency. Nonlinear distortion of the waveforms due to the harmonics and relative phase shift is highest at the points that exhibit the maximum THD in the domain.³⁴ Figure 3(b) shows the distribution of THD at the target plane. The points with the maximum THD have harmonics with a total power that is 45% that of the fundamental. When compared to the PPP, the heat deposition has a more localized distribution due to the higher harmonics. It can also be seen that the points with the highest THD also have the highest heat deposition. To quantify the deviation of the spatial structure of the acoustic parameters as the source pressure is increased, two metrics were used. Equations (2) and (3) define the peak signal-to-distortion ratio (PSDR) and the nonuniformity index (NUI), respectively,

Here, m is the data matrix of the field, and g is the binary matrix of the star-shaped pattern (ground truth), while n_X and i are, respectively, the number and index of the rows in the data matrix, and n_Y and j are, respectively, the number and index of the columns

Here, the mean and the max of the data matrix were calculated for the region of interest (ROI), which is the area enclosed by the star pattern. NUI was selected to study the effect of nonlinear distortion on the spatial localization of the acoustic parameters, while PSDR is used to assess the effect of the distortion on the reconstruction quality. PSDR and NUI were calculated as the source pressure gradually increased from $P0=0.1$ MPa to $P0=1.8$ MPa. Figure 3(c) shows the PSDR drops by 1 dB for the PPP and by 2 dB for the heat deposition. The PSDR for the intensity is essentially constant, while it increases slightly for the PNP. Figure 3(d) depicts a similar story where the NUI for the PNP and intensity change only slightly. However, the NUI increases for PPP and more rapidly for heat deposition with a maximum value of 700%.

Cavitation can cause significant desired or unintended bio-effects in the tissue in medical acoustics.²³ The mechanical index (MI) was used to calculate the inertial cavitation zones within the target plane. The MI is calculated by dividing the PNP in MPa by the square root of the excitation frequency in MHz.³⁶ As the MI involves a trade-off between PNP and excitation frequency, the MI can be a more comprehensive and reliable measure than PNP alone.³⁷ The cavitation zones were calculated for the simulated nonlinear field and a linear prediction with a source pressure of $P0=1.8$ MPa. Figure 4 depicts the predicted saturated pixels with a contour line encapsulating the cavitation zones. The structural similarity index (SSIM) between the star pattern and the obtained cavitation zones was also calculated. The predicted cavitation zones for the linear (top row) and nonlinear (bottom row) fields can be seen. The focusing gain of the PNP decreases due to nonlinear distortion. As a result, in the linear case, larger cavitation zones are predicted. The MI saturation levels were selected according to experimentally tested cavitation thresholds of ex vivo tissue.³⁸ Figures 4(a) and 4(b) depict the MI saturation level of 5, Figs. 4(c) and 4(d) depict the saturation level of 5.5, and Figs. 4(e) and 4(f) depict the saturation level of 6. When nonlinear effects in the sound field were taken into account, the SSIM decreased by an average of 0.32. As a result, nonlinear modeling is critical to avoid overestimation of the PNP and, as a result, cavitation effects. Therefore, nonlinear modeling is essential to avoid overestimation of the PNP and, consequently, the cavitation effects.

As opposed to a spherically focused ultrasound field which has a clear focal point and focal length, a holographically modulated ultrasound field has a complex spatial distribution with intricate regions of high and low pressure. As a result, it is critical to consider the acoustic parameters of the entire field rather than just the target plane. Figure 5 shows two cross-sections of the sound field, the YZ plane (Top row) and the XZ plane (Bottom row) both extend from the hologram plane at Z = 0 mm to Z = 30 mm (5 mm beyond the target plane). Two contour lines are also plotted at a pressure drop of −3 and −6 dB. When compared to the linear distributions in (a) and (b), the PPP in (c) and (d) show higher localization, while the PNP in (e) and (f) show significant spreading. The plots also highlight focal areas that occur before the target plane, as shown in the YZ plane, or after, as shown in the XZ plane. Because nonlinear effects are inherently cumulative, it is critical to examine the entire sound field. In computing acoustic holograms, the pressure field is only constrained at the target plane, thus, significant bio-effects or shock wavefronts may form before or after the target plane at undesirable locations.

Finally, as a special case of holographic modulation, we will consider an axisymmetric ring pattern. Previously reported experimental and numerical studies on spherically focused transducers have revealed a spatial focal shift of the PPP and PNP at high intensities.^39–41 It would be interesting to look into axisymmetric holographic ultrasound fields and investigate whether they exhibit similar behavior. An axisymmetric simulation for a ring pattern, as shown in Fig. 6, was also implemented to this end. IASA was used to construct the desired pressure by imposing a binary amplitude ring pattern with a radius of 20 mm at the target plane. The simulations were more memory efficient due to the axisymmetric geometry, and up to 60 harmonics were included. The spatial domain had dimensions of 25 mm in the radial direction and 30 mm in the axial direction with the target plane located at Z = 25 mm in the axial direction. The axial cross section of the linear sound field with a snapshot of the ring pattern at the target plane is shown in Fig. 6(a). Figures 6(b) and 6(c) show the PPP and PNP with a black contour line at −1 dB. A couple of nonlinear effects can be deducted from observing Figs. 6(a)–6(c). The PPP focal region becomes narrower and spatially shifts further from the source (right), while the PNP slightly shifts toward the direction of the source (left). The curves in Fig. 6(d) represent the spatial location of the maxima of the peak pressures as the source pressure increases. At a source pressure of 3 MPa, the maximum PPP shifts by 2 mm away from the source, while the maximum PNP shifts by 0.8 mm toward the source. The behavior of the maximum pressures presented in Fig. 6(d) is due to the effects of nonlinear refraction. Compressive phases travel faster and self-defocus, while rarefaction phases travel slower and self-focus.⁴⁰ As nonlinear effects increase at the target plane, the beam width of the focal decreases by 50% for the PPP and remains constant for the PNP.

In conclusion, we have numerically examined the effects of nonlinear distortion on holographically modulated ultrasound generated by HALs. Using the k-space pseudospectral method, we have performed homogeneous three-dimensional nonlinear acoustic simulations with power law absorption. Due to the computationally intensive nature of the problem, the model was implemented in C++ on a high-performance computing cluster. The 3D spatial distribution of the PPP, PNP, intensity, heat deposition, and cavitation zones were then analyzed for a range of excitation levels. It was shown that the nonlinear waveform distortion leads to pronounced effects on the spatial structure of the sound field. Due to diffraction and focusing effects, the waveform becomes asymmetric with different PPP and PNP patterns. The areas of maximum nonlinear distortion become increasingly focused in terms of the PPP distribution, the NUI rises to 190%, and the PSDR drops by 1 dB. Harmonic generation also causes the heat deposition distribution to become extremely confined with a NUI of 700%. Nonlinear effects have also been shown to flatten the peak negative pressure distribution with minimal impact on the PSDR or NUI. In addition, nonlinear effects were shown to be vital for correctly estimating cavitation zones. When designing HALs for high-power therapeutic ultrasound techniques, it is crucial to consider the dynamic behavior of the spatial structure of various acoustic parameters in a nonlinear field. Our findings contribute to the informed use of HALs in high-intensity applications. We also envision that this investigation can be a precursor for encompassing nonlinear effects in the realm of computer-generated holography, allowing for the creation of specific nonlinear acoustic fields to improve the efficacy and precision of therapeutic procedures. Furthermore, manipulating the amplitudes or phases of higher harmonics could improve applications such as harmonic imaging and enhancing peak positive or peak negative pressures. The present investigation considered a homogeneous domain with the physical parameters of water. Our findings are also applicable to other nonlinear domains; however, according to the power law absorption, larger acoustic absorption will reduce the amplitude of the harmonics while a higher nonlinearity parameter will lead to faster harmonic generation.⁴² In addition, in heterogeneous domains such as biological tissue, transmission, reflection, and refraction at the interfaces should also be considered.

See the supplementary material for verification of the numerical approach.

This work was supported by the U.S. National Science Foundation (NSF) under EAGER and CAREER grants, Award Nos. CMMI 2121933 and CMMI 2143788, which are gratefully acknowledged.