Generating arbitrary ultrasound fields with tailored optoacoustic surface profiles

Recently, the generation of arbitrary acoustic fields has been the subject of increasing interest in physical acoustics. Several works have proposed the use of multi-foci or continuously patterned fields for applications in particle manipulation,¹ ultrasound therapy,² neural stimulation,³ and non-destructive testing.⁴ Generating these fields with conventional piezoelectric arrays poses a significant challenge, particularly at high frequencies. This has led to a number of works exploring novel alternatives including sparse 2-D piezoelectric arrays,³ acoustic kinoforms created using metamaterials,⁵ or profiles fabricated using three-dimensional (3-D) printing techniques.⁶ 

The optical generation of ultrasound (OGUS) using the photoacoustic effect is particularly well suited for the generation of these patterned fields. The lack of electrical connections allows a “transducer” to be deposited on any surface that can be illuminated. The pattern of the absorbing surface can be tailored to generate almost any chosen wavefield when it is excited by a short optical pulse (Fig. 1).

Several previous works have investigated tailored profiles for the generation of specific fields. Concave lenses have been used to generate single high pressure foci for cell cleaving and histotripsy.⁷ Axicon surfaces have been used to create non-diffracting acoustic fields for imaging,⁸ and circular ramps have been used to generate rotating acoustic fields for use in particle manipulation.⁹ Recently, Chan et al. showed that a stepped surface profile can be used to generate two pulses at a single focal point with a single optical pulse. The same work demonstrated that 3-D printing with transparent materials is well suited to the creation of these surface profiles.¹⁰ 

The question of how to generate a desired acoustic field using different configurations of piezoelectric arrays has been addressed via optimisation in a number of previous works.^11,12 However, an approach has yet to be developed for the design of optoacoustic surface profiles. In this work, an optimisation approach using random downhill binary search¹³ is presented for the flexible design of surface profiles for the generation of arbitrary OGUS fields.

In free space, the acoustic signal p(r, t) at a point $r=(x,y,z)∈ℝ3$ and time t due to a photoacoustic source distribution p₀(r)δ(t) can be written as¹⁴ 

where G is the free-space Green's function for the wave equation, c₀ is the sound speed, and p₀ is compactly supported in $Ω⊂ℝ3$⁠. In this paper, Ω is restricted to a thin region of thickness δ around a surface defined as follows. Consider a square region of the xy plane, [−L/2, L/2] × [−L/2, L/2], tiled by N² square pixels of side Δ = L/N centred on the coordinates (x_i, y_i). Each square (i, j) is now displaced in the z-direction to a height h_ij. The surface around which p₀ is defined is the discontinuous surface formed by the set of squares centred at (x_i, y_i, h_ij). For practical reasons, the surface is restricted to those squares, which lie in a circle of radius L/2.

The inverse problem of finding a source p₀ that will generate any arbitrarily defined time-varying field p(r, t) is severely ill-posed, especially when p₀ is restricted to a surface. To overcome this, here the inverse problem was restricted to finding a set of heights h_ij that generates a large amplitude pressure field at just M discrete focal points r_m. The time at which the large amplitude occurred was allowed to vary. The surface height was also discretised such that $hij=nΔ, n∈ℕ$⁠, for ease of fabrication. The inverse problem was solved using an optimisation approach, with the cost function defined as

Here, $av(pmax(r))$ and $sd(pmax(r))$ are the average and standard deviation of the temporal peak pressures at the focal points r_m, α and β are weighting terms, and S is a term penalising large variations in height between neighbouring pixels. The value of p_max (r) was calculated from Eq. (1) using a ray-based approach, with the surface source approximated by a set of point sources of amplitude A located at $(xi,yj,hij),p0≈∑i,jAijδ(x−xi,y−yi,z−hij)$⁠. The amplitudes A_ij were used to account for the spatial profile of the illumination.

The first term in Eq. (2) rewards increases in peak pressure at each of the target points. The second penalises variations in peak pressure between different focal points. The third term penalises large variations in h_ij between neighbouring squares, thereby acting as an additional constraint on the source. This was introduced as large local variations in height might lead to structural instabilities, uneven absorber deposition, and unwanted reflections. The smoothness term S is given by a second order accurate central difference approximation of the Laplacian

The two weights α and β in Eq. (2) are assigned heuristically by calculating profiles for a range of values of α and β and simulating the fields they generate, with the backing material, using the k-Wave toolbox.¹⁵ In this work, α was set to 1 and β was assigned, after initialising the surface, such that $βS=av(pmax(r))$⁠.

The optimisation was initialised by choosing $hij(0)$ in the following way. First, a set of spherical sections (bowls) centered on the focal points r_m were defined as the surfaces $bm={(x−xm)2+(y−ym)2+(z−zm)2=zm2 | x2+y2≤L2/4}$⁠. For this step, if the foci were distributed over multiple planes, the area occupied by each bowl was adjusted starting from the edge to ensure that the initial pressure at each target point was approximately the same. Second, the union of these surfaces, $∪{bm}$⁠, was sampled at the source points (x_i, y_j) and rounded to the nearest nΔ to give $hijmulti$⁠, which is multivalued with up to M possible surface heights for each source pixel (i, j). In the third step of the initialisation, one of these values must be chosen for each source pixel, as the absorbing surface must be single-valued due to the fact the light is all absorbed by the first absorbing surface it encounters. The selection process involved taking each pixel (i, j) in turn, at random but not repeating any pixel, and considering a 5 × 5 pixel region around it. If any pixels in this region have already been assigned a single-valued height, the value of $hijmulti$ chosen was that which minimised the average difference to these surrounding heights, with one condition. It was found, in simulation, that large planar regions with the same value of h_ij generated a high peak pressure at locations away from the targeted focal points. To ameliorate this, only height differences with the neighbouring pixels of >3Δ were considered. When multiple values of $hijmulti$ had equal difference to the surrounding heights, or when no pixels in the surrounding region had yet been assigned single-valued heights, the value of $hijmulti$ that minimised the sum of the first two terms in Eq. (2) was chosen. The result was a single-valued initial surface profile, $hij(0)$⁠.

Following the initialisation, iterative refinements to the surface heights h_ij were obtained by minimising Eq. (2), within a specified range of heights, for each pixel sequentially in random order. The range evaluated for each position was the set of heights for which p(r_m, t_m) ≠ 0 for at least one m, where t_m is the time at which the peak occurs. This process gave the surface profile $hij(1)$⁠; it was then iterated to give the converging series $hij(k)$⁠. The iterations were terminated when the number of pixels changed from the previous iteration fell below a threshold, which was set at 0.5% of the total pixels on the surface. Figures 2(a) and 2(b) illustrate the initial (Fig. 2(a)) and final (Fig. 2(b)) states of a surface profile designed to generate 4 foci at different depths. The acoustic field generated by this profile when illuminated by a single optical pulse may be seen in Figs. 2(c) and 2(d). This was simulated using the k-Wave toolbox¹⁵ and clearly demonstrates the generation of four discrete foci.

One alternative way to choose the surface profile would be to form an array of M bowl transmitters each focussed on one focal point r_m. However, such an approach would suffer from considerably worse lateral focussing due to the solid angle subtended by each bowl at its focal point being smaller, by a factor of about M, than the solid angle available in the optimisation approach. In the optimisation approach, any pixel can contribute to the focus at any point. Additionally, it allows the relative peak pressure generated at each focal point to be controlled. The relative time at which the peak pressure occurs can also be controlled by adjusting the cost function to optimise the pressure at select time points. This comes at the expense of the absolute peak pressure generated at each point. This optimisation approach therefore has several useful features for the design of surface profiles for the generation of arbitrary fields.

To demonstrate that the optimisation approach is effective, samples were fabricated via additive manufacture. This was done using a high-resolution polyjet printer (Objet350 Connex, Stratasys, Eden Prairie, MN, USA). The substrate was printed using a transparent material (VeroClear-RGD810) that is rigid and nearly colourless. The absorbing surface was formed through aerosol deposition of an optically absorbing polymer composite (Gloss Super, Plastikote, Valspar, USA).

An example surface profile is shown in Fig. 3(b). This was calculated for a circular aperture of diameter 23 mm discretised using a spacing Δ of 200 μm. The target pattern was a continuous distribution of high pressure resembling the numeral “7.” This was chosen as it is asymmetric and continuous and therefore represents a challenging target distribution. This makes it a good test of the limits of the inversion.^5,16 This “7” had a height of 12 mm and a width of 6 mm and was placed in a plane 25 mm from the bottom of the surface profile. The “7” was discretised using 61 evenly spaced points for input into the optimisation. It was assumed that the laser pulse incident on the surface had a top-hat profile. The forward model in Eq. (1) was solved accounting for frequencies up to 6 MHz. To prepare the sample from the calculated surface profile, the backing material was added to each position up to the calculated height. A further 2 mm was added across the bottom of the surface to provide rigidity. The backing material is optically transparent and acoustically absorbing so this layer can be increased in thickness arbitrarily to eliminate the backward travelling wave. The surface mesh and the 3-D printed transparent sample can be seen in Figs. 3(a) and 3(b).

The field generated by this surface profile was characterised experimentally. A fibre coupled Q-switched Nd:YAG laser (Minilite, Continuum, San Jose, CA, USA) was used to excite ultrasound from the fabricated sample. The optical pulses had a wavelength of 1064 nm, a pulse length of 6 ns, pulse energies of 14 mJ, and an approximately Gaussian profile. Measurements of the acoustic field were made using a previously described optical scanner.¹⁷ This uses a planar Fabry-Perot thin film interferometer for acoustic detection over a 2-D plane.

A schematic of the set-up used is shown in Fig. 4. The fibre was attached to a circular tube with a length of 10 cm and a diameter of 2.5 cm. The natural divergence of the laser beam expanded the beam to 2.5 cm by the end of the tube. The samples were attached to a retaining ring, which was connected to the other end of the tube. The front face of the samples was then immersed in water approximately 1 cm above the optical scanner. Signals were recorded over a 22.20 × 22.05 mm area using a step size of 0.15 mm. Four averages were taken at each position. The 3-D wavefield was calculated from the experimental measurements using linear acoustic holography. Prior to this, the data were band-pass filtered between 30 kHz and 6 MHz and were spatially upsampled by a factor of 2. The filter was applied to remove low frequency noise introduced by the detection system and to limit the measurement data to the design bandwidth.

The experimental measurements were compared against a simulation carried out using the k-Wave toolbox. The simulation grid was 288 × 288 × 384 grid points with a spacing of 0.1 mm. The calculated h_ij was upsampled by a factor of 2 in (x, y) such that each position (i, j) was represented by 4 pixels and then inserted as a source mask in the simulation. The medium for the simulation was inhomogeneous. The backing material VeroClear was inserted as a region with a sound speed of 2495 ms⁻¹ and a density of 1190 kg m⁻³. The rest of the medium was set to a sound speed of 1500 ms⁻¹ and a density of 1000 kg m⁻³. A bipolar pulse filtered at 6 MHz was used to represent the acoustic signal generated from the sample. The spatial profile of the laser was approximated by adding a 2-D Gaussian apodisation to the source. This had a standard deviation of 5 mm.

The peak pressure in the target plane of the surface for both the experimental data and the simulation may be seen in Fig. 5. In both cases, the numeral 7 has been clearly realised and the two show excellent agreement in the overall shape of the pattern and in the low background pressure. The focal gain of the surface profile in the experimental data is 2. This gave pressures of approximately 60 kPa across the target pattern. While not a focus of this work, it has been shown previously that the pressure can be increased up to MPa levels with appropriate design of the optical absorber.¹⁰ The signal to noise ratio (SNR) in the target plane of the experimental data is 4.1. This was quantified by creating a mask encompassing the target pattern and calculating the ratio of the average peak pressure within the mask to outside it. The time window between the first and last occurrence of a peak across the target pattern was 2.6 μs. There are however noticeable differences in the amplitude generated across the 7 in both datasets. These could be corrected by measurement of the apodisation A_ij in the experimental set-up for input into both the optimisation and the simulation.

Several factors determine the complexity of the acoustic fields that can be generated and the focal gain and SNR that can be achieved. The surface profile consists of a set of partial bowl transducers each centred on a target point and taking up a fraction of the aperture. The focal gain is therefore fundamentally limited by the aperture size, the centre frequency, and the total area occupied by the target pattern. However, the optimisation allows certain positions on the aperture to contribute to the peak pressure at multiple points so the precise scaling varies with the distribution of target points and the bandwidth of the acoustic signal. In addition to these factors, the performance also improves with smaller pixel size and the sizes of features that can be realised in the pattern are diffraction limited.

This paper has demonstrated that complex patterned distributions of pressure can be generated using tailored optoacoustic surface profiles and a single optical pulse. An optimisation approach for the design of these profiles for the generation of arbitrary distributions of pressure has been introduced. It has also been shown that the calculated profiles can be cheaply and accurately manufactured using widely available techniques such as 3-D printing. In the future, these arbitrary wavefronts could enable applications in several areas of physical acoustics. Promising uses include the precise correction for aberrations introduced by the skull⁵ and the manipulation of particles.⁶ 

This work was funded by the Engineering and Physical Sciences Research Council and the UCL EPSRC Centre for Doctoral Training in Medical Imaging.