Modelization and optimized speckle detection scheme in photorefractive self-referenced acousto-optic imaging.

A photorefractive BSO single crystal can be used for axially resolved acousto-optic imaging of thick scattering media in absence of a reference beam. This configuration renders the experimental setup easier to realize for imaging through thick scattering media with an improved optical etendue. We present here a model and simulations that explains these results. It is based on the spatial heterogeneity of the speckle pattern incident on the crystal. Optimization of the detector position and of the speckle grain size is confirmed by the model.


INTRODUCTION
Optical imaging through thick biological tissues of several centimeters requires indirect methods in order to retrieve a millimeter resolution over the entire volume. The major drawback comes from scattering that generates multiple optical paths, so that the origin of the photons that propagate through the medium is lost. Many attempts have been made over the past few decades based on flux techniques such as diffuse optical tomography [1,2], but the extraction of the signal is not straightforward, since it requires solving an inverse scattering problem. A more recent concept couples light with an ultrasonic beam within the medium [3]. Dealing with the acousto-optic (AO) effect, the ultrasound (US) locally modulates the optical path, and thus the phase of light crossing the US field. This phase modulation mainly creates two weak optical sidebands (< 1%), shifted from the career frequency ω L of the US frequency +/ − ω U S . The magnitude of these so-called tagged photons depends on the optical properties of the medium, such as absorption and scattering. A 3D image of the optical properties of the medium is thus performed with a scan of the US beam. Mainly two approaches are implemented to detect the optical sidebands. The first one consists in a highly spectrally resolved flux measurement with a confocal Fabry-Perot interferometer [4], or in spectral holeburning filtering at liquid helium temperature [5]. A second approach uses a coherent detection of the tagged photons, initially performed by digital off-axis holography with camera detection (CCD, CMOS) [6,7]. This coherent detection requires accurate sampling of the speckle and the use of a fast camera, since the speckle decorrelates in living tissues on a millisecond time scale or less [8]. For these reasons, wavefront adaptive setups have been implemented with photorefractive (PR) crystals and large-area single detectors, operating much faster than a camera [9]. In dynamic PR holography, a refractive index grating is recorded within the volume of the crystal, resulting from the interference between the tagged photons and a plane wave off-axis reference beam. Consequently, the reference diffracts a spatial replica of the tagged photons' wavefront, which is time averaged according to the relative response time τ P R of the crystal (1-20ms), depending on light irradiance (a few 100mW/cm 2 ) [10]. The use of a reference beam to record the hologram generally brings an important coherent background noise (beam fanning), which limits the signal-to-noise ratio (SNR) of detection, since the AO signal is less than two orders of magnitude compared to the DC-flux components when dynamic axial resolution along the ultrasonic beam is considered. In this paper, we will show that AO imaging with PR detection can also be performed without applying a reference beam onto the crystal. This type of interaction was demonstrated by Kamshilin et al. [11] for application of PR correlation filtering to the detection of vibrating rough surfaces or vibrating multimode fibers. The proof of principle for AO imaging is illustrated by Benoità la Guillaume et al. [12] with a dynamic 2D image of ink inclusions embedded within 40mm of a multiply scattering phantom. The PR crystal used is the sillenite type Bi 12 SiO 20 (BSO), known for its high photosensitivity in conventional holographic experiments in the blue-green spectral region. We outline that this configuration is much simpler to implement after having optimized the light-collection efficiency and the speckle size incident on the crystal. The experimental setup detailed in Ref. [12] is shown in  We apply an excitation burst (carrier at 2.23 M Hz) on the US transducer to get a good axial resolution of 1.35 mm while the lateral resolution, given by the tranverse profile of the US, is 2 mm. The used gel phantom has a reduced scattering coefficient µ ′ s of about 10 cm −1 with two included china ink cylinders to simulate local absorbers. The detection of the AO signal is performed with a large area photodiode of 3.6×3.6mm 2 placed behind the optical system (L 3 L 4 ) that improves the collection of the speckle light. As analyzed in [12], the self referenced interaction can be interpreted as follows: in absence of US the incident speckled wavefront creates in the PR BSO a photoinduced space charge field which modulates the refractive index through the diffusion of carriers at no applied electric field. This complex phase volume hologram results from the interference of the incident speckle field with self-diffracted beam fanning noise fields. The speckle pattern photo-induces micrometer size random refractive index dots within the BSO volume. The operating conditions of the self-referenced PR holography for acousto-optic imaging are explained as follows: when the US transducer is off, the incident speckle field is perfectly phase matched with the recorded phase volume hologram. It thus results an optimum coherent diffraction on the detector. When the US transducer is now operating over short pulse duration the incident speckle structure is changed and it diffracts a wavefront which does not perfectly match with the complex phase hologram recorded in the volume of the BSO. In such conditions the diffracted intensity on the detector is reduced or enhanced as shown in [12] and permits an acousto-optic imaging of a phantom embedded in the scattering medium. It was also shown experimentally in [12] that in the case of self referenced holographic configuration an optimum signal detection is obtained when the two following conditions are fulfilled: first, the speckle grain size recorded in the PR crystal must be of the order of the Debye screening length. Second, the position of the detector in the transverse plane must be optimized due to the spatial frequency response of the photorefractive effect when speckle is involved.
The objective of the paper is to develop hereafter an analytic model and numerical simulations which confirm these optimum conditions for signal detection in the conditions of self referenced speckle beam holography with a PR crystal operation at zero electric field.

THEORETICAL MODEL AND NUMERICAL SIMULATIONS
According to the experiment described in [12], we first define the incident speckle field (i.e. the reference speckle field) and its intensity inside the BSO crystal respectively as: S r (r) and I r (r) = |S r (r)| 2 , where r is the position vector inside the PR medium. This random distribution of coherent light intensity I r photoinduces an internal space charge electric field in the whole volume of the PR crystal noted E sc (r). When the classical approximations are made in the PR model where electrons are considered as the single specy of charge carriers [13], the Ohm's law governing the current density and the rate equations of the PR charge carriers lead to the following differential equation [14,15] : where ∇ is the nabla operator, k B is the Boltzmann constant, T is the temperature, e is the elementary charge, k 0 is the Debye screening wave number (typically for BSO k 0 ≃ 4.65µm −1 [16]) and I d is a uniform background illumination responsible of dark current photoconductivity. Using the Fourier transform properties, Eq. (1) becomes: where FT denotes the 3D Fourier transform and k is the wavevector. Assuming that the variation of the PR variables along the propagation axis can be neglected and a two- (2) may be rewritten as: whereẼ scx (k t ),Ẽ scy (k t ) andĨ r (k t ) are respectively the 2D Fourier transforms, calculated in the (x, y) plane, of the transverse components of the space charge field and of the speckle intensity. k t = k x + k y is the transverse component of the wavevector and α = k B T I d e is a constant. These equations show that each internal space charge field component corresponds to the product of an odd function, presenting extreme values for ±k 0 , withĨ r (k t ) that corresponds to the autocorrelation function of the speckle intensity.
Following the analytical expression of the internal space charge field given by Eq. 3, numerical calculations are performed to illustrate the electric field components. Figures 2(a) and 2(b) represent respectively the intensity of a typical speckle pattern in the (x, y) spatial domain and its spectrum in the (k x , k y ) spatial frequency domain. Figure 2(c) depicts the modulus of the Fourier transform of the speckle intensity |Ĩ r (k t )|. The calculations are performed at the wavelength λ = 0.532µm and with the sampling steps ∆x = ∆y = λ 4 in the (x, y) plane. The speckles are simulated using a complex field amplitude with a real part and an imaginary part that obey to a gaussian random distribution. The 1 e 2 half width σ s of the speckle grains is adjusted by filtering, in the Fourier domain, the spectrum of the speckle with a gaussian filter. We verified that the statistical property of the intensity fluctuations are in agreement with the thermal statistic of speckles [17] where the probability density function for intensity verifies the distribution law p(I) = 1 I e − I I with I the mean intensity. We also verified that standard deviation of the spatial fluctuations of intensity is equal to the mean intensity : are considered here. In fig. 2(a), the mean size of the speckle grains is σ s = 0.21µm that corresponds to the Debye screening length such that σ s = l D = 1 k 0 . In the (k x , k y ) spatial 3, we assume that the pedestal of the autocorrelation function of the speckle intensity has a gaussian shape with a 1 e 2 half width σ k related to the size σ s of the speckle grains. We can easily demonstrate that the amplitude of the space charge field components exibit maxima for the spatial frequencies k x = ±σ k and k x = ±k 0 when the size of the speckle grains is respectively greater and smaller (or equal) than the Debye screening length.
The crystallographic orientation of the BSO crystal used in [12] is shown in Fig. 1  light is incident on the plane (x, y) where the x axis is perpendicular to the plane (0,0,1) and the y axis is perpendicular to the plane (1,-1,0). In [12], measurements of the angular dependency of the PR signal were performed with a photodiode placed in the far field of the PR crystal and that scanned the PR signal along the k x spatial frequencies where k y ≈ 0.
For this spatial frequency range, the internal space charge field has only a component along the x axis (i.e. theẼ scx (k t ) component) and two birefringence axes appear, one parallel to the laboratory x axis, and one parallel to the y axis. The refractive indices found by the light polarized along these axes and propagating in the directions such that k y = 0, are n 0 + ∆n x and n 0 + ∆n y where n 0 is the linear refractive index of the BSO crystal and [13]    ∆n x (k x ) = 0 As scattered light is strongly depolarized, we emphasize that only the component of the speckle field polarized along the y axis reads the modulation of the refractive index ∆n y induced by the componentẼ scx of the internal space charge field. These results suggest that the x component of the speckle field provides no contribution to the PR signal. Then, in the experimental setup, the use of a polarizer oriented along the y axis at the output of the BSO crystal could improve the signal to noise ratio of the PR signal. For the y axis polarization considered in this paper, BSO exhibits both photoinduced birefringence and large optical activity at the green wavelength. In such conditions wave progation of an incident depolarized speckle beam is complex and its detailed analysis may be out of the scope of this paper. When taking account of the optical activity for linearly polarized plane waves, it leads to a reduction of the BSO holographic efficiency [16].
Because the modulation of the refractive index is weak (10 −5 to 10 −6 ), the PR phase hologram induced by the speckle and read by the y polarized light can be expressed by : where L is the thickness of the hologram. Then, the y component of the output speckle field is given by : S out ry (x, y) = S in ry (x, y)H y (x, y).
Using the Fourier transform expression of the space charge field component given by the Eq. (3) and the expression of the phase function of the hologram given by the Eq. 5, the (k x , k y = 0) spatial frequency spectrum of the output speckle is given by: where ⊗ denotes the convolution product and δ(k x ) is the Dirac function centered at the zero spatial frequency and β a new constant including all the constants of the previous equations. We can see that the amplitude of the output spectrum of the speckle diffracted by the photorefractive hologram is made of the incident speckle field (i.e. the zero order of diffraction), and a second term corresponding to the diffraction of the incident speckle by its own hologram. From an algebric point of view, this term is of the same sign as k x and explains the origin of the reduction or the enhancement of the PR signal with respect to the position of the detector along the k x axis. We stress out that this term is null in the mean direction of the incident light (i.e. the zero spatial frequency) so that no PR signal can be recorded in this direction. This property can be physically explained if considering the speckle field within the crystal as a superposition of plane waves whose incidences are included within the aperture of the lenses that shine the cristal : a multi two-wave mixing process is performed along the mean propagation axis z, leading by symmetry to a photorefractive gain which can be positive or negative, depending on the sign of the angle of the plane wave with respect to z. As a consequence, the superposition of all these contributions cancels the effect along z.
Because of the small phase modulation induced by the hologram, the diffracted speckle is very similar to the reference speckle and the contribution of the PR signal in the output speckle is very weak. Then, the PR signal is obtained by subtracting the background level of light such that: In order to study the diffraction efficiency of the PR hologram, numerical simulations are performed with different sizes of the speckle grains varying in the range [0.1µm − 3µm] and speckles with a mean intensity kept constant. Figures 4(a) and 4(b) show, for different sizes of the speckle grains, the variation of the modulus of the x component of the space charge field and of the PR signal with respect to the normalized spatial frequency kx k 0 . In agreement with Eq. 3, the amplitude and the position of the maximum of the space charge field depend on the size of speckle grains. The best efficiency is obtained at the spatial frequencies k x = ±k 0 when the size of the speckle grains corresponds approximately to the Debye screening length (σ s ≈ 1 k 0 = 0.21µm). The experimental results reported in Fig. 2 of [12] show the angular dependency of the PR signal measured with a large area photodiode of 3.6×3.6 mm 2 placed at a distance of 25mm of the PR crystal (in this case the optical system L 3 L 4 in the experimental setup is removed). In order to fit with these experimental conditions, the PR signal curves presented in fig. 4(b) are the result of the convolution product between the PR signal calculated with Eq. 8 and a square function, with a width ∆k D = 1.7 µm −1 , that simulates the 0.14 rad angular aperture of the photodiode. Moreover, to smooth out the intensity fluctuations of the PR signal due to the small number of speckle grains in the area of the simulated photodiode, we average the calculated PR signal over more than 10 numerical realizations to obtain these curves. Finally, the PR signal is plotted simultaneously with respect to the normalized spatial frequency kx k 0 and with θ = asin kxλ 2π , the external angle (outside the crystal). In fig. 4(b), the shape of the PR signal is in a good agreement with the results presented in [12] where amplitude of the PR signal was maximized for a position of the detector corresponding to an external angle θ = ±17 • when the size of speckle grains was on the order of l D . These results also suggest that experimental parameters used in [12] were not optimized. Indeed, fig. 4(b) shows that the PR signal is maximized for larger speckle grains (0.4µm ≤ σ s ≤ 0.9µm) while the amplitude of the space charge field is maximized when the size of the speckle grains fits with the Debye screening length. Secondly, the efficiency of the PR signal is optimized when the detection is performed in a direction corresponding to a spatial frequency smaller than k 0 . Physically, this result differs from the two-wave mixing experiments where the photorefractive gain is optimized when the grating wave number K equals k 0 [16]. Indeed, in this case the terms related to the speckle in Eq. 7 will be replaced by Dirac functions δ(k x ±K) that leads to a diffraction efficiency which is governed only by the impulse response of the PR crystal βkx and maximized when K = k 0 . With speckles, convolution product in Eq. 7 implies that diffraction efficiency depends on the impulse response of the PR crystal as well as the spatial frequency width of the speckle σ k = 1 σs . Finally, finite size of the detector is also a factor to be taken into account to optimize measurement of the PR signal. Figure 4c shows the amplitude of the PR signal convolved with a square function with a width ∆k D and calculated as a function of the normalized spatial frequency kx k 0 and the normalized half width of the gaussian envelop of the speckle spectrum σ k k 0 . The black dotted curve depicts the position of the detector that gives the higher PR signal as a fonction of σ k . Vertical and horizontal black dotted lines make visible the couple of parameters (σ k = 0.20k 0 , k x = 0.35k 0 ) optimizing the PR signal. We can stress out that this optimized couple of parameters depends on the angular aperture of the detector. Now, in order to modelize the ultrasound modulated optical tomography by selfreferenced photorefractive holography reported in [12], let us consider the case where an incident speckle tagged by a short US pulse, noted S in t , is diffracted by the hologram induced by the reference speckle S r . Because the interaction between the short US pulse with the light propagating in the scattering medium leads to changes of the speckle pattern, the tagged speckle is transiently not perfectly matched with the complex phase hologram recorded in the whole volume of the BSO crystal. As a consequence, the level of the PR signal varies during the time of the US pulse.
In our simulations, the tagged speckle is modelised by the coherent superposition of the reference speckle S r with a "noise" noted S n (i.e. a fully decorrelated speckle) so that the amplitude of the tagged speckle may be written in the (x, y) domain as: where γ represents the ratio of the incident light interacting with the US pulse in the scatering medium. Then, the intensity of the tagged speckle is given by: Because the noise speckle S n is fully decorrelated with the reference speckle S r , the mean value of the interference term between the two speckles is null. Figure 5 shows the typical result of the interference between two speckles fully decorrelated. Therefore, this term can be neglected when γ ≪ 1 so that the intensity of the tagged speckle is assumed to be the incoherent sum of the intensities of the two speckles. According to this assumption, the detected intensity of the tagged speckle that reads the PR hologram after US interaction can be written, in the spatial frequency domain, as : Experimentally, the position of the photodiode is adjusted along the k x dimension to maximize the PR signal. Then, the measured AO signal is directly the intensity of the spectrum of the tagged speckle at this spatial frequency and the DC background corresponds to the intensity of the untagged speckle measured at the same spatial frequency. Hence, we will define the relative AO signal as : where ∆k D denotes the integration of the intensity of the spectrum of the diffracted speckle field over the spatial frequency bandwidth ∆k D of the photodiode. For numerical calculation, the mean radius of the speckle grains in the (x, y) is adjusted to 0.21 µm (i.e. the Debye screening length). Figure 6(a) represents the typical shape of the intensity of the spatial frequency spectrum of a tagged speckle and the two yellow squares symbolise the optimized position and the area of the detector centered on the spatial frequencies k x = ± k 0 2 (see Fig.  4b). For these two positions of the photodiode and using Eq. 12, we calculate the variation of the relative AO signals as a function of the level γ of the "noise" in the tagged speckle.
Assuming that the scattered light totally filled the volume of the scattering medium, the maximum value taken by γ is determined very approximately to 10 −3 by calculating the ratio between the focusing volume of the US pulse (∼ 5mm 3 ) and the total volume of 13 the scattering medium (< 40 3 mm 3 ). Figure 6(b) gives the variation of the relative AO signals, respectively noted R AO+ and R AO− when the detector is centered on + k 0 2 or − k 0 2 , as a function of the level of noise in the tagged speckle. Because of the small number of speckle grains in the area of the simulated detector, the curves correspond to the AO signals averaged over few tens of numerical realizations. In agreement with the shape of the PR signal given by the fig. 4(b) and the Eq. 11, the relative AO signal measured around + k 0 2 and − k 0 2 respectively decreases and increases linearly, when the level of noise increases. All these results give a satisfactory quantitative interpretation of the experimental results reported in [12]. However, the developed model mainly outlines the optimal values of the speckle grain size and of the detector position and, for the moment, does not permit to quantify the holographic efficiency and the AO signal dynamic range. In future work, several experimental parameters, such as fraction of tagged photons, beams polarization, light collection efficiency on the detector need to be carefully considered in order to refine amplitude consistency between the model and experiment.

CONCLUSION
The self referenced wavefront acousto-optic imaging technique has been analysed in the condition where a single speckle field is projected on a PR crystal. In these experiments, no reference beam is applied and the incident speckle creates a complex space charge field inducing a phase volume hologram in the BSO. The model leads to determine the parameters controling the efficiency of the hologram. We show that the PR signal efficiency is maximized for a spatial frequency that depends on the crystal Debye screening wave number, the width of the spatial frequency spectrum of the speckle but also on the angular aperture of the detector. It results that there is an optimum position of the photodetector and an optimum size of the speckle grains to detect the speckle pattern in presence of the tagged photons in an acousto-optic imaging experiment. The numerical simulations agree well with the experimental results previously reported. The method is robust and easy to implement in comparison with conventional holography. It can be applied with other types of nonlinear media in particular operating in the near IR region for acousto-optic imaging in highly scattering biological tissues and in turbid media and vibrometry. This work was supported by the LABEX WIFI (Laboratory of Excellence within the