Numerical Simulation of Enhanced Photoacoustic Generation and Wavefront Shaping by a Distributed Laser Array

: In photoacoustic imaging, the use of arrayed laser sources brings several advantages. Acoustic waves can be generated with flexible control of wavefronts, bringing functionality such as ultrasonic beam steering and focusing. The use of arrays reduces the optical intensity while increas ‐ ing the strength of the ultrasonic wave, bringing the advantages of improved signal ‐ to ‐ noise ratio (SNR) while avoiding laser ‐ induced damage. In this paper, we report a numerical model for study ‐ ing the generation and shaping of acoustic wavefronts with laser arrays. The propagation of me ‐ chanical waves, photoacoustically generated by thermal expansion, is simulated and discussed in detail. In addition, a partially delayed distributed array is studied both theoretically and quantita ‐ tively. The developed model for wavefront control through time ‐ delayed laser pulses is shown to be highly suited for the optimization of laser array generation schemes.


Introduction
Photoacoustic imaging as a non-invasive detection technology has been widely used in biomedical imaging [1], nondestructive testing [2], and other fields.Due to the thermoelastic effect [3], when a laser beam irradiates a sample, optical absorption in the material causes the light energy to initiate thermal expansion, resulting in deformation of the illuminated area and its surroundings.If the illumination is a fast transient pulse or if the light is modulated periodically, then the material and its surroundings deform correspondingly, thereby generating an acoustic wave.
In photoacoustic imaging, structural features or defects in the object can be characterized by detecting reflections or other perturbations in the ultrasonic waves.The spatiotemporal features of the ultrasonic wave are related to those of the incident laser light, but although an intense beam can generate strong acoustic waves, it may damage the material.However, laser intensities far below the damage threshold generate acoustic waves that are relatively weak.The signal-to-noise ratio (SNR) of the reconstructed image can be quantitatively expressed as SNR ~ S × γ × A × P × M /N (1) where S is the photo-generated ultrasound intensity, γ is the sensitivity of the detector, A is the imaging area, P is the number of projections, M is the number of signal averages, and N is the noise floor of the system [4].Thus, if the excited ultrasound signal intensity, which is proportional to the optical energy, is reduced by a factor of K, to recover the SNR, averaging over K 2 times is needed.Therefore, in order to obtain a higher SNR, researchers tend to increase the magnitude of the generated photoacoustic signal.In current biomedical photoacoustic imaging research, to enhance the imaging quality or increase the imaging depth, the laser pulse energy density used in the experiment is usually higher than the ANSI specified threshold of 20 mJ/cm 2 [5], which can result in ablation of the tissue [6,7].It is vital to explore methods for generating strong acoustic waves using energy densities below the laser ablation threshold.To achieve efficient photoacoustic wave generation, Murray et al. used a continuous wave laser with intensity modulation to excite ultrasound waves [8], which greatly improved the SNR and simultaneously kept the sample surface intact.However, this method resulted in the generation of multiple acoustic modes and reflected waves, and interpretation of the signals was complicated.Huang et al. proposed a photoacoustic method that uses a laser array instead of a single light source to excite acoustic waves [9], which increases the intensity of the ultrasonic wave without increasing the energy density illuminating the sample surface.
Research of Matt Clark's group demonstrated that waves can be focused and steered using a spatially distributed laser source, with a spatial light modulator or specially designed computer-generated holograms [10][11][12][13].They also developed a method using laser-induced phased arrays (LIPAs) to detect defects in metals and provide in-process monitoring during additive manufacturing (AM) [13,14].It has also been shown that temporal control of the elements in a linearly arrayed laser source is an efficient way to enhance the intensity of the generated acoustic wave or steer the propagation direction [15][16][17][18][19][20][21].Although wave steering and signal enhancement by controlling the spatial or temporal distribution of the generated laser source have been demonstrated, the method for manipulating the wavefront of a laser-generated acoustic wave is still uninvestigated.
The use of acoustic metamaterials to manipulate the wavefront of propagating acoustic waves from various sources has been proposed and developed.However, it is difficult to apply these to laser-generated acoustic waves due to the operation bandwidth and the requirement of the wave to be transmitted through or reflected from the meta-surface.
Beamforming of acoustic waves through array sensors has been widely used to improve the contrast, spatial resolution, and signal-to-noise ratio of ultrasound images.Based on our understanding of previous research, two categories of techniques are conceptualized as the beamformer approaches, either at the excitation-side and the probingside of acoustic waves.
The traditional probe-side beamforming method achieves the purpose of signal enhancement by processing the received signal from the perspective of signal acquisition rather than directly generating stronger acoustic signals from the perspective of excitation.This is a fatal flaw in extremely absorbent biological tissues.
Traditional probe-side beamformer methods include the delay and sum (DAS) method [22], the minimum variance (MV) method [23,24], etc.The principle of these methods is to align the target beams according to the delay time and add them sequentially to enhance the signal.DAS is independent of echo data, which leads to the main lobe width being too wide and the side lobes intensity being too high [24].The minimum variance (MV) method, e.g., adaptive beam synthesis, calculates a dynamic weight value based on the echo data, makes full use of the characteristics of the echo signal, reduces the side lobe signal, and improves the resolution of the image.The basic principle is to minimize the array output by keeping the gain in the desired direction constant.
For traditional excitation-side forming, Von Ramm et al. achieved beam steering and enhancement by applying a time delay to each piezoelectric plate in an 8-element line array [25].The distance between the line sources is about one wavelength.The delay time is controlled by the line source spacing and rotation angle.
Shi-Chang Wooh [26] proposed that a larger inter-element spacing can help improve the directivity of the steered wave, yet only to a certain extent; the grating lobes will be introduced (or the beam is not steerable) when the inter-element spacing exceeds its critical value, which is about half the wavelength of the sound wave and controlled by the desired maximum steering angle.
The traditional excitation-side beamformer can achieve controllable superposition between different sources, thereby achieving high-degree-of-freedom beam direction control.In principle, it can also achieve the wavefront formation and control discussed in this article.However, a systematic study of the process of exciting, shaping, and steering photoacoustic waves using laser array sources is therefore still essential.The purpose of this article is to explore whether it is feasible to realize the formation and control of the wavefront in the process of photoacoustic excitation, as well as the mechanism of its formation and propagation.
In this paper, the characteristics of surface and bulk waves excited by a laser array with partial time delay are numerically simulated using the finite element method (FEM).We verify whether, by applying a time delay between different elements of the laser array, the photoacoustic wave can be enhanced and whether this enhanced wave still propagates as a wavefront maintaining its complete shape.

Numerical Simulation
Both the finite element method (FEM) and finite-difference time-domain (FDTD) method can be used to model the process of light-excited acoustic waves.In terms of computational speed, FDTD has a great advantage.However, considering the complexity of photoacoustic imaging, the FEM method was selected as it is more versatile and rigorous.
Three configurations for exciting photoacoustic waves were simulated, namely single-pulse single-point excitation (Model 1), single-pulse laser array excitation (Model 2), and laser array with time-delayed excitation between spots (Model 3).Except for the number of excitation light sources and the time delay between the spots, other parameters in the simulation remained unchanged, including the sample size, material, energy density, and pulse duration.
The modeled sample was a homogeneous, isotropic cylindrical aluminum block with a height of 4 mm and a radius of 4 mm.The laser array was incident perpendicularly to the top surface of the cylinder.Material parameters of aluminum from reference [27] were used for the simulation, with a surface wave sound velocity of 2940 m/s, a transverse wave sound velocity of 3080 m/s, and a longitudinal wave sound velocity of 6320 m/s.A laser source was used to generate the ultrasound, with a pulse width of 20 ns, pulse energy of 160 μJ, and energy density of 510 mJ/cm 2 .The laser spot radius on the sample surface was 100 μm, and it was assumed that 80% of the light was absorbed.Figure 1 shows the aluminum cylinder in the Cartesian coordinate system, with three laser spots A, B, and C focused on the top surface.The simulation used the solid heat transfer module and the solid mechanics module of COMSOL Multiphysics to simulate in the transient field Model 1 and Model 2 were used as reference simulations to compare with Model 3. In Model 1, a single excitation pulse was used to illuminate the single spot "B".Model 2 was based on Model 1, with excitation spots "A" and "C" added on either side of the single laser spot B. The size and energy density of these three laser spots were the same as for Model 1.As shown in Figure 1, through pre-simulation results, the distance between adjacent laser sources was selected as 0.9 mm, which is about one wavelength of the excited acoustic wave.Model 3 was based on Model 2 with a time delay imposed on the central spot B. It can be seen from the geometric properties of the propagation of acoustic waves in Figure 2a that as the waves propagate, the difference between d1 and d2 becomes smaller.When the distance between adjacent sources is much smaller than d1, and the wavefront has traveled sufficiently far, δd = d2 − d1 < ε × v. ε is defined as the allowable time error, which is chosen to be as same as the simulation time step, i.e., 0.01 μs, v is the velocity of the excited acoustic waves.During this time, the acoustic wave travels about one-tenth of the acoustic wavelength.The new merged wavefront due to superposition can be considered as a continuous wavefront that can propagate in an infinite medium.The distance Y at which the superposed acoustic waves add constructively on the yaxis can be calculated from where L is the distance between adjacent sources, and ε × v is the maximum acceptable error.In Model 2, the distance Y for acoustic waves without a time delay was calculated to be 13.5 mm, which is larger than the size of the sample.The corresponding traveling time is about 4.5 μs, calculated from the velocity of the surface acoustic wave.Furthermore, using this arrangement in practical applications would make it difficult to image features such as defects or impurities within the sample.Figure 2b shows the propagation of acoustic waves when excitation of the central source B is delayed.If the delay in exciting the central source B is δt, the delayed acoustic wave from B adds constructively to the acoustic waves from the symmetrical sources A and C at a time t, which can be calculated from when the delay time is chosen to be 0.1 μs, constructive wave superposition occurs 0.5 μs after the first excitation at a distance on the y-axis of about 1.2 mm, which is less than 10% of the distance without any time delay.When the delay time is increased to 0.2 μs, the time for constructive superposition is shortened to 0.325 μs with a propagation distance of only 0.375 mm.
To ensure that the superposed propagating waves do not separate again, δd is required to be always less than the acceptable error value.
where t is the propagation time after constructive interference, t0 is the time that constructive interference first occurs, and δt is the delay time.For all times t > 0, δd < ε × v, that is, there is no solution when δd ≥ ε × v.We can see that when the delay time is 0.1 or 0.2 μs, the solution of the inequality is an empty set.Therefore, after constructive superposition occurs, the enhanced wavefront propagates as a complete new wavefront without separating again.In the results and discussions section, we will confirm that the enhanced wavefront propagates as a complete wavefront by fitting and reconstructing the wavefront.

Solid Heat Transfer
The energy distribution of the laser can be described as where I0 is the maximum power density, and f(r) and g(t) are the spatial and temporal distributions of the laser pulse, respectively.The thermal conduction equation can be written as where T(r, z, t) represents the temperature distribution at time t; ρ is the material density; Cp is the thermal capacity of the material; and k is the thermal conductivity.
We consider the excited spot as a heat source loaded at the center of the cylinder surface.To simplify the simulation and reduce the computation time, we consider other boundaries of the cylinder as thermally insulated boundaries.The boundary conditions of the other surfaces can therefore be written as And the boundary conditions of the absorption layer can be written as where A(T) is the optical absorptivity of the focus plane [28].

Solid Mechanics
When the laser pulse is focused on the surface of the sample, a transient displacement field is excited due to thermoelastic expansion.In a homogeneous medium, the displacement satisfies: where u is the transient displacement, ρ is the density of the material, α is the thermoelastic expansion coefficient, and λ and μ are the first-order and second-order Lamé constants, which are determined by the material-related quantities represented by the strain-stress relationship.The first-order Lamé constant λ represents the compressibility of the material, which is equivalent to the bulk elastic modulus or Young's modulus, and the secondorder Lamé constant μ represents the shear modulus of the material.The values of λ and μ for aluminum are 6.10 and 2.49, respectively [29].

Results and Discussions
Figure 3 shows the z-direction displacement due to the surface wave propagating in the xy-plane.Comparing Model 1, Model 2, and Model 3, we can see that a wavefront showing constructive superposition appears first in Model 3. The color bar used in each subfigure in Figure 3 is the same, ranging from −1 × 10 −6 (dark blue) to 0 (light blue) to 3 × 10 −6 mm (red).It can be seen that due to the time delay, the wavefront in Model 3 lags behind Model 1 and Model 2, but its intensity is significantly enhanced.In Figure 4, we plot the displacement for the three models at points at which the distance to the x-axis and y-axis are both 1, 1.5, 2, and 2.5 mm, referred to as point 1, point 2, point 3, and point 4, respectively.It can be seen that the peak displacement for Model 3 is much larger than that for Model 1 and Model 2 because of the constructive superposition of the acoustic wavefronts.The simulation results agree closely with the theoretical model, but the noise of Model 3 is larger than that of Model 1 and Model 2, and this problem needs to be considered in the subsequent data processing.MATLAB was used to fit the wavefront data to confirm that the enhanced wavefront continues to propagate in the form of a constructively superposed wave.A surface composed of points of equal amplitude is considered a wavefront surface.For the purpose of choosing a fitting method with a simple form and high accuracy, we used second to fifthorder polynomial fitting with the Zernike polynomial.Rayleigh waves propagate in the form of cylindrical wavefronts along the free surface of the medium [30], so only the first three orders were considered when using Zernike polynomials for fitting.The first three orders of the Zernike polynomial are shown in Table 1.Due to symmetry, only the region with x > 0 was considered when fitting.Due to acoustic absorption, the surface wave decays exponentially inside the material.The surface wave on the z-axis far from the heat source is very weak and can barely be discriminated from noise, so only data within 0.1 mm of the heat source was used for fitting.
Figure 5 depicts the wavefronts of Model 1, Model 2, and Model 3 fitted by polynomials of different orders.The fitting error is displayed in Table 2.For higher-order polynomials, the fitting error is reduced.The results using third-order Zernike polynomial fitting are similar to third-order polynomial fitting.In both Model 1 and Model 3, fitting using fourth-order and fifth-order polynomials offers significantly better accuracy compared with second-order and third-order fitting, and it is recommended to adopt the fourthorder polynomial fitting as the difference between fourth-and fifth-order polynomial fitting is marginal.Regardless of the order of fitting for Model 2, the fitting error is larger than that for Model 1 and Model 3, and the fitted results are also not ideal.It can be seen from the original data of Model 2 in Figure 5 that weaker wavefronts appear distinct from the main wavefront.This implies that the use of multiple sources without time delay for simultaneous excitation not only cannot form a single enhanced wavefront, but the weaker wavefronts also lead to difficulties in interpreting ultrasonic images and extracting data.For further verification, we conducted numerical calculations based on the parametric indirect microscopic imaging (PIMI) method [31], which is a method that can optically visualize the sound field wavefront by measuring the phase retardation and polarization ellipse orientation angle of the probe light.Linearly polarized light with rotation of its polarization angle was used as the probe beam incident on the surface of the sample.The polarization states of the reflected light were perturbed by the surface acoustic wave.By measuring the reflected light at four different polarization angles, I(0°), I(90°), I(45°), and I(135°), the Stokes parameters, the phase delay δ, and the polarization orientation angle ψ were calculated.The calculated polarization parametric images of the probe beam can be related to the ultrasonic wave field via the photo-elastic effect [31].
The visualization of the acoustic field is shown in Figure 6.The Stokes parameters S0, S1, S2, S3 phase delay δ and polarization orientation angle ψ can be obtained from I(0°), I(90°), I(45°), I(135°).I(0°) is the intensity of the linearly horizontal polarized light component, I(90°) is the intensity of the linearly vertical polarized light component, I(45°) is the intensity of the linearly +45° polarized (L+45P) light component (L for linearly and P for polarization), I(135°) is the intensity of the linear −45° polarized (L−45P) light components.From Figure 6, we found that S3 and ψ show high sensitivity to the laser-induced surface acoustic wave, while the wavefront can barely be seen in the phase retardation δ.As a comparison, the propagation of acoustic waves in the xy-plane induced by three acoustic emitters at the same locations as the laser array source was simulated using MATLAB's k-wave toolbox [32].The heat sources in Model 1, Model 2, and Model 3 were replaced by sound sources with the same size, shape, distributions, and delay time.The simulation area was a square with a side length of 8 mm, and an 800 × 800 mesh was used.The speed of sound in the medium was taken to be 2940 m/s, the same as the surface wave sound velocity of aluminum.The uniformly distributed source was defined by the acoustic pressure, which was assigned a value of −0.1 Pa because thermally induced SAWs have a negative initial acoustic pressure.The initial frequency of the acoustic signal was specified as 3.    Figure 10 shows normalized FEM simulation results and k-wave simulation results separately displayed on the same graph.The simulation results of the two methods have the same trend, but there are slight differences in details, k-wave simulations are simply based on a superposition of sound waves, while FEM simulates the complete process of sound waves generated by thermal expansion.There are also some differences between the different methods.In Model 1, Model 2, and Model 3, the arrival time of sound waves at detection points in FEM is earlier than that in k-wave.This is because the starting point of sound wave propagation in FEM is at the boundary of the heat source (strain caused by temperature gradient), while the starting point of sound wave propagation in the kwave simulation is at the center of the sound source.The propagation distance of photogenerated acoustic waves in FEM is reduced by half the heat source radius compared to k-wave, which leads to the sound waves in the FEM simulation at the detection point arriving earlier.In addition to the main wave, there are waves that are generated by multi-point excitation in Models 2 and 3, but these fail to fully participate in the superposition.The kwave approach was only used as a comparative simulation to verify the suitability of the model.The actual photoacoustic process is closer to FEM simulation.Therefore, only the signal differentiation problem in the FEM model is considered here.Figure 10 shows that the excited sound wave has only a single peak in the negative direction, which can be used as the main peak of detection.
In the FEM simulation, there is slight vibration after the sound wave passes.This is because of the coarsening of edge mesh and the inhomogeneity of mesh division due to the compromise of calculation time.This kind of slight vibration can be eliminated by refining the mesh and decreasing the time step.
Simulations of the bulk wave are shown in Figure 11.Clearly, the laser array enhances the intensity of the bulk wave.As expected, the amplitude of the acoustic wave generated by the laser array is larger than that generated by a single laser source with the same energy density.In Figure 11a-c, the difference between the arrival times of the peak of the bulk wave in Model 1 and Model 2 at Points 5, 6, and 7 (The points on the z-axis that are 1, 2, and 3 mm away from the central heat source B) are 0.22, 0.10, and 0.06 μs, respectively, after the excitation time.This is because, for detection points further from the excitation source, the propagation paths of bulk waves originating from spots A and C (symmetrically located on either side of spot B) gradually approach the propagation path of the wave originating from spot B. Because there are two sources, the dominant bulk wave propagating along the z-axis results from a superposition of the waves excited at spots A and C. In Figure 11a, the smaller peak is the acoustic wave excited by source B. In Figure 11b,c, this peak adds constructively to the waves excited by A and C, which increases the width of the wave.Comparing the peak amplitudes of waves at locations 1, 2, and 3 mm from spot B on the z-axis, it can be seen that the amplitudes of the bulk waves generated in Model 2 are, respectively, 1.1, 1.9, and 2.4 times of those excited by Model 1.In the far-field, the bulk wave excited by a distributed laser array is enhanced by a factor of more than 2 compared to a single spot.This phenomenon is analogous to the mechanism of the widely used synthetic aperture method for the enhancement of the detection SNR [33]. Figure 12 shows the attenuation curve, obtaining the amplitude of the z-direction displacement at six points along the z-axis: 1, 1.5, 2, 2.5, 3, and 3.5 mm away from source B. It can be seen that the slope of the attenuation curve for Model 1 is steeper than that for Model 2. This phenomenon arises because the wave has a plane component.

Conclusions
This study compares surface and bulk acoustic waves excited by a single laser source, a laser array composed of three sources without time delay, and a laser array of three sources with a time delay imposed on the central source.It is shown that for the same excitation energy density, arrayed photoacoustic excitation leads to surface acoustic waves that are around four times stronger in particular directions than that of a single source.Moreover, by applying an appropriate time delay to the center source, the time and distance required for constructive superposition to form a complete wavefront can be effectively shortened.The use of the laser array has a significant enhancement on the amplitude of bulk acoustic waves propagating into the object, analogous to the synthetic aperture method in ultrasonic detection.In a future study, we will systematically investigate the characteristics of wave generation by laser array sources with different energy spatial distributions, temporal profiles, and delays between pulses in order to obtain a better understanding of the mechanisms of wave generation using laser arrays.Using a laser array with an appropriate delay time can enhance the intensity of the excited acoustic signal without damaging the sample, and the excited acoustic wavefront still propagates in a regular shape with little clutter.This method provides a new way to manipulate the wavefront and improve the generation efficiency in nondestructive photoacoustic imaging applications.

Figure 1 .
Figure 1.Cylindrical aluminum block and the three laser illumination spots A, B, and C on its top surface.

Figure 2 .
Figure 2. (a) Schematic diagram of laser array excitation with no time delay between pulses.(b) Schematic diagram of laser array excitation with a central pulse (B) delayed.

Figure 3 .
Figure 3. Displacement in the z-direction due to the surface wave propagating in the xy-plane at different times.(a) Model 1; (b) Model 2; (c) Model 3. The corresponding times are 0.6, 0.8, 1, 1.2, and 1.4 μs (from top to bottom).

Figure 4 .
Figure 4. Displacement in the z-direction as a function of time of Model 1, Model 2, and Model 3 for (a) point 1, (b) point 2, (c) point 3, and (d) point 4.

Figure 5 .
Figure 5. Processed data fitted with 2-5 order polynomials and Zernike polynomials at 1.2 μs.From top to bottom: original data interpolation results, second, third, fourth, and fifth-order polynomial

3
MHz, approximately equal to the frequency of the sound wave excited in the FEM simulation.The input signals of the acoustic sources A, B, and C in Model 3 are shown in Figure 7.The input signal for source B has a time delay of 200 ns relative to sources A and C.

Figure 7 .
Figure 7.The input signals of the acoustic sources A, B, and C. The wavefronts for Model 1, Model 2, and Model 3 at different times, shown in Figure 8, are generally consistent with the sound field distributions simulated by FEM.Four probes were placed at Points 1, 2, 3, and 4 (at 1, 1.5, 2, and 2.5 mm to both x-axis and y-

Figure 11 .
Figure 11.Displacement in the z-direction as a function of time for Model 1 (single laser source) and Model 2 (laser array) at (a) Point 5, (b) Point 6, and (c) Point 7.

Figure 12 .
Figure 12.The peak z-direction displacement of the wave for different locations for the cases of single laser and laser array sources.

Table 1 .
First 3 orders of the Zernike polynomial in Cartesian coordinate.

Table 2 .
Fitting error: SSE (sum of squares for error), R-square, and RMSE (root mean squared error).