High-sensitivity multispeckle diffuse correlation spectroscopy

Abstract. Significance: Cerebral blood flow is an important biomarker of brain health and function as it regulates the delivery of oxygen and substrates to tissue and the removal of metabolic waste products. Moreover, blood flow changes in specific areas of the brain are correlated with neuronal activity in those areas. Diffuse correlation spectroscopy (DCS) is a promising noninvasive optical technique for monitoring cerebral blood flow and for measuring cortex functional activation tasks. However, the current state-of-the-art DCS adoption is hindered by a trade-off between sensitivity to the cortex and signal-to-noise ratio (SNR). Aim: We aim to develop a scalable method that increases the sensitivity of DCS instruments. Approach: We report on a multispeckle DCS (mDCS) approach that is based on a 1024-pixel single-photon avalanche diode (SPAD) camera. Our approach is scalable to >100,000 independent speckle measurements since large-pixel-count SPAD cameras are becoming available, owing to the investments in LiDAR technology for automotive and augmented reality applications. Results: We demonstrated a 32-fold increase in SNR with respect to traditional single-speckle DCS. Conclusion: A mDCS system that is based on a SPAD camera serves as a scalable method toward high-sensitivity DCS measurements, thus enabling both high sensitivity to the cortex and high SNR.


Introduction
Measuring cerebral blood flow noninvasively and with high sensitivity is critical for clinical applications such as measuring the oxygen metabolic rate 1,2 and monitoring intracranial pressure. 3,4 Furthermore, although neuroscience applications such as functional activation mapping 5,6 and noninvasive brain-computer interface 7,8 have been pursued primarily using functional magnetic resonance imaging and near-infrared spectroscopy (fNIRS), such applications could in principle benefit from functional cerebral blood flow measurements. [9][10][11] Diffuse correlation spectroscopy (DCS) 12 is a promising noninvasive optical technique for monitoring cerebral blood flow 13,14 and for measuring cortex functional activation during finger tapping 9 and visual stimulation 10,11 tasks. DCS measures deep-tissue dynamics by coupling coherent light into the subject and measuring the fluctuations in the speckle field created by the light diffusing out of the subject. 12,15,16 Increasing the source-detector separation (ρ) of DCS optodes increases the proportion of detected photons that travel beneath the scalp and skull, deep into the brain cortex. [17][18][19] However, the increase in sensitivity to deep tissue comes at the expense of a decrease *Address all correspondence to Edbert J. Sie, edbertsie@fb.com; Francesco Marsili, fmarsili@fb.com in the signal-to-noise ratio (SNR) of the measured signal. Indeed, fewer diffuse photons reach the detector when the detector is farther away from the source due to absorption, scattering, and radial spread from the source. To increase the brain sensitivity of DCS, we developed multispeckle DCS (mDCS), a method that can extend ρ without sacrificing SNR by performing thousands of independent speckle measurements in parallel.
When we inject coherent light into a dynamic scattering medium such as those shown in Figs. 1(a) and 1(e), a dynamic speckle pattern emerges as shown in Fig. 1(b). DCS estimates Fig. 1 Experimental setup for multispeckle DCS. (a) Experimental setup consisting of a highcoherence-length (∼9 m) CW laser emitting at 785 nm wavelength, a diffuser plate, and a 32 × 32 pixel SPAD camera detector. The laser output was coupled to the phantom using an SMF. The diffused photons were collected by an MMF in a transmission geometry and coupled to the SPAD camera. We could change the fiber-SPAD distance (z) with a translating lens tube. The setup was used to calibrate the diameter of the speckles. The SPAD camera consisted of 1024 SPADs arranged in a 32 × 32 array, with a pixel pitch of 50 × 50 μm and a pixel active area of 6.95 × 6.95 μm. To measure the speckle turnover time trace shown in Fig. 1(c), we mounted the diffuser plate on a motorized rotation stage and rotated the plate at an angular speed of 15 deg ∕s. (b) Snapshots of speckle patterns recorded with the SPAD camera at varying z. These snapshots are recorded using the diffuser plate setup shown in Fig. 1(a). By adjusting the fiber-SPAD distance (z ¼ 106, 59, 30 mm), we can tune the average diameter of the speckles (d ¼ 8, 4, 2 pixels). The input laser power was tuned to avoid saturating the SPAD camera. Replacing the diffuser plate with a milk phantom did not change the diameter of the speckles, as long as we used the same fiber core diameter (D) and fiber-SPAD distance (z). The white circles denote the average diameter of the speckles that we estimated with 2D spatial autocorrelation (Appendix B). The frame exposure time was T bin ¼ 10 μs. (c) Detected photon counts versus time from a single pixel as shown in (b), recorded by rotating the diffuser plate around its axis. This time trace is the raw signal for calculating the intensity autocorrelation function g 2 ðτÞ. To resolve the dynamics of the speckles in this setup, we used T bin ¼ 10 μs, which is shorter than the correlation time of the speckles τ c ¼ 557 μs. (d) Normalized intensity autocorrelation function g 2 ðτÞ of the photon counts shown in (c). The g 2 curve was normalized by the value of g 2 at τ ¼ 10 μs. The correlation time constant τ c (¼557 μs) was estimated by fitting the data with the function g 2 ðτÞ ¼ 1 þ β expð−τ 2 ∕τ 2 c Þ for a ballistic motion, where β is the coherence factor 15,20 [see Appendix E for more details on the functional form of g 2 ðτÞ]. We normalized g 2 ðτÞ − 1 and used τ c as a fitting parameter. (e) Liquid phantom setup consisting of a source (SMF) and detector (MMF) optodes arranged in a reflection geometry. The source-detector separation (ρ) was 11 mm. This setup was used to measure the SNR gain and validate our mDCS model. the speed of the scatterers by measuring the correlation time constant of the speckles (τ c ), which depends on the speed of the scatterers. 16 We can estimate τ c by measuring the intensity of a speckle versus time and by calculating the corresponding intensity autocorrelation function (g 2 ), as shown in Figs. 1(c) and 1(d).
Based on the correlation noise model, 21,22 we can increase the SNR of the g 2 curve [defined as the ratio between mean and standard deviation of g 2 ðτÞ calculated over several integration periods] by increasing the detected photon count rate (N ph ) and the integration time (the time over which the intensity of the speckle field was recorded for each g 2 calculation, T int ). However, T int is limited by the time scales of the dynamics we need to measure (for example, T int should be ≲100 ms to measure pulsatility), and N ph is limited by the skin and eye maximum permissible exposure (MPE) to laser radiation. 23 Therefore, to increase SNR we need to consider parameters other than T int and N ph .
Traditional DCS systems employ single-mode fibers (SMFs) to couple one speckle onto a single-photon detector. There have been attempts to increase N ph using a multimode fiber (MMF) to couple multiple speckles onto the same detector. 21,24,25 Measuring multiple speckles on a single detector indeed increases N ph linearly with the number of detected speckles, but at the same time it decreases the magnitude of g 2 − 1 (¼β) by the same factor. Under typical experimental conditions for DCS, the increase in N ph and the decrease in β compensate each other, so measuring multiple speckles with the same detector does not improve the SNR of g 2 (see Sec. 3.2 for more details on the SNR of DCS measurements). Therefore, we will henceforth refer to the above systems (one speckle/one detector and multiple speckles/one detector) as single-speckle DCS systems. mDCS is a new promising method to improve on the SNR of DCS through multiple independent measurements of the speckle field. The mDCS technique has been pursued in various ways using 28 stand-alone single-photon avalanche diodes (SPADs), 10,26 8 stand-alone SPADs, 27 4 pixels of a 5 × 5 SPAD array, 28 and interferometric near-infrared spectroscopy detecting 20 speckles with a CMOS camera. 29 While these prior approaches are important proofs of principle of mDCS, the largest improvement in SNR reported so far is an SNR gain (defined as the ratio between the SNR achieved with mDCS and with single-speckle DCS using detectors with the same performance) of 5. 26 Here, we report on a method to perform 1024 independent DCS measurements using a kilopixel SPAD camera, which provides an SNR gain of 32. Our technique allows for DCS measurements with 36 times more speckles and approximately six times higher SNR gain with respect to the current state-of-the-art for mDCS. Our technique is ∼3 to 32 times more sensitive than single-speckle DCS using state-of-the-art high efficiency (∼70% at 780 nm) single photon detectors, depending on the count rate. With the ongoing development for megapixel SPAD technology, 30,31 our approach offers a scalable implementation toward mDCS instruments with >10;000 times more speckles and >180 times higher SNR gain than the current state-of-the-art and >100 times higher SNR gain than single-speckle DCS.

Experimental Setup
Figures 1(a) and 1(e) show our experimental setup. We coupled a 785-nm wavelength CW laser to the scattering media with an SMF. We collected the light diffusing out of the scattering media with an MMF coupled to a SPAD camera (PF32, Photon-Force). The scattering media were a diffuser plate to calibrate the diameter of the speckles [ Fig. 1(a)] and a liquid phantom to measure the SNR gain [ Fig. 1(e)] (see Appendix A for details on the phantom setup).
To ensure that each pixel of the SPAD camera provided an independent measurement of the speckle field, we adjusted the average diameter of the speckles on the SPAD camera [d, see Fig. 1(b)] to be equal to or smaller than the size of the pixels of the SPAD camera (a ¼ 6.95 μm, see Appendix A for details on the SPAD camera). The average diameter of a speckle obeys the following relationship: 32 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 9 9 d ¼ where λ is the wavelength of the light (785 nm), z is the distance between the detector fiber end and the SPAD camera, and D is the fiber core diameter (910 μm). We calibrated our system by measuring d at different z with the setup shown in Fig. 1(a) (see Appendixes B and C for details on the speckle diameter calibration). Based on Eq. (1), the diameter d does not depend on the particular medium used, as long as D and z are fixed. We measured the SNR gain of our mDCS system using a liquid phantom in a reflection geometry, as shown in Fig. 1(e), to better mimic the dynamics and optical properties of human tissue. The fiber-SPAD distance was z ¼ 8.1 mm, which corresponds to d ¼ a. For the SNR gain measurements, the SPAD camera recorded the dynamics of the speckle field with a frame exposure time of T bin ¼ 4 μs for up to 2 × 10 6 frames (8 s).

Calculating the Intensity Autocorrelation Functions
We calculated the intensity autocorrelation function g 2 as a function of time lag τ for each pixel (i) of the camera: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 5 5 4 where nðtÞ is the number of detected photons in time bin t, τ is the time lag, and the square bracket h: : : i denotes the average over the integration time T int . 33,34 We calculated the ensemble average of the g 2 curves measured over M pixels as E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 6 ; 4 7 3 3 Results and Discussion 3.1 SNR Gain using a SPAD Camera Figure 2(a) shows a representative g 2 i curve, and Fig. 2(b) shows the ensemble average measured over all of the pixels, g 2 j M¼1024 . We measured g 2 with T int ¼ 50 ms and repeated the measurement for up to 8 s, for a total of 160 g 2 measurements per pixel. We calculated the time average fMEAN½g 2 ðτÞg and standard deviation fSTD½g 2 ðτÞg of g 2 over all of the 160 integration periods.  Figure 2(c) shows the measured SNR½g 2 ðτÞ ¼ MEAN½g 2 ðτÞ − 1∕STD½g 2 ðτÞ for increasing M (solid circles) and the SNR calculated using our mDCS noise model (to be discussed below) with no fitting parameters (solid lines). The SNR increased with increasing ensemble size, which is consistent with the prediction of our model. We estimated the SNR gain by calculating the ratio SNR½g 2 ðτÞ∕SNR½g 2 i ðτÞ for the first bin (τ ¼ 4 μs). Figure 2(d) shows that the measured SNR gain (solid circles) increased with the number of pixels as ffiffiffiffi ffi M p (solid line) and reached a maximum value of 32 when we used all of the 1024 pixels of the camera, in agreement with our prediction.

Multispeckle DCS Noise Model
Zhou et al. 21,22 reported on a noise model for single-speckle DCS that allows for estimating the SNR of g 2 from several experimental parameters that can be measured independently. To evaluate our findings quantitatively, we extended the noise model to mDCS by assuming that each speckle of the speckle field is an independent realization of the same ergodic random process (see below for details of the model). For ergodic random processes, ensemble and time statistics are the same, so observing one speckle for T int is equivalent to observing M speckles for T int ∕M. Our model predicts that the SNR of the first bin in the shot-noise limited regime, in which most in vivo DCS measurements operate (see below for details of the noise components of the speckle field).
Under the assumption of speckle ergodicity, we extended the single-speckle DCS noise model 21,22 to mDCS by including the number of independent speckle field measurements (M) alongside the integration time (T int ) in the expression for STD½g 2 ðτÞ: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 4 ; 1 1 6 ; 1 6 5 where T bin is the frame exposure time (which is equal to the correlator time bin interval), T int is the integration time, M is the number of independent speckle field measurements, β is the Fig. 2 Achieving SNR gain of 32 using a 32 × 32 pixel SPAD camera. For this measurement, we used a CW laser at 785-nm wavelength coupled to the liquid phantom shown in Fig. 1(e). We used input power of 33 mW and ρ ¼ 11 mm to demonstrate the SNR gain using a SPAD camera. Under operating conditions using the MPE limit for skin (28.5 mW at 785 nm) and larger ρ, the SNR gain remains the same but the SNR magnitude depends on the measured count rate as we discussed in the next section. The fiber-SPAD distance was adjusted to z ¼ 8.1 mm to satisfy the condition d ¼ a (see Appendix B for the estimation of d and Appendix C for the effect of d on the SNR 35,36 ). (a) Time statistics of the pixel intensity autocorrelation function g 2 i ðτÞ. The solid line corresponds to the mean value of g 2 i ðτÞ over 160 integration periods (T int ¼ 50 ms) and the shade corresponds to AE2 standard deviations. The measured β ¼ 0.2 is lower than the value expected from the case of unpolarized light (β ¼ 0.5) because the speckle diameter matched the pixel active area (d ¼ a) 35 (see Appendix C). (b) Time statistics of the ensemble intensity autocorrelation function g 2 ðτÞ.
The solid line corresponds to the mean value of g 2 ðτÞ over 160 integration periods (T int ¼ 50 ms). The standard deviation of g 2 is 32 times smaller than that of g 2 i and not as apparent in the plot. (c) SNR of g 2 ðτÞ for increasing size of the pixel ensemble M ¼ 1, 16, 64, 256, 512, and 1024. Averaging over a larger ensemble of pixels leads to higher SNR. The SNR measured over 160 integration periods (solid circles) is compared with the SNR calculated using the mDCS noise model with no fitting parameters (solid lines). All of the parameters of the noise model were obtained from experiments. (d) SNR gain of g 2 ðτ ¼ 4 μsÞ as a function of the size of the pixel ensemble M ¼ 1 to 1024. The SNR gain measured over 160 integration periods (solid circles) increases as the square root of the ensemble size (solid line), as predicted by our model. coherence factor, τ c is the speckle correlation time, m is the delay time bin index such that τ ¼ mT bin , and hni (¼N ph × T bin ) is the number of detected photons within T bin per pixel. To derive the expression for STD½g 2 ðτÞ shown in Eq. (4), we simplified the intensity autocorrelation function with g 2 ðτÞ ≈ 1 þ β expð−τ∕τ c Þ, which is a good approximation under the operating conditions of our liquid phantom measurements [see Appendix E for more details on the functional forms of g 2 ðτÞ].
We write Eq. (4) in the following form: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 6 ; 6 5 1 where A, B, and C are functions of τ and other correlation parameters. In the limit of low photon q .
Therefore, for low count rates, the SNR is proportional to , and it is limited by the shot noise of the speckle field (fluctuation in the number of photons in a speckle). At higher photon count rates (hni ≫ ffiffiffiffiffiffiffiffiffi A∕C p ), C dominates, which leads to . Therefore, for high count rates, the SNR is limited by the speckle noise (fluctuation in the coherence time of a speckle) and independent of the N ph . The SNR transitions between the shot-noise and speckle-noise limited regimes at the threshold count rate, , which is apparent from the gradual changes in the slope of the SNR versus N ph curves as shown in Fig. 3(b).
The mDCS noise model is based on the assumption that the dominant source of noise is the intrinsic statistical temporal fluctuation of the speckle field, which is due to shot noise and speckle noise. Thus, the validity of the noise model is limited to measurements with negligible extrinsic noise sources, such as dark counts from the detector, after-pulsing, pixel cross-talk, and background photons. Depending on the relative magnitudes of extrinsic and intrinsic noises, increasing the number of speckles coupled to a pixel (d < a) may have different effects on the SNR of g 2 . If shot noise is dominant, measuring multiple speckles on a single detector increases N ph linearly with the number of detected speckles, but at the same time it decreases β by the same factor. Since these two effects approximately compensate each other, measuring multiple speckles with the same detector does not improve the SNR of g 2 . If speckle noise is dominant, the number of speckles coupled to a pixel does not affect the SNR because the SNR does not depend on N ph or β. Finally, if extrinsic noise is dominant, coupling multiple speckles to a pixel may increase the SNR by increasing N ph .
Since τ c in brain tissue is in the range of tens of μs, which becomes shorter for larger sourcedetector distances ρ, 26 most in vivo DCS measurements are likely to operate in the shot-noise limited regime. For instance, for T bin ¼ 1 μs, τ ¼ 1 μs, τ c ¼ 10 μs, and β ¼ 0.5, the threshold count rate is N th ∼ 10 5 cps per pixel, which is over an order of magnitude larger than the count rates reported for DCS measurements on an adult human head at ρ ¼ 29 mm (∼10 4 cps per speckle). 26 To support our ergodicity assumption and validate our model experimentally, we compared the SNR of g 2 and g 2 i at τ ¼ 4 μs for increasing T int or N ph . Figure 3(a) shows the SNR versus T int with different M ¼ 1, 32, and 1024. Increasing the ensemble size by 32 times from 1 to 32 at T int ¼ 10 ms (see arrow from point A to point B) resulted in the same SNR gain that would be achieved with one pixel and 32 times longer T int (see arrow from point A to point C), as predicted by our model. This result supports our ergodicity hypothesis that pixel ensemble averaging and time averaging result in the same SNR gain. Similarly, the SNR versus N ph curves in Fig. 3(b) show that increasing the ensemble size by 878 times is equivalent to increasing the count rate of one pixel by ∼30 times. The solid lines in Fig. 3 are predictions based on our model with no fitting parameters.
The measurements shown in Fig. 3(b) also illustrate how our ensemble averaging technique is applicable to photon-starved applications, such as the detection of deep-tissue dynamics. Figure 3(b) shows that the SNR measured with only one pixel and ∼100 detected photons per T int (point D) is comparable to the SNR measured by averaging over 878 pixels and as low as ∼4 detected photons per pixel per T int (point E). mDCS reached the predicted SNR gain even in the case of an unprecedentedly small number of detected photons available to calculate g 2 for each pixel.
In summary, Figs. 2(c), 2(d), 3(a), and 3(b) show that our mDCS noise model (solid lines) could predict the experimental results (solid circles) well without any fitting parameters. In particular, we confirmed experimentally the mDCS noise model prediction that the SNR is propor- in the shot-noise limit. Our results confirm that we can increase the SNR beyond the limits imposed by the integration time T int and the photon count rate N ph by increasing the ensemble size M.

Conclusions and Outlook
We have demonstrated a scalable method for mDCS measurements to enhance the SNR gain by a factor of 32 with respect to single-speckle DCS using a kilopixel SPAD camera. Our technique can be used to achieve higher sensitivity to the cortex than traditional DCS using a larger sourcedetector distance (ρ), while maintaining higher SNR. We have also extended the DCS noise model to the case of multiple independent speckle field measurements and validated our model experimentally. Our results indicate that our mDCS model works well even in extremely photonstarved conditions (down to four detected photons per pixel available to calculate g 2 ), thus enabling measurements at a large ρ.
Due to the recent investments in LiDAR technology for automotive and consumer electronics applications, high-performance large-pixel-count SPAD cameras with improved detection efficiency are rapidly becoming more available and less expensive. 30,37 While the first megapixel SPAD array was recently reported, 31 for mDCS to take full advantage of the extraordinary advances of SPAD camera technology, we envision the need for real-time data compression schemes implemented in the read-out FPGA of the SPAD camera or directly on chip. 30 Fig. 3 mDCS noise model. The fiber-SPAD distance was adjusted to z ¼ 8.1 mm to satisfy the condition d ¼ a. The SNR is the ratio between the mean and the standard deviation of g 2 measured over 160 integration periods. (a) SNR (τ ¼ 4 μs) of g 2 for increasing integration times T int and fixed N ph ¼ 1.5 × 10 5 cps, measured (solid circles) and calculated (solid lines) using different pixel ensemble sizes M ¼ 1, 32, and 1024. The SNR gain achieved by increasing the ensemble size to M ¼ 32 (see arrow from points A to B) can also be achieved by increasing T int by 32 times (see arrow from points A to C). (b) SNR (τ ¼ 4 μs) of g 2 for increasing detected photon count rates N ph and fixed T int ¼ 50 ms, measured (solid circles) and calculated (solid lines) using different pixel ensemble sizes M ¼ 1, 32, and 878. The SNR measured with only one pixel and ∼100 detected photons per T int (point D) is comparable to the SNR measured by averaging over 878 pixels and as low as ∼4 detected photons per pixel per T int (point E). Here, we used only up to M ¼ 878 pixels to avoid the DCR contribution from hot pixels (see Appendix D for details on DCR characterizations on the SPAD camera).
Finally, our mDCS technique can also be implemented in the time domain to enhance sensitivity to deep tissue. [38][39][40][41]

Appendix A: Experimental Setup
In our experiments, we used a high-coherence-length (∼9 m), 785-nm wavelength CW laser (Thorlabs, DBR785P) that was coupled to dynamic scattering media with an SMF (5-μm core diameter, 730-nm cutoff wavelength). The source power was adjusted by an attenuator to produce up to 33 mW at the output of the SMF. We collected the light diffusing out of the dynamic scattering media with an MMF (910-μm core diameter, 0.22 NA) coupled to a 32 × 32 pixel SPAD camera (Photon-Force, PF32). The SPAD camera consisted of 1024 SPADs arranged in a 32 × 32 array, with a pixel pitch of 50 × 50 μm and a pixel active area of 6.95 × 6.95 μm. Each pixel was surrounded by electronics necessary to bias and quench the SPAD, as well as to record photon detection events, with a detection efficiency of ∼8.2% at an excess bias of 1.7 V and wavelength of 785 nm. About ∼85% of the pixels had an average dark count rate (DCR) of ∼24 Hz. Each pixel was configured to record between 0 and 127 detection events per frame. We ran the measurements in the photon-counting mode using frame exposure times in the range of 4 μs < T bin < 10 μs, with the corresponding frame rate between 250 and 100 kfps. As shown in Fig. 1, our scattering media were a diffuser plate to calibrate the diameter of the speckles [ Fig. 1(a)] and a liquid phantom to measure the SNR gain [ Fig. 1(e)].
The diffuser plate consisted of a 1-in.-diameter ground glass diffuser (Thorlabs, DG10-120). The source fiber (SMF), the diffuser phantom, and the detector fiber (MMF) were arranged in a transmission geometry with a 10-mm gap between the phantom and the detector optode. The other end of the detector fiber was coupled to the SPAD camera with a translating lens tube that allowed us to adjust the MMF-SPAD distance in the range of z ¼ 3.5 to 200 mm. The image of the speckles was recorded using the SPAD camera at T bin ¼ 10 μs. The diameter of the speckles was determined using two-dimensional (2D) spatial autocorrelation of the speckle pattern as described in Appendix B. To measure the speckle turnover time trace, we mounted the diffuser plate on a motorized rotation stage and rotated the plate at an angular speed of 15 deg ∕s.
The liquid phantom was a mixture of milk and water with a volume ratio of 1:8 in a 19 × 19 × 21 cm 3 black Noryl plastic container at room temperature. The temperature of the room was 21 AE 0.5°C. We did not actively stabilize the temperature of the milk phantom. We brought the connectors of the source (SMF) and detector (MMF) fibers into contact with the surface of the mixture. The source-detector distance was ρ ¼ 11 mm. The other end of the MMF optode was coupled to the SPAD camera with an MMF-SPAD distance of z ¼ 8.1 mm to achieve d ¼ a. The frame exposure time was T bin ¼ 4 μs.

Appendix B: Calibrating the Diameter of the Speckles Using 2D Spatial Autocorrelation
In this experiment, we used a diffuser plate in transmission geometry with a 785-nm CW laser source, an SM source fiber optode (4.4-μm core diameter, 0.13 NA), and an MM detector fiber optode (400-μm core diameter, 0.50 NA). The size of the speckles could be adjusted by tuning the fiber-SPAD distance via the well-known relationship, 32 d ¼ λz∕D, where d is the average diameter of the speckles, λ is the laser wavelength, z is the fiber-SPAD distance, and D is the MMF detector fiber core diameter. Figure 4(a) shows three images of the speckle pattern recorded by the SPAD camera. As predicted, the average size of the speckles increased when we increased z. We quantified the average diameter of the speckles by calculating the 2D spatial autocorrelation function of each image, as shown in Fig. 4(b). Figure 4(c) shows a linecut from panel (b) at z ¼ 106 mm and a Gaussian curve fit that gives the estimated d as the 1∕e value of the Gaussian. We repeated the measurement and fitting procedure at increasing z for three different MMF detector fiber core diameters (D ¼ 200, 400, and 600 μm). As shown in Fig. 4(d), the measured diameters (solid circles) were in good agreement with the calculated diameters (solid lines) based on the d ¼ λz∕D equation. The use of the 2D spatial autocorrelation was limited to speckle diameters larger than the size of a pixel pitch (50 × 50 μm). We used the calibrated z from the above equation to estimate smaller speckle diameters. We applied this calibration to the milk phantom experiments, for which we used an MM detector fiber with 910-μm core diameter and adjusted z ¼ 8.1 mm to obtain d ¼ 6.95 μm, matching the diameter of the speckles with the length of the pixel active area.

Appendix C: Effects of Speckle Size per Pixel on SNR
Since the fiber core diameter determines the diameter of the speckles on the SPAD pixel, it is important to investigate how using different fiber core diameters, thus the diameter of the speckles, affects the ensemble SNR of g 2 . We used the liquid phantom setup to perform the SNR versus fiber core diameter study. Figure 5(a) shows that larger fiber core diameters resulted in higher SNR than that of the smaller ones for the same z. The SNR increased with decreasing z for all of the fibers we tested. For the largest core fiber that we tested (which had 910-μm core diameter), the SNR saturated at ∼250 for z < 10 mm. The maximum SNR was achieved when we tuned the diameter of the speckles to be comparable to the length of the pixel active area (6.95 × 6.95 μm). We examined this by plotting the Fig. 4 Measuring the diameter of the speckles using 2D spatial autocorrelation. (a) Snapshots of the speckle pattern at increasing fiber-SPAD distances z ¼ 59, 106, and 132 mm. These snapshots show that the average size of the speckles gets larger at increasing z. (b) 2D spatial autocorrelation images of the speckle pattern from (a). These images provide more quantitative information about the diameter of the speckles at respective z. (c) A linecut of the 2D spatial autocorrelation images from the z ¼ 106 mm in (b). We estimated the measured average diameter of the speckle (solid circles) through Gaussian fitting (solid curve). (d) Measured diameter of the speckles (solid circles) at varying detector MMF fiber core diameters. The straight lines are comparison based on the d ¼ λz∕D calculation. ensemble SNR of g 2 (4 μs) (from all fiber core diameters) against the number of speckles per pixel active area, as shown in Fig. 5(b). The number of speckles per pixel active area was calculated as the ratio between the length of the pixel active area and the diameter of the speckles. Figure 5(b) shows that all of the data points fall on a single curve. This suggests that, for a given count rate (N ph ) and integration period (T int ), the g 2 SNR in mDCS depends solely on the number of speckles per pixel active area. The SNR gradually increased until it reached one or more speckles per pixel active area [ Fig. 5(b)]. Here, S denotes the number of speckles per pixel active area. While β decreased with increasing S, STD also decreased with increasing S due to larger N ph from having more detected speckles, thus keeping the SNR the same at S ≥ 1 per pixel active area. We observed that β increased if we increased the diameter of the speckles, as predicted, and saturated beyond the Nyquist spatial sampling rate, which we reached when the diameter of the speckles was larger than twice the length of the pixel active area. 35,36 This is consistent with the observed onset of saturated SNR for the 910-μm fiber core diameter at z ≤ 10 mm [ Fig. 5(a)]. The results presented in this work were achieved using a 910-μm fiber core diameter at z ¼ 8.1 mm at which the diameter of the speckles matches the length of the pixel active area (d ¼ a).

Appendix D: Effects of Hot Pixels in SPAD Camera on SNR
Physical defects in SPAD pixels can effectively increase the DCR due to trapped carriers and after-pulsing. We refer to pixels that have high DCR (>100 Hz) as hot pixels and those with low DCR (<100 Hz) as cool pixels. We identified hot pixels in the kilopixel SPAD camera and investigated their effects on the intensity autocorrelation function g 2 and SNR of g 2 .
We characterized the DCR of the SPAD camera by recording the count rate during which the camera's aperture was covered with a blank. Figures 6(a)-6(c) show the distribution of DCR across 146 hot pixels (white) and 878 cool pixels (blue). The hot pixels make up <15% of the total pixels. Figure 6(d) shows that the mean DCR of hot pixels was 26.6 kcps and that of cool pixels was 24 cps, with a mean DCR of 3.8 kcps across all pixels.
We measured the ensemble average of g 2 − 1 (τ ¼ 4 μs), which we referred to as the coherence factor β [ Fig. 6(e)], and SNR of g 2 (τ ¼ 4 μs) [ Fig. 6(f)] from cool pixels (blue), all pixels (orange), and hot pixels (yellow) at decreasing count rate per pixel (N ph ). The integration time was T int ¼ 50 ms, and the ensemble average was performed over 160 integration periods. Figure 6(e) shows that the ensemble average of the cool pixels (blue circles) resulted in β values that remained constant over a wide range of N ph . However, the ensemble average of the hot pixels (yellow circles) resulted in β values that abruptly decreased at decreasing N ph down to N ph ∼ 2.6 × 10 4 cps. The ensemble average of all pixels (orange) resulted in β values that Fig. 5 Tuning the diameter of the speckle across the condition of one speckle per pixel. (a) SNR of g 2 increased with decreasing fiber-SPAD distance. Larger fiber core diameters resulted in higher SNR, particularly prominent for 910-μm fiber core diameter at fiber-SPAD distance of z ≤ 10 mm. (b) SNR of g 2 versus number of speckles per pixel active area (6.95 × 6.95 μm). The SNR increased upon reaching the condition of one speckle per pixel active area and achieved a maximum value between S ¼ 1 to 10 speckles per pixel active area.
were slightly lower than those of the cool pixels down to N ph ∼ 3.8 × 10 3 cps. The minimum N ph were different for the three different ensemble averages, which were given by the mean DCR values shown in Fig. 6(d). The corresponding SNR of g 2 values [ Fig. 6(f)] showed a clear decrease in magnitude when the count rates approached the minimum N ph . The SNR versus N ph plot shows the importance of DCR characterization of the SPAD camera, in particular for DCS measurements at very low N ph .
The decreasing β values with decreasing N ph were due to DCR contribution from each pixel, which is expressed as We calculated the predicted β values for each pixel using Eq. (6); they are plotted in Fig. 6(e) (solid lines) for corresponding ensemble averaging of cool pixels, hot pixels, and all pixels. The calculated β values are consistent with the measurements. In particular, the DCR contribution to the decreasing β values is apparent for the ensemble average of hot pixels.  4 μs), or coherence factor β, across selected pixels at varying N ph (T int ¼ 50 ms). Here, we used a liquid phantom with source-detector separation of 11 mm and a 785-nm CW laser. The standard deviations are shown by the widths of the colorshaded regions. Solid circles show the measured g 2 − 1, and solid lines show the predicted g 2 − 1 as calculated using Eq. (6). (f) Measured SNR of g 2 (τ ¼ 4 μs) across selected pixels at varying N ph (over 160 integration periods).
9 Appendix E: Fitting Functions of g 2 τ The functional form of the intensity autocorrelation function g 2 ðτÞ is determined by the dynamic and optical properties of the scattering media and the geometry surrounding the light source and the detector. In principle, we can predict the functional forms of g 2 either by directly estimating the resulting phase factor for all of the scattered electric fields for simple systems or by solving the correlation diffusion equation for more elaborate systems. The most commonly used dynamical scattering media in DCS experiments involve scatterers that are undergoing random flow motion (ballistic) or Brownian motion (diffusive). It has been shown 15,16,20,21 that, in the case of ballistic motion, the functional form can be simplified to g 2 ðτÞ ¼ 1 þ β expð−τ 2 ∕τ 2 c Þ at small time lags, τ ≪ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3μ a ∕ðμ 0 s k 2 0 αv 2 Þ p , where μ a and μ s 0 are the absorption and reduced scattering coefficients of the medium, k 0 is the wavenumber of the light, α is the percentage of light scattering events from moving scatterers, and v is the standard deviation of the speed of the scatterers. Similarly, in the case of diffusive motion, g 2 ðτÞ ¼ 1 þ β expð−τ∕τ c Þ at small time lags, τ ≪ μ a ∕ð2μ 0 s k 2 0 αD b Þ, where D b is the diffusion coefficient of the moving scatterers. In our work, the scattering media were a rotating diffusive plate [ Fig. 1(a)] and a liquid phantom [ Fig. 1(e)], and the details of the experimental setups were discussed in Appendix A. Figure 7 shows the normalized g 2 ðτÞ (solid circles) measured from the rotating diffuser plate [ Fig. 7(a)] and the liquid phantom [ Fig. 7(b)]. We also plotted the functional forms of g 2 ðτÞ that were used to fit the measured data: g 2 ðτÞ − 1 ¼ β expð−τ 2 ∕τ 2 c Þ (red curve) and g 2 ðτÞ − 1 ¼ β expð−τ∕τ c Þ (blue curve). As shown in Fig. 7, the intensity autocorrelation function of the rotating diffuser plate can be fitted well with the ballistic motion functional form expð−τ 2 ∕τ 2 c Þ and that of the liquid phantom with the diffusive motion functional form expð−τ∕τ c Þ over a wide range of time lags. The fitting results are consistent with the fact that the scatterers in the rotating diffusive plate behave like ballistic scatterers and those in the liquid phantom behave like diffusive scatterers. (a) Normalized g 2 ðτÞ − 1 (solid circles) measured from the rotating diffuser plate and two different functional forms of g 2 ðτÞ that were used to fit the measured data: c Þ (red curve) and g 2 ðτÞ − 1 ¼ β expð−τ∕τ c Þ (blue curve). The g 2 curve was normalized by the value of g 2 at τ ¼ 10 μs. The measured g 2 ðτÞ can be fitted well with g 2 ðτÞ − 1 ¼ β expð−τ 2 ∕τ 2 c Þ, where τ c ¼ 557 μs, consistent with the ballistic motion dynamics of the rotating diffuser plate. We measured β ¼ 0.7 with the rotating diffuser plate setup. The resulting β was higher than the value expected from the case of unpolarized light (β ¼ 0.5) likely due to a polarization bias in the detected light. The speckle diameter was tuned to be three times larger than the pixel size (d > 3a). (b) g 2 ðτÞ − 1 (solid circles) measured from the liquid phantom and two functional forms of g 2 ðτÞ (solid curves). The measured g 2 ðτÞ can be fitted well with g 2 ðτÞ − 1 ¼ β expð−τ∕τ c Þ, where τ c ¼ 197 μs, consistent with the diffusive motion dynamics of the liquid phantom. We measured β ¼ 0.2 with the liquid phantom setup. The resulting β was lower than the value expected from the case of unpolarized light (β ¼ 0.5) and results from having the speckle diameter matched with the pixel active area (d ¼ a). The maximum value for unpolarized light (β ¼ 0.5) can be obtained by increasing the speckle diameter beyond the Nyquist rate (d > 2a).
Edbert J. Sie is a research scientist at the Brain-Computer Interface (BCI), Facebook Reality Labs (FRL), developing high-sensitivity optical systems to measure neural activity. He received his PhD in physics from Massachusetts Institute of Technology (MIT) in 2017 and a postdoctoral fellowship from Stanford University, where he led research in ultrafast optics and imaging techniques, time-of-flight detection, quantum materials, and atomically thin semiconductors.
Hui Chen is an optical engineer in the BCI team at FRL Research at Facebook. He received his PhD in physics from the University of Connecticut in 2014. After that, he worked as an optical engineer in Calmar Laser and then as a development scientist at Corning Incorporated, before joining Facebook.
E-Fann Saung is an optomechanical engineer at FRL Research. He received his BS and MS degrees in mechanical engineering from Illinois Institute of Technology and San Diego State University, respectively, and is a licensed professional engineer in machine design and materials.
Ryan Catoen is a software engineer at FRL Research. Previously, he studied computer science at the University of Waterloo. His research interests focus on accelerating the development of next-generation noninvasive brain-computer interfaces.
Tobias Tiecke received his PhD in experimental quantum physics from the University of Amsterdam, the Netherlands, in 2009 and was a postdoctoral researcher in quantum optics at Harvard University. He currently works as a research scientist manager at Facebook, Inc.
Mark A. Chevillet is a research director at FRL leading research into a wearable silent speech interface featured by The Economist, IEEE Spectrum, and The Wall Street Journal. Previously, he directed research programs in BCI, artificial intelligence, and cognitive performance at Johns Hopkins University Applied Physics Laboratory, with joint appointments in neuroscience and cognitive science, the Science of Learning Institute, and the Kavli Neuroscience Discovery Institute.
Francesco Marsili is a research scientist manager at FRL developing diffuse optical techniques for the measurement of biosignals. Previously, he led research in single photon detectors, optical communication, quantum optics, superconductivity, and nanofabrication at the Jet Propulsion Laboratory, the National Institute of Standards and Technology, and MIT. He received his PhD in physics from the École Polytechnique Fédérale de Lausanne, Switzerland, in 2009.