Spectral power density of the random excitation for the photoacoustic wave equation

The photoacoustic phenomena was described by Alexander Graham Bell in 1881 but its application as an imaging method has attracted much attention after the improvement of the laser, computer and transducer technologies.^1–17 Photoacoustic imaging (PAI) uses absorption of the electromagnetic energy as a contrast mechanism to analyze the system.¹³ In PAI, a short pulsed laser is sent to tissue to obtain an acoustic impulse response. The tissue absorbs the laser energy and then the absorbed light leads to thermal expansion generating a pressure wave. The propagating photoacoustic wave is received by a detector called as a transducer. Due to the transducer and the signal processing, the absorption profile of the tissue composition is obtained.^1,2 Applications of PAI range from tissue and molecular imaging to clinical medicine and biology.¹⁸ 

It is well known that in acoustics, the cross-correlation of observed waves due to the excited random sources can yield the extraction of the Green's function.^19,20 For example, Snieder obtained the Green's function from the diffusive and the acoustic waves due to sources described by the Dirac delta functions.^20,21 This method is also utilized when the sources are placed inside a domain or on its surface.²⁰ 

Obtaining the Green's function by the correlation of the waves resulting from the randomly distributed sources has attracted much intention. The representation theorems of the correlation and convolution types, and time reversal invariance provide the utilization of the method. On the other hand, the formalism does not require time reversal invariance necessarily. For example, Snieder obtained the Green's function for the diffusive medium.²¹ There are a lot of applications of this technique in the fields of ultrasound,^21–28 exploration seismology,^21,29,30 structural engineering,^31,32 numerical modeling,³³ and crustal seismology.^{21,23–25,34}

In this work, the superposition of the Green's function and its complex conjugate (time reversal counterpart) is retrieved from the correlation of photoacoustic waves excited by random point sources inside the volume. The formalism presented here makes it possible to obtain the photoacoustic pressure waves at the positions of the random sources and the power density spectrum of the photoacoustic system in the frequency domain.

The value of this study roots from the fact that the cross-correlation measurements can be used to retrieve functional information related to the flow dynamics inside the tissue. For example, J. Brunker and P. Beard³⁵ quantify Doppler time shifts using the cross-correlation of photoacoustic signals. Then, using the measured time shifts, they estimate the velocity of the red blood cells. The motivation behind the use of the cross-correlations is that the motion of the red blood cells cause fluctuations in the photoacoustic signal. It is also very important to note that red blood cells inside the blood can be modeled as randomly distributed (spatially uncorrelated) point sources.³⁶ 

The photoacoustic wave equation is given by

Here p(r, t) is the solution of the equation representing a photoacoustic wave. β, κ, v_S and T(r, t) stand for the thermal coefficient of volume expansion, the isothermal compressibility the speed of sound, and the rise in temperature at position r and time t, respectively.^1,37

If the acoustic confinement and thermal confinement conditions are satisfied, then the photoacoustic wave equation takes the following form

Defining

Eq. (2) becomes

The convolution type of the representation theorem^20,38,39 can be obtained from the integral

Here p_A and p_B represent the photoacoustic waves resulting from source terms

Using the Fourier transformation, f(r, t) = ∫f(r, w)e^−iwtdω, for p and s gives

in the frequency domain, ω.

The adjoint of Eq. (6) is

where ^* represents complex conjugation.

Substituting the Fourier transformation of the photoacoustic wave equation in Eq. (5) leads to obtaining an expression for the representation theorem of the convolution type

Evaluation of the following integral provides an expression for the representation theorem of the correlation type

where E_i and Q_i are the Fourier transform and the adjoint of the Fourier transform of the photoacoustic wave equation, respectively. Substituting Eqs. (6) and (7) in the correlation type^20,38,39 integral leads to

For a laser pulse, the temporal and space parts of the source can be considered the Dirac delta function.¹ Calasso et al. studied the photoacoustic effect resulted from heat deposition at a point in an inviscid fluid describing the source term in the wave equation by the Dirac Delta function to determine the fluid motion.⁴⁰ For this reason, setting

and writing s_i into

gives

Inserting these s_i and p_i states into Eq. (8) leads to

Using the product rule for differentiation

and applying Gauss's theorem, Eq. (14) takes the following form

Following similar steps, the correlation type integral is obtained

Assuming that S is a sphere with a sufficiently large radius,^41,42 the gradient of G is approximated by

where the Green's function for the Helmholtz equation is

and n is the unit vector pointing in the positive r direction. The surface integral on the right hand side of Eq. (16) is zero²⁰ so that the reciprocity provides

Writing Eq. (18) Eq. (17) and using the reciprocity, the representation theorem of the correlation type leads to

In order to express the right hand side of Eq. (21) in three dimensions, we integrate Eq. (21) over the radius of the sphere r, which is now treated as a variable, between 0 and R. Here, R is greater than r_A and r_B (R > r_{A, B} > 0). In this way, the integration of the surface integral over r leads to a volume integral. The integral on the left hand side of Eq. (21) becomes

and the integral on the right hand side of Eq. (21) becomes

Note that the radius of the sphere, r, varies from 0 to R. In fact, the surface integral in Eq. (21) is obtained by the far field approximation. On the other hand, in the near field, r ∼ 0, the contributions of the surface integrals in Eq. (17) are very small because of the fact that the terms ∇G(r, r_A) and ∇G^*(r, r_B) fluctuate randomly and the total contribution from these integrals can be neglected. For this reason, the volume integral in Eq. (23) contributes significantly in the far field and it can be omitted in the near field. Therefore, we assume that the radius of the sphere ranges from zero to R.

Hence, in three dimensions, Eq. (21) can be generalized to the following expression

Consider spatially uncorrelated random point sources at different location inside the sphere of radius R with the following relation²¹ 

in order to obtain a relation between the casual and the acasual Green's functions from the correlation of the solutions of the photoacoustic wave equation generated by the random point sources where ∣s(ω)∣² and brackets represent the power spectrum of the point sources and the expectation value, respectively.

Multiplying Eq. (24) with the power spectrum ∣s(ω)∣² gives

and then

Here p(r_A, ω) and p(r_B, ω) are the pressure waves at r_A and r_B resulting from point source B at r_B and point source A at r_A, respectively. Eq. (29) states that the multiplication of the superposition of the Green's function and its adjoint with the power spectrum of the excitation can be written in terms of the correlation of the photoacoustic waves at r_A and r_B, respectively.

The Green's functions for Eqs. (6) and (7) are given by the following expressions

and

respectively. Writing expressions (30) and (31) into Eq. (29), the power density spectrum takes its final form

Note that r_A and r_B are the point source locations. It seems that it is impossible to measure p(r_A, ω) and p^*(r_B, ω) corresponding to the photoacoustic pressure waves at the locations of the point sources A and B, respectively, via detectors. In order to obtain these pressure waves, a method is presented in the next section. Eq. (32) states that the spectral power density of the system can be obtained correlating the calculated photoacoustic pressure waves at different locations due to the random excitation. Therefore, the impulse response of photoacoustic systems can be obtained utilizing the method presented in the next section and the representation theorems of the convolution type and the correlation type.

M. A. Anastasio et al. studied the photoacoustic tomography (PAT) reconstruction problem in the frequency domain with inverse source concepts.⁴³ They also obtained a mathematical expression between the photoacoustic data and the optical absorption distribution in the frequency domain. In this work, the inverse problem dealt with is obtaining the pressure waves at the sources' locations from the measurements of the photoacoustic pressure waves on a spherical surface surrounding sources different from Anastasio's work.⁴³ In this work, point sources are considered volume sources not at the bounding surface. Reconstruction of the absorption map of the tissue is the inverse problem in PAT.⁴³ 

A pressure wave, p(r), can be measured via a transducer at a position of r. For two point sources at r₁ and r₂ lying in the same direction in which polar coordinate θ = 0 (or x = 0 and y = 0 line), the Fourier transformation of the photoacoustic wave equation in frequency domain leads to

When r ≠ r_i, Eq. (33) reduces to Helmholtz equation. For azimuthal symmetry, the general solution of Helmholtz equation in spherical coordinates is

where μ = cos θ, P_l(μ) is the Legendre function and

At the position of the trancducer r = R₁, we take the general solution of Eq. (33) with r ∈ [r_i, R₁]. Hence, the experimental p(R₁, t) in ω,

In order to find a_l and b_l, another pressure wave needs to be measured at a different location, r = R₂, so that

Multiplying Eqs. (35) and (36) with

and

where the orthogonality relation of the Legendre polynomials

Solving Eqs. (37) and (38) gives

and

where

and the detector is assumed to be capable of measuring all possible azimuthal angles at a fixed radius.

The solution between the point sources, r ∈ [r₁, r₂], is

The continuity of the pressure wave at the locations of Dirac delta point sources, r = r_i, requires

For azimuthal symmetry, the Laplacian takes the following form

where the Laplace-Beltrami operator,^44–46 Ω_LB, is given by

for three dimensions.

Writing the Laplacian and

where

Setting

ε → 0⁺ and making use of the continuity of the pressure wave, Eq. (47) becomes

Here, we define

Writing

and

where prime represent differentiation with respect to r.

Solving Eqs. (51) and (52) for c_l and d_l in terms of a_l and b_l gives

and

where

Finally, the pressure wave at the position of the source, r = r₂, is

The pressure wave at the location of the other point source (r = r₁),

can also be determined.

It is also possible to find the solution of the innermost region by using the continuity and the discontinuity of the photoacoustic wave at r = r₁. For this region, since

Performing the similar calculations, e_l is obtained

This methodology can be generalized to n finite number of point sources lying in the same θ direction (See Fig. 1). In other words, all the point sources are assumed to be co-linear.

The pressure waves are

and

for the (1)^th (leftmost) interval and the (n + 1)^th (rightmost) interval, respectively.

For the intermediate regions, the pressure waves in terms of linear combinations of

Constants a_il and b_il can be obtained by applying the boundary conditions at the boundaries of the intervals. By using the continuity of the pressure wave and the jump of its derivative, the following equations at the boundaries of the intervals (at r_i's) are obtained

and

where

where primes stand for differentiation with respect to r and

Multiplying Eq. (66) by the matrix

where the Wronskian,

By using these transfer matrices for a finitely many number of the Dirac delta functions, the coefficients of the rightmost region can be related to the coefficients of the leftmost region as

where we set b₂ = 0 since

Hence, when there exists n finitely many number of random sources, the pressure wave in the leftmost region that cannot be measured via a detector can be determined by calculating total transfer matrix,

Eq. (70) shows that when there exists finitely many number of point sources, the pressure wave at the innermost region can be calculated by using only the detector readings at the outermost region regardless of the pressure waves between the innermost and the outermost regions.

Therefore, pressure waves for diffusive media can be measured at the location of sources which are at the bounding surface. On the other hand, the pressure waves at the locations of the point sources not on the boundary but throughout a volume bounded by a surface cannot be measured directly for photoacoustic media. For this reason, by detecting pressure waves at any positions that are greater than the positions of the point sources and then using the continuity and the discontinuity of the pressure wave, the pressure waves p_A and p_B at the source positions r_A and r_B, respectively, emitted from the sources can be determined. Writing the calculated p_A and p_B waves in Eq. (32), the spectral power density of the random sources can be obtained. The spectral power density contains information about the source's optical absorption and excitation.⁴⁸ The signal's behavior with the help of the spectral power density can also be estimated. Hence, our results can be very useful for photoacoustic systems in order to analyze the tissue composition.

In this work, the representation theorems of the convolution type and the correlation type for the photoacoustic wave equation are recovered for lossless medium. The source term of the equation can be approximated by a Dirac delta function so that correlating measurements of the photoacoustic waves due to the excitation of random point sources can lead to the superposition of the Green's function and its conjugate (time reversal) inside a volume bounded by a surface. Even though it seems that direct measurement of the pressure waves at the locations of the point sources is an obstacle, in this work, it is shown that these pressure waves can be determined from the solutions of the photoacoustic wave equation by applying boundary conditions at the source locations in the frequency domain.

For photoacoustic systems, it is not possible to measure pressure waves at source locations throughout a volume by transducers. Even though the pressure waves cannot be measured directly at the source locations, they can be calculated utilizing the pressure waves measured at any locations different from the positions of the point sources and then using the continuity and the discontinuity of the pressure waves at the locations of sources. Inserting the calculated pressure waves into the expression obtained from the representation theorem of the correlation type, the spectral power density of the system excited by point sources is obtained. The spectral power density of a system has information about optical properties and excitation of the system and hence it is very useful for photoacoustic systems in order to analyze tissue composition. The spectral power density can also give an estimation about the behavior of the signal of the excited system.

The methodology given here can be used to obtain the spectral power density of the excitation due to two random sources and may be generalized to any finitely many number of point sources case. For any finite number of point sources case, by using the overall transfer matrix, the pressure wave at the innermost region is obtained from only the pressure waves measured at the outermost region without obtaining the waves between the outermost and innermost regions.

It has been shown that the velocity of a random distribution of micron-scale absorbers can be measured by the correlation of photoacoustic waves.³⁵ In other words, time shift measurements obtained by cross-correlation of photoacoustic waves make it possible to measure the blood flow in the microvasculature. Fluctuations in the photoacoustic signal resulted from the motion of the red blood cells enable the utilization of the cross-correlations of the photoacoustic waves. As a physical motivation, the red blood cells can be treated as spatially uncorrelated or random point sources.³⁶ 

The methodology can also be used to analyze a photoacoustic system excited by random sources throughout a volume bounded by a closed spherical surface and hence it can be useful for biomedical optics.

We would like to thank Osman Teoman Turgut for his valuable comments on the manuscript. This research is supported in part by Marie Curie Reintegration Grant 268287, Bogazici University Research funding BAP 7126, BAP 6033, BAP 6190 and TUBITAK Grant 112T253.