Interferometric near-infrared spectroscopy (iNIRS) reveals that blood flow index depends on wavelength

Blood flow index (BFI) is an optically accessible parameter, with unit distance-squared-over-time, that is widely used as a proxy for tissue perfusion. BFI is defined as the dynamic scattering probability (i.e. the ratio of dynamic to overall reduced scattering coefficients) times an effective Brownian diffusion coefficient that describes red blood cell (RBC) motion. Here, using a wavelength division multiplexed, time-of-flight- (TOF) - resolved iNIRS system, we obtain TOF-resolved field autocorrelations at 773 nm and 855 nm via the same source and collector. We measure the human forearm, comprising biological tissues with mixed static and dynamic scattering, as well as a purely dynamic scattering phantom. Our primary finding is that forearm BFI increases from 773 nm to 855 nm, though the magnitude of this increase varies across subjects (23% ± 19% for N = 3). However, BFI is wavelength-independent in the purely dynamic scattering phantom. From these data, we infer that the wavelength-dependence of BFI arises from the wavelength-dependence of the dynamic scattering probability. This inference is further supported by RBC scattering literature. Our secondary finding is that the higher-order cumulant terms of the mean squared displacement (MSD) of RBCs are significant, but decrease with wavelength. Thus, laser speckle and related modalities should exercise caution when interpreting field autocorrelations.


Introduction
Diffuse Correlation Spectroscopy (DCS) [1][2][3] is an established optical technique for non-invasive estimation of biological tissue blood flow index (BFI) [1].BFI is defined as the product of α, the ratio of the dynamic reduced scattering coefficient (µ ′ s,dyn ) to the total tissue reduced scattering coefficient (µ ′ s ), and D B , an effective red blood cell (RBC) Brownian diffusion coefficient.BFI is a proxy for nutritive blood flow, and correlates with conventional perfusion metrics [4][5][6][7][8][9][10][11][12][13].BFI is typically measured in the wavelength range from 767-855 nm [14][15][16][17][18][19][20][21][22][23][24], though there is a growing push in the DCS field to measure around 1064 nm, a wavelength range which confers several benefits [25][26][27][28].To the extent that BFI is a valid proxy for tissue perfusion, one expects that it should not depend on the wavelength of light used to probe tissue.Indeed, most studies in the field have treated BFI as wavelength-independent, and at least one study explicitly states this assumption [23].However, a wavelength-dependent α could introduce inconsistencies in BFI measurements [23].Directly testing wavelength dependence of α experimentally is challenging as there are numerous confounding factors.For instance, in multi-wavelength DCS, different sources or collectors for the two wavelengths may introduce variability in the tissue volume probed.Also, DCS autocorrelations depend not only on BFI, but also the time-of-flight distribution (DTOF), which is also wavelength-dependent. Thus, differences in BFI can be masked by other experimental factors in multi-wavelength DCS [23].
Interferometric near-infrared spectroscopy (iNIRS) was recently introduced to provide autocorrelations from biological tissue with TOF resolution [21,[29][30][31][32][33].Here, to address the above questions, we further advance the iNIRS approach with wavelength-division multiplexing [33].The two system wavelengths, 773 and 855 nm, span the most frequently-used DCS wavelengths.The dual-wavelength iNIRS setup was meticulously designed to encompass key pivotal features: shared source and collector for both wavelengths, simultaneous autocorrelation measurements at both wavelengths without time multiplexing, and TOF-resolved autocorrelations with matched TOF resolution at both wavelengths.These elements contribute to a robust framework for recovering BFI experimentally, avoiding many issues that may confound BFI comparisons in DCS.
Our primary finding is a wavelength-dependence of α in biological tissue.This finding, while unexpected for the DCS community, is in fact already implied by the RBC light scattering literature [34][35][36].This result should inform interpretation of DCS in its numerous applications.Our secondary finding; that that higher-order cumulant terms are significant but decrease with wavelength, should inform interpretation of laser speckle and related modalities.

Theory
iNIRS recovers a TOF-resolved field autocorrelation, g 1 (τ s ,τ d ), where τ s is TOF and τ d is time lag [21,29,31,32].Allowing for a finite TOF resolution, the iNIRS experimental data closely relates to DWS theory [37], which is based on simple fundamental principles.According to DWS, the field autocorrelation, g DWS 1 (τ s ,τ d ), is given by where k is the medium wavenumber, and v c is the medium speed of light.Specifically, incorporating a dynamic scattering probability, α= µ ′ s,dyn /µ ′ s , DWS yields Though well-known and widely used, Eq. ( 2) includes a number of assumptions: approximation of a distribution of the number dynamic scattering events by a single mean value, neglect of higher order MSD cumulants, and neglect of random flow.Therefore, Eq. ( 2) may not describe g 1 (τ s ,τ d ) well at early τ s and late τ d [38], as is required for this work.On the other hand, Bonner and Nossal's seminal dynamic light scattering theory [39], if extended to include TOF [38], incorporates a Poisson-distributed number of dynamic scattering events and implicitly includes higher MSD cumulants by integrating over the dynamic phase function [38].Their theory can be further modified to allow both Brownian motion and random flow, a feature missing from Eq. ( 2), to drive the mean squared displacement (MSD) of dynamic scatterers [38].The error arising from exclusion of the second cumulant of MSD increases with higher dynamic scattering anisotropy, g dyn , and fewer dynamic scattering events [38].Since RBCs have g dyn ∼0.98 and we need to analyze g 1 (τ s ,τ d ) at early τ s , we use an expression from [38] which extends Bonner and Nossal's theory and incorporates fewer assumptions than Eq.(2).
In Eq. ( 3), m(τ s ) is the average number of dynamic scattering events at τ s , q 2 (θ) is the momentum transfer as a function of dynamic scattering angle θ, <∆r 2 (τ d )⟩ is the MSD of the RBCs, and angle brackets ⟨•⟩ θ denotes an angular weighted average over θ.In contrast to Eq. ( 2), the more general Eq. ( 3) can describe g 1 (τ s ,τ d ) even with few dynamic scattering events and high g dyn [38].Hereafter, we will drop the subscript θ in the angular average for sake of readability.
In simulation, Eq. ( 3) was shown to better describe TOF-resolved autocorrelations, particularly in the slowly-decaying "tails" where DWS theory of Eq. ( 2) fails [38].Taylor series expansion of the exponential term in Eq. (3) up to three MSD cumulants yields where and Note that the fourth and the higher order terms are contained in O[⟨∆r 2 (τ d )⟩ 4 ].Incorporating hybrid motion [40,41], , where v 2 is the second moment of the Gaussian velocity distribution and arranging different powers of τ d (see full details in Appendix A), we have O(τ d 4 ) contains fourth and higher orders of τ d .It is also worthwhile noting that MSDcumulants [Eq.(4)] do not necessarily correspond with τ d -cumulants [Eq.(8)] when hybrid motion is assumed.Importantly, in Eq. ( 8), different order terms for Brownian motion and random flow contribute to both the second-and third-order τ d terms, with opposite signs.Moreover, setting v to 0 and neglecting second-and higher orders of τ d in Eq. ( 8) yields conventional DWS [Eq.( 2)], given m(τ s )=µ s,dyn v c τ s , where µ s,dyn is the dynamic scattering coefficient.

Experimental setup
Dual-wavelength iNIRS [Fig.1] comprises two distributed feedback (DFB) lasers with center wavelengths of 773 nm (Eagleyard EYP-DFB-0773-00075-1500-TOC03-0002) and 855 nm (Eagleyard EYP-DFB-0855-00150-1500-TOC03-0000).We achieved synchronous wavelength tuning by sinusoidally modulating the currents via two independent current controllers (SRS LDC-501), with modulation input voltages from a common function generator (TTi TGF4042).Wavelength tuning with 2.5 V pp from the function generator provided a TOF resolution of 75 ps FWHM for 773 nm and 67 ps FWHM for 855 nm.To equalize the TOF resolution, we introduced a series resistor to create a voltage divider at the input of the current controller driving the 855 nm laser.This enabled us to achieve a TOF resolution of 58 ps for both wavelengths with 3.3 V pp from the function generator.Both configurations were used in this study.The sinusoidal drive frequency was 100 kHz, yielding 5 µs field autocorrelation resolution [30].
The lasers emitted light with elliptical spatial beam profiles.To couple this light into a single-mode (SM) fiber, for each wavelength, we collimated each beam using a convex aspheric lens and adjusted the beam shape to be circular with an anamorphic prism pair.Subsequently, for each wavelength, the collimated beam passed through a free space isolator before being focused into a SM fiber-based 90:10 splitter, with 90% of the light directed to the sample arm and the remaining 10% to the reference arm.
The two 90% ports were connected to a custom wavelength-division multiplexer (WDM) characterized by 98% transmission and ∼20 dB isolation.The WDM combined both wavelengths before illuminating the sample.The remitted light was collected through a similar WDM, this time to separate the two wavelengths.This approach ensured that both wavelengths illuminated and collected light from the same spot on the tissue, in accordance with our experimental requirements.Slight differences in the probed tissue regions due to variations in optical properties at the two wavelengths were unavoidable, however.
Subsequently, the light from each arm of the splitter WDM was recombined with the appropriate reference arm light of the same wavelength via a 50:50 combiner.Finally, both arms of both 50:50 combiners were connected to New Focus 1807 dual-balanced detectors (DBD) which provided the interference signal.At this stage, any cross-talk in the WDMs was further suppressed by interferometric gating whereby only the sample light of the correct wavelength was amplified with reference light.The DBD output voltages were sampled by an Alazar Technologies 9416 fast data acquisition card (DAQ) at 100 MHz.
The DAQ card derived its clock from the reference clock output of the function generator.All measurements presented in this paper were continuously acquired following a trigger signal, generated from the function generator, and sent to the DAQ to commence data acquisition.The trigger signal to the DAQ is a square wave, synchronous to the sinusoidal signal used to tune the DFB lasers.This mode of triggered and synchronized acquisition, coupled with low timing jitter of the function generator, allowed us to omit recording of the trigger signal, thereby simplifying our data processing and conserving hard drive space.

Experimental measurements
With the dual wavelength iNIRS system, measurements were performed on an Intralipid mixture, prepared by mixing 21 mL of 20% Intralipid in 500 mL water.No absorbing agent was added.Measurements were performed at room temperature.The source-detector separation was 11 mm for the Intralipid measurements.
Measurements were also performed on the forearms of three healthy males at rest.The procedures and protocols were approved by the NYU Langone and UC Davis Institutional Review Boards.The source-detector separation was 6 mm for the human measurements.The probe was placed on the dominant forearm right above the brachioradialis muscle.The combined skin and adipose tissue thickness over the brachioradialis muscle was measured by a skinfold caliper for each subject.Measurements at both wavelengths were acquired simultaneously for 40 seconds.Though major findings were ultimately verified across subjects (age 31.67 yrs.± 4.04 yrs.for N = 3), we mainly present measurements and analysis from one subject (age: 28 yrs).

Data acquisition and processing
iNIRS provides time-of-flight (TOF) resolution through wavelength tuning of the source [21,[29][30][31][32][33].To ensure accurate tracking of instantaneous wavelength variations, a calibration measurement, hereafter referred to as the instrument response function (IRF) measurement, was conducted.During the IRF measurement, the sample in Fig. 1 is replaced with a single-mode (SM) patch cord of known length.This process generates a fringe pattern at the output of the dual-balanced detectors (DBD) for each complete wavelength sweep.
We initiate the IRF processing by identifying stationary points, which correspond to time instances where the instantaneous wavelength reaches either the maximum or minimum of the sweep.The phase information derived from the Hilbert transform of the recorded fringes serves as a proxy for the optical frequency.A stationary point occurs where the time derivative of the phase becomes zero.It is important to note that for accurate wavelength tracking using phase information, the time derivative of the phase must change sign at every alternate sweep.We designate a sweep as a "forward sweep" when the optical frequency increases from minimum to maximum and a "backward sweep" when it decreases from maximum to minimum.Since our DFB lasers exhibited non-identical forward and backward sweeps, separate phase-versus-time calibration curves were generated for each sweep type by averaging the phase over many sweeps.The fringe envelope for each sweep type was also determined from the IRF measurement, and later used for Gaussian shaping of the measurement data.
With the envelope and phase-versus-time calibration curves for the sweep types in hand, we are prepared for measurements using the setup shown in Fig. 1.As previously mentioned, triggered acquisition ensures that the first stationary point always occurs a fixed time after acquisition starts, and always for the same sweep transition (e.g.forward-to-backward sweep).A low phase noise function generator ensures subsequent stationary points occur after each sweep cycle.Once stationary points are identified, the time axis of each sweep is transformed into phase (or by proxy, optical frequency) using the corresponding phase-versus-time curve from the IRF measurement.
Thereafter, we combined the measured data from the two sweep types, followed by resampling, Gaussian shaping and finally, Discrete Fourier Transform, to obtain the complex mutual coherence function Γ rs (τ s ,t d ) as described in Ref. [30].Here t d refers to delay time.Finally, the unnormalized iNIRS field autocorrelation function, G 1 (τ s , τ d ), is obtained from Eq. ( 9) [29,30]: Background measurements are acquired in the same way as sample interference measurements, except while blocking light transmission through the sample arm.Data processing steps for the background measurement mirror those of sample interference measurements.Subtracting the background autocorrelation G BG 1 (τ s , τ d = 0) from the measurement G 1 (τ s , τ d = 0) yields the desired temporal point spread function (TPSF), used to recover optical properties of the probed region.

Data analysis
In order to compare the scaling approaches described in the previous section, we need to quantify the degree of agreement between the predicted g λ 2 1 (τ s ,τ d ), obtained from scaling the the τ d -axis of g λ 1 1 (τ s ,τ d ), and the measured g λ 2 1 (τ s ,τ d ), obtained from our dual wavelength iNIRS setup.To facilitate this comparison, we utilize a Bland-Altman plot [42] to visually evaluate agreement between predicted and measured values.This method involves plotting the mean of the predicted and measured g λ 2 1 (τ s ,τ d ) values on the X-axis and their difference on the Y-axis for all τ d at a fixed τ s .In this analysis, a smaller difference in the Bland-Altman plot signifies a better agreement between the predicted and measured values.Here, the Bland-Altman plots offer insights into the degree of agreement for every τ d .Note that the difference at τ d = 0 will invariably be zero due to normalization.
The Bland-Altman plots offer detailed analysis at individual τ s .To comprehensively assess agreement across all TOFs (τ s ), we employed the coefficient of determination (R 2 ) [43].R 2 quantifies the goodness of agreement between the predicted and measured g λ 2 1 (τ s ,τ d ).We used 1-R 2 to compare scaling approaches, where a lower value signifies better agreement between predicted and measured quantities.

Results
We begin by fitting G 1 (τ s ,τ d ) from Intralipid (g dyn ∼0.6) to the simple DWS exponential model [unnormalized form of Eq. ( 2)], routinely used in the DCS field.For Intralipid data, we fit the entire autocorrelation curve i.e., g 1 (τ s ,τ d ) ranging from 1 to 0. TOF-resolved autocorrelations decay exponentially for both wavelengths [Fig.3, Fig. 4], consistent with DWS theory.Here the absolute medium "EARLY", "MID" and "LATE" τ s are 60, 330 and 500 ps, respectively.This definition of time gates is consistent with previous DCS literature [28,44].Unlike the Intralipid phantom, G 1 (τ s ,τ d ) from in-vivo measurements (g dyn ∼0.98) do not decay exponentially.Thus, G 1 (τ s ,τ d ) from both wavelengths were fitted to the unnormalized form of Eq. ( 13) [Fig.5, Fig. 6].For in-vivo data, we fit g 1 (τ s ,τ d ) ranging from 1 to 0.1.Since a single exponential model does not fit in-vivo G 1 (τ s ,τ d ), we will not use this model for fitting in-vivo measurements further.Instead we will use the autocorrelation functions directly to test both proposed scaling approaches (Section 2.5) on the Intralipid and in-vivo data.
For the in-vivo forearm measurements, constant µ ′ s,dyn -scaling performed best in predicting g 773 1 (τ s , τ d ) from g 855 1 (τ s , τ d ) [Figs.[10][11][12].The τ d -axis of 855 nm was scaled using Eq. ( 12) with recovered µ ′ s of 1.23 and 1.44 mm −1 for 855 and 773 nm respectively.Finally, to comprehensively compare the scaling techniques, we calculated 1-R 2 (where R 2 is the coefficient of determination) across different τ s to show how well they predicted g 773 1 (τ s , τ d ) from g 855 1 (τ s , τ d ).For the Intralipid mixture, we find that the constant α-scaling approach is better than constant µ ′ s,dyn -scaling for every τ s [Fig.13(a)].In the human forearm in-vivo [Fig.13(b)], we find that constant µ ′ s,dyn -scaling is better.The preceding analysis scaled the τ d -axis of the autocorrelations at one wavelength to predict the other wavelength.However, theory [Eqs.( 14)-( 17)] predicts a more complex relationship where higher-order cumulants scale more strongly with wavelength.So, we proceeded to consider scaling of individual terms in the τ d -cumulant expansion.As the first τ d -cumulant is sufficient for Intralipid, the results are very straightforward.Namely, the constant α-scaling approach applied to Eq. ( 11) successfully predicted ξ(τ s ) at 773 nm from ξ(τ s ) at 855 nm [Fig.14].These results on the exponential fit parameter are unsurprising given the results already presented, which showed that an exponential describes Intralipid G 1 at all TOFs, and that the constant α-scaling approach works best for Intralipid.
Next, we proceeded to consider scaling of terms in the τ d -cumulant expansion for in vivo measurements.The in-vivo G 1 (τ s ,τ d ) were fitted to Eq. ( 13) for the coefficients ξ 1 (τ s ), ξ 2 (τ s ) and ξ 3 (τ s ) [Fig.15(a), Fig. 16(a) and Fig. 17(a), respectively], for both wavelengths.The prediction of these three coefficients at 773 nm from their counterparts at 855 nm was accomplished through scaling based on Eqs. ( 15), ( 16) and ( 17     Note that 1-R 2 is shown, where R 2 is the coefficient of determination.The constant α-scaling approach performs best for the purely dynamic Intralipid phantom whereas the constant µ ′ s,dyn -scaling approach performs best in-vivo, with mixed static and dynamic scattering.At late τ s , differences between scaling approaches are less apparent due to low SNR.13) fitting of G 1 (τ s ,τ d ) from human forearm for both the wavelengths.(b) Same as (a) but ξ 1 (τ s ) at 855 nm was scaled using Eq. ( 15).13) fitting of G 1 (τ s ,τ d ) from human forearm for both the wavelengths.(b) Same as (a) but ξ 2 (τ s ) of 855 nm wavelength was scaled using Eq. ( 16).13) fitting of G 1 (τ s ,τ d ) from human forearm for both the wavelengths.(b) Same as (a) but ξ 3 (τ s ) at 855 nm was scaled using Eq. ( 17).
Next, we investigated the ability to predict 773 nm autocorrelations from the τ d -cumulant expansion order at 855 nm.Constant µ ′ s,dyn -scaling [Eqs.( 15)-( 17), also referred here as theory-based scaling] and a simpler approach where all the constants were scaled by the ratio of the wavelengths squared, here (855/773) 2 , referred here as λ 2 -scaling, were investigated.For both approaches, collectively referred to as higher cumulant-scaling, scaled coefficients were plugged into Eq.( 13) to generate predictions.A direct Eq.( 13) fit of g 773 1 (τ s , τ d ) serves as a benchmark.We found that theory-based scaling of the three τ d -cumulant coefficients at 855 nm predicted g 773 1 (τ s , τ d ) best [Fig.18, Fig. 19, and Fig. 20].Finally, we summarized the error (1-R 2 ) of the g 773 1 (τ s , τ d ) prediction from g 855 1 (τ s , τ d ) using the higher cumulant-scaling methods across τ s [Fig.21].The 1-R 2 values obtained from Eq. ( 13) fit of g 773 1 (τ s , τ d ) serves as a benchmark.We see that the approach where the τ d n coefficient is scaled by λ 2n outperforms the approach where all coefficients are scaled by λ 2 .This suggests that higher-order τ d terms scale more strongly with wavelength than the first-order (exponential decay) term.13); (b) Predicting g 773 1 at early τ s by scaling the recovered 855 nm ξ 1 (τ s ), ξ 2 (τ s ) and ξ 3 (τ s ) based on theory [Eqs.(15), ( 16) and ( 17)]; (c) Same as (b) but scaling each constant with (855/773) 2 ; and (d) Bland-Altman plot showing difference between fitting and measurement of (a), between theory-based prediction and measurement of (b), and between λ 2 scaling-based prediction and measurement of (c).(EARLY τ s : 60 ps, MID τ s : 330 ps and LATE τ s : 500 ps) Fig. 19.Same as Fig. 18 but for mid τ s of 330 ps.Comparison of the prediction error (1-R 2 ) of higher cumulant-scaling (based on theory, and λ 2 -scaling).The error of a direct g 773 1 (τ s , τ d ) fit to Eq. ( 13) serves as a benchmark for comparison.At late τ s , differences between scaling approaches are less apparent due to low SNR.

Discussion
This study employs a new dual-wavelength iNIRS setup to compare TOF-resolved autocorrelations across wavelengths for the first time.Our investigation uncovers a previously-overlooked wavelength-dependence of α, the probability of dynamic scattering.The immediate consequence of this primary finding is that BFI, the optical index of blood flow, is wavelength-dependent. Secondary findings highlight the importance of incorporating higher MSD cumulants in modeling the later time lags of field autocorrelations ("tails"), particularly at early and even moderate TOFs.

Wavelength dependence of α
We applied two scaling methods to estimate 773 nm autocorrelations from 855 nm autocorrelations: the constant α-scaling approach and the constant µ ′ s,dyn -scaling approach.
We observed that the constant α-scaling approach accurately predicted 773 nm autocorrelations from 855 nm autocorrelations in an Intralipid mixture [Fig.7, Fig. 8, Fig. 9, and Fig. 13(a)], consistent with α=1.As expected, the constant µ ′ s,dyn approach performed worse in this purely dynamic scattering phantom, as Intralipid reduced scattering does vary with wavelength [45].Yet in-vivo, the constant µ ′ s,dyn approach better predicted 773 nm autocorrelations from 855 nm autocorrelations [Fig.10, Fig. 11, Fig. 12, and Fig. 13(b)].Though surprising, this result in fact aligns with RBC measurements in the literature [34][35][36], where reduced scattering is similar at 773 nm and 855 nm.To further underscore this point, α was derived [Fig.22] from prior RBC scattering literature [34], under various assumptions for wavelength-dependent tissue scattering [46].Though derived purely from the literature, Fig. 22 recapitulates the major experimental result of our work: the wavelength-dependence of α.This result implies that BFI=αD B is also wavelength-dependent (as clearly, D B is wavelength-independent).BFI was derived from the autocorrelation decay rate at τ s 400 ps [31].Consistent with the discussion, in-vivo forearm BFI ratios between wavelengths differed from 1, and also varied between subjects [Table 1].This variation may occur because different skin types yield different wavelength dependencies of forearm reduced scattering [47,48]; yet dynamic (RBC) scattering should not depend on skin type.The mean and standard deviation of BFI for the in-vivo measurements were 1.86 × 10 −7 (±3.14 × 10 −8 ) mm 2 /sec for 855 nm and 1.57 × 10 −7 (±4.72 × 10 −8 ) mm 2 /sec for 773 nm, across three subjects who participated in the study.
A wavelength-dependent α implies that static and dynamic scattering need not co-vary with wavelength.Fortuitously, the constant µ ′ s,dyn approach worked well for the two chosen wavelengths in our in-vivo forearm experiments.Yet, from a physical standpoint, as refractive index variations leading to scattering are different for static and dynamic tissue, we can envision situations (such as different tissues or wavelengths) where neither α nor µ ′ s,dyn are invariant with wavelength.To enable meaningful comparisons, one solution is to evaluate µ ′ s BFI=µ ′ s,dyn D B , as a substitute for BFI alone.Such an approach mitigates concerns arising from differential behavior of µ ′ s,dyn and µ ′ s with wavelengths, but is still subject to variability in µ ′ s,dyn with wavelength, if any.

Better autocorrelation modelling by including higher MSD cumulants
Building on theoretical predictions [38], in this work we experimentally demonstrated the importance of considering higher-order MSD cumulants in autocorrelation modeling.We observed that the common single exponential model of DWS works well for Intralipid (g dyn ∼0.6 [49][50][51]), but inadequately captures the slowly decaying TOF-resolved field autocorrelation tails in vivo (g dyn ∼0.98 for RBCs [52][53][54]).This limitation is particularly pronounced at early TOF and is visible even at moderate TOFs of 330 ps [Fig.5 and Fig. 6].Our findings underscore the need for considering higher MSD cumulants to describe autocorrelation tails in laser speckle, when determining camera exposure times, and DCS, especially with short source-collector separations, when selecting fitting ranges for autocorrelation functions.Interestingly, higher-order MSD cumulants are almost never considered in the field.Admittedly, their impact would be lessened in g 2 = 1+|g 1 | 2 (intensity autocorrelation) as compared to g 1 (field autocorrelation).Addressing these considerations is crucial for achieving more accurate and comprehensive interpretations of tissue dynamics.
The multiple wavelengths in this study shed further light on the nature of deviations from the single exponential DWS model [Eq.( 2)].As anticipated by theory [Eq.( 8)], we found more prominent autocorrelation tails for the shorter wavelengths (773 nm) versus longer wavelengths (855 nm).Interestingly, we observed a fortuitous alignment of the in-vivo autocorrelation tails (late τ d ) between the two wavelengths in Fig. 10(a), Fig. 11(a) and Fig. 12(a).This alignment appears to result from two cancelling errors: first, the faster initial decay slope at the shorter wavelength and second, the more prominent tail at the shorter wavelength.
Last, to further test the importance of higher-order cumulants, we scaled the 855 nm τ dcumulant coefficients to predict autocorrelation functions at 773 nm (Fig. 18, Fig. 19, Fig. 20 and Fig. 21).Results reinforce the theoretical models [e.g.Eqs. ( 4), (8), and (14)] that predict stronger scaling of higher cumulants with wavelength.It is noteworthy that the scale factors were derived from Eq. ( 14) while neglecting random flow; however, random flow may still impact the tails of the autocorrelation function.Likewise, while the signs of the τ d -cumulant coefficients [ξ 1 (τ s ), ξ 2 (τ s ), and ξ 3 (τ s )] align with expectations from a pure Brownian diffusion model [Fig.15, Fig. 16 and Fig. 17], we cannot conclusively exclude the impact of random on the autocorrelations [38].We also cannot exclude impact of multiple layers (skin, fat, muscle) with different dynamics which may also affect autocorrelation tails [31].
Finally, it is important to note that the results on the effects of α are based mainly on the first cumulant (DWS theory) which describes the initial decay at early τ d of the autocorrelation function.This is distinguished from the higher-order cumulants which affect the autocorrelation "tails" at late τ d .Hence the two results are discussed separately and independently.

Conclusions
Our experimental findings substantiate an increase in the dynamic scattering probability, α, by 23 ± 19% from 773 nm to 855 nm in the forearms of 3 human subjects.This implies that BFI is wavelength-dependent, and depending on scattering, may be subject dependent as well.Therefore, we advocate for a meticulous approach to interpreting DCS-derived BFI that accounts for a possibly wavelength-dependent α, among numerous other factors that can confound comparisons between DCS systems.Additionally, our experiments fortify prior theoretical predictions underscoring the importance of higher MSD cumulants in accurate autocorrelation modeling at early to intermediate TOFs.

Appendix A
Here we show how we arrive at Eq. ( 8) from Eq. ( 4).We begin by taking natural logarithm on both sides of Eq. ( 4) and multiply them with -1/ m(τ s ) to obtain:

Fig. 1 .
Fig. 1.Dual wavelength iNIRS setup.All fibers are single mode with angle polished connectors.I/C controller: current controller; DFB: distributed feedback laser; L,L': lenses; APP: anamorphic prism pair; WDM: wavelength division multiplexer; DBD: dual balanced detector.Bold lines denote optical paths, while dotted lines denote electrical paths.Path 1 (red) is for 855 nm, while path 2 (blue) is for 773 nm.

Fig. 2 .
Fig. 2. (a) Light paths in an emulsion such as Intralipid include only dynamic scattering events.(b) Light paths in vascularized biological tissues include both static and dynamic scattering events.Therefore, for an emulsion, µ ′ s,dyn =µ ′ s and α=1, whereas in biological tissue with blood flow, µ ′ s,dyn ≠µ ′ s and α≠1, in general.

Fig. 4 .
Fig. 4. Same as above but for 773 nm.The IRF and corresponding TPSF are shown in the inset of Fig. 4(a).

Fig. 6 .
Fig. 6.Same as Fig. 5 but for 773 nm.The IRF and corresponding TPSF are shown in the inset of Fig. 6(a).

Fig. 13 .
Fig. 13.Ability to predict g 773 1 (τ s , τ d ) from g 855 1 (τ s , τ d ) by scaling the τ d -axis for (a) Intralipid phantom and (b) human forearm.Note that 1-R 2 is shown, where R 2 is the coefficient of determination.The constant α-scaling approach performs best for the purely dynamic Intralipid phantom whereas the constant µ ′ s,dyn -scaling approach performs best in-vivo, with mixed static and dynamic scattering.At late τ s , differences between scaling approaches are less apparent due to low SNR.

Fig. 20 .
Fig. 20.Same as above but for late τ s of 500 ps.Differences in scaling approaches are less apparent due to low SNR.

Fig. 21 .
Fig. 21.Comparison of the prediction error (1-R 2 ) of higher cumulant-scaling (based on theory, and λ 2 -scaling).The error of a direct g 773 1 (τ s , τ d ) fit to Eq. (13) serves as a benchmark for comparison.At late τ s , differences between scaling approaches are less apparent due to low SNR.