Comparative Characterization of Nonlinear Ultrasound Fields Generated by Sonalleve V1 and V2 MR-HIFU Systems.

Performance characteristics of the V1 and V2 Sonalleve arrays determined from the hydrophone measurements combined with simulations are presented and analyzed in the next subsections. The arrays are compared in terms of the boundary conditions obtained from the holography measurements, the dimensions of the focal regions in the linear propagation regime, the manifestation of nonlinear effects, and the output levels needed for formation of developed shocks. The last subsection shows the results for an axially symmetric equivalent source matched to each array. The importance of including the central opening in the equivalent source model of the V2 array is demonstrated. The first step in the acoustic field characterization of Sonalleve V1 and V2 arrays was the collection of hydrophone measurements at a low output level associated with linear acoustic propagation. These measurements included 2-D holography scans acquired in a prefocal plane. Magnitude and phase distributions of the measured pressure field holograms are shown in Fig. 2 . Notably, the V1 magnitude plot ( Fig. 2a ) has comparable pressures on and off axis whereas the V2 array ( Fig. 2c ) has a region near the beam axis with lower pressures. This difference in the pressure pattern is caused by the larger central opening of the V2 array. The phase distributions of the pressure field holograms ( Fig. 2b , d ) are axisymmetric, showing the sphericity of the wavefronts. In addition, the phase axial symmetry indicates good alignment between the holographic measurement planes and spatial orientation of therapeutic arrays because phase is sensitive to the angular positioning inaccuracy. To align the holography plane with modeling coordinates and have peak pressures at the focus remaining on the z axis, the hologram plane was rotated 0.8° around the x axis and shifted by 2.75 and 0.45 mm along the x and y axes, respectively. Alignment procedure was based on the calculation of two rotation angles from parameters of the line drawn through pressure field maxima in several planes parallel to the holography plane in the vicinity of the focal maximum [ 30 ]. Then, based on these angles, transformation of the Cartesian coordinates was performed and used to numerically calculate the angular spectrum representation of the hologram written for the new coordinates, where hologram was perfectly aligned. After alignment, the measured holograms from Fig. 2 were then backpropagated to the plane z = 0 mm at the transducer apex using the angular spectrum method. The resulting source holograms ( Fig. 3 ) were employed as boundary conditions to the 3-D Westervelt model. Even though these source holograms were calculated on a plane and not on the spherical surface of each transducer, some relevant features are still identifiable. Comparing Fig. 1 and Fig. 3(a , c) , it is seen that both arrays had non-functional elements. For the V1 array, a missing element in the upper-left region is apparent. For the V2 array, there were three missing elements in two sectors of the lower-left quadrant as well as one absent element in the upper-right sector close to the central opening. The boundary conditions also allow for estimating the approximate size of the central opening for each transducer array (20 mm for V1 and 44 mm for V2). Phase distributions shown in Fig. 3(b , d) are approximately axisymmetric. The reconstructed source holograms were used to evaluate the structure of linear acoustic fields for each array. Forward propagation from each source hologram yielded the fields depicted in Figs. 4 and 5 . Validation of these numerical simulations was performed by direct comparison of holography-based calculations against independent hydrophone measurements ( Fig. 4 ). In the Fig. 4 plots, pressure amplitudes were normalized relative to corresponding peak pressure values in the measurements. The coordinates in Fig. 4 denote distances relative to the focal maximum for each array. These comparisons show that for both arrays the field structure obtained by linear simulation from the holography-based boundary condition is in good agreement with direct hydrophone measurements. As previously assumed, the V2 array turned out to be less focused and its focal region is larger compared to the V1 array ( Figs. 4 , 5 ). Notably, the dimensions of the focal regions differ primarily in the axial direction while remaining virtually identical in the transverse directions ( Fig. 4 ). The −6 dB dimensions of the linear focal lobes at 1.2 MHz are 1.5 mm × 1.5 mm × 9.2 mm for the V1 array and 1.6 mm × 1.6 mm × 12.5 mm for the V2 array. Two-dimensional pressure amplitude distributions ( Fig. 5 ) in the axial planes x z and y z have identical structures for each array ( Fig.5 a , b for V1 and Fig.5 d , e for V2). In addition, focal regions in the transverse x y plane have a circular shape, demonstrating that both arrays generate focused linear fields with axial symmetry. To determine the focusing gain of each array, a nominal source intensity I 0 was defined as the ratio of power W 0 of the measured hologram to the surface area S of 256 transducer elements: I 0 = W 0 / S . From I 0 , a nominal source pressure p 0 at the array elements was then defined in the plane wave approximation: p 0 = 2 ρ 0 c 0 I 0 . Finally, the pressure focusing gain G was defined as a pressure amplitude at the focus p F divided by nominal source pressure amplitude p 0 at the array elements: G = p F / p 0 . In this way, the pressure focusing gain of the V1 array ( G = 67 ) is 1.3 times higher than that of the V2 array ( G = 51 ) . The lower V2 focusing gain is consistent with its longer focal region ( Fig. 5 ). Nonlinear acoustic field characterization included the results of 3-D nonlinear acoustic modeling with boundary conditions reconstructed from the measured holograms and direct fiber optic hydrophone measurements of focal waveforms, both performed for a wide range of acoustic powers. The output levels used in simulations for the V1 and V2 arrays are provided in Tables I and II , respectively. For each output level, the first two columns in each table list the system setting for ampvals and acoustic power and the corresponding acoustic output power based on hydrophone measurements. To determine actual acoustic powers (column 2 of Tables I and II ) for different settings, we first calculated the power of each measured hologram using the angular spectrum approach (formula (83) from Ref. 21 ). The holograms were measured at 259 ampvals and 430 ampvals for V1 and V2 arrays, respectively. Then, the power for another value of ampvals was calculated by multiplying the hologram power by a scaling factor determined as the square of the ratio of pressure values from the near-source measurements. Note that everywhere below, the power will be indicated as the acoustic power determined in this way (column 2 of Tables I and II ). Also, the nominal source pressure p 0 at the array elements along with the nominal source intensity I 0 calculated as described in the previous subsection are given in the third and fourth columns, respectively. Comparisons of peak positive pressure p + and peak negative pressure p - obtained from both modeling and measurements at the focus of the V1 and V2 arrays are presented in Fig. 6 over the range of all measured output levels. Measured waveforms were processed by considering averages of peak values over 8 consecutive cycles, with mean values plotted as circles ( Fig. 6 ). In experiment, multiple waveforms at powers less than 150 W were averaged during acquisition to minimize noise inherent to the FOPH. At higher output powers, where the signal-to-noise ratio was naturally improved, averaging was avoided to decrease the exposure time given concerns about cavitation at the tip of the FOPH. Numerical simulations show good agreement with the FOPH data for both V1 and V2 arrays ( Fig. 6 ). Across the entire range of output levels, simulation results for the peak positive pressure remain within 7 MPa of the corresponding measurement data with the largest relative discrepancy of 10%. For the peak negative pressure, simulated and measured data differ less than 2.5 MPa with the largest discrepancy of about 20%. The standard deviation in measured peak pressures over 8 acoustic cycles was in the range of 3-10% with lower deviation values at low powers. All mentioned above discrepancy values are consistent with those presented in [ 7 ]. As illustrated in Fig. 7 , waveforms at the focus were measured and simulated over a range of output levels. The waveforms in the first row ( Fig. 7 a , e ) correspond to the threshold of quasilinear waveform distortion, at which 10% of the total wave intensity is distributed over harmonics of the fundamental frequency [ 31 ]. Presented in the third row ( Fig.7 c , g ) are waveforms with a fully developed shock that is characterized by zero-pressure level of the bottom edge of the shock front, i.e. the shock amplitude is equal to the peak positive pressure [ 11 ]. In this case, the shock amplitude normalized to the source pressure p 0 reaches a maximum. The quasilinear case ( Fig. 7 a , e ) and the case of developed shock formation ( Fig.7 c , g ) are conventional thresholds for characterizing the strength of nonlinear effects. Below the quasilinearity threshold, propagation can be considered to be linear. Beyond the quasilinear threshold, nonlinear effects become prominent, with a shock appearing near the positive peak of the waveform ( Fig.7 b , f ). With further increases in the output level, the bottom edge of the shock moves toward the zero pressure level until the shock is fully developed. Beyond the level of developed shock formation in the focal waveform, the growth rate of the shock amplitude and peak pressures slow down and their values gradually saturate ( Fig.7 d , h ). Interestingly, the acoustic powers corresponding to these threshold characteristic cases differ less than 10% for V1 and V2 arrays: The quasilinearity threshold is reached at about 60 W for both arrays. Developed shocks form at 525 W and 515 W for the V1 and V2 arrays, respectively. As observed for the linear case, the nonlinear field of the V1 array has higher peak pressure values than those of the V2 array at the same power ( Figs. 6 and 7 ). Nonlinear simulations allow tracking the shock characteristics more easily than the FOPH measurements because of better temporal resolution and the absence of uncertainties associated with deconvolution of the measured waveforms with significant components at higher harmonic frequencies. Shock amplitudes A sh in the simulated acoustic waveforms, shown in Fig. 6 , were calculated by determining the beginning and end of the shock front from time points at which the time derivative of pressure decreases to 2.5% of its peak value. This method for determining the shock amplitude A sh has been proposed earlier and described in detail in [ 10 ], [ 32 ], and [ 33 ]. Focal values of the peak positive pressure p + , peak negative pressure p - , and the shock amplitude A sh are listed in the three last columns of Tables I and II for each array. Shock fronts in the focal waveform starts to form at lower acoustic power (about 175 W) for the V1 array than for the V2 array (about 220 W). In addition, the Sonalleve V2 system produces shocks with amplitudes that are about 10-15 MPa lower than the V1 version when considering the same power output level above 350 W ( Fig.6 ). The amplitude of the developed shock is 83 MPa for the V2 system (power of 515 W) and is 96 MPa for the V1 one (power of 525 W). Thus, substantially higher power levels are required for the V2 system to reach the same shock-wave exposure conditions in treatments like boiling histotripsy [ 5 , 6 , 15 ]. For completeness, a final comparison of the shapes and dimensions of the nonlinear focal regions for V1 and V2 arrays was made. Such a comparison of simulations along focal axes is depicted in Fig.8 for the peak positive pressure p + and peak negative pressure p - . Dimensions of the corresponding focal regions and focusing gains p F + / p 0 and p F - / p 0 are listed in Table III for increasing values of the power. Nonlinear propagation effects lead to smaller peak-positive and larger peak-negative focal regions than in the linearly focused beam ( Fig. 8 , Table III ). In addition, axial and transverse dimensions of the focal region for the peak positive pressure p + change nonmonotonically [ 34 ]. Initially, increasing the source power leads to a decrease in the p + focal area reaching a minimum when shock formation occurs at the focus. Further increase in power leads to the formation of shocks in a larger area around the focus and saturation effect begins close to the focus, which causes p + focal dimensions to grow. Unlike the peak positive pressure p + , the dimensions of the focal area of peak negative pressure p - change monotonically with array output power. More specifically, dimensions in all directions x , y , and z increase with output power within the nulls of the pressure and the focal maximum slightly moves towards the transducer ( Fig. 8 , Table III ). Focusing gains of the peak positive pressure p F + / p 0 and the peak negative pressure p F - / p 0 change in a similar way as the dimensions of their focal regions. Initially, the peak positive pressure focusing efficiency p F + / p 0 increases with the array power due to more efficient focusing of higher harmonics generated in the beam and differences in their relative diffraction phase shifts [ 35 , 36 ]. The maximum value of the focusing gain p F + / p 0 is 251 at about 300 W for the V1 array and 209 at about 350 W for the V2 array. Note that these maximum focusing gains for p + are 3.7 and 4.1 times higher than the linear focusing gains, respectively. After formation of the shock front, focusing efficiency p F + / p 0 drops due to the absorption of the wave energy at the shocks that occurs prefocally. Focusing gain for the peak negative pressure p F - / p 0 changes monotonically and decreases by a factor of 1.8 relative to the linear focusing gain over the operating range of output levels for each array. A visual representation of the shapes of p + and p - focal regions is provided in Figs. 9 and 10 for the power level at which developed shocks form at the focus. The p + distributions are very narrow in the transverse directions x and y ( Fig. 9 ). At the −6 dB level, the size of the p + focal area along x , y , and z axes is only 0.6 mm × 0.6 mm × 6.1 mm for the V1 array and 0.6 mm × 0.6 mm × 7.7 mm for the V2 array. In addition, the drop in p + to 90% of the maximum occurs in a region of about only 0.2 mm in the transverse x and y directions. Note that such a small width of the focal area in the transverse directions is comparable to the size of the FOPH hydrophone tip in experiments. Thus, the peak positive pressure p + measured by hydrophone can be influenced by an averaging effect and is very sensitive to accurate positioning of the tip. Consequently, peak positive pressures are sometimes underestimated in hydrophone measurements of nonlinear fields at very high source output levels [ 16 ]. Focal regions of the peak negative p - pressures are tear-drop shaped and are significantly larger than the corresponding sizes of the peak positive focal regions ( Figs. 9 , 10 ). At the −6 dB level, the size of the p - focal area along x , y , and z axes is 1.9 mm×1.9 mm × 11.0 mm for the V1 array and 2.2 mm × 2.2 mm × 15.7 mm for the V2 array. For both V1 and V2 arrays, the locus of the peak value is slightly shifted (about 1 mm) toward the array ( Fig.8 c , g ; Fig.10 ). Nominal parameters of the Sonalleve V1 and V2 arrays and their equivalent sources with and without the central opening are listed in Table IV . For each equivalent source represented by an annular spherical segment, the diameter of the central opening was chosen to match that of the physical array: 20 mm for V1 and 44 mm for V2. Accordingly, only two geometric parameters (focal distance and outer diameter) were fitted using the analytical solution ( 2 ) applied to normalized pressure levels. Similar to the previous studies, the focal distances of equivalent sources appeared to be slightly longer than the focal distances of the arrays to better match the focal lobe of their fields [ 9 ]. The outer diameters of equivalent sources differed from those of the physical arrays by no more than 3%. Active surfaces of both equivalent sources turned out to be about 1.7 times greater than the nominal surface of the 256 active elements of the arrays, inversely proportional to their filling factors (62.8% for V1 and 56.5% for V2). Acoustic power of each equivalent source was chosen to match the same linear pressure amplitude at the focus and characterized by the power coefficients listed in Table IV . Each coefficient is the ratio of the array’s acoustic power calculated and scaled from the holography measurements to the power of the corresponding equivalent source. For example, an acoustic power of 100 W of the equivalent source for the V1 array matches the output of the actual V1 array at 134 W. As shown in Fig. 11 , the equivalent sources accurately reproduce corresponding actual axial pressure distributions over the main focal lobe under linear propagation conditions. In this regard, equivalent sources with and without a central opening provide the same performance. Simulations of nonlinear fields generated by each equivalent source were performed using the open software “HIFU beam” over the entire range of output levels. Comparisons of the simulated peak positive, p + , and peak negative, p - , focal pressures are presented in Fig. 12 . The acoustic power displayed in these plots refers to the corresponding power of the actual array ( i.e ., the displayed power uses the relevant coefficient from Table IV to recalculate power from that used for the equivalent source). For the Sonalleve V1 system, inclusion of the central opening in the equivalent source has no significant effect on focal peak pressures in the nonlinear field ( Fig. 12a ). In contrast, inclusion of a central opening for the V2 system is crucial in accurately capturing nonlinear behaviors with an equivalent source. More specifically, the presence of a central opening leads to later saturation of the focal peak positive pressure p + and a 15% higher saturation level ( Fig. 12b ). This example demonstrates that equivalent sources may need to include a central opening when representing transducers with a relatively large one. In this case, the equivalent source can be readily adapted by simply including a central opening of the same size as the original transducer. Note that the equivalent sources with central openings (in the form of an annular spherical segment) provided accuracy of 1% for focal peak positive pressure p + ( Fig. 12 ). Such good agreement allows specialists working with Sonalleve MR-HIFU systems to use the equivalent source parameters identified here to accurately simulate the expected nonlinear acoustic fields in water or layered tissue using a tool such as the “HIFU beam” software described previously. Specific instructions for conducting such simulations are provided in the Appendix .

The transducer arrays in both the Sonalleve V1 and V2 systems comprise 256 circular elements arranged on a spherical surface. The elements are 6.6 mm in diameter with a 1.2 MHz operating frequency. Based on the design of each array, different element locations for each system are illustrated as projections onto a flat surface in the plots in the top row of Fig. 1 . The V1 and V2 geometries are described respectively by 127.8 and 135.9 mm apertures; 120 and 140 mm radii of curvature; and formally calculated F-numbers F # (i.e., ratio of radius of curvature to the aperture) of 0.94 and 1.03. Thus, the V2 array is less focused than the V1 version. Additionally, the newer V2 array has elements located in eight symmetrical sectors with a significantly larger central opening (about 44 mm in diameter, compared to about 20 mm for the V1). The same Sonalleve driving electronics was used for both arrays. The system software allows detailed control for the magnitude and phase of each element [ 13 ]. For the present study, we only consider the driving configuration in which all elements are driven in phase with the same amplitude. Accordingly, the acoustic beam is not electronically steered so that the focus remains on the geometric axis of the transducer. Although the system typically delivers clinical treatments with intensity correction to account for electronic steering and power feedback control to ensure consistent output levels, these features were disabled for the present characterization efforts. The output level for each array was controlled in the system software using an “ampval” setting that roughly corresponds to the voltage applied to the elements. Each array was calibrated at the factory such that every ampval setting corresponds to an electric and acoustic power level. Although this calibration is helpful, it was not used here. Instead, each ampval setting was directly calibrated by measuring a corresponding near-source pressure level (in the linear propagation regime). This approach inherently calibrates the output levels of the relatively short measurement pulses (on the order of tens of cycles) used for source characterization in the present study (see Section II-C ). In contrast, the factory calibration is based on much longer pulses. The measurements used to characterize the V1 array were presented in detail in Ref. [ 7 ]. A fully analogous approach with the patient table again outside the magnet room was used for the V2 array. To accommodate the different patient table associated with the V2 system, a slightly modified experimental arrangement was utilized as depicted in the bottom row of Fig. 1 . In both cases, the array is positioned in an oil bath inside the patient table with an acoustic window above the transducer provided by a 50 μ m plastic membrane. For the V2 array, the top of the table is further separated from the transducer by a layer of degassed water (i.e., a “water disc”). In both cases, a cylindrical acrylic tank was mounted to the tabletop to hold water in which hydrophone measurements could be performed. For the V1 arrangement, the tank was threaded into the tabletop and sealed with an O-ring. For the V2 arrangement, the bottom of the tank was sealed with its own 50 μ m membrane, and this entire assembly was placed on top of a ledge surrounding the tabletop membrane. Altogether, the propagation path into the measurement tank for the V1 system involved a single membrane separating the transducer in its oil bath from the water bath used for measurements. For the V2 system, the path involved three membranes between the array and the measurement bath: one at the bottom of the water disc, one at the top of the water disc (this membrane serves as the exterior of the patient table in normal use), and a final membrane used to seal the bottom of the acrylic tank. Note that a very thin layer of water was used to provide coupling between the tabletop membrane and the membrane at the bottom of the tank. Despite the more complicated path for the V2 system, the plastic membranes are designed to have minimal impact on the acoustic propagation. As for conditions inside the tank, both tanks held a similar volume of water with an inside diameter of 184 mm and a water depth of at least 230 mm. For the V1 array, the geometric focus was located about 100 mm above the bounding membrane with the array moved vertically up 17.5 mm from the “home” position defined by system software. Similarly, the focus of the V2 array was about 70 mm above the membrane with the array moved vertically up 15 mm from the “home” position. This tank geometry provided hydrophone access for measuring the ultrasound field proximal to the focus. Finally, we note that in both cases the water was degassed to about 10% of saturation and maintained at temperatures from 21–25°C. For the V2 array, there is potential for additional refraction to occur (beyond that at the oil-water interface) if a temperature mismatch exists between the tank and the water disc. To minimize this potential, active cooling of the water disc was disabled so that the disc temperature remained at about 20°C during measurements. In comparing measurements for the V2 array, any changes in refraction at this interface were neglected. Using the test configuration depicted in Fig. 1 , new hydrophone measurements for the V2 array were acquired using a custom LabVIEW program (National Instruments Corp., Austin, TX). Key aspects of the acquisition include synchronized movement of the hydrophone using a 3-D positioner (VXM stepper motor controllers and Unislide linear positioners, Velmex Inc., Bloomfield, NY); triggering of the driving electronics using a function generator (Model HP33120A, Keysight Technologies, Santa Rosa, CA); and capturing of the hydrophone signal using a PC-based digitizer (Gage Compuscope CSE1422, Vitrek, LLC, Poway, CA). Measured waveforms were later processed in Matlab (The MathWorks Inc., Natick, MA). Source characterization measurements followed the same approach described for the V1 array in Ref. [ 7 ]. With this approach, hydrophone measurements at a low output level are first conducted to characterize the linear acoustic field. These measurements include a 2-D holography scan in a prefocal plane as well as independent measurements in the focal region to validate the recorded hologram. This hologram defines the pattern of vibrations of the source as a boundary condition for modeling. To complete the boundary condition, the amplitude of this vibration pattern must be scaled as a function of the source output level. Accordingly, additional measurements at a single near-source point (where the field remains nearly linear) are made over a range of output levels. For both the V1 and V2 arrays, the low-amplitude measurements were acquired using a capsule hydrophone with a nominal aperture of 200 μ m (Model HGL-0200 with AH-2020 preamplifier; Onda Corp., Sunnyvale, CA). For the V2 array, the holography scan was made in a plane transverse to the beam axis at a distance 40 mm proximal to the focus. The scan aperture was 80.4 × 80.4 mm with a step size of 0.6 mm. When triggered at each scan point, the Sonalleve driving electronics were programmed to deliver an 80-cycle pulse at an amplitude of 430 ampvals (50 W nominal acoustic power based on the factory calibration). This scan was designed to provide a time window (beginning at 102.6 μ s and lasting 10 cycles) over which the recorded waveforms could be analyzed to define a steady-state hologram in terms of the pressure magnitude and phase. To complement the measured hologram by calibrating a full range of output levels, near-source measurements were made at a point on-axis, 40 mm proximal to the focus. These measurements utilized 80-cycle pulses at output levels from 87 to 2859 ampvals (i.e., nominal acoustic powers from 5 to 900 W based on the factory calibration). In addition to the holography and near-source measurements used to define boundary conditions to the modeling, additional measurements at the focus were conducted over a full range of output levels to validate the results of nonlinear simulations. Measurements of focal waveforms were conducted with a fiber optic probe hydrophone (FOPH) (Model FOPH 2000; RP Acoustics, Leutenbach, Germany), which utilizes a 100 μ m diameter optical fiber and has a nominal 100 MHz bandwidth. To minimize deflection of the fiber tip, all waveforms measured with the FOPH were acquired with the fiber approximately parallel to the ultrasound beam and later deconvolved based on impulse-response data provided by the manufacturer. The overall source characterization approach described in Ref. [ 7 ] and above in Section II-C uses both measurement-based simulations and independent validation measurements. In order to compare such simulations with independent measurements, it is instructive to note two areas that can pose challenges: (1) misalignment of the experimental and theoretical coordinate systems, and (2) amplitude calibration of the simulation boundary conditions relative to the independent validation measurements. Regarding the alignment of measurement and modeling coordinates, we note that the measured hologram is designed to capture the entire 3-D field. Hence, even if the beam axis is not exactly perpendicular to the scan plane, no error is introduced. Here we accounted for any misalignment during the step in which the measured hologram is backprojected to define a source hologram as a boundary condition for modeling (see Section II-E ). More specifically, the backprojection reconstructs the field in a plane at an axial distance that corresponds to the apex of the physical transducer. Then, this plane is adapted to modeling coordinates by rotating and centering it to ensure that it is perpendicular to the true axis of the beam, with peak pressures at the focus remaining on axis. Simulated axial pressures are then compared to FOPH measurements of focal waveforms. Because nonlinear beams can be focused to a very small spot and the location of this spot generally shifts along the beam axis at different amplitudes, care must be taken in identifying the measurement location in theoretical coordinates. Here we selected the measurement position for each array by finding the location of peak positive pressure at an elevated output level with nonlinear focusing. For the V1 array, this output level was 820 ampvals (152 W acoustic by factory calibration); for the V2 array, the selected output level was 873 ampvals (150 W acoustic by factory calibration). A second challenge in comparing measurement-based simulations with independent FOPH measurements pertains to the calibration of output levels. As described in Ref. [ 7 ], fiber optic hydrophones can be calibrated to absolute pressures through well-known relations: (a) the Gladstone–Dale equation describing the optical index of refraction in water as a function of density, and (b) the Tait equation of state to relate density and pressure in water. Accordingly, we accept the FOPH measurements to be calibrated (with some associated uncertainty). In contrast, we take the holography measurements to be initially uncalibrated. Even though the 200 μ m capsule hydrophone used for these measurements has a calibrated sensitivity value at 1.2 MHz, this value neglects the impact of directivity. As reported in separate work [ 14 ], a hologram measured with a similar 200 μ m hydrophone for a comparable focused transducer operating at 1.5 MHz underestimated the beam’s true power by about 25%. Although a 25% underestimate in power is nontrivial and cannot be fully corrected without relevant directivity data for the hydrophone, we have found that these directivity effects do not have much impact on the structure of the field near the focus. Consequently, the approach used for both V1 and V2 arrays is to determine an effective hydrophone sensitivity such that the model boundary conditions based on holography accurately represent the true acoustic power as measured under quasilinear conditions by the FOPH. This sensitivity is then used to consistently scale all simulation boundary conditions based on the near-source measurements made across all output levels. In this way, we establish consistent boundary conditions for simulations to evaluate the ability of the model to quantify nonlinear waveform distortion and shock formation as output levels increase. This approach meets the goals of the current study. It would be possible to take a different approach in which directivity measurements are made to characterize the capsule hydrophone and more accurately define the amplitudes of measured holograms. Then, independent uncertainties (both for simulation boundary conditions derived from holograms and FOPH measurements) could be considered in comparing simulation results with validation measurements. The nonlinear acoustic fields generated in water by V1 and V2 arrays at increasing output levels were modeled based on the one-directional version of the Westervelt equation. The 3-D model includes diffraction, nonlinearity, and thermoviscous absorption and has been shown to accurately represent the nonlinear acoustic fields generated by different types of HIFU transducers [ 15 , 16 ]. Further details of the model are described in Refs. 7 , 15 , 17 . For completeness, we include here a brief description of the model and the numerical algorithm used for its implementation. Using a retarded time coordinate τ , we write the Westervelt equation to describe forward propagation: (1) ∂ 2 p ∂ z ∂ τ = c 0 2 Δ p + β 2 ρ 0 c 0 3 ∂ 2 p 2 ∂ τ 2 + δ 2 c 0 3 ∂ 3 p ∂ τ 3 . Here p is the acoustic pressure and Δ p denotes the Laplace operator acting on p over three spatial coordinates x , y , and z . As shown in Fig. 1 , z is parallel to the beam axis while x and y denote transverse coordinates. The retarded time is defined relative to time t as τ = t - z / c 0 , where c 0 is the speed of sound. Other acoustic parameters ρ 0 , β , and δ are the density, nonlinear parameter for the propagation medium, and sound diffusivity, respectively. In the modeling, the values of these parameters ( c 0 = 1485 m / s , ρ 0 = 997 kg / m 3 , β = 3.5 , δ = 4.33 × 10 - 6 m 2 / s ) were set to correspond to the experimental conditions in water at room temperature. A model boundary condition was defined at the apex plane of the array ( z = 0 ) as a pressure distribution based on holography measurements [ 18 ]. With this approach, the angular spectrum method was used to linearly backpropagate the field represented by the aligned hologram [ 19 , 20 , 21 ]. Simulations based on the Westervelt equation (1) were performed at increasing output levels and the results were compared to direct pressure measurements at the focus made with a fiber optic hydrophone. The presence of oil surrounding the array (see Fig. 1 ) was not explicitly included in simulations of nonlinear forward propagation. This approach is reasonable because nonlinear propagation effects occurred almost entirely in water (in or near the focal region). Moreover, because the simulation boundary conditions were defined from a hologram measured in water, refraction at the oil-membrane-water interface was implicitly accounted for. A method of fractional steps with an operator splitting procedure of second-order accuracy was used to solve the Westervelt equation (1) [ 22 ]. For each propagation step along the beam axis, the splitting procedure was implemented by dividing Eq. (1) into several simpler equations that separately govern diffraction, nonlinearity, and absorption behaviors. Both time-domain and frequency-domain representations of the pressure field were used in the numerical solution. At shorter distances, in the near field of the array, where the shock fronts are not yet formed, a frequency-domain approach was employed. As the degree of nonlinear waveform distortion increased and more harmonics were required, the numerical algorithm automatically switched to a shock-capturing time-domain Godunov-type scheme [ 23 ]. The switch to the time-domain scheme was performed when waveform steepness reached a threshold such that the amplitude of the 10th harmonic exceeded 1% of the harmonic amplitude at the fundamental frequency. Parameters of the numerical scheme were set as follows: the axial step Δ z varied from 0.4 mm in the near field to 0.1 mm in the focal region of the beam; the transverse step sizes were set at Δ x = Δ y = 0.02 mm ; the maximum number of harmonics included in the calculations was N max = 1000 . An equivalent source method is based on the idea that a single-element piston source with simple, axisymmetric geometry (either flat or spherically curved) can generate the same nonlinear acoustic field in the focal region as more complicated real sources. The basis for this approach relies on the condition that nonlinear effects related to the real source are most pronounced in the focal region; consequently, an equivalent source with matching behavior under linear propagation conditions will accurately describe the corresponding nonlinear field at higher output levels [ 11 ]. This approach is appealing because the computational burden for simulating nonlinear fields is much smaller for the equivalent axially symmetric source. In addition, the equivalent source model makes it possible to estimate nonlinear acoustic field parameters not only when focusing at the geometrical focus but also when steering off-axis. It has been shown that the peak positive and peak negative pressures, as well as the shock amplitude would be the same in the steered off-axis focus as in the geometrical focus. Higher power of the array is required in this case to compensate for the effect of steering and it should be scaled the same way as in the linear focusing conditions [ 24 ]. Here, for each of the V1 and V2 arrays we consider two possible equivalent source models: a single-element, spherically curved bowl that vibrates uniformly and an annulus of such bowl with a central opening. The bowl-shaped equivalent source with uniform vibrational velocity amplitude u 0 on its surface is represented by a boundary condition defined by its aperture D , radius of curvature F , frequency f , and nominal pressure amplitude p 0 char = ρ 0 c 0 u 0 . For the annular equivalent source, the diameter d of the central opening is added to parameters listed above in defining the representative boundary condition. The simplest way to identify suitable parameters of an equivalent source is to use the on-axis analytical solution of the Rayleigh integral for a uniformly vibrating, spherically curved source [ 18 ]: (2) p ( z , t ) = p 0 char e - i ω t 1 - z / F e ikz - e ikR max , where i is the imaginary unit, z is the axial distance from the source apex to the observation point, and R max is the distance from the observation point to the edge of the source: (3) R max = F 1 + ( 1 - z / F ) 2 - 2 ( 1 - z / F ) cos arcsin D 2 F . The real part of expression ( 2 ) gives the distribution of the pressure amplitude on the beam axis, which can be compared with that of the real source. To obtain the solution for the equivalent source with a central opening, the equation (2) can be used twice: first to calculate the field with no opening as above and then again to subtract the field of a source with aperture D equal to the diameter d of the central opening. A method that implements this approach for defining equivalent sources has been proposed and validated first for the case of a strongly focused ( F # = 0.9 ) single-element spherical bowl transducer without a central opening ( d = 0 ) [ 25 , 26 ]. The method then was proven to be applicable for accurate simulations of nonlinear fields generated by multi-element focused transducers, including a V1 system that has small central opening [ 11 ]. The effect of the central opening on linear focused fields has been studied analytically [ 27 , 28 ]. It was demonstrated that increasing the size of the central opening leads to a shift of the field maximum toward the transducer, a decrease in the number of on-axis lobes, and an elongation of the focal region along the axis in conjunction with transverse narrowing. However, the influence of including a central opening on the nonlinear fields generated by equivalent sources has not yet been studied. The utility in defining accurate equivalent sources for the Sonalleve V1 and V2 arrays lies in the potential for a wide range of users to conduct relatively simple simulations of nonlinear acoustic fields generated by the arrays at different output levels. The equivalent source model is axially symmetric (i.e., two-dimensional), and can be used as a boundary condition in a freely available software tool “HIFU beam” [ 12 ] (link is given in reference [ 29 ]). This software tool is designed for simulating high-intensity focused ultrasound fields generated by single-element transducers and annular arrays with propagation in water or in flat-layered media that mimic biological tissues [ 12 ]. The software uses shock-capturing methods that allow for simulating strongly nonlinear acoustic fields with high-amplitude shocks. Here, the “HIFU beam” software with equivalent-source boundary conditions corresponding to the V1 and V2 arrays was used to simulate nonlinear acoustic fields in water at different output levels. The software was run in wide-angle parabolic approximation mode (“WAPE”) of solving the Westervelt equation for a one-layered propagation medium (water), using the physical properties for water as indicated in the previous section. Accordingly, simulations included nonlinear effects and thermoviscous absorption, while power-law absorption effects typical for biological tissues were disabled. For this problem statement, the simulator solves the one-way Westervelt equation with radial symmetry, which can be written in the retarded time coordinate system as follows: (4) ∂ 2 p ∂ τ ∂ z = c 0 2 ∂ 2 p ∂ z 2 + ∂ 2 p ∂ r 2 + 1 r ∂ p ∂ r + β 2 ρ 0 c 0 3 ∂ 2 p 2 ∂ τ 2 + δ 2 c 0 3 ∂ 3 p ∂ τ 3 . Eq. (4) is a 2-D version of Eq. (1) and it takes into account the same physical effects. Parameters of the numerical grid set in the “HIFU beam” simulator were as follows: radial step Δ r = 0.025 mm , axial step Δ z = 0.025 mm , maximal number of harmonics N max = 1000 .

A comparative characterization of nonlinear acoustic fields generated by the V1 and V2 Sonalleve therapeutic arrays is presented. The characterization of both arrays was performed using a combination of hydrophone measurements and numerical modeling. This approach uses acoustic holography measurements of the linear field in order to set a boundary condition to the 3-D nonlinear numerical model based on the Westervelt equation. Nonlinear simulations were carried out for a wide range of acoustic powers, and results were validated by comparison with independent FOPH measurements. Comparative calibration analysis demonstrates that at the same acoustic power, the V2 array generates 10-15 MPa lower shock amplitudes at the focus compared to the V1 array. This difference is caused mainly by a smaller focusing angle of the V2 array. Note, however, that the V2 array has a larger central opening than the V1 array and the presence of a larger opening leads to higher shock amplitudes. If the V2 array had smaller central opening comparable to V1, the resulting shock amplitude would be even smaller. Formation of a developed shock at the focus occurs at approximately the same acoustic power (about 520 W) for both arrays. At the same acoustic power, the size of focal area for the peak pressures the V2 array is larger in the axial direction than for the V1 array; however, it is almost identical in transverse directions. Weaker focusing of the V2 array leads to widening of the focal lobe, but the presence of large central opening, as has been shown previously, leads to its narrowing [ 28 ], all together resulting in the same transverse size of the focal area of the V1 and V2 arrays. As one can see from the Table III , in nonlinear propagation the tendency remains the same: the transverse dimensions of the focal areas for the peak pressures are fairly similar with a difference of less than 10% for both arrays at the same power. Summing up everything above, the V2 array produces an acoustic field with lower peak pressures and longer focal region than the V1 array. Axisymmetric equivalent source models in a form of a spherical segment were constructed for both arrays. The importance of defining equivalent sources with a central opening that matches the actual array was demonstrated. Based on features of the acoustic fields of the V1 and V2 arrays, the V1 version may be preferable in applications requiring very high-amplitude shocks such as boiling histotripsy. The V2 version may be preferable in applications where the size of the focal region is important (e.g., thermal heating in clinical applications with quasilinear fields). Note that the acoustic powers reported here are based on measured holograms that utilize short pulses; these powers are somewhat different from nominal power values based on factory calibrations. Based on our experience, the difference between acoustic powers based on measured holograms and nominal system powers can be up to 20%. Reconciling these different values would depend on consideration of various measurement uncertainties and is beyond the scope of the current study, which is focused on understanding the basic capabilities of two different therapeutic arrays.

MR-HIFU is being used for various noninvasive therapeutic applications [ 1 , 2 ]. The main bioeffect in the clinical use of MR-HIFU technology is thermal ablation caused by absorption of acoustic energy and, as a result, heating and thermally coagulating targeted tissue. The Sonalleve MR-HIFU system (Profound Medical, Canada) bears the CE marking for treating uterine fibroids, adenomyosis, desmoid tumors, osteoid osteoma, and bone metastases as well as FDA approval for treating osteoid osteoma [ 3 , 4 ]. The standard treatment procedures utilize acoustic powers from about 100-300 W, which corresponds to quasilinear or weakly nonlinear ultrasound wave distortion at the focus. However, the technical characteristics of the Sonalleve system make it possible to deliver much higher acoustic pulses with powers up to 1000 W. At such power levels, shocks fronts of up to 100 MPa develop in the focal waveform, which allows use of the Sonalleve MR-HIFU platform for therapies requiring the presence of high-amplitude shocks. One such therapy is boiling histotripsy that uses millisecond-long HIFU pulses with shocks to mechanically emulsify tissue [ 5 , 6 ]. The Sonalleve platform has two versions of the therapeutic transducer array. Though both V1 and V2 versions comprise 256 array elements, each version has a different element pattern arranged around a central opening of a different size [ 7 , 8 ]. Nonlinear acoustic field characterization of the Sonalleve system is needed for the design of shock-wave treatments like boiling histotripsy and is also required for accompanying regulatory approvals. Characterization of the V1 Sonalleve system has been performed previously in water for the complete range of available output levels using a combination of acoustic holography measurements and nonlinear modeling [ 7 ]. Hydrophone field measurements of the V2 array have been reported for output levels up to half (500 W) of the available system power [ 8 , 9 ]. In this study, a comparative characterization of the nonlinear acoustic fields generated by the V1 and V2 Sonalleve systems is presented. For the V1 array, the data from previous calibration work are used [ 7 ]. For the V2 array, new measurements and modeling were performed to comprehensively calibrate the system in the same way that was done for the V1. A comparative analysis of the V1 and V2 arrays reveals important ways in which the different array geometries impact shock amplitudes at the focus. The influence of the central opening on nonlinear effects is also considered. Nonlinear effects are most prominent when waves propagate in the same or nearly the same direction, i.e. with small angles relative to the transducer axis in focused beams. It has been shown that shocks form at lower focal pressures and therefore have smaller amplitudes for less focused transducers [ 10 ]. The presence of the central opening excludes these coherently interacting waves to yield a somewhat higher effective focusing angle, higher pressures at shock formation, and thus higher peak positive pressures and shock amplitudes at the focus. To complete the study, we define axisymmetric equivalent sources for both arrays. The equivalent source model is based on the premise that the nonlinear effects are mostly concentrated in the focal region of the HIFU beam [ 11 ]. Thus, matching the focal acoustic fields for real and equivalent sources at low power, in the linear propagation regime, leads to the same focal acoustic fields at proportionally higher power levels, in nonlinear propagation regimes [ 11 ]. Here we also demonstrate that the central opening must be taken into account when constructing an equivalent source model for the V2 system with larger opening. The accuracy of equivalent source models for the V1 and V2 systems is validated. Accordingly, various specialists who work with the Sonalleve MR-HIFU systems can use these simplified source models to accurately and efficiently simulate the nonlinear acoustic fields generated with different system settings. In particular, the equivalent source boundary conditions are readily implemented with the freely available “HIFU beam” software, which allows simulation of nonlinear acoustic fields in water or in a flat-layered medium imitating biological tissues [ 12 ].