Tandem-pulsed acousto-optics: an analytical framework of modulated high-contrast speckle patterns

Recently we presented acousto-optic (AO) probing of scattering media using addition or subtraction of speckle patterns due to tandem nanosecond pulses. Here we present a theoretical framework for ideal (polarized, noise-free) speckle patterns with unity contrast that links ultrasound-induced optical phase modulation, the fraction of light that is tagged by ultrasound, speckle contrast, mean square difference of speckle patterns and the contrast of the summation of speckle patterns acquired at different ultrasound phases. We derive the important relations from basic assumptions and definitions, and then validate them with simulations. For ultrasound-generated phase modulation angles below 0.7 rad (assuming uniform modulation), we are now able to relate speckle pattern statistics to the acousto-optic phase modulation. Hence our theory allows quantifying speckle observations in terms of ultrasonically tagged fractions of light for near-unity-contrast speckle patterns.


Introduction
Both photoacoustics and acousto-optic tomography (Resink et al 2012) combine the use of light and sound in turbid media. Daoudi et al (2012) and Hussain et al (2012) demonstrated that the combination of both techniques makes it possible to perform fluence-compensated photoacoustic measurements, opening up the possibility of quantitative measurements of the optical absorption coefficient. In the underlying algorithm, the detected amount of ultrasonically tagged light is an important quantity. It would be ideal if both measurements can be done in the same setup, with the same laser system. This would give a more reliable method at lower cost, while increasing acquisition speed. However, the two techniques have different laser requirements. Photoacoustics (PA) needs short (typically ~5 ns), high-energy (>mJ ) pulses, and acousto-optics (AO) needs a long coherence length (>1 m) and a temporal resolution that allows for recording the dynamic behavior of speckles under the influence of ultrasound. Traditionally acousto-optic measurements use a quasi-CW (continuous-wave) laser so that a slow optical detector, often a CCD (charge-coupled device) camera, can integrate over multiple ultrasound cycles. A camera is used for its high number of optical detector elements so that acquiring information for a great number of speckles is possible, thus increasing the SNR (signal-to-noise ratio) (Leveque et al 1999).
Recently we have shown (Resink et al 2014a, b) that it is possible to combine the important properties for PA and AO in one laser system, thus reducing the complexity of the system. One of these methods is the speckle contrast method (Li et al 2002). A short laser pulse will give an instantaneous speckle pattern on the CCD and will thus have ideally a contrast of unity and will not give us any information on speckle dynamics. One option to perform a speckle contrast measurement is to integrate over a great number of short optical pulses, each recorded at a different phase of the ultrasound. However in a dynamic medium like biological tissue this requires a high pulse rate. Recently we have shown that AO probing can be done with two nanosecond laser pulses, each addressing a different phase of the ultrasound, where the speckle patterns due to each individual laser pulse are either subtracted or added. Here we investigate this new method and describe a way to quantify the signal such that the estimated amount of tagged light is comparable with the more traditional methods like speckle contrast that are well investigated (Li et al 2002, Zemp et al 2006. For clarity we show that the speckle contrast can be derived from the same equations as the AO signal in the sum-and-difference method that we presented in Resink et al (2014b).

The relation between phase modulation and the tagged fraction
By applying ultrasound to a small tissue volume the light will be locally modulated by small index of refraction changes and scatterer displacements. The light is injected into the turbid medium and propagates along a great number of different paths through the medium. A fraction of these paths will overlap with the ultrasound and will be phase modulated. Here we restrict ourselves to a linearly polarized speckle pattern. The electric field for a single polarization of path n can be written as where |E n | is the amplitude of the electric field, ω l the oscillation frequency of the unmodulated electric field, t the time, φ n is the phase of the light over path n at t = 0, and δ n is the phase modulation amplitude. The phase modulation has the frequency of the ultrasound ω us , which for brevity will be denoted ω. This phase oscillation itself has a phase ϕ n and is randomly distributed over the interval (0, 2π). When we investigate the influence of the phase modulation on the speckle pattern we can omit the fast oscillation of the electric field factor. The electric field then is described as: When δ n << 1 we can write equation (2) as a sum of phasors as shown in figure 1: one static phasor that represents nontagged light and two phasors that rotate with the ultrasound frequency, one clockwise (c) and one anticlockwise (a), leading to where E A n n is the amplitude of the stationary phasor. The phases φ c,n and φ a,n of the rotation of these phasors are chosen such that at t = τ n the two rotating phasors are in line with each other and perpendicular to the stationary phasor; hence c n n a n n n , , 0 The fraction of light r n that is tagged for this light path is then defined as the ratio of the energy of the tagged paths over the total amount of energy. Let us calculate the tagged fraction as a function of the phase modulation amplitude δ n . The modulus of one rotating phasor (clockwise or anticlockwise) is the square root of half the energy of the light that is tagged. The length of the stationary phasor A n is the square root of the energy that is not tagged, which is related to the energy of tagged light by energy conservation. The phase modulation amplitude can then be approximated as: The fraction of light that is considered tagged is the average value of r n , since some light paths might be strongly phase modulated, while other phasors remain stationary. Thus the average fraction R of tagged light is expressed as  We now have the relation between phase modulation and tagged fraction which will be used to derive the equations that will set the phase modulations for a given R in our simulations.

Modulation of speckle intensity by tagged light
Speckles are the result of the interference of n randomly phased electric fields E n (t). The intensity of a single speckle is given as where E tot (t) is the total electric field for the speckle which, with the help of equation (3) can be written as where E 0 is the sum of all stationary phasors, E c and E a are the results of the clockwise and anticlockwise rotating phasors respectively. The terms E 0 , E c , and E a have statistics resulting from a random walk in the complex plane, where the step sizes are given by δ E A E , , n n n n 1 2 and δ E n n 1 2 respectively. The phases and magnitudes of the resulting phasors are random because the phases of each of the initial phasors are also random due to the random paths the light travels in the sample. In the complex plane E tot (t) describes a ellipsoid with its principal axis under a random angle resulting from this process. The center of this ellipsoid is at E 0 and the length of the major and minor axis are |E c + E a | and |E c − E a |.
By substituting equation (9) in equation (8) we obtain The second harmonic term appears as cross terms from E a and E c . The magnitude of I 2ω is small compared with I 1ω but has comparable amplitude to a a and will not be neglected yet. We now have an expression that gives us the intensity of a single speckle as a function of time.

The effect of tagged light on the speckle contrast
The contrast of a speckle pattern is defined as where σ is the standard deviation of the speckle pattern with intensity distribution ⃗ I x ( ) and ⃗ x is the spatial average operator, hereafter denoted as . We recognize in I 0ω as given by equation (11) the summation of three statistically independent speckle patterns that are stationary in time. The phase relation between E 0 , E c , and E a cancels out for I 0ω . Further, I 1ω and I 2ω are periodic over 2π and π respectively and have an average value of 0.
In the quasi-CW AO speckle contrast method the intensity is averaged over an integer number of ultrasound cycles, and the contributions of the two dynamic terms in I(t) are therefore negligible. The speckle contrast becomes I I I  I I I   I  I  I  II  II  II   I  I  I  I I  I I  I The relative contributions of I 0 , I c , and I a are determined by the tagged fraction R. The amount of tagged light is equally divided over I c and I a , and the untagged light I 0 is the remaining part; thus, when I 0ω is normalized, For speckle patterns with a contrast of unity we are allowed to use and for uncorrelated speckle patterns we can use Combining equations (13)-(16) we obtain an expression that links the tagged fraction to the contrast for ideal speckle patterns integrated over an integer number of ultrasound cycles for small fractions of tagged light: For R << 1 we obtain Equation (19) is similar to the result from the calculations described in Li et al (2002) for small tagged fractions and close-to-unity contrast situations. So we have an equation describing the speckle contrast as a function of the tagged fraction.

The AO difference method
In the difference method, instead of integrating the speckle pattern over an integer number of ultrasound cycles we use two short light pulses at opposite ultrasound phases. We obtain speckle patterns ⃗ I x ( ,0) at time t = 0 and ⃗ π ω I x ( , / ) at t = π/ω, the difference between which is induced by the ultrasound. Both speckle patterns are normalized. In equations (10) and (11) there is a background term I 0ω that, being static, does not change between these two speckle patterns. I 2ω is small and periodic with half the ultrasound period and does not contribute to the speckle differences. I 1ω has an opposite sign between the two speckle patterns and is thus maximally different.
Using equation (11) the speckle difference can be written as We want to express the influence of the ultrasound as a single number instead of a complete spatial pattern. Therefore we want to perform spatial averaging, however the spatial average of equation (20) is zero. Let us investigate the mean square difference, x cx Further, the cross term has a spatial average of zero, so by substituting equation (22) in (21) we get the relation between the tagged fraction R and the mean square difference of two speckle patterns: (1  ) 8 (1  ) 8 .
This expression provides an estimate of the tagged fraction R from the difference between two instantaneous speckle patterns for opposite ultrasound phases.

Materials and methods
We test the relations given by equations (7), (18), and (24) on simulated speckle patterns with modulation. On the output optode of a sample we define 2500 wiggling phasors on random positions and with a random phase. The phase modulation of these phasors is simulated using equations (2) and (8). For each pixel of the simulated camera the distance and thus the phase relation between the pixel and the origin of all 2500 output positions of the optode is calculated. Then for each pixel we calculate the intensity associated with the sum of all the electric field components. This gives us one instantaneous speckle pattern. For the speckle contrast simulations we take 20 snapshots of the speckle pattern equally spaced over one ultrasound cycle. We vary the tagged fraction R from 0 to 0.8 and use random phase variation amplitudes δ n chosen such that equation (7) is satisfied for several phase modulation distribution functions of δ n . We test eight scenarios which are plotted in figure 2. For each value of R and each distribution function a histogram is generated with 51 bins of δ n for values between -π and π. We fill each histogram with 2500 randomly chosen values for δ n within the constraints of the distribution function we test. The color scale represents the number of values of δ n within that bin. Four classes of amplitude modulation distribution functions make up these eight distributions. The first class has one phase modulation distribution and assumes all values for δ n are equal, and is shown in the first row of figure 2. The second class assumes a phase modulation distribution where all the modulation amplitudes are randomly chosen between a lower and an upper limit. These distributions are designed such that either the probability of finding a phase modulation is equal over the whole range or its square is equally distributed. These are shown in rows 2 and 3 of figure 2 respectively.
The third class assumes one part of all the light is not modulated and thus has cases in which δ n = 0. The remaining light is modulated with values for δ n chosen equal. The modulated fractions of 10%, 30%, 50%, and 70% are shown in rows 4, 5, 6, and 7 respectively.
For the fourth class of distributions we assume that the values of δ n have a Gaussian distribution. The distribution is shown in row 8 of figure 2. We expect the distribution to be more similar to that of class 3 when the tagging is close to the optodes. A fraction of the light can escape the medium without interacting with the ultrasound when the optode size is larger than the ultrasound focus. The Gaussian distribution might resemble a more realistic phase modulation distribution and we assume it most closely resembles reality for most experiments. In the equations of figure 2, ξ denotes a random number: ξ −1…1 is a random number from a evenly distributed set of numbers between −1 and 1, ξ 0…1 is evenly distributed between 0 and 1, and ξ Gauss is a Gaussian distribution with a standard deviation of 1 and an average of 0.
Using equation (2) we calculate the electric fields that form the speckle pattern. The light interferes on a virtual CCD array of 300 × 300 pixels that pointwise-samples the speckle pattern. The average speckle size is 3 pixels resulting in approximately 10 4 speckles. To calculate the speckle contrast we add 20 instantaneous speckle patterns taken equally spread over one full ultrasound cycle to simulate integration on the camera over a finite time. The simulated speckle contrast is used to estimate the tagged fraction using equation (18).
The instantaneous speckle patterns at ω = t 0 and ω π = t are stored for analysis of the two speckle patterns. For the difference method we plot equation (24) as a function of the set tagged fraction R that also determines δ n via equation (7) for the simulated speckle patterns, and compare this against equation (24).

Results
To illustrate the simulations we show speckle patterns for the 0 and π phase shifts in figures 3(A) and (B) respectively. The difference between the speckle patterns is plotted in figure 3(C) while the blurred speckle pattern that is integrated over one ultrasound cycle is shown in figure 3(D). All the wiggling phasors have an identical phase modulation amplitude of 0.35 rad, leading to a tagged fraction R = 0.058.
For each phase modulation function the simulated tagged fraction was determined using equations (18) and (24) for the contrast and difference methods for a range of set tagged fractions. For the Gaussian phase modulation distribution the results are plotted in figure 4.
In figure 4 the simulated tagged fraction is plotted versus the set tagged fraction. An ideal method would follow the line y = x. Both methods correctly predict the tagged fraction for values between 0 and 0.2 while the simplified methods represented by equations (19) and (23) only predict R correctly for a smaller range. For tagged fractions above 0.2, we see that both methods start to fail to reproduce the set tagged fraction. The contrast method performs better Here the simulated tagged fraction is shown for the contrast method (R c , red +) and difference method (R diff , blue ×). Data obtained with the approximated versions of equations (19) and (23) 2 ) are shown in black. Phys. Med. Biol. 60 (2015) 4371 for high tagged fractions than the difference method. For low tagged fractions the difference method is superior to the contrast method and shows much smaller statistical scatter. We test both methods for several distributions of phase modulation amplitudes that satisfy equation (7). The result is shown in figure 5 as the fractional difference from the ideal behavior for the speckle contrast method (left) and the mean square difference method (right).

S G Resink and W Steenbergen
The error bars depict the standard deviation over 15 simulations. The red line is calculated by substituting the result of equations (17)

Discussion and conclusion
We have theoretically explored the relationship between phase modulations of light and the resulting speckle pattern modulations. For a certain light trajectory, we modeled tagged light as two counter-rotating phasors that are necessary to approximate a periodic phase modulation. The light represented by these rotating phasors is a fraction of the total light that forms the speckle pattern on the detector. This is the fraction of tagged light. The relation between the fraction of tagged light and the phase modulation is given in equation (7). We explored the effect of phase modulations on the difference between speckle patterns for a π rad ultrasound difference (difference method), and on the contrast of the speckle pattern integrated over a complete ultrasound cycle (contrast method). We showed the relationship between the signal obtained with the difference method (equation (24)) and the speckle contrast method (equation (18)) for ideal simulated speckle patterns. These relationships were tested for simulated speckle patterns, with several phase modulation angle distributions. The analytical expressions allow for extraction of tagged fractions from speckle observations, and for simulated speckle patterns these predictions were shown in agreement with the input values for fractions of tagged light up to 0.2. For very small fractions the error bars of the contrast method in particular are large because of the stochastic properties of speckle patterns. The natural spread of speckle contrast is on the order of the contrast reduction caused by the phase modulation due to the random nature of the speckle contrast (Duncan et al 2008). Therefore the difference method is expected to perform better for small acousto-optic signals in experiments. The trend of the contrast method for low tagged fractions does however obey the model. For larger fractions R of tagged light we often see an over-or underestimation compared to the ideal behavior. This is because the assumption in equation (3) of small values for δ n does not hold. A good example of this is given in the fourth row of figure 5. In that case only a small fraction of all the light is modulated. That small fraction of light must be highly phase modulated to achieve the same amount of tagged light on the camera. When the phase modulation is large the counter-rotating phasor model represented by equation (3) no longer holds. The simulation results depend on the phase modulation distribution, and we conclude that the model holds a tagged fraction of up to 20% for most tested distributions. The valid range of the model can be extended by making more accurate approximations and adding more terms to equation (3). However, considering the simplicity of the counter-rotating phasor model, the range of validity is remarkably large. When aiming for high signal strengths, tagged fractions R in the range 0.1-0.2 are associated with contrast values in the range 0.8-0.9 (according to equation (18)), and hence reductions in contrast of approximately 10-20%. Such relative ΔC values have been observed in Resink et al (2014b), which suggests that the theoretical framework presented here will provide a valuable quantification framework for situations encountered in practice.