Ultrasound Axicon: Systematic Approach to Optimize Focusing Resolution through Human Skull Bone

The use of axicon lenses is useful in many high-resolution-focused ultrasound applications, such as mapping, detection, and have recently been extended to ultrasonic brain therapies. However, in order to achieve high spatial resolution with an axicon lens, it is necessary to adjust the separation, called stand-off (δ), between a conventional transducer and the lens attached to it. Comprehensive ultrasound simulations, using the open-source k-Wave toolbox, were performed for an axicon lens attached to a piezo-disc type transducer with a radius of 14 mm, and a frequency of about 0.5 MHz, that is within the range of optimal frequencies for transcranial transmission. The materials properties were measured, and the lens geometry was modelled. Hydrophone measurements were performed through a human skull phantom. We obtained an initial easygoing design model for the lens angle and optimal stand-off using relatively simple formulas. The skull is not an obstacle for focusing of ultrasound with optimized axicon lenses that achieve an identical resolution to spherical transducers, but with the advantage that the focusing distance is shortened. An adequate stand-off improves the lateral resolution of the acoustic beam by approximately 50%. The approach proposed provides an effective way of designing polydimethylsiloxane (PDMS)-based axicon lenses equipped transducers.


Introduction
Focused ultrasound (FUS) in pulsed mode is unique among transcranial brain stimulation methods in combining exceptional spatial resolution (on the millimeter scale) [1][2][3] with the potential to target subcortical structures (deeper than 10 cm) [4] through the intact skull. It also has potential for inducing neuronal excitation or suppression without evidence of tissue damage [5].
Recently, we demonstrated the advantages of focusing ultrasound (US) through polydimethylsiloxane (PDMS)-based axicon lenses to selectively drive brain activity [6]. The ultrasound axicon is shown in Figure 1. It has the shape of a cone. As the cone angle (φ) decreases, the focus moves closer to the lens. Developing low-intensity applications include the opening of the blood-brain barrier and ultrasonic neuromodulation. Both techniques have recently been extended to human subjects and are under active research. Spherically FUS transducer has been, until now, the most commonly used for transcranial focusing ultrasound, but it has a large focal length (several centimeters) that may hinder its coupling with the head, for example, when the focal plane is too close to skull inner ivory layer. With spherical segment ultrasound transducers, FUS is commonly delivered through a big plastic bag containing degassed water placed over the scalp. This is due to the fact that, for cortical stimulation, the fact that, for cortical stimulation, the acoustic beam should be focused a few millimeters deep from the skull surface [3]. In this sense, the great advantage of the axicon lens is its ability to suppress the near-field and maintain a very near focus from the lens face outward. PDMS-based axicon lens affixed on the face of conventional transducers makes it possible to build more compact devices without liquids that can regasify or leak, with a greater spatial resolution.
However, for an axicon lens to offer high spatial resolution and depth, control is necessary to adjust the separation between the transducer and axicon lens, called stand-off (δ), for each particular case [7]. Given the geometry of the lens, if δ is not properly adjusted, internal reflections can occur, making the configuration useless as a focused transducer. The lens reduces the focal length (F) from the near-field distance (N) of the attached US transducer. The relation between axicon lens angle and the value of F/N produced was described in [8] for acrylic plastic/oil combination, but the effect of the stand-off was not modeled and the adequate value of δ was not described elsewhere. In [7], the stand-off of the different settings was adjusted experimentally to obtain the best signal-to-noise ratio (SNR) for ultrasonic contact inspection. In our approach, time domain ultrasound simulations, based on the k-space pseudo-spectral methods, play a key role in enabling the modeling and systematic design of ultrasonic axicon lenses with optimum stand-off for high-resolution ultrasonic applications, as transcranial stimulation, in high-resolution mapping, and in the evaluation of a wide variety of defects. Finally, the effects of human skull on axicon fields were tested. The polydimethylsiloxane (PDMS) fills the conical cavity of the lens and the stand-off. An ultrasonic lens is formed by the PDMS/plastic interfaces. The window is a thin material layer between the conical cavity and the face of the lens; (b) schematic ultrasonic beam pattern for transducers with axicon lenses. These are described by its total included cone angle (φ). Depth of focus (DOF) is the focal region and F is the desired focal length. Both DOF and F depend on the sound velocity into the material used for the interface.
Skull is generally constituted of three relatively homogeneous layers: The outer and inner tables of solid ivory bone, and the central layer of diploe of cancellous bone, with a blood-and fat-filled porous structure. The dimensions of the blood-and fat-filled inclusions are random, and an average thickness in the direction normal to the surface is about 0.6 mm. [5], so the fundamental frequency used was chosen to help alleviate the concerns for acoustic energy absorption or refraction by the skull. For frequencies lower than about 0.5 MHz, reflection of sound is the principal cause of insertion loss. At frequencies between about 0.6 and 0.9 MHz, the absorption loss in the diploe layer begins to limit sound transmission, so the oscillations are damped out and the insertion loss increases linearly with frequency. At about 0.9 MHz, the scattering loss in the diploe layer begins to limit sound transmission, so the insertion loss begins to increase as the fourth power of frequency [5].
A full-wave nonlinear ultrasound model based on the k-space pseudo-spectral method was developed and released as part of the open-source k-Wave Acoustics Toolbox [9,10]. This model can account for the propagation of nonlinear ultrasound waves in homogeneous or heterogeneous media The window is a thin material layer between the conical cavity and the face of the lens; (b) schematic ultrasonic beam pattern for transducers with axicon lenses. These are described by its total included cone angle (φ). Depth of focus (DOF) is the focal region and F is the desired focal length. Both DOF and F depend on the sound velocity into the material used for the interface.
Skull is generally constituted of three relatively homogeneous layers: The outer and inner tables of solid ivory bone, and the central layer of diploe of cancellous bone, with a blood-and fat-filled porous structure. The dimensions of the blood-and fat-filled inclusions are random, and an average thickness in the direction normal to the surface is about 0.6 mm. [5], so the fundamental frequency used was chosen to help alleviate the concerns for acoustic energy absorption or refraction by the skull. For frequencies lower than about 0.5 MHz, reflection of sound is the principal cause of insertion loss. At frequencies between about 0.6 and 0.9 MHz, the absorption loss in the diploe layer begins to limit sound transmission, so the oscillations are damped out and the insertion loss increases linearly with frequency. At about 0.9 MHz, the scattering loss in the diploe layer begins to limit sound transmission, so the insertion loss begins to increase as the fourth power of frequency [5].
A full-wave nonlinear ultrasound model based on the k-space pseudo-spectral method was developed and released as part of the open-source k-Wave Acoustics Toolbox [9,10]. This model can account for the propagation of nonlinear ultrasound waves in homogeneous or heterogeneous media with power acoustic absorption and without restrictions on the directionality of the waves. The accuracy of the implementation of nonlinear ultrasound model in the k-Wave Toolbox was validated using experimental measurements of the ultrasound made with a linear diagnostic ultrasound probe and a membrane hydrophone [11].
A computational model for elastic wave propagation in heterogeneous media can be constructed based on the solution of coupled first-order acoustic equations given in Equations (1)-(3) using the Fourier pseudo-spectral method. This uses the Fourier collocation spectral method to compute spatial derivatives, and a leapfrog finite-difference scheme to integrate forwards in time. Using a temporally and spatially staggered grid, the field variables are updated in a time stepping [12,13]. To simulate free-field conditions, a perfectly matched layer (PML) is also applied to absorb the waves at the edge of the computational domain [14]. Without this boundary layer, the computation of the spatial derivate via the FFT causes waves leaving one side of the domain to reappear at the opposite side. The use of the PML thus facilitates infinite domain simulations without the need to increase the size of the computational grid.
where u is the acoustic particle velocity, ρ 0 is the medium density, ρ is the acoustic density, c is the thermodynamic sound speed, p is the acoustic pressure, and p 0 = p(t = 0) is the initial pressure distribution. There are two main stages in this work: The definition of general design equations for PDMS-based axicon lenses with optimal stand-off, through comprehensive ultrasound simulations, for an optimal transcranial focusing; and the measurements and simulations performed to determine focusing performance of the proposed lenses, through a human skull phantom.

Materials and Methods
In order to develop a design model, axicon lenses with φ angles between 80 • and 170 • in steps of 5 • were simulated with increments of δ in steps of λ/4 with 8 grid points per wavelength (λ) in the stand-off medium. The final pressure field along with the RMS beam pattern was calculated. The normalized cross-section profile area at the focus F was determined, for each δ. The minimum area represents the greatest improvement with the optimum stand-off, which minimizes transmitted energy outside of the main beam and improves lateral spatial resolution with lower possible sidelobes.
The domain was discretized using a grid point spacing of 250 µm (giving a maximum supported frequency of 2.06 MHz), and a grid size of 512 × 512 grid points (corresponding to a domain size of 128 × 128 mm). Simulations were run on an NVIDIA ® GTX 950 graphics processing unit (Santa Clara, CA, United States) using the MATLAB ® Parallel Computing Toolbox (Natick, MA, United States). The simulation of all angles/stand-off combinations can be completed in approximately 60 h. By default, numbers in MATLAB ® are stored in double precision. However, in almost all cases, k-Wave does not require this level of precision. In particular, the performance of the PML generally limits the accuracy to around 4 or 5 decimal places. We use a PML thickness of 20 grid points that gives a transmission coefficient of −100 dB. This corresponds to a reduction in the signal level of 1 × 10 −5 , which is significantly less than double precision. Further, there will also be uncertainties in the definition of the materials properties. A list of the main simulation inputs is given in Table A1. In this work, a heterogeneous medium was defined as a layered interface on both sides of the conical cavity of the axicon lens as shown in Figure A1. The convective nonlinear effects from the convection of mass was considered. However, at low frequencies and amplitudes, nonlinearity will only have a small effect on the wave field. At higher frequencies and amplitudes, this effect become more important.
The accuracy of the implementation of ultrasound model with the k-Wave Toolbox was validated in our previous work [6] using experimental measurements of FUS made with the same axicon lens attached transducer and a needle hydrophone (Force Technology MH28) within a 6-L anechoic test tank. Our previous study has already shown that there is a good agreement between the simulated and experimental beam patterns. In this work, we also characterized the acoustic pressure amplitude of the beam pattern of the axicon lens when FUS was transmitted through a human skull phantom for experimental validation of simulated transcranial ultrasound propagation. There is negligible conversion to shear waves in the layers of skull when the incidence angle is within about 20 • of normal [5]. The ability of bone to support shear waves can affect transcranial transmission, although the changes to the intracranial field are typically negligible for ultrasound applied at normal or near-normal incidence. Therefore, we will model only longitudinal waves. The phantom was created from the parietal portion of a mesh segmented from MRI head image data. Clear Med610 3D printing resin was used to create the skull bone phantom. The acoustic properties of Stratasys™ materials were recently reported in [15]; thus, these measurements were not repeated as part of the current work. The reported and measured property values, and estimated uncertainty in those measurements, are shown in Table 1.
To test the effects of a human skull on FUS fields, we inserted a 5 mm thick fragment of parietal bone phantom between the transducer and the hydrophone, as shown in Figure 2. The transducer described in our previous work [6] has an ultrasonic piezo-disc-type element of 28 mm diameter (SMD28T21F1000R, Steminc Steiner & Martins, Inc., Davenport, FL, USA) of PZT-4 mounted on stainless-steel housing operating in thickness mode vibration at 445 kHz. Epoxy (Resoltech 1040, Resoltech, Rousset, France) resin was used in order to build the conical cavity of the lens with an angle φ of 144 • . Degassed PDMS (Sylgard 184, Dow Corning, Midland, MI, USA) was used to fill the conical cavity of the lens and for the lens-transducer interface with a stand-off δ of 30 mm. The ultrasonic lens is formed by the epoxy/PDMS interface. The specifications are summarized in Table 2. small effect on the wave field. At higher frequencies and amplitudes, this effect become more important.
The accuracy of the implementation of ultrasound model with the k-Wave Toolbox was validated in our previous work [6] using experimental measurements of FUS made with the same axicon lens attached transducer and a needle hydrophone (Force Technology MH28) within a 6-L anechoic test tank. Our previous study has already shown that there is a good agreement between the simulated and experimental beam patterns. In this work, we also characterized the acoustic pressure amplitude of the beam pattern of the axicon lens when FUS was transmitted through a human skull phantom for experimental validation of simulated transcranial ultrasound propagation. There is negligible conversion to shear waves in the layers of skull when the incidence angle is within about 20° of normal [5]. The ability of bone to support shear waves can affect transcranial transmission, although the changes to the intracranial field are typically negligible for ultrasound applied at normal or near-normal incidence. Therefore, we will model only longitudinal waves. The phantom was created from the parietal portion of a mesh segmented from MRI head image data. Clear Med610 3D printing resin was used to create the skull bone phantom. The acoustic properties of Stratasys™ materials were recently reported in [15]; thus, these measurements were not repeated as part of the current work. The reported and measured property values, and estimated uncertainty in those measurements, are shown in Table 1.
To test the effects of a human skull on FUS fields, we inserted a 5 mm thick fragment of parietal bone phantom between the transducer and the hydrophone, as shown in Figure 2. The transducer described in our previous work [6] has an ultrasonic piezo-disc-type element of 28 mm diameter (SMD28T21F1000R, Steminc Steiner & Martins, Inc., Davenport, FL, USA) of PZT-4 mounted on stainless-steel housing operating in thickness mode vibration at 445 kHz. Epoxy (Resoltech 1040, Resoltech, Rousset, France) resin was used in order to build the conical cavity of the lens with an angle φ of 144°. Degassed PDMS (Sylgard 184, Dow Corning, Midland, MI, USA) was used to fill the conical cavity of the lens and for the lens-transducer interface with a stand-off δ of 30 mm. The ultrasonic lens is formed by the epoxy/PDMS interface. The specifications are summarized in Table  2.

Estimation of the Axicon Lens Angle
The relation between axicon lens angle φ and the value of F/N produced (see Appendix A) is illustrated graphically in Figure 3. The transducer near field length N is given by N = D 2 f /4c. The coefficient of determination R 2 is 0.97 with p-value significance level of < 0.00001. This relation is for δ = δ Optimum :

Estimation of the Axicon Lens Angle
The relation between axicon lens angle φ and the value of F/N produced (see Appendix A) is illustrated graphically in Figure 3. The transducer near field length N is given by N = D 2 f/4c. The coefficient of determination R 2 is 0.97 with p-value significance level of < 0.00001. This relation is for δ = δOptimum: The following relation, based on our study summarized in Table A2 (see Appendix A), was found experimentally valid for the lens described. The optimum ratio of F/N appears to be between 0.1 and 0.3. This produces focal beam diameters (dF with the axicon lens equipped transducers and dN of a conventional transducer), and depths of focus (DOF) of similar ratio according to: The following relation, based on our study summarized in Table A2 (see Appendix A), was found experimentally valid for the lens described. The optimum ratio of F/N appears to be between 0.1 and 0.3. This produces focal beam diameters (d F with the axicon lens equipped transducers and d N of a conventional transducer), and depths of focus (DOF) of similar ratio according to: For values of F/N > 0.4, some evidence of the original near field still remains. As the value F/N decreases below 0.4, all evidence of the original near field is rapidly suppressed in the lens system. This contrasts with the behavior of spherical lenses where some original near field is always present.
We observe that there is an inversely proportional relationship between F/N and κδ values, as shown in Figure 4 for different values of angle φ. The relative percentage loss of lateral resolution (LLR) as a function of κδ is shown in the same figure, where κ is the wave number. LLR is relative to the best lateral resolution that can be achieved with κδ value that minimizes energy transmitted outside of the narrowest possible main beam. Since the speed of sound in PDMS is lower than in water, as the value of δ increases, the wavelengths-weighted average sound speed through the heterogeneous medium with step discontinuity of velocity decreases, which effectively reduces the ratio of F/N.

Reflectivity Effect on the Internal Walls of the Lens Housing
To illustrate the effect of outer case and inner isolation ( Figure 5) on the lateral resolution of the lens, the relative percentage LLR as a function of κδ, and the normalized lateral beam profile are shown in Figure 6 for two different housing. With a reflecting housing, without inner sleeve, there will be more internal reflections in the lens system, since there are abrupt transitions of acoustic impedance with the inner walls. As a result, the lateral resolution of the lens decreases, compared to non-reflective housing with inner isolation. For example, with a PDMS-based lens at 445 kHz and a cone angle of 130 • , assembled with the optimum stand-off inside a housing of Inconel-625 with internal reflectivity of 0.8 dB down, the relative LLR value is reduced by half and energy transmitted outside the main beam is 10% higher, compared to inner sleeve reflectivity of 40 dB down.

Estimation of the Optimum Stand-Off
We found a numerical relationship for δ based on the simulation results of 2420 angle/stand-off combinations for axicon lenses (see Appendix A). The relation between the optimum value of δ (which improves spatial resolution) and the value of F/N is shown in Figure 7. The linear regression equation is given by: The coefficient of determination R 2 is 0.96 with a p-value significance level of < 0.00001.
For values of F/N > 0.4, some evidence of the original near field still remains. As the value F/N decreases below 0.4, all evidence of the original near field is rapidly suppressed in the lens system. This contrasts with the behavior of spherical lenses where some original near field is always present.
We observe that there is an inversely proportional relationship between F/N and κδ values, as shown in Figure 4 for different values of angle φ. The relative percentage loss of lateral resolution (LLR) as a function of κδ is shown in the same figure, where κ is the wave number. LLR is relative to the best lateral resolution that can be achieved with κδ value that minimizes energy transmitted outside of the narrowest possible main beam. Since the speed of sound in PDMS is lower than in water, as the value of δ increases, the wavelengths-weighted average sound speed through the heterogeneous medium with step discontinuity of velocity decreases, which effectively reduces the ratio of F/N.

Reflectivity Effect on the Internal Walls of the Lens Housing
To illustrate the effect of outer case and inner isolation ( Figure 5) on the lateral resolution of the lens, the relative percentage LLR as a function of κδ, and the normalized lateral beam profile are shown in Figure 6 for two different housing. With a reflecting housing, without inner sleeve, there will be more internal reflections in the lens system, since there are abrupt transitions of acoustic impedance with the inner walls. As a result, the lateral resolution of the lens decreases, compared to non-reflective housing with inner isolation. For example, with a PDMS-based lens at 445 kHz and a cone angle of 130°, assembled with the optimum stand-off inside a housing of Inconel-625 with internal reflectivity of 0.8 dB down, the relative LLR value is reduced by half and energy transmitted outside the main beam is 10% higher, compared to inner sleeve reflectivity of 40 dB down.

Estimation of the Optimum Stand-Off
We found a numerical relationship for δ based on the simulation results of 2420 angle/stand-off combinations for axicon lenses (see Appendix A). The relation between the optimum value of δ (which improves spatial resolution) and the value of F/N is shown in Figure 7. The linear regression equation is given by: The coefficient of determination R 2 is 0.96 with a p-value significance level of < 0.00001.         As an example, Figure 8 shows the focusing behavior of PDMS-based 144° axicon lens with four different values of δ (a: 0 mm, b: 20.5 mm, c: 34 mm, and d: 45.25 mm). These different settings are indicated in Figure 9 showing the relative LLR as a function of κδ. With a stand-off of 34 mm, the lateral spatial resolution improves by up to 40%, compared to the same lens without stand-off. The optimum stand-off predicted by Equation (6) was checked for different frequencies. Figure  10    As an example, Figure 8 shows the focusing behavior of PDMS-based 144 • axicon lens with four different values of δ (a: 0 mm, b: 20.5 mm, c: 34 mm, and d: 45.25 mm). These different settings are indicated in Figure 9 showing the relative LLR as a function of κδ. With a stand-off of 34 mm, the lateral spatial resolution improves by up to 40%, compared to the same lens without stand-off.
The optimum stand-off predicted by Equation (6) was checked for different frequencies. Figure  As an example, Figure 8 shows the focusing behavior of PDMS-based 144° axicon lens with four different values of δ (a: 0 mm, b: 20.5 mm, c: 34 mm, and d: 45.25 mm). These different settings are indicated in Figure 9 showing the relative LLR as a function of κδ. With a stand-off of 34 mm, the lateral spatial resolution improves by up to 40%, compared to the same lens without stand-off. The optimum stand-off predicted by Equation (6) was checked for different frequencies. Figure  10

Acoustic-Field Experimental Scan
The hydrophone scans were performed both without and through the phantom. For the scan with the skull, the starting distance to the transducer was increased to 10 mm to avoid collision between the skull and hydrophone. Experimental beam patterns produced are shown in Figure 11. The lateral dimension of FUS beam cross-profiles measured at −6 dB drop of the pressure at the focus was 3.5 mm in the free space condition and 5 mm after transcranial transmission. We also characterized the acoustic field in the axial direction, perpendicular to the lens face and skull. The FUS pressure half width of the half maximum was 22 mm in the free space condition and 18 mm after transcranial transmission. Under these conditions, transmission of 445 kHz FUS-axicon lens through the skull led to an approximately 40% loss in lateral resolution of the acoustic beam, and on the other hand, an approximately 18% increase in the axial resolution. When FUS was transmitted through the skull, the acoustic pressure dropped by half. The insertion loss of our skull phantom was approximately −6 dB.
Intracranial focal characteristics obtained by simulation using the same configuration of Table 2 for different thicknesses of the skull are indicated in Table 3. There is a good coincidence between the simulation and the values obtained experimentally in water for the skull thickness of 5 mm, separated 2 mm from the lens, as shown in Table 4.

Acoustic-Field Experimental Scan
The hydrophone scans were performed both without and through the phantom. For the scan with the skull, the starting distance to the transducer was increased to 10 mm to avoid collision between the skull and hydrophone. Experimental beam patterns produced are shown in Figure 11. The lateral dimension of FUS beam cross-profiles measured at −6 dB drop of the pressure at the focus was 3.5 mm in the free space condition and 5 mm after transcranial transmission. We also characterized the acoustic field in the axial direction, perpendicular to the lens face and skull. The FUS pressure half width of the half maximum was 22 mm in the free space condition and 18 mm after transcranial transmission. Under these conditions, transmission of 445 kHz FUS-axicon lens through the skull led to an approximately 40% loss in lateral resolution of the acoustic beam, and on the other hand, an approximately 18% increase in the axial resolution. When FUS was transmitted through the skull, the acoustic pressure dropped by half. The insertion loss of our skull phantom was approximately −6 dB.
Intracranial focal characteristics obtained by simulation using the same configuration of Table 2 for different thicknesses of the skull are indicated in Table 3. There is a good coincidence between the simulation and the values obtained experimentally in water for the skull thickness of 5 mm, separated 2 mm from the lens, as shown in Table 4.

Acoustic-Field Experimental Scan
The hydrophone scans were performed both without and through the phantom. For the scan with the skull, the starting distance to the transducer was increased to 10 mm to avoid collision between the skull and hydrophone. Experimental beam patterns produced are shown in Figure 11. The lateral dimension of FUS beam cross-profiles measured at −6 dB drop of the pressure at the focus was 3.5 mm in the free space condition and 5 mm after transcranial transmission. We also characterized the acoustic field in the axial direction, perpendicular to the lens face and skull. The FUS pressure half width of the half maximum was 22 mm in the free space condition and 18 mm after transcranial transmission. Under these conditions, transmission of 445 kHz FUS-axicon lens through the skull led to an approximately 40% loss in lateral resolution of the acoustic beam, and on the other hand, an approximately 18% increase in the axial resolution. When FUS was transmitted through the skull, the acoustic pressure dropped by half. The insertion loss of our skull phantom was approximately −6 dB.
Intracranial focal characteristics obtained by simulation using the same configuration of Table 2 for different thicknesses of the skull are indicated in Table 3. There is a good coincidence between the simulation and the values obtained experimentally in water for the skull thickness of 5 mm, separated 2 mm from the lens, as shown in Table 4.

Discussion
Although the analysis was carried out with a single element transducer with no aberration correction, and a specific skull geometry designed to approximate the varied shaped of the skull, it is expected that the relative influence of different medium properties and aspects of medium geometry will be maintained. Phantom geometry is the major material influence on the intracranial field and sound speed is shown to be the most influential acoustic property in focus pressure, position, and volume. From the experimental beam patterns shown in Figure 11, the unexpected focusing property of the skull in axial axis may be described as a nonlinear effect that causes the beam to rotate back toward the skull insertion point, creating a more compact pressure cigar-shaped acoustic field. This was also observed using segmented-sphere transducers [16,17]. Thus, the skull is not an obstacle for transcranial focusing of US and may exert an additional acoustic lensing effect to enhance spatial resolution under certain conditions. Out of focus, the sound pressure decreases with a very steep slope. From the comparison in Table 4, there is a good coincidence between the simulation and the values obtained experimentally in water for skull thickness of 5 mm, separated 2 mm from the lens. Regarding the influence of the thickness skull in the focus of axicon lenses, since on both sides of the skull there is a discontinuity in the acoustic impedance, with different velocities of sound propagation, and the average sound speed is the acoustic property that most influences the focal distance, we observe as expected, small variations in the position, diameter, and depth of focus, for different thicknesses of the cranial bone, resulting in the average deviations of the focus less than 1 mm.
On the other hand, transmitting FUS through human cranial bone caused an approximately 40% loss in lateral resolution of the acoustic beam, estimated by the intensity full width at half maximum. However, this loss of resolution is compensated by adequate stand-off of axicon lens. In addition, we

Discussion
Although the analysis was carried out with a single element transducer with no aberration correction, and a specific skull geometry designed to approximate the varied shaped of the skull, it is expected that the relative influence of different medium properties and aspects of medium geometry will be maintained. Phantom geometry is the major material influence on the intracranial field and sound speed is shown to be the most influential acoustic property in focus pressure, position, and volume. From the experimental beam patterns shown in Figure 11, the unexpected focusing property of the skull in axial axis may be described as a nonlinear effect that causes the beam to rotate back toward the skull insertion point, creating a more compact pressure cigar-shaped acoustic field. This was also observed using segmented-sphere transducers [16,17]. Thus, the skull is not an obstacle for transcranial focusing of US and may exert an additional acoustic lensing effect to enhance spatial resolution under certain conditions. Out of focus, the sound pressure decreases with a very steep slope. From the comparison in Table 4, there is a good coincidence between the simulation and the values obtained experimentally in water for skull thickness of 5 mm, separated 2 mm from the lens. Regarding the influence of the thickness skull in the focus of axicon lenses, since on both sides of the skull there is a discontinuity in the acoustic impedance, with different velocities of sound propagation, and the average sound speed is the acoustic property that most influences the focal distance, we observe as expected, small variations in the position, diameter, and depth of focus, for different thicknesses of the cranial bone, resulting in the average deviations of the focus less than 1 mm.
On the other hand, transmitting FUS through human cranial bone caused an approximately 40% loss in lateral resolution of the acoustic beam, estimated by the intensity full width at half maximum. However, this loss of resolution is compensated by adequate stand-off of axicon lens. In addition, we find that the PDMS provides a smaller focal zone, which is desirable for neurostimulation [6]. All this allows a higher resolution, comparable to spherical transducers (the most commonly used for ultrasonic brain therapy), but with the advantage that the near field is eliminated, and the focus distance is shortened. For other applications, a wider focal zone and a line focus such as that obtained with glycerin or ethylene glycol may be desirable [7]. From Equation (5), the elimination of the near field, by means of an appropriate axicon lens, enables transducers featuring the same wavelength/diameter ratio to produce the same focal spot size. Thus, large diameter, low-frequency transducers may be used. This is useful for brain stimulation where low frequencies are required for penetration of the skull.
One problem of devices with axicon lenses are the relatively high sidelobes [18]. How much this will affect will depend on the proposed applications. With the calculation of the optimum value of δ, by Equation (6), a better lateral resolution is achieved. This relation between axicon lens stand-off and the value of F/N is applicable to high-resolution epoxy resin/PDMS lenses or other similar combination.

Conclusions
The numerical approach proposed in this paper provides a complete and effective way of designing axicon lenses for many high-resolution applications, such as mapping or detection. This is also suitable for focused ultrasound through human skull bone. In view of providing better lateral resolution with lower sidelobes, the use of design programs for this task is not as straightforward. The choice of a good starting point is an important factor for successful optimization. It is easier to obtain a starting design using relatively simple formulas and then use it in a lens design program for future analysis and optimization. We believe that this will be an effective way of designing axicon lenses, for example, to build focused windows to the brain for clinically viable transparent cranial implant for chronic ultrasonic therapy and stimulation of the brain.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A k-Wave Simulations
The simulations are based on the k-space pseudo-spectral method and implemented in MATLAB ® with the open-source k-Wave toolbox. The functions based on the coupled first-order acoustic equations (named kspaceFirstOrder2D) are named with four input structures (kgrid, medium, source, and sensor). These structures define the properties of the computational grid, the distribution of medium properties, stress and velocity source terms, and the locations of the sensor points used to record the evolution of the wave field over time. The propagation of the wave field is computed step by step in a 2-D layered medium, with the pressure values at the sensor elements stored after each iteration. These values are returned when the time loop has completed. A list of the main simulation inputs is given in Table A1. One of the advantages of k-Wave is that the spatial gradients are calculated by using FFTs rather than using a finite-difference stencil. This means that for linear simulations, only two points per wavelength are required (Nyquist). For nonlinear simulations, the number of points per wavelength at the fundamental frequency will depend on the highest frequency of interest. The amplitudes of the harmonics should decay smoothly. Initially, the computational grid is defined using the function makeGrid. The time steps used in the simulation are defined by the object property kgrid.t_array. The time array was defined using the function makeTime, which calculate within the simulation functions using the time taken to travel across the longest grid diagonal at the slowest sound speed, and a Courant After the computational grid, the properties of the propagation medium shown in Figure A1 are defined by the objects medium.sound_speed and medium.density. Also, the object medium.BonA represents the nonlinear properties of the medium. With medium.BonA = 0, k-Wave will include convective nonlinear effects in the model, but not include the nonlinear relationship between the acoustic pressure and acoustic density. If medium.BonA is undefined, k-Wave instead solves linearized equations.
The parameters medium.alpha_coeff (a 0 ) and medium.alpha_power (y) describe the power law acoustic attenuation in the medium, where the attenuation is of the form a = a 0 × fˆy. The power law absorption and acoustic nonlinearity of water is specified by: display_mask = source.p_mask; % assign the input options input_args = {'DisplayMask', display_mask, 'PMLInside', false,… 'DataCast', 'gpuArray-single', 'PlotPML', true, 'PlotLayout', true}; % run the simulation sensor_data = kspaceFirstOrder2D(kgrid, medium, source, sensor,… input_args{:}); The computed focal length (F) and optimum stand-off (δ) values for different lens angles (φ) of epoxy and PDMS axicon, with a transducer near-field distance in the water of 68.2 mm, are shown in Table A2. Table A2. Values, obtained by running a 512 × 512 simulation, of focal length (F), ratio of F/N, optimum stand-off (δOptimum), focus diameter (dF), and maximum improvement of lateral resolution (MILR), for epoxy/PDMS-515 kHz ∅28mm axicon lens with different angles. The computed focal length (F) and optimum stand-off (δ) values for different lens angles (φ) of epoxy and PDMS axicon, with a transducer near-field distance in the water of 68.2 mm, are shown in Table A2. Table A2. Values, obtained by running a 512 × 512 simulation, of focal length (F), ratio of F/N, optimum stand-off (δ Optimum ), focus diameter (d F ), and maximum improvement of lateral resolution (MILR), for epoxy/PDMS-515 kHz ∅28mm axicon lens with different angles. to ultrasound frequency used for 445 kHz is 3.37 mm in water, so that for this frequency, the dimensions of the skull inclusions are smaller than one-half the wavelength and, therefore, ultrasound will not be severely scattered by these inclusions.

Appendix C Material Characterization
Clear Med610 material was used to create the skull bone phantom described in Section 2 and Appendix B. The acoustic properties of Clear Med610 are equal to those of the VeroBlack that were reported in [15]; thus, these measurements were not repeated. However, independent measurements of the sound speed and attenuation of epoxy resin and PDMS were conducted and are shown in Table A3. To obtain reliable simulations, it is very important to accurately know the propagation velocity of ultrasound through the materials used for the lens. For materials characterization, both were poured into cylindrical molds and left to set at 25 • C for 48 h. The velocity of sound in epoxy resin and PDMS samples was determined by time of flight technique using an ultrasonic echoscope Digital-Echograph 1090 of Karl Deutsch. This measurement was performed in reflection mode with 2 MHz probe Karl Deutsch S6WB2.25. The sound speed in these materials has a weak dependence on frequency, less than 1% measured from 500 kHz to 2 MHz. The Attenuation of ultrasound was determined in reflection at the applied frequency of 2 MHz. The dimensions and weight of the test samples were measured by a micrometer and a digital scale, respectively.

Appendix D Guideline for Axicon Lens Design
This guide is based on design Equations (4)-(6) for Epoxy/PDMS materials combination.
(1) Calculate transducer near field, as: N = D 2 f 4c where D is transducer diameter, f is transducer frequency, and c is sound velocity of the material under inspection (2) Select the desired lens focus F and check that: Values of 0.1 are rarely used because it gives a very near focus. Values or 0.4 give a profile similar to the transducer but remove the N zone.