**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.


**Table A1.** Summary of the simulation inputs for relationships estimation of the axicon lens angle and the optimum stand-off.

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–Friedrichs–Lewy (CFL) number of 0.1, where CFL = *c*0∆*t*/∆*x*. The computational grid is defined by:

% size of the computational grid

Nx = 512; % number of grid points in the x (row) direction

x = 128e-3; % size of the domain in the x direction [m]

dx = x/Nx; % grid point spacing in the x direction [m]

% create the computational grid

kgrid = makeGrid(Nx, dx, Nx, dx);

% create the time array

[kgrid.t\_array, dt] = makeTime(kgrid, medium.sound\_speed);

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*<sup>0</sup> × *fˆy*. The power law absorption and acoustic nonlinearity of water is specified by:

medium.alpha\_power = 2; %[dB/(MHzˆy cm)] medium.alpha\_coeff = 2.17e-3; %[dB/(MHzˆy cm)] medium.BonA = 4.96;

Next, a time varying pressure source is defined by assigning a binary source mask to source.p\_mask (which defines the position of the source points) along with a time varying source input to source.p. Here, a single sinusoidal time series is used to drive the transducer element. The pressure source is defined by:

```
% define a transducer element
source.p_mask = makeLine(Nx, Ny, startpoint, pi/2, 112);
% define a time varying sinusoidal source
source_freq = 0.515e6; % [Hz]
source_mag = 1.0; % [Pa]
source.p = source_mag*sin(2*pi*source_freq*kgrid.t_array);
% filter the source to remove any high frequencies not supported %by the grid
source.p = filterTimeSeries(kgrid, medium, source.p);
```
Finally, to visualize the acoustic beam produced by the axicon lens, a sensor mask covering the entire computational domain was defined. Only the total beam pattern is required, thus at each time step, k-Wave only updates the maximum and r.m.s. values of the pressure at each sensor point by setting sensor.record to {'p\_final', 'p\_max', 'p\_rms'}. The sensor mask is given here:

% create a sensor mask covering the entire computational domain % using the opposing corners of a rectangle sensor.mask = ones(Nx, Ny); % set the record mode to capture the final wave-field and the % statistics at each sensor point sensor.record = {'p\_final', 'p\_max', 'p\_rms'}; When the input structures have been defined, the simulation is started by passing them

to kspaceFirstOrder2D.

The simulation is invoked by:

% create a display mask to display the transducer

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{:});

**Figure A1.** (**a**) Sound speed mask for the different mediums; (**b**) density mask for the different mediums.

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.

