Laser photothermoacoustic heterodyned lock-in depth profilometry in turbid tissue phantoms

Ying Fan, Andreas Mandelis, Gloria Spirou, I. Alex Vitkin, and William M. Whelan Center for Advanced Diffusion-Wave Technologies, Department of Mechanical and Industrial Engineering, 5 King’s College Road, University of Toronto, Toronto, Canada M5S 3G8 Department of Medical Biophysics, University of Toronto and Ontario Cancer Institute/Princess Margaret Hospital, 610 University Ave, Toronto, Canada M5G 2M9 Department of Radiation Oncology, University of Toronto and Ontario Cancer Institute/Princess Margaret Hospital/University Health Network, 610 University Ave, Toronto, Canada M5G 2M9 Ryerson University, Dept. of Mathematics, Physics and Computer Science, 350 Victoria Street, Toronto, ON, Canada M5B 2K3 Received 26 January 2005; revised manuscript received 5 July 2005; published 4 November 2005


I. INTRODUCTION
In recent years, the field of photoacoustic ͑PA͒ applications to soft tissue imaging, cancerous lesion detection, and subdermal depth profilometry has enjoyed very rapid development ͓1-3͔.This is so because PA detection has shown concrete promise of depth profilometric imaging in turbid media at depths significantly larger than accessible by purely optical methodologies ͓4͔.
In state-of-the-art photoacoustic imaging systems, pulsed lasers have always been the choice of signal source ͓4-6͔.Some continuous-wave experimental system configurations, however, have been reported with the reverse effect, acoustooptic imaging ͓7͔.The advantages of using pulsed laser sources include the following.͑a͒ Considerable signal strength, which yields acceptable signal-to-noise ratios ͑SNR͒, can be obtained right after the short pulse.͑b͒ Following optical absorption of a short laser pulse by turbid tissue, optical-to-thermal energy conversion and localized photothermoacoustic ͑PTA͒ volume expansion, the acoustic signals received within approximately 1 s after the end of the laser pulse are essentially thermally adiabatic.Therefore, the signal carries information about the thermal shape of the absorber, which substantially coincides with its geometric shape before any significant heat conduction can deform the image ͓8͔.Pulsed PTA detection, however, presents disadvantages in terms of laser jitter noise and thermal noise within the wide bandwidth of the transducer and hard-tocontrol depth localization of the contrast-generating subsurface features.Furthermore, in order to avoid any detrimental effects to human tissue, the pulse energies must be limited to below 5 mJ/ pulse at the expense of SNR.
Frequency-domain PTA methodologies can offer alternative detection and imaging schemes with concrete advantages over pulsed laser photoacoustics.This paper reports the theoretical and experimental development of a new frequency-domain PTA analytical technique and instrumentation featuring frequency sweep ͑chirp͒ and heterodyne modulation and lock-in detection of a continuous-wave laser source at 1064 nm wavelength.The advantages of the system include the following.͑a͒ Low fluence of harmonicallymodulated cw laser, with the concomitant advantage that a much higher optical energy density can be deposited during a given length of time without damaging the tissue, the thermal effects being further mitigated by thermal diffusion.͑b͒ The superior SNR of the lock-in amplifier compared to pulsed laser averaged transients can offset much of the SNR deterioration at high frequency ͑MHz͒ range.Frequency chirps may also recover the strength of the high-frequency Fourier components through fast-Fourier transformation of the frequency-domain impulse response, thus matching the major advantage of pulsed-laser excitation.͑c͒ Possible multichannel lock-in signal processing and image generation in quasireal time ͓9,10͔.

II. THEORY OF PHOTO-THERMO-ACOUSTIC WAVE GENERATION IN A TURBID LAYER
In this section, a mathematical model is developed to study laser-induced photothermoacoustic waves in turbid media.This model includes both the scattering and the absorption properties of the PTA medium.The contributions of the diffuse photon-density-wave replacing the purely absorbing thermoacoustic signal source are also considered.
Figure 1 shows the geometry used for the onedimensional mathematical model.The configuration closely corresponds to the experimental geometry.It contains three coupled layers: the top and bottom layers are composed of a fluid with the middle layer composed of a solid.The top layer is assumed to be semi-infinite fluid and occupies the spatial region −ϱϽz ഛ −L.It has density f and speed of sound c f .The solid layer has thickness L, density s , speed of sound c s , specific heat at constant pressure C Ps , optical absorption coefficient at the laser wavelength a , optical scattering coefficient s , bulk modulus K s and isobaric volume expansion coefficient ␤ s .The bottom layer extends from 0 ഛ z Ͻϱ.The reason for not considering the finite thickness of the bottom layer is that in our experiments no reflections from the fluid-container interface are observed.
An analytical solution of the coupled PTA problem in the form of spectral integrals can be obtained by converting the time-domain equations to their frequency-domain counterparts using Fourier transformations ͑FT͒ ͓11,12͔.For a harmonic optical source, the Fourier transform of the radiative transfer equation yields Eq. ͑1͒, which is satisfied by the diffuse photon density wave ͑DPDW, or diffuse radiant energy fluence rate͒ field ͓13͔, d ͓W m −2 ͔:

͑1͒
Here a source strength depth distribution is assumed that decreases exponentially into the turbid medium ͑Bouguet's law͒ with total attenuation ͑extinction͒ coefficient: Also, where I 0 is the laser fluence, g is the mean cosine of the scattering function of the photon field over all spatial directions described by the solid angle.In view of the almost entirely forward scattering of photons in tissue, g values range between 0.6 and 0.98 ͓14͔.
The complex diffuse-photon wave number is defined as where Here is the speed of light ͑Ϸ10 10 cm/ s for light propagating in turbid media͒; D is the optical diffusion coefficient, in units of length.The general solution of Eq. ͑1͒ is d ͑z,͒ = A 1 e p ͑z+L͒ + A 2 e − p ͑z+L͒ + Be − t ͑z+L͒ .͑5͒ Constant B can be determined as Constants A 1 and A 2 can be solved using the boundary conditions for the DPDW ͓15͔ where A =2D͑1+r 21 ͒ / ͑1−r 21 ͒ϵ2D.r 21 is the internal reflectance, defined as the ratio of the upward-to-downward hemispherical diffuse optical fluxes at the boundary.Therefore,

͑8͒
To complete the solution, the coherent photon-density field, c = I 0 e − t ͑z+L͒ , must be added to the diffuse-photon-density distribution, d .The total photon density field t = d + c .
In frequency domain, the thermal-wave equation can be written as where s ͑z , ͒ is the thermoelastic temperature rise above ambient.␣ s and s are, respectively, the thermal diffusivity and conductivity of the solid medium.The general solution of Eq. ͑9͒ is + C 4 e − t ͑z+L͒ .͑10͒ Constants C 2 , C 3 , C 4 can be solved as In the fluid z ഛ −L and z ജ 0, the temperature field can be written, respectively, as , C b can be solved using the boundary conditions of thermal continuity at the fluid/solid interfaces: By introducing in the solid a particle or molecule displacement potential, s ͑z , ͒, the coupled wave equations in the solid and fluid can be easily solved.The displacement potential is related to the magnitude of the one-dimensional displacement vector, U s ͑z , ͒, as Due to laser PTA excitation by a large spot-size laser beam, further expanded by intra-solid optical scattering, only longitudinal waves are assumed to propagate in an isotropic solid.This assumption allows the use of the Helmholtz equation which is satisfied by the displacement potential, s : where k s = / c s is the acoustic wavenumber in the solid for small-amplitude acoustic perturbations.The general solution to this equation is + G 5 e − p ͑z+L͒ + G 6 e − t ͑z+L͒ + G 9 e s ͑z+L͒ .͑17͒ Constants G 3 , G 4 , G 5 and G 9 are found to be Inside the fluid, since wave sources are of a potential nature, liquid motion will be potential-driven motion.By introducing a scalar potential of the velocity field,

͑19͒
where the subscript i =1,2 indicates the top and bottom fluid, respectively, one can obtain the photothermoacoustic wave equation ͓Eq.͑20͔͒ for a nonviscous fluid: where k f = / c f is the wave number for small-amplitude acoustic perturbations in the fluid.The small-amplitude pressure change in the fluid is related to the velocity potential, fi by The general solutions to Eq. ͑20͒ can be written as The constants ͑G 1 , G 2 , G 7 , G 8 ͒ in Eqs.͑17͒ and ͑22͒ can be determined through the boundary conditions of stress and velocity continuity at the two interfaces z =0,−L: Substituting the displacement potential, temperature field, and velocity potentials into the boundary conditions, Eq. ͑23͒ can be written as where and G 7 , the constant of interest for the fluid pressure field determination, can be solved as The pressure field can then be calculated as Theoretical simulations were performed for the simple case of a solid turbid layer immersed in water.Three input parameters, the optical absorption coefficient, optical scattering coefficient and thickness of the solid were changed independently for each simulation to illustrate the time-domain PTA signal generation through the developed theory.Table I ͓16,17͔ presents the optical and elastic properties used as input parameters for the mathematical model.Equation ͑27͒ was used to calculate the laser-induced acoustic field within a user-selected frequency range.The time-domain results were obtained from their frequency-domain counterparts using inverse Fourier transformation ͑IFT͒.Figures 2 and 3 illustrate the effects of varying optical absorption ͑penetra-tion͒ depth on the PTA signal from a 5-mm-thick turbid layer.The two peaks in each plot indicate the acoustic waves generated at the top and bottom surfaces of the turbid layer.It can be observed from these figures that increasing the optical absorption coefficient always results in an increase of the signal amplitudes corresponding to the first peaks.The magnitude of the second peak, however, is affected by the amount of energy transmitted to the bottom of the turbid layer and the amount of energy absorbed by that area.For example, Fig. 3, which features an optical absorption coefficient of 4 cm −1 , shows a degraded ratio of the peak magnitudes ͑second peak/first peak͒ compared to that of Fig. 2.This indicates that although the optical absorption is higher in Fig. 3, significant amount of laser energy is absorbed during the light transmission process, resulting in a smaller amount of optical fluence being available to reach the back surface of the solid.
Figure 4 shows the simulated PTA signal obtained from a 10-mm-thick turbid solid with an optical absorption and a scattering coefficient of 100 m −1 .The time delay between the two acoustic peaks is equal to the thickness of the solid divided by the speed of sound in the medium.Due to the increase in thickness, the energy loss during light transmission is more significant in Fig. 4, which results in a smaller amplitude of the second acoustic peak.
Figures 5͑a͒ and 5͑b͒ feature the same material properties, but with an increased optical scattering coefficient, s , from 0.5 cm −1 to 8 cm −1 .This results in a significant decrease of the acoustic signals ͑peaks at around 46 s͒ at the bottom surface of the turbid layer, which is due to the combined effects of energy absorption and scattering during light transmission.The PTA signals at the top surfaces ͑peaks at around 39 s͒, however, show slight increases as the scattering coefficient increases.This phenomenon, which is due to the localization of the optical source closer to the surface, is also evident in the experimental results shown in Sec.V.

IV. EXPERIMENTAL CONFIGURATION: PTA FREQUENCY SWEEPS AND HETERODYNED SUBSURFACE IMAGING
A rapid frequency-scanning ͑chirp generation͒ system was designed and implemented to allow for fast PTA depth profilometry from subsurface inhomogeneities through a wide range of depths at a fixed spatial coordinate.Unlike the stateof-the-art PTA imaging techniques, which use pulsed laser as an excitation source, this frequency-scanning system used a near IR cw laser as an optical source.The advantage of using a low fluence, harmonically-modulated cw laser is to gain a much higher optical energy density deposition tolerance than pulsed laser excitation without damaging the probed tissue for biomedical applications.Frequency chirps may also recover the strength of the high-frequency Fourier components through fast-Fourier transformation of the frequency-domain impulse response, thus matching the major advantage of pulsed-laser excitation.Besides their advantage of speed for full frequency-spectrum acquisition compared to point-bypoint frequency scans, another major feature of frequencyswept photothermoacoustic signals is their ability to be Fourier transformed into a time-delay domain inverse spectrum which carries the depth profilometric information in a series of time sequences equivalent to the impulse response of the acoustic system ͓17͔.Therefore, frequency-swept PTA signals are hybrid between time-and frequency-domain, a unique property which can be used to enhance signal-tonoise ratio and delay-time control significantly through heterodyning and lock-in noise filtering.The block diagram of the frequency-sweep heterodyne PTA imaging system is shown in Fig. 6.The laser used to generate PTA pressure waves is an ytterbium fiber laser ͑IPG Photonics, 1064 nm͒.A frequency-swept ͑chirp͒ signal is generated by a function generator ͑FG1, Stanford Research Systems, DS345͒ to drive the acousto-optic modulator ͑Neos Technologies, N15180-1.06-Gap͒ and modulate the intensity of the laser beam.The chirp signal of FG1 is triggered by a delay-pulse generator ͑Stanford Research Systems, DG 535͒, which is also used to trigger the second function generator ͑FG2͒.The laser beam is directed onto the specimen using an optical mirror.An acoustic transducer ͑Panametrics, V382͒ with a central frequency of 0.5 MHz and a −8 dB frequency bandwidth of 0.1-1 MHz is used to receive the acoustic signal.The received signal is amplified by a preamplifier ͑Panametrics, 5676͒.At the mixer 1 ͑Mini-circuits, ZAD-3͒, the output of FG2 is mixed with the output of the PTA ultrasonic transducer.The output of the mixer is further sent to a low-pass filter ͑LPF1, Stanford Research Systems, SR 640͒.This signal is then mixed with a single harmonic frequency, 0 , using a second mixer.The single-frequency signal is generated by the internal oscillator of a high-frequency lock-in amplifier ͑LIA; Stanford Research Systems SR 844͒.The output signal of mixer 2 is filtered using the second low-pass filter ͑LPF2͒ and then sent to the LIA as the input signal.The amplitude and phase of the LIA output are stored in a computer for display and analysis.
The PTA signal generation flow chart associated with the circuit of Fig. 6 is described in Fig. 7.The linear frequency swept ͑chirp͒ signal generated by FG1 can be written as cos͓͑a + bt͒t͔, where a = 0.1 MHz is the starting frequency and b = 0.9 MHz/ ms is the sweep rate.This chirp signal is triggered by the DG 535 delay-pulse generator, which also triggers FG2.The output of FG2 is also a chirp signal, delayed by a controlled time,, through the delay pulse generator.This signal can be written as cos͕͓a + b͑t − ͔͒t͖.The intensity of the laser beam is modulated by the acousto-optic modulator according to the chirp signal generated by FG1.At the acoustic transducer, the received signal can be written as: cos͕͓a + b͑t − z / c͔͒t͖, where z represents the depth at which the acoustic signal originates, and c is the speed of sound in the probed medium.Due to the linear relationship between the depth and the delay time when the transducer receives the signal, this expression shows that the information at a certain depth can be related to the frequency components of the chirp signal.At mixer 1, the frequency components from the two input channels ͓a + b͑t − z / c͒ and ͓a + b͑t − ͔͔͒ undergo addition and subtraction.The resulting high frequency part ͓2a + b͑2t − z / c − ͔͒ is then removed by low-pass filter 1.The remaining low frequency part ͓b͑ − z / c͔͒ of the signal is down-shifted and contains a wide spectrum of frequencies, each with time-delayed information.Since both the input signals of mixer 1 contain the same starting frequency, a, this component is canceled out after subtraction.The output signal from LPF1 can be written as cos͕͓b͑ − z / c͔͒t͖.To perform PTA depth profilometry, mixing this output ͑second mixer͒, with a single harmonic fixed frequency, 0 , generated by the internal oscillator of the LIA, and low-pass filtering the two sidebands of the mixed signal, yields an output which can be represented by cos͕͓ 0 − b͑ − z / c͔͒t͖.The frequency down-shifted output is detected by the LIA.By scanning the chirp delay time, , a nonzero LIA signal output is expected at 0 only when = z / c.Therefore, scans at a fixed spatial coordinate are equivalent to depth coordinate scans and can yield information from different probe depths in the sample at a fixed lateral coordinate point.Scanning over a predetermined two-dimensional area of the sample will generate a subsurface 3D image.

V. EXPERIMENTAL RESULTS, COMPUTATIONAL FITS TO DATA, AND DISCUSSION
The new PTA imaging system has been characterized using phantoms with tissuelike optical properties.For biomedical imaging, substantial contrast is expected to arise from differences in optical absorptions between healthy and malignant tissue owing to tumor angiogenesis, which gives rise to the presence of increased blood flow in the latter ͓4͔.The optical absorption coefficients at 1064 nm are around 0.1 cm −1 and 10 cm −1 for blood-deficient dermis and oxygenated blood, respectively ͓18͔.The effective scattering coefficient for breast tissue is around 1.2 cm −1 ͓4͔.Solid phantoms were made of plastisol, mixed with different percentages of titanium dioxide and plastic color to closely mimic the scattering and absorption properties of human tissue.Three types of phantom specimens were tested: ͑A͒ ͑7.3± 1͒-mm-thick single-layer solid phantoms with varying optical absorption coefficients, a , ranging from 0.25 cm −1 to 1 cm −1 .͑B͒ ͑3 ± 0.5͒-mm-thick single-layer solid phantoms with a fixed optical absorption coefficient, a =1 cm −1 , and varying optical scattering coefficients, s , ranging from 1 cm −1 to 5 cm −1 .͑C͒ An absorbing phantom ͑ a =3 cm −1 ͒ embedded ͑4 mm deep͒ inside a scattering medium ͑ s = 1.3 cm −1 ͒.The optical properties of each phantom specimen were obtained from the literature ͓19͔.
The solid lines in Figs. 8 and 9 are the experimental results of single point scans on type A specimens, obtained using the PTA imaging system, while the dashed lines are the simulated results.The chirp signal covered a frequency range from 0.1 MHz to 1 MHz and the step size of the delay time used for the scan was 0.1 s.The frequency-domain simulated results were calculated by substituting the thickness, observation distance ͑the distance from the acoustic transducer to the top surface of the specimen͒, material properties, and the chirp frequency range into Eq.͑27͒.The corresponding time-domain pressure fields were obtained by applying inverse Fourier transformation to the frequency-domain results.Good agreement was obtained between the experimental and numerical results.Due to the large size of the laser beam ͑ϳ4 mm͒ and the short distance from the imaging layer to the interface, the PTA behavior was expected to be very similar to the 1D situation and to be adequately interpreted by our 1D theory, as observed.
To obtain the best fits to the entire frequency record of the pressure responses, the exact values of the bulk modulus, K s , isobaric volume expansion coefficient, ␤ s , thermal diffusivities, ␣ f,s , and thermal conductivities, f,s , were found not to be as important as the speed of sound, sample thickness, observation distance ͑the distance from the top surface of the turbid layer to the acoustic receiver͒ and the optical properties of the sample.For the secondary parameters, the listed values in Table II are used for the numerical simulations.An important parameter used for the theoretical fits is the speed of sound, c s = 1390 m / s, which was obtained using time-offlight measurements.The other primary parameters, including the optical coefficients, the sample thickness, and the observation distance are listed in the caption of each figure.
Some important features can be observed from the experimental and numerical results.These features include ͑A͒ A decrease in optical absorption coefficient always results in diminished signal amplitude.͑B͒ Because the optical absorption depths of all the samples are commensurate with their thickness, the second echo observed in the amplitude and phase plots corresponds to direct absorption at the back interface and the launching of an acoustic pulse traveling toward the front surface ͑transducer location͒.The magnitude of the second echo ͑peak͒ increases relative to that of the first peak, as a decreases.This suggests that PTA depth profilometry of low-absorption-coefficient regions may yield full three-dimensional images of sub-surface lesions and result in better contrast than more absorbing features, provided the laser fluence is sufficient to retain a satisfactory signal-tonoise ratio.
The solid lines in Figs. 10 and 11 are the experimental results of single point scans on type B specimens, obtained using the PTA imaging system, while the dashed lines are the simulated results.To obtain the best fits to the entire frequency record of the pressure responses, the most important parameters are the speed of sound, sample thickness, observation distance and the optical properties of the sample.Since both type A and type B materials were manufactured using the same material, plastisol, the speed of sound, c s = 1390 m / s, was used for the numerical curve fitting.The other primary parameters, including the optical coefficients, the sample thickness, and the observation distance are listed in the caption of each figure.The absorption coefficients were obtained using optical measurements and the scattering coefficients were obtained from numerical curve fitting.
Figure 12 shows the depth profilometric image of the cross-section of an absorber ͑4 mm by 3 mm͒ embedded in a scattering medium ͑type C specimen͒.The right-hand-side vertical axis indicates the time delay and the left-hand side vertical axis indicates the equivalent depth.The horizontal axis is the spatial scan coordinate.The range of the delay time is precisely controlled to be 13-23 s to cover an area of interest ͑object area͒.This level of image depth control is a major advantage of this PTA technique compared to conventional pulse-laser diagnostics.The step size of the delay time is 0.1 s, which corresponds to a distance of around 150 m in water.The horizontal pixel size is 0.5 mm.The front surface of the absorber is clearly visible with a sharp increase of the PTA signal.The bottom surface of the embedded object features lower signal amplitude, which is due to the attenuation of optical and acoustical energy.The depth resolution, ⌬s, of the system is around 0.75 mm.This value   is calculated, using ⌬s = ⌬ ϫ c, where c Ϸ 1500 m / s is the sound velocity inside the probed medium and ⌬ Ϸ 0.5 s is the uncertainty of the delay time, which is caused by system jitter of the function generator.The horizontal spatial resolution is observed to be approximately 1 mm.The vertical and horizontal resolutions are comparable to those obtained from a pulsed laser photoacoustic system ͓4͔, which uses a 32element transducer array as a sound receiver.One of our future tasks is to develop and adopt a transducer array for the cw PTA imaging system to improve system resolution.

VI. CONCLUSION
A linear frequency-domain PTA theory has been developed for a composite liquid-solid-liquid one-dimensional geometry, which includes both the optical scattering and absorption effects and natural mixed ͑rigid and free͒ boundary conditions at the solid-liquid interface.The resulting predic-tions were compared to experimental data obtained using a newly developed PTA imaging system.The PTA imaging system features linear frequency sweep, heterodyne modulation and lock-in detection of a continuous-wave laser source at 1064 nm wavelength.This system can be used to perform subsurface slice-by-slice imaging from operator-determined, precisely-controlled depths.Good agreement was obtained between the numerical and experimental results on solid tissue phantoms with different optical absorption and scattering coefficients, indicating that the PTA depth profilometry system using a continuous-wave laser source and frequencyswept, heterodyned detection, may be applicable to biomedical imaging of blood rich tissue as in the case of subsurface cancerous tumors.

FIG. 1 .
FIG. 1. Geometry used for formulating the PTA problem.Symbols are defined in the text.

TABLE I .
Elastic properties used as input parameters for the numerical simulation.

TABLE II .
Elastic properties used as input parameters for the mathematical model.