Beyond diffuse correlations: deciphering random flow in time-of-flight resolved light dynamics

: Diﬀusing wave spectroscopy (DWS) and diﬀuse correlation spectroscopy (DCS) can assess blood ﬂow index (BFI) of biological tissue with multiply scattered light. Though the main biological function of red blood cells (RBCs) is advection, in DWS/DCS, RBCs are assumed to undergo Brownian motion. To explain this discrepancy, we critically examine the cumulant approximation, a major assumption in DWS/DCS. We present a precise criterion for validity of the cumulant approximation, and in realistic tissue models, identify conditions that invalidate it. We show that, in physiologically relevant scenarios, the ﬁrst cumulant term for random ﬂow and second cumulant term for Brownian motion alone can cancel each other. In such circumstances, assuming pure Brownian motion of RBCs and the ﬁrst cumulant approximation, a routine practice in DWS/DCS of BFI, can yield good agreement with data, but only because errors due to two incorrect assumptions cancel out. We conclude that correctly assessing random ﬂow from scattered light dynamics requires going beyond the cumulant approximation and propose a more accurate model to do so.


Introduction
Diffusing Wave Spectroscopy (DWS) was formulated in the late 1980s to probe dynamics of scatterers in soft matter based on intensity fluctuations of multiply scattered light. DWS models the field autocorrelation as an integral over photon paths, assuming uncorrelated motion at different scattering events [1][2][3], in the light diffusion approximation. In biological tissues, DWS and its differential formulation of diffuse correlation spectroscopy (DCS) can measure red blood cell (RBC) dynamics, which correlate empirically with blood flow [4][5][6]. The assumptions underlying DWS theory are powerful and widely used, even in distant fields such as laser speckle flowmetry [7,8]. However, given recent experimental advances in measuring time-of-flight-(TOF-) resolved light dynamics [9][10][11] in detail for the first time, some assumptions in DWS/DCS of blood flow merit further scrutiny. Here we focus on two important and inter-related assumptions; 1) the nature of RBC motion, and 2) the first cumulant approximation of the average over dynamic scattering angles.
Particle motion can be ordered, as in advective, randomly oriented, flow, or disordered, as in Brownian motion, or a combination of both [7]. DWS/DCS measurements of particle motion in biological tissues such as the rat brain [12,13], piglet brain [14], human brain [6,[15][16][17][18], breast [19,20], head and neck [21], and skeletal muscle [22] have shown that experimental data is best fitted by ignoring advection of RBCs, and only considering Brownian motion. Thus, the blood flow index (BFI) in DWS/DCS has units of a Brownian diffusion coefficient.
The starting point of dynamic light scattering theories [23,24] is the normalized field autocorrelation of m th order paths: exp i m j=1 q j · ∆r j , where denotes an average over particle displacement vectors ∆r j and dynamic momentum transfer vectors q j , which are assumed to be independent, and i = √ −1. In the multiple scattering regime, assuming uncorrelated dynamic scattering events while neglecting both random flow and higher order terms in the cumulant expansion [23], DWS expresses the normalized TOF-resolved field autocorrelation function, g 1 , as: In Eq. (1),m is the average number of dynamic scattering events at time-of-flight (TOF), τ s , q 2 (θ) θ is the angle-averaged momentum transfer for dynamic scattering, D B is the Brownian diffusion coefficient of dynamic particles, and τ d is autocorrelation time lag.
It is notable that most studies of blood flow in biological tissues utilize Eq. (1), although this expression was first derived for uniform colloidal suspensions, foams, and gels with unordered (Brownian) motion. A major biological function of RBCs is the transport of gases in the circulation; thus while hydrodynamic diffusion of RBCs is well-documented [25,26], RBCs should also possess characteristics of ballistic motion [24]. Therefore, it seems incongruous that blood flow, in this conventional sense, is typically assumed not be observable through light scattering techniques.
In this work, we present a plausible explanation to help resolve these discrepancies. We simulate field autocorrelations with a custom extended Monte Carlo method that explicitly allows for a mixture of static and dynamic particles with different scattering phase functions. We attempt to accurately model both the high scattering anisotropy and low volume fraction of RBCs in tissue. We assess the results to critically examine the first cumulant approximation, or first cumulant expansion, which is common in DWS. We show that under certain circumstances, errors due to the cumulant approximation, on one hand, and neglecting random flow, on the other, can cancel, enabling the conventional theory in Eq. (1) to describe autocorrelations well, even while ignoring advection. We then propose a more accurate approach that avoids these assumptions and might enable direct quantification of RBC advection in dynamic light scattering measurements.

Field autocorrelation models
In this section, we develop candidate models for TOF-resolved field autocorrelations, highlighting key assumptions in each.

Nature of scatterer motion
In DWS/DCS, the mean-squared displacement at time lag τ d , ∆r 2 (τ d ) , of RBCs is modeled as either random flow, given by where v 2 is the second moment of the Gaussian velocity distribution, or as Brownian motion, given by ∆r 2 (τ d ) = 6D B τ d , where D B is the Brownian diffusion coefficient. In this work, we allow for RBCs to exhibit a "hybrid" of both random flow and Brownian motion, given by where the two types of motion are assumed to be independent.

Adapting Bonner and Nossal's theory
To assess the validity of the cumulant approximation (or first cumulant expansion), we start from Bonner and Nossal's theory [24], which implicitly includes all cumulant terms by directly integrating over the dynamic scattering phase function. This theory envisions tissue as a mixture of static (tissue matrix) and dynamic scatterers (i.e. RBCs and other scattering blood components). We extend this theory as follows: 1) the lower limit of number of dynamic scattering events, m, is set to zero to include purely static scattering paths, 2) the Henyey-Greenstein (and later, the more general Gegenbauer kernel) scattering phase function for dynamic scattering is assumed in deriving the normalized autocorrelation for a single dynamic scattering event (g 1 ss ), 3) the average number of dynamic scattering events,m(τ s ), is assumed to bem(τ s ) = p dyn µ s τ s c/n, where p dyn is the probability of dynamic scattering, µ s is the static scattering coefficient, c is speed of light in vacuum, and n is the index of refraction of medium, and 4) as described above, we explicitly allow for RBC Brownian motion, in addition to random flow considered Bonner and Nossal's original theory.
It is appropriate to remark on p dyn , the probability that a scattering event is dynamic. Note that p dyn is given by p dyn = V B µ s,blood /µ s , where V B is blood volume fraction and µ s,blood is blood scattering coefficient, and does not depend on anisotropy factor. Given µ s = 100 cm −1 , at a hematocrit of 33%, µ s,blood /µ s ≈ 7.8 at 850 nm [27], and p dyn = 10% corresponds to V B = 1.3%. On the other hand, the DCS probability of dynamic scattering, α [14,17], is given by , where g dyn and g stat are anisotropy factors of dynamic and static scatterers, respectively. If p dyn (1-g dyn ) << (1-g stat ), α p dyn (1-g dyn )/ (1-g stat ). For example, α 2% when p dyn = 10%, g dyn = 0.98 (typical for RBCs [27]) and g stat = 0.9.
Bonner and Nossal [24] assumed that multiple static scattering in tissue is sufficient to achieve random illumination of dynamic particles and uncorrelated dynamic scattering events. This enabled them to construct the field autocorrelation in the multiple dynamic scattering regime from the single dynamic scattering field autocorrelation, g 1 ss . Our extension of their expression is: In Eq. (2), P m is the probability that a photon will experience m scattering events with dynamic particles before leaving the medium, g 1 ss is the autocorrelation for a single scattering event from a moving particle with random illumination, as a function of time lag, τ d : In Eq. (3), the angle brackets θ denote an angular weighted average over dynamic scattering angle θ, and p HG (θ) is the Henyey-Greenstein scattering phase function for dynamic scatterers with an anisotropy of g dyn [28][29][30][31]: and π ∫ 0 p HG (θ) sin(θ)dθ = 1 / 2π (5) Figure 1(a) compares the Henyey-Greenstein (HG) phase function for g dyn = 0.6 (typical for Intralipid [32]) to g dyn = 0.98 (typical for RBCs [27]), with polar plots in the inset. Forward scattering dominates for g dyn = 0.98. Assuming that the number of dynamic scattering events P m has a Poisson distribution [24], as detailed in Appendix A, Eq. (2) becomes:

Critical examination of the cumulant approximation
The cumulant approximation conveniently brings the angular average of Eq. (3) into the exponent, providing a very simple and intuitive expression for the autocorrelation function. In this section, we revisit the cumulant expansion, paying special attention to errors incurred by the high scattering anisotropy of RBCs. Commonly, in the derivation of DWS [23], the normalized field autocorrelation is succinctly expressed as: In contrast to the Bonner and Nossal summation [Eq.
(2)], which posits an integer distribution of dynamic scattering events, Eq. (7) assumes that the number of dynamic scattering events is equal to the (possibly non-integral) average value. Note that while both Eqs. (6) and (7) yield the same first cumulant term, higher order cumulant terms are different. Here, we examine the cumulant expansion of the extended Bonner and Nossal expression in Eq. (6), which incorporates fewer assumptions. In this case, Taylor expansion can be applied to the natural logarithm of Eq. (6), as detailed in Appendix B, to yield the first and second cumulant terms: (8) Note that third and higher cumulant terms are contained in O( ∆r 2 (τ d ) 3 ), and the subscript θ in angular average has been dropped for readability. We examine the second cumulant to determine the error in the first cumulant approximation. Then, Eq. (8) can be simplified as: where R is the magnitude ratio of the second cumulant to first cumulant-squared, given by: where and q 4 (θ) = 2 4 k 4 sin 4 (θ/2) .
As shown in Fig. 1(b), R slowly increases with g dyn when g dyn < 0.9, and rapidly increases with g dyn when g dyn > 0.9. For example, whenm = 214.3 (τ s = 1000 ps and µ s = 100 cm −1 ), R increases 3.2 times from g dyn = 0.6 to g dyn = 0.9 and increases 4.8 times from g dyn = 0.9 to g dyn = 0.98, or overall 15.4 times from g dyn = 0.6 to g dyn = 0.98. Furthermore, the second cumulant term vanishes as R approaches zero [Eq. (9)], after many dynamic scattering events. This happens more rapidly for smaller g dyn ( Fig. 1(b)). Thus, the error in the first cumulant approximation merits particular attention for large g dyn . Ultimately, the functional form, or "shape" of the autocorrelation decay with respect to time lag, τ d , depends on which polynomial orders of τ d contribute in the exponent. With hybrid motion, different cumulant orders can contribute to the same polynomial order of τ d . Thus, another way to simplify Eq. (8) is to group terms according to the order of τ d : The term O(τ d 3 ) in Eq. (13) includes contributions from higher order τ d terms. Although written in different ways, Eqs. (9) and (13) both represent expressions for the cumulant expansion. Importantly, Eq. (13) clearly shows that both Brownian motion and random flow contribute to the τ d 2 term, representing the lowest order deviation from the pure exponential decay of DWS.
Note that the lag time range for neglecting higher order cumulants (τ d << 141 µs) coincides with the lag time range for neglecting random flow (τ d << 220 µs), thus this special case has little practical significance, given our assumed parameters. Nevertheless, it is included for completeness to isolate each individual assumption.

Special case 3: first order τ d term dominates (A1&2)
While special case 1 compares magnitudes of ordered and unordered motion, in special case 3, we examine τ d terms which determine the ultimate shape of the autocorrelation. From Eq. (13), the general condition for the first order τ d term to dominate the τ d 2 term is given by: Note that Eq. (15) can be more or less stringent than the condition in special case 1, depending on the relative magnitudes of ordered and unordered motion in the denominator. When Eq. (15) holds, Eq. (13) can be simplified, up to O(τ d 3 ), to Importantly, Eq. (16) is equivalent to conventional DWS [Eq. (1)], which is derived by neglecting both higher order cumulant terms and random flow (A1&2). Table 1 summarizes all three special cases.

First cumulant for random flow cancels second cumulant for Brownian diffusion
From Eq. (13), the coefficient of τ d 2 includes contributions from the first cumulant term for random flow and second cumulant term for Brownian motion. The velocity distribution may be selected such that these two contributions cancel, minimizing deviations from an exponential decay at small τ d . The distribution standard deviation, v*, is then given by: Note that this value of v* ensures that Eq. (16) is valid across lags. For an HG phase function with g dyn = 0.98 and D B = 6 × 10 −11 m 2 .s −1 [26,33,34], v* has a value of 1.25 mm/s. This is a reasonable value of RBC speed [35,36]. Thus, while our velocity distribution is chosen to eliminate the second cumulant, it remains well within physiological norms.

Monte Carlo simulations
Because experimental systems with Brownian motion and random flow are challenging to contrive, in this study, Monte Carlo is taken as the gold standard. A previous Monte Carlo algorithm simulating the time resolved reflectance of a semi-infinite medium [37] was modified to calculate the normalized field autocorrelation. In the modified version, both dynamic scattering (moving particles, with a customized anisotropy g dyn ) and static scattering (background tissue scattering, with a fixed anisotropy, g stat = 0.9) are explicitly simulated, with a probability of dynamic scattering p dyn . Only dynamic scattering events contribute to the recorded total accumulated momentum transfer Y dyn = Σ j (1-cosθ j ) [34,37,38], where θ j is the dynamic scattering angle at dynamic scattering site j (Fig. 2). By explicitly distinguishing dynamic and static scattering events, this approach obviates assumptions, including the cumulant approximation and similarity relation that are often invoked in standard Monte Carlo approaches [39]. Even though their momentum transfer is not stored, static scattering events are important in determining the Y dyn distribution. Unless otherwise stated, the Henyey-Greenstein phase function [Eq. (4)] is assumed for dynamic scattering. Monte Carlo simulations stochastically estimate the probability density P(Y dyn ,τ s ). When a photon reaches the detector, Y dyn associated with that photon is scored in a P(Y dyn ,τ s ) histogram. Then, the overall normalized field correlation is determined as the weighted average of field correlations across all dimensionless momentum transfers experienced by detected photon paths [3,37,40]: Equation (18) assumes uncorrelated motion at different scattering events, as the dynamic momentum transfer is summed in the exponent, which is equivalent to multiplying autocorrelations for individual dynamic scattering events. Note that the displacements at all scattering events are assumed to be Gaussian, independent, and identically distributed with variance ∆r 2 (τ d ) . Table  2 summarizes the simulation parameters, selected from previous studies of human brain [41][42][43][44] at a wavelength of 850 nm. To study the dependency of DWS/DCS validity on anisotropy of dynamic scatterers (g dyn ), a range of g dyn values (between 0.6 and 0.98) and corresponding speeds (v*), selected to eliminate the τ d 2 term in Eq. (17), were simulated. The range of values for g dyn encompasses particles of interest such as Intralipid (g dyn ≈ 0.6) [32,45,46], 1 µm diameter polystyrene microspheres (g dyn ≈ 0.9) [47][48][49][50][51], and red blood cells (g dyn ≈ 0.98) [27,52,53]. Unless mentioned, other parameters are kept constant as shown in Table 2. The propagation of a photon in the medium continues until the photon escapes from the medium, or until TOF exceeds a pre-defined 1 ns threshold. Photons that escape the medium within 5 cm from the source are collected. Absorption in the medium is included by multiplying the photon weight by the albedo at every scattering event [37,54]. The TOF bin width for the simulation is 5 ps, and the autocorrelation is assigned to the center point of each time bin. The normalized TOF-resolved field autocorrelation should not depend on absorption for our narrow TOF bins. Each simulation launched 300 million photons and took an average of 35 hours to run on an Intel i9 Model 9900K 3.6-GHz processor.

First cumulant approximation (A2)
First the consequence of ignoring higher order cumulants on the field autocorrelation is examined by comparing g 1 derived from Monte Carlo simulations (MC hybrid) to those derived from the first cumulant approximation (A2) [Eq. (14)], with hybrid motion, at two different TOF values, τ s = 47.5 ps (Fig. 3(a)) and 197.5 ps (Fig. 3(b)), and various g dyn values. These TOF values correspond to photon path lengths of 1.01 cm and 4.23 cm in the medium, respectively. In general, both methods agree well at small τ d whereas the long time lag tail of A2 is lower than that of MC hybrid, with more disagreement at early TOF ( Fig. 3(a)) and higher g dyn (Fig.  3(b)-(c)). While inaccuracy of the cumulant approximation at early TOF (few dynamic scattering events) is well-known [23], its dependency on anisotropy of dynamic scatterers has not been rigorously investigated. Surprisingly, discrepancies in the tails of the autocorrelation extend beyond τ s ∼ 500 ps for g dyn = 0.98 (Fig. 3(c)). A zoomed-in linear scale of Fig. 3(c) shows that, for high g dyn , neglecting higher order cumulants significantly underestimates the tails even at large TOFs (Fig. 3(d)). The role of g dyn in lowest order error term in the cumulant approximation is described by R, the ratio of the second cumulant term to the square of the first cumulant term [Eq. (10) and Fig.  1(b)]. Larger R values at early TOF and larger g dyn ( Fig. 1(b)) predict larger errors in the first cumulant approximation.

Neglecting random flow (A1)
Next the consequence of ignoring random flow on the field autocorrelation is examined by comparing Bonner and Nossal's theory for Brownian motion alone (A1) to MC hybrid (Fig. 4). As described earlier, A1 accounts for all cumulant terms, but ignores random flow. A1 agrees well with MC hybrid at small τ d , but decorrelates slower than MC hybrid at large τ d . In other words, ignoring random flow increases the autocorrelation tails. This phenomenon is more obvious at shorter τ s (Fig. 4(a)) and at higher g dyn (Fig. 4(b)-(d)) because for fixed particle dynamics, the autocorrelation decays slower for high g dyn (e.g. g dyn = 0.98) than for low g dyn (e.g. g dyn = 0.6).

Conventional DWS (A1&2)
As mentioned earlier, conventional DWS ignores both random flow and higher cumulant terms (A1&2). Even though A1 and A2 are individually incorrect, Eq. (16) is still approximately satisfied as the velocity distribution width was chosen [Eq. (17)] so that it is valid up to O(τ d 3 ). Thus A1&2 surprisingly agrees well with the gold standard MC hybrid, especially when g dyn < 0.8 (Fig. 5(a)), or when g dyn = 0.98 and τ s > 200 ps (Fig. 5(c)-(d)). The fact that A1&2 performs better than A1 and A2 is surprising because A1 and A2 include only one incorrect assumption  each, whereas A1&2 includes two incorrect assumptions. This result can be explained by the fact that the autocorrelation tails can be reduced either by neglecting higher order cumulant terms (observed in A2 in Fig. 3) or by including random flow (suggested in A1 in Fig. 4). Thus, for A1&2, ignoring flow increases the tail whereas ignoring higher order cumulants reduces the tail. Therefore, in this case, two incorrect assumptions (Fig. 5) provide better predictions than just one ( Fig. 3 and Fig. 4). Residual minor errors in Fig. 5 are due to O(τ d 3 ) terms in Eq. (13). To illustrate the role of the first cumulant approximation in DWS theory for pure Brownian motion or diffusion, in Fig. 6, Bonner and Nossal's theory (A1) and Monte Carlo simulations (MC diff.) are compared to DWS theory (A1&2) In this case, all cumulant terms are accounted for in MC diff. and A1, while only the first order cumulant term is accounted for in A1&2. As expected, A1 and MC diff. agree well while the tail of A1&2 is much lower, even at τ s = 997.5 ps (Fig. 6(c)-(d)). This further confirms that neglecting higher order cumulants reduces the autocorrelation tails.

Full cumulant expansion (Bonner and Nossal) with hybrid motion
To avoid assumptions in DWS/DCS theory, Bonner and Nossal's approach of directly integrating over the scattering phase function can be used. Thus, a B&N hybrid model accounts for all cumulant terms, and for both Brownian diffusion and random flow of dynamic scatterers. As shown in Fig. 7, B&N hybrid is an improvement over A1&2, agreeing well with MC hybrid regardless of anisotropy of dynamic scatterers (Fig. 7(a)-b), and the agreement is observed at short (Fig. 7(a)) and long τ s (Fig. 7(c)-(d)). Thus, eliminating incorrect assumptions altogether achieves excellent agreement with the gold standard.

Discussion
Using a customized Monte Carlo of light scattering in realistic tissue models as a gold standard, we highlight that errors due to two common, incorrect assumptions, one regarding the nature of RBC motion and the other regarding the form of the field autocorrelation, can cancel, yielding much better predictions of TOF-resolved field autocorrelations than just one incorrect assumption. This suggest that goodness-of-fit alone is a flawed measure by which to assess assumptions underlying a theory. By carefully examining each assumption, we arrive at a more accurate approach to assess random flow and further improve goodness-of-fit.

Accuracy of the cumulant approximation
The cumulant approximation, or first cumulant expansion, is widely used in light scattering theory, because it yields simple analytic expressions [1,7,8,23,55]. Here we argue that errors in the cumulant approximation require special attention, particularly at late time lags for highly anisotropic dynamic particles such as RBCs and, as discussed further below, in tissues with low probabilities of dynamic scattering. As such, we encourage the view that the cumulant approximation, which is central in the diffusion approximation for correlation transport, be viewed and evaluated essentially independently of the diffusion approximation for radiative transport. While the cumulant approximation can be improved by restricting fitting to early time lags (τ d ), our results show significant disagreements at late time lags when τ s ≈ 200 ps ( Fig.  3(b) and Fig. 6(b)). Importantly, as described in the next section, assessing RBC random flow requires analyzing longer time lags where the accuracy of the cumulant approximation becomes questionable.

Revealing random flow
The appropriate description of RBC motion (ordered vs. unordered) has been the topic of continued debate in the field of laser speckle [8,55]. Certainly, as discussed in Section 2.3.1, at sufficiently short time lags (i.e. τ d << 220 µs), Brownian motion dominates random flow. Meanwhile, random flow just impacts the tails of TOF-resolved field autocorrelations (for instance, where A1 is higher than MC hybrid in Fig. 4). Therefore, assessing random flow is thus only feasible at later time lags.
As discussed in Section 2.3.2, the cumulant approximation is invalid if τ d ≈ 141 µs or larger, precisely the lag time range needed to assess random flow. This revelation adds a new wrinkle to the debate on the nature of RBC dynamics, as the cumulant approximation is assumed in standard models for ordered and unordered flow in the literature [7]. Given our assumptions (g dyn = 0.98, v = 1.25 mm/s and D B = 6 × 10 −11 m 2 .s −1 ) we hypothesize that the first cumulant approximation cannot accurately assess random flow.
To test this hypothesis, we treated the MC hybrid model as experimental data and fit either the A2 or B&N hybrid model to this data, with D B and v as free parameters. Figure 8 compares the fitted D B and v values to the expected (actual) values. Two different fitting ranges of the MC field autocorrelation are examined: g 1 MC > 0.01 ( Fig. 8(a)-(c)), and g 1 MC > 0.5 (Fig. 8(d)-(f)). B&N hybrid accurately recovers D B and v while A2 underestimates D B and v. As A2 ignores higher order cumulants, the fitting procedure reduces D B to approximate the long time lag tail in g 1 MC . With a restricted fitting range (to exclude the tail of g 1 MC ), A2 more accurately estimates D B (Fig. 8(d)). However, notably, A2 always yields v∼0 mm/s, suggesting that A2 is a poor model for quantifying random flow.
Thus, particular caution is warranted in applying the first cumulant approximation for small numbers of dynamic scattering events. Thus, for laser speckle flowmetry in parenchymal regions with sparse blood vessels, the range of autocorrelation lags (or the range of exposure times) must be chosen judiciously to avoid model errors.

Dynamic scattering probability dependence
In the above sections, we show that one confound in applying DWS/DCS theory to tissue measurements is the exceptionally high g dyn of RBCs. Here, we highlight another important confound; the low dynamic scattering probability (p dyn ) of tissues. As mentioned in Section 2.4, a smaller V B implies a smaller p dyn . To illustrate the effect of p dyn on DWS/DCS accuracy, p dyn = 30% (V B = 3.9%, Fig. 9(b)) was compared to p dyn = 10% (V B = 1.3%, Fig. 9(a)), used heretofore in this work. Obviously, smaller p dyn values result in a slower autocorrelation decay due to the smaller number of dynamic scattering events. For both p dyn , A1&2, representing conventional DWS/DCS, underestimates the autocorrelation tails of MC hybrid, which are well-represented by B&N hybrid. Importantly, disagreement between A1&2 and MC hybrid persists to later TOFs for smaller p dyn , due to the smaller number of dynamic scattering events. For example, A1&2 agrees well with MC hybrid when p dyn = 30% at 497.5 ps (Fig. 9(b)), but still decays faster when p dyn = 10% ( Fig. 9(a)).
To quantify agreement between various models as p dyn changes, the coefficient of determination (r-squared) versus τ s is shown (Fig. 10). The MC hybrid model is taken as the gold standard, and hence has an r-squared of unity. Different ranges of the field autocorrelation function are examined, corresponding to g 1 MC > 0.01, 0.2, and 0.5. As expected, B&N hybrid agrees best, with an r-squared approaching unity in all cases, while A1&2 (conventional DWS/DCS) is the second best, followed by A1 (B&N diffusion). Limiting the fitting range greatly improves the   9. Accuracy of DWS/DCS assumptions depend on the probability of dynamic scattering, p dyn : g 1 from less accurate models (the first cumulant approximation for Brownian motion (A1&2) and B&N diffusion (A1)), disagrees with more accurate models (B&N hybrid and MC hybrid, the gold standard) up to longer τ s for p dyn = 10% (a) than for p dyn = 30% (b) when g dyn = 0.98. r-squared for A1, especially at early TOF. For example, the lowest r-squared value for A1 is well below 0.4 (0.9) if the tail of autocorrelation function is included (Fig. 10(a),(d)), and above 0.4 (0.9) if the tail is excluded (Fig. 10(c),(f)) when p dyn = 10% (30%). These results are expected since the first cumulant term for Brownian motion dominates early time lags. Fig. 10. Quantification of results in Fig. 9. Coefficient of determination (r-squared) based on comparing the gold standard MC hybrid model to different models: A1, B&N hybrid, and A1&2 when p dyn = 10% (a-c), and p dyn = 30% (d-f). Note that by definition, the MC hybrid model has an r-squared of unity across TOFs.
Increasing p dyn (the number of dynamic scattering events) improves the r-squared for A1 and A1&2. For instance, A1 yields r-squared below 0.6 at early TOFs when p dyn = 10% ( Fig. 10(a)), and well above 0.6 at early TOFs when p dyn = 30% (Fig. 10(d)). The improvement in A1 is expected since with more dynamic scattering events, the first cumulant term for Brownian motion dominates more of the autocorrelation decay. The improvement in A1&2 is also expected since higher order cumulants are less important with more dynamic scattering events.

Proposed improved expression
A common approach (COM) to determine field autocorrelations is a weighted average over total dimensionless momentum transfer Y tot [4,40,56]: In homogenous media, note that this approach is equivalent to DWS, if Y tot = µ s (1-g stat )s, where s=τ s c/n is the photon path length in the medium. In Eq. (19), Y tot is the total dimensionless momentum accumulated from all scattering events and can be obtained from a standard Monte Carlo method that does not explicitly simulate dynamic scattering [38]. To account for dynamic scattering events, the momentum from simulations is multiplied by α (note that α ≈ 0.2p dyn when g dyn = 0.98 and g stat = 0.9). This approach implicitly assumes the cumulant approximation and a simple similarity relationship between dynamic and total momentum transfer: Y dyn = αY tot .
Here, we argue that the most accurate approach is to use a Monte Carlo simulation that includes both static and dynamic scattering events (Section 2.4), as this approach correctly performs angular averaging over the scattering phase function, implicitly including all cumulant orders. To achieve accurate results with a standard Monte Carlo code that does not explicitly include separate dynamic and static scattering events, we propose an improved expression (IMP): Rewriting Eqs. (19) and (20), respectively, to take a weighted average of field autocorrelations over all detected photon paths, each with weight w i [37,57], we arrive at and In Eqs. (21) and (22), w i is the weight and τ si is the TOF for the detected photon i, g 1 ss is the normalized single scattering autocorrelation function [Eq. (2)],m i = αY tot,i / (1-g dyn ) is the average number of dynamic scattering events (Alternatively, similar results are obtained with m i = p dyn µ s s i ). To compare the common method [Eq. (21)] to our more accurate customized Monte Carlo method (MC hybrid), and to validate the proposed improvement [Eq. (22)], we record Y tot using the GPU-based code Monte Carlo eXtreme (MCX) developed by Fang et al. [38]. The MCX simulation, which accumulates Y tot from all scattering events [39], is used to determine the field autocorrelation for both the common approach [Eq. (21)] and the improved approach [Eq. (22)]. The simulation launches 10 8 photons at the source, which took approximately 5 minutes on an NVIDIA Dual GeForce RTX 2080 GPU. The simulation returns a list of all detected photons and their individual momentum transfers and partial paths. Absorption is included by applying the Beer-Lambert law to determine the photon weight, w i . To approximate the 5 cm detection radius in the gold standard MC simulations (Section 2.4), in the MCX simulations, 11 detectors with radii of 0.5 mm were evenly spaced from 0 to 5 cm from the source. All detected photons were included when applying in Eqs. (21) and (22).
As shown in Fig. 11, the proposed method (IMP), with a more accurate autocorrelation (adapted from Bonner and Nossal), closely approximates our customized Monte Carlo method (MC hybrid) which explicitly simulates dynamic and static scattering, resulting in the best agreement with near unity r-squared both when p dyn = 10% (α = 2%) [ Fig. 11(e)-(f)], and when p dyn = 30% (α = 6%) [ Fig. 11(g)-(h)]. On the other hand, the common approach (COM hybrid), which assumes the cumulant approximation, underestimates the autocorrelation tails when τ s is less than 300 ps. While agreement improves for τ s > 300 ps, these discrepancies are significant because the cumulant approximation is widely assumed to be valid in DCS/DWS across a range of source-detector separations. In summary, our proposed method [Eq. (22)]: 1) improves the underestimation of autocorrelation tails observed in the common method (and DWS/DCS theory), and 2) enables the use of highly optimized existing Monte Carlo tools. Another potential advantage of the proposed method is that it can be extended to multi-layered tissues, by tracking momentum transfer or partial path length in each tissue type, as is possible in MCX.

Gegenbauer phase function
Bonner and Nossal's theory and the MC method incorporate all cumulants and depend on the entire scattering phase function, unlike DWS theory (A1&2), which just depends on cosθ . In this section, we examine the effect of phase function choice on the autocorrelation. Early goniophotometric measurement of RBC suspensions showed that light scattering by RBCs is well-described by a Gegenbauer (GE) phase function [58][59][60][61]. The GE phase function which is given by: Equation (23) requires |g GE | ≤ 1 and a >-0.5. Note that g GE cosθ in Eq. (23), and that only if a = 0.5, Eq. (23) reduces to Eq. (4), with g GE = g dyn = cosθ . Thus, the GE phase function is a generalization of the HG phase function. If cosθ = 0.98 and a = 1.2 [58], we obtain that g GE = 0.893. As shown in Fig. 12(a), HG and GE phase functions are approximately the same when cosθ = 0.6 but are significantly different when cosθ = 0.98, especially at large scattering angles. The cumulant approximation error, R, is 2.35 times lower at cosθ = 0.98 for GE ( Fig.  12(b)). The inset in (b) shows the cumulant approximation error ratio between HG and GE, which is independent ofm. The upper x-axis in (b) shows the g GE value corresponding to cosθ on the lower x-axis.
As above, the velocity distribution v* is selected so that the first cumulant term for random flow cancels second cumulant term for Brownian diffusion. Like R, v* depends on phase function [Eq. (17)], and v* = 0.817 mm/s when a = 1.2, g GE = 0.893 and D B = 6 × 10 −11 m 2 .s −1 [26,33,34]. Figure 13 shows that while the autocorrelation tails are greater for the HG phase function (B&N hybrid, HG) than the GE phase function (MC hybrid, GE and B&N hybrid, GE), all have longer tails than DWS (A2), which ignores higher order cumulants. The smaller cumulant approximation error for the GE phase function is explained by smaller R values in Fig. 12(b).  Figure 14 also confirms that ignoring random flow (A1, GE) increases the autocorrelation tail. Thus assuming a GE rather than HG phase function does not change the essential conclusions regarding autocorrelation tail; namely, that neglecting higher order cumulants reduces the tail, neglecting random flow increases the tails, and that both assumptions in tandem lead to better agreement with the gold standard MC simulation. Moreover, Fig. 14 shows that conventional DWS/DCS (A1&2) with two incorrect assumptions agrees well with the more accurate models that use the GE phase function for both p dyn = 10% and p dyn = 30%.  Fig. 9 with the GE phase function (a = 1.2, g GE = 0.893): once again g 1 from less accurate models (A1&2 and A1), disagrees with more accurate models (B&N hybrid and MC hybrid, the gold standard) up to longer τ s for p dyn = 10% (a) than for p dyn = 30% (b) when cosθ = 0.98 Figure 15 shows r-squared versus τ s for different ranges of the field autocorrelation. As before, the MC model is taken as the gold standard. Like Fig. 10, Fig. 15 shows that MC hybrid, GE agrees well with A1&2 (two incorrect assumptions) at all TOFs but disagrees with A1, GE (one incorrect assumption) at earlier TOFs, and that r-squared for A1, GE improves in tissues with higher p dyn , and when the range is limited to early time lags.

Simplifying assumptions
To arrive at our proposed improved expression for field autocorrelations in biological tissue [Eq. (22)], several assumptions were adapted from the original Bonner and Nossal theory [24]. Specifically, it was assumed that RBCs move independently, that averaging over photon paths and scattering events is equivalent to ensemble averaging, that the distances between static tissue source points and moving cells is large compared to displacements for decorrelation, and that static scattering randomizes the direction of incident light between dynamic scattering events. In incorporating a hybrid model of RBC motion, we have also assumed that advective and Brownian components of the motion are independent.
Importantly, vascular geometry was not explicitly simulated in the current study. RBCs were assumed to be uniformly distributed throughout the tissue volume, parameterized by the probability of dynamic scattering (p dyn ). In a previous study with a proscribed parallel vessel geometry, p dyn was empirically found to be proportional to vessel radius-squared and inversely proportional to vessel spacing-squared [34]. While our simplified tissue model enabled us to investigate the cumulant approximation without assuming a specific geometry, future studies should investigate realistic vascular geometries, accounting for the issues highlighted in the present work.

Conclusion
We have shown that field autocorrelations in realistic biological tissue models with hybrid motion are accurately predicted by our adaptation of Bonner and Nossal's original theory for RBC scattering. Importantly, neglecting higher order cumulants reduces the field autocorrelation tails whereas neglecting random flow increases them. Thus, for particular choices of random flow distribution and Brownian diffusion coefficient, DWS theory with pure Brownian motion might agree with experimental data, leading to the incorrect conclusion that random flow is absent. Furthermore, if experimental data are fit using the DWS cumulant approximation, allowing for hybrid motion, dynamics are underestimated. Here, we argue that dynamics can be accurately assessed by analyzing field autocorrelations with a model that includes higher order cumulants. Altogether, this work supports that random flow may be assessable in TOF-resolved field autocorrelations, but higher order cumulants should be incorporated into models to ensure its accurate assessment.

Appendix A
Here, we derive Eq. (6). In Bonner and Nossal's original theory [24], the normalized field autocorrelation function can be expressed as where I m is the autocorrelation function for a path, and P m is the probability of such a path. When the lower limit in the summation is set to zero to include purely static scattering paths, we re-write Eq.(A.1):