Interferometric near-infrared spectroscopy directly quantifies optical field dynamics in turbid media

DAWID BORYCKI, OYBEK KHOLIQOV, AND VIVEK J. SRINIVASAN* Department of Biomedical Engineering, University of California Davis, California 95616, USA Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun, Poland Present address: Institute of Physical Chemistry, Polish Academy of Sciences, Kasprzaka 44/52, 01-224 Warsaw, Poland *Corresponding author: vjsriniv@ucdavis.edu


INTRODUCTION
Coherent light scattered from a turbid medium generates a random intensity distribution or speckle pattern [1,2], which fluctuates in time as the scatterers move.This effect serves as the basis for sensing and imaging methods based on dynamic scattering of coherent light (DSCL).In dynamic light scattering (DLS), the time scale of fluctuations of single scattered light, quantified using the intensity autocorrelation [3,4] or power spectrum [5,6], is used to determine scatterer dynamics.Assuming the timedependent field U s t d obeys zero-mean Gaussian statistics, the field autocorrelation g 1 τ d , which encodes scatterer dynamics, can be obtained from the intensity autocorrelation g 2 τ d via the Siegert relationship [7][8][9][10][11], (1) where I s t d jU s t d j 2 denotes the instantaneous intensity, and where τ d denotes the time lag for correlation, the brackets indicate ensemble averaging, and β ∈ 0; 1 is a parameter accounting for the number of measured speckles.If coherent detection is used (as in this work), β 1.
DSCL techniques have two main limitations.Firstly, the photon time-of-flight (TOF) is not measured.Accurate determination of particle dynamics from multiply scattered light is challenging if the TOF distribution is unknown.Secondly, conventional intensity-based DSCL methods are limited to systems in which the Siegert relation holds [3].If the number of scattering paths are few or correlated, Gaussian statistics do not hold, invalidating the Siegert relation [23].Additionally, the Siegert relationship cannot be applied to non-ergodic media [24], including samples with static components [10,25].
When the field reemitted by a turbid medium contains a static, time-independent, additive component U c , U s t d is written as where U f t d is a zero-mean Gaussian process.The interpretation of Eq. ( 4) is that U c does not change over experimental time scales, whereas U f t d does.After substituting Eq. ( 4) into Eq.( 2), g 1 τ d takes the form [10,11,26] (5) is the autocorrelation of the dynamic term, while I c U c U c , Īf G 1;f 0 and Īs hI s t d i t d I c Ī f .Assuming ergodicity of U f , ensemble averaging to determine G 1;f is equivalent to temporal averaging with respect to t d .
If β 1, the intensity autocorrelation is [10,11,26] By solving Eq. ( 6) for γ 1 τ d and assuming it is real and nonnegative, one obtains which yields a different estimation for the field autocorrelation than the Siegert relationship with β 1; nevertheless, both agree for η 0. Assuming coherent detection (β 1), η ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 − g 2 0 p , since γ 1 0 1.However, unless β is known a priori or fixed at a known value (e.g., 1) through careful experimental design [27], η cannot usually be estimated independently of β.
Approaches to address non-ergodicity include the use of higher-order intensity autocorrelations [11,28], ensemble averaging of independent speckles [24,26] by inducing sample motion [29], spatial diversity [10,16,30], insertion of an additional ergodic medium [31], or fitting based on an assumed functional form of γ 1 τ d [25].However, these approaches entail either added complexity or extra measurements.Yet another approach is to deliberately introduce a strong static component to increase η, such that the central term in Eq. ( 6) dominates and Reγ 1 τ d can be assessed [24].However, in this "heterodyne" approach, the relative fluctuation size is reduced [11] and η is artificially modified.Here, we present a straightforward approach to directly measure g 1 τ d from a single time series and obtain both η and γ 1 τ d from Eq. (5) (see Eq. (S9) of Supplement 1).
To accomplish this, we demonstrate interferometric nearinfrared spectroscopy (iNIRS) [32] for direct determination of TOF-resolved field and intensity autocorrelations.This is achieved through the analysis of the spectral interference between light traversing the sample and reference paths.This work extends beyond the intensity autocorrelations [32] demonstrated previously and determines field autocorrelations directly from the phase of the spectral interference.

MEASUREMENT OF FIELD AUTOCORRELATIONS
The spectral density is measured as a function of optical frequency, ν, over time, t d , by a Mach-Zehnder interferometer (MZI) with a narrow-linewidth tunable laser.The time for one laser sweep is assumed to be rapid compared to the time scale of the field fluctuations.The registered signal is then given by Sν; where S DC ν; t d is the DC component, assumed to be dominated by the reference arm, and W rs ν; t d is the crossspectral density between complex spectral amplitudes of the two paths [33].By the Wiener-Khintchine theorem, W rs ν; t d can be used to recover information about the sample field.Namely, after suppressing S DC ν; t d , and inverse Fourier transforming the resulting signal, one obtains where F −1 denotes the inverse Fourier transform, Γ rs is the mutual coherence function, and τ s is the conjugate variable to ν, representing the TOF delay between interferometer arms.U r and U s are hypothetical, not measured, reference and sample fields whose complex mutual coherence function, Γ rs τ s ; t d [inverse Fourier transform of W rs ν; t d ], relates to the complex transmitted sample field.We assume that the complex conjugate term Γ rs −τ s ; t d can be excluded.This, in practice, is achieved by adjusting the optical path mismatch between interferometer arms so the two terms in the last line of Eq. ( 8) do not overlap.
By repeatedly measuring Γ rs τ s ; t d over time t d (distinct and separable from τ s ), the TOF-resolved first-order iNIRS autocorrelation is determined using where by a convolution integral, where I 0 τ s denotes the instrument response function (IRF), which depends on the spectral range over which W rs is measured [32]: Δν ≈ cΔλ∕λ 2 0 (c is the speed of light, Δλ denotes the wavelength bandwidth, and λ 0 stands for the center wavelength in free space).Hence, iNIRS measures the convolution of the TOFresolved autocorrelation G 1 with the IRF.For large Δλ, the IRF approaches a delta function.Therefore, G iNIRS 1 extracts TOF-resolved dynamics.In addition, where I s τ s ; t d jΓ rs τ s ; t d j 2 and denotes convolution with respect to τ s .Thus, Ī iNIRS s τ s , the measured temporal point spread function (TPSF), is equal to the true photon TOF distribution, hI s τ s ; t d i t d , convolved with the IRF.
In the presence of a static component, I c , Eq. ( 11) can be written as After convolving with the IRF, one obtains [cf.Eq. ( 12)] where G iNIRS that decorrelates over time.Accordingly, Eq. ( 13) now reads as τ s ; 0. Therefore, by using temporal field autocorrelations, iNIRS distinguishes static and dynamic contributions.Finally, the TOF-resolved intensity autocorrelation can also be determined experimentally using iNIRS: Consequently, iNIRS measures both the TPSF and TOFresolved field and intensity autocorrelations, overcoming ambiguities inherent in intensity-based DSCL techniques.This allows verification of the Siegert relation, and potentially, the study of non-Gaussian and non-ergodic light scattering where Eq. ( 1) does not hold.
The iNIRS experimental setup is depicted in Fig. 1(a).The frequency-swept light is divided into reference and sample arms.The collimated beam from the sample arm irradiates the turbid medium, which attenuates, broadens, and delays the incident light distribution.The light paths are finally combined by a fiber coupler and detected by a differential detector.The digitized electronic signal approximates 2ReW rs ν; t d .Since the laser is swept in frequency, short photon paths through the sample produce smaller electronic beat frequencies than longer paths, provided that the reference arm TOF shorter than the ballistic sample arm TOF [Fig.1(b)].By 2ReW rs ν; t d to Γ rs τ s ; t d through Fourier analysis [Eq. ( 8)], the beat frequency, amplitude, and phase of the electronic signal [Fig.1(b)] yield the sample TOF, field magnitude, and field phase, respectively.This novel feature enables plotting TOF-resolved analytic field time courses [Fig.1(c)] based on Γ rs .As the field is complex, consisting of real and imaginary parts, these plots are naturally threedimensional.
The coherent addition of multiple light paths within the TOF range defined by the IRF [Eq.( 13)] leads to fluctuations in jΓ rs τ s ; t d j ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I s τ s ; t d p and φ s τ s ; t d argΓ rs τ s ; t d , as the phases of these light paths change over time [Fig.1(c)].The TPSF can be determined by incoherent averaging over t d at each value of τ s [Fig.1(d) and Eq. ( 12)].The field of ballistic photons does not decorrelate (although not present in this work, the field of photons with purely static scattering also does not decorrelate).The field of late photons decorrelates faster than that of early photons, as more scattering events lead to increased momentum transfer [Fig.1(e)].In the DWS limit [13], the field decorrelation rate linearly depends on the number of reduced scattering events, μ 0 s υτ s (with υ denoting the speed of light in the medium and μ 0 s standing for the reduced scattering coefficient).
The fluctuating complex field time courses (second column of Fig. 2) were used to estimate g iNIRS 1 [Eq.( 10)] and g iNIRS intensity-based estimate is noisier, particularly for long-photon TOFs [see Fig. 2(c)].As suggested by our simulations, the additional phase information in the optical field reduces the autocorrelation estimation error and bias (see Supplement 1).The non-ergodicity, which invalidates the Siegert relationship, was confirmed through a statistical analysis.The intensity and phase distributions for a complex Gaussian field with a static background of zero phase angle [Eq.( 4)] are given by [2] p I I s ξ exp−ξI s I c I 0 α; where ξ Ī −1 f , α 2ξ ffiffiffiffiffiffiffiffi I s I c p , I 0 α is a modified Bessel function of the first kind, and I c is the static intensity.Equation ( 16) is the modified Rician distribution.For I c 0, Eq. ( 16) reduces to the negative exponential distribution, p I I s ξ exp−ξI s , while Eq. ( 17) becomes p φ φ s 1 2π .Figure 5 depicts the intensity and phase distributions obtained for the c p 5.2% IL20 phantom at three fixed values of τ s .The intensity distributions were also fitted with modified Rician and negative exponential distributions using a maximum likelihood approach.For short paths, the intensity distribution obeys modified Rician statistics, and the negative exponential fit is poor.However, as the static contribution decreases with the increasing path length, the phases become more uniform [Fig.3(d)], the intensities approach the negative exponential distribution, and zero-mean Gaussian statistics describe the field.Thus, the statistical distributions in Fig. 3 confirm the autocorrelations in Fig. 2.
To further confirm that the static field component invalidates the Siegert relationship, the detection path was cross-polarized, suppressing ballistic paths.Figure 4 depicts the TPSF, complex field time courses, and autocorrelation estimates for c p 5.2%, where the invalidity of the Siegert relation was previously most apparent.Comparing cross-to co-polarization results [Fig.2(a)], the TPSF now comprises only a dynamic component [Fig.4(a)].Consequently, the transmitted field does not contain a constant term [Fig.4(b)].As a result, the validity of Eq. ( 1) is restored for all paths [see Figs.4(c) and 4(d)].Exclusion of the static component by polarization gating is also confirmed by statistical distributions in Fig. 5. Here, the measured light intensities, even for short paths, are exponentially distributed and the phases are uniform.

DISCUSSION
In iNIRS, as in other DSCL methods, multiple speckles must be observed in order to estimate the autocorrelations, with improved precision being achieved with more temporal speckle evolutions or "lifetimes" [36].However, our simulations (Supplement 1) suggest that fewer speckle lifetimes are required with field-based measurements to achieve the same precision as intensity-based measurements.Though results were obtained with a measurement time of 0.8 s in the main text and 0.2 s in Supplement 1, these times could be further reduced by employing improved estimators or utilizing both forward and backward laser sweeps.
It is important to note that the interferometer must be more stable for field-based than intensity-based measurements.As a rule of thumb, the phase of the interferometer must be stable over the decorrelation time scale of interest.In the present study, this criterion was satisfied without taking special measures to stabilize the interferometer.The advantages of field-based DSCL, including statistical efficiency and freedom from constraints of the Siegert relationship and ergodicity, are expected even if alternative methods (e.g., phase shifting) are used to determine the optical field.Importantly, iNIRS yields TOF-resolved field autocorrelations, from which sample dynamics can be derived through DWS [13] without explicit knowledge of the absorption.
In summary, we demonstrated iNIRS for the measurement of the intensity and field autocorrelations simultaneously with the TOF resolution.This approach to address non-ergodicity uses a single set of measurements on the light naturally reemitted from an intact sample.We used this unique capability to directly examine the Siegert relationship, assumed in many DSCL techniques.We found a breakdown of this relationship for small scatterer concentrations and short photon TOFs.Methods to measure the optical field will improve the efficiency and quantitative capabilities of DSCL techniques.
(a)].For each sample, N 40; 000 consecutive interference signals were acquired with a laser sweep duration of ∼10 μs at a sweep repetition rate of 50 kHz and an incident power of 25 mW.Each signal was processed to estimate Γ rs τ s ; t d , and temporal averages of I s τ s ; t d were calculated to obtain the TPSFs [first column of Fig. 2].As c p increases, I iNIRS c τ s decreases rapidly due to ballistic attenuation, while Ī iNIRS f τ s decreases more gradually.Therefore, with increasing c p , the relative contribution of the dynamic intensity, Ī iNIRS f τ s , increases.Note that I iNIRS c τ s has the same width as the IRF due to the convolution in Eq. (

Fig. 2 . 1 ,
Fig. 2. Transmitted field and intensity dynamics as a function of TOF and scatterer concentration for parallel polarization: (a) c p 5.2%, (b) c p 5.5%, and (c) c p 5.8%.The first column shows the TPSFs, with constituent static and dynamic components, and the second column shows complex field time courses for three τ s values, marked on the TPSF plots.The third column shows that the directly estimated field autocorrelations, g iNIRS 1