Precision and bias in dynamic light scattering optical coherence tomography measurements of diffusion and flow

: We quantify the precision and bias of dynamic light scattering optical coherence tomography (DLS-OCT) measurements of the diffusion coefficient and flow speed for first and second-order normalized autocovariance functions. For both diffusion and flow, the measurement precision and accuracy are severely limited by correlations between the errors in the normalized autocovariance function. We demonstrate a method of mixing statistically independent normalized autocovariance functions at every time delay for removing these correlations. The mixing method reduces the uncertainty in the obtained parameters by a factor of two but has no effect on the standard error of the mean. We find that the precision in DLS-OCT is identical for different averaging techniques but that the lowest bias is obtained by averaging the measured correlation functions before fitting the model parameters. With our correlation mixing method, it is possible to quantify the precision in DLS-OCT and verify whether the Cramer-Rao bound is reached.


Introduction
Dynamic light scattering optical coherence tomography (DLS-OCT) relies on the measurement of fluctuations of scattered light and coherence gating to obtain simultaneous depth-resolved information about diffusive and translational motion of particles.Initially, DLS-OCT was used for measuring diffusion coefficients [1] and flow speeds of particle suspensions [2][3][4][5].Later, several improvements have been suggested for increasing the DLS-OCT flow velocity dynamic range [6,7].
DLS-OCT has the advantage over phase-resolved Doppler OCT that a flow can be measured for zero Doppler angle.It can also be used for particle sizing where the particle size is determined from the estimated diffusion coefficient using the Stokes-Einstein relation.The combined flow and diffusion estimation is particularly interesting for in-line particle sizing during process control [8].Sensitivity and precision of phase-resolved Doppler OCT has been widely studied and reported in literature [9][10][11][12].However, there is very little information available about the precision and bias of DLS-OCT diffusion and flow measurements.These are crucial for reliable measurements, especially in medical and pharmaceutical applications.Effects of noise and bias in OCT on the measured autocorrelation function have been reported [13], but their influence on the underlying parameters remains unclear.
In this work we perform simulations, measurements, and a theoretical analysis to quantify the precision and bias of diffusion and flow measurements using DLS-OCT.In our analysis, we consider both the first and second-order normalized autocovariance functions, g 1 (z, τ) and g 2 (z, τ), diffusive particle motion, and, for both dilute and non-dilute suspensions, translational particle motion.We derive analytical expressions for the highest attainable precision using the Cramer-Rao bound and compare them with the obtained results.We also assess the bias for different averaging techniques.

Theory
The geometry for OCT diffusion and flow measurements is described in [6,7].The propagation of the optical beam is described by a Gaussian beam along the z-direction.We assume that the scattering process is stationary.

Correlation functions for particle diffusion
For a non-flowing particle suspension, the normalized depth-dependent autocovariance of the OCT complex-valued signal in a backscattering geometry, including the effect of SNR, is given by [3,4,[13][14][15] where E(z, t) is the depth and time-dependent complex-valued OCT signal, I(z, t) is the OCT signal intensity, τ is the autocovariance time lag, A 1 (z) is the autocovariance amplitude containing the effect of a diminishing signal-to-noise with depth [13], SNR(z) is the depth-dependent experimental signal-to-noise ratio [7,13], D is the particle diffusion coefficient given by Stokes-Einstein equation, and q = 2nk 0 is the scattering wavenumber for the OCT backscattering probe configuration with the incident light wavenumber in vacuum k 0 and the medium refractive index n.Equation (1) for g 1 (z, τ) is also known as the first-order autocorrelation function where ⟨E(z, t)⟩ t = ⟨E * (z, t)⟩ t = 0.The signal-to-noise correction describes the non-zero time lag autocorrelation.At τ = 0, the normalized correlation coefficient is unity.The autocorrelation function of the mean-subtracted OCT signal intensity, also known as the Pearson correlation or autocovariance, decays twice as fast than g 1 (z, τ) [1,15] and can be expressed with the normalized second-order autocovariance using the Siegert relation [16,17] g 2 (z, τ) = ⟨(I(z, t) − ⟨I(z, t)⟩ t ) (I(z, t + τ) − ⟨I(z, t)⟩ t )⟩ t ⟨I(z, t)⟩ )︂ 2 = A 2 (z)e −2Dq 2 τ , (2) where A 2 (z) is a depth-dependent amplitude factor for the intensity autocorrelation function.For a mean-subtracted intensity autocovariance, the Siegert relation states that g 2 (z, τ) is the square of |︁ |︁ g 1 (z, τ) |︁ |︁ .In Eq. ( 2) we have assumed that the average number of particles in the scattering volume, N s , is sufficiently large (N s ≫ 1) [16][17][18][19].This ensures that the particle probability distribution in the scattering volume and the scattered light fluctuations follow Gaussian statistics.For a non-flowing particle suspension this requirement is almost always satisfied whenever the backscattered OCT signal intensity from every voxel is high [7].

Correlation functions for particle flow
For diffusing and flowing particle suspensions the first-order normalized autocovariance function magnitude is where v(z) is the depth-dependent total flow speed, θ is the Doppler angle, w z is the coherence function waist, w 0 is the Gaussian beam waist in focus, and A 1 (z) is the same SNR correction factor as given in Eq. (1).We take the absolute value to get rid of the phase component in g 1 (z, τ) originating from the particle translational motion along the optical axis [2][3][4].
For flow measurements we focus on the second-order normalized autocovariance function, g 2 (z, τ), that does not depend on the phase, is easier to implement, and can also be implemented in phase-unstable OCT systems.Here, we differentiate between very dilute and non-dilute sample regimes.For the flowing non-dilute particle suspensions, where the number of particles in the scattering volume, N s , is much greater than one, g 2 (z, τ) is again obtained using the Siegert relation [3,4,6,14,15] where A 2 (z) is the same SNR correction factor as given in Eq. ( 2).Here, we have neglected the effect of a gradient of the axial velocity on the autocovariance function [20].In the non-dilute regime the scattering process is strictly Gaussian [7].
In dilute suspensions, the expected number of particles in the scattering volume can be significantly lower than 1.Therefore, for stationary particles certain depth voxels would give zero signal.However, due to the translational particle motion, the scattering signal is obtained from every voxel during the acquisition time.When the number of particles in the scattering volume is very low, the scattering process becomes non-Gaussian, and the Siegert relationship does not apply anymore [18,21,22].In the dilute case with N ≲ 1, the second-order normalized autocovariance is [7,13,19] in which g 1 (z, τ) can be incorporated as follows Here, A 3 (z) is given by where w(z) is the radius of the local beam waist, N s (z) is the average depth-dependent number of particles in the scattering volume [7,19], and κ(z) is the kurtosis of the noise-subtracted complex field distribution.The kurtosis can also be expressed as a ratio of the average squared noise-subtracted intensity to the squared mean noise-subtracted intensity.It can be obtained using the measured OCT signal intensity I(z, t) and the signal-to-noise ratio, which simplifies to κ(z) = 2 for the Gaussian scattering process with )︂ −2 .The kurtosis can be depth-dependent for the non-Gaussian process due to different signal intensities and particle flow speeds at different depths.In our previous study [7] we assumed κ(z) = 2, even though the scattering process was non-Gaussian.This made our analysis simpler considering that κ(z) minimally affects A 3 (z) when the signal-to-noise ratio is sufficient.Here Eqs. ( 5)-( 8) use no assumptions on the scattering process, κ(z), and do not impose any restrictions on the number of particles.

Precision of diffusion and flow estimation
Precision is defined as the spread of the measurement around the mean and can be quantified by the square root of the measurement variance.For DLS-OCT the precision is governed by random errors, i.e., unpredictable changes in the measured intensity.Since the particle diffusion coefficient and the flow speed are determined by fitting the models from Sec. 2.1 and Sec.2.2 to the measured g 1 (z, τ) and g 2 (z, τ), the precision of the fit parameters is determined by the random errors of the normalized autocovariance function at different time delays.The maximum obtainable precision in the fit parameters is determined by the Cramer-Rao lower bound (CRLB).It is computed by inverting the Fisher information matrix [23].For calculating the Cramer-Rao lower bound the probability distribution functions of g 1 (z, τ) and g 2 (z, τ) must be known.The correlation coefficients from a long time series data are approximately normally distributed [24][25][26].Typical diffusion or flow measurements contain at least thousands of temporal sampling points, which is more than enough to assume that g 1 (z, τ) and g 2 (z, τ) are normally distributed at every time delay.In general, errors (residuals) at different time delays can be correlated.Calculating the CRLB becomes challenging because we don't know the error correlations and their impact on the model bias.To address this, we assume uncorrelated residuals.In this case, for M observations g(z, τ 1 ), g(z, τ 2 ), . . .g(z, τ M ) with N model (fit) parameters p 1 , p 2 , . . .p N , the Fisher matrix F is an N × N symmetric matrix given by where τ m is the discretized τ for index m.The observations depend on the model parameters and have a variance σ 2 (z, τ m ).In our analysis the observations are g 1 (z, τ m ) or g 2 (z, τ m ) and the summation is performed over all considered time delays.The covariance matrix of the model parameters is obtained by inverting the Fisher information matrix.Consequently, the variances of our model parameters, σ 2 p j , are represented by the diagonal entries of the covariance matrix.This relationship is described by the equation

Precision of estimating the diffusion coefficient
For a non-flowing particle suspension the diffusion coefficient is determined by fitting Eq. ( 1) to the real part of g 1 (z, τ m ), with A 1 (z) and D being the fit (model) parameters.The lowest attainable variance CRLB in the fitted diffusion coefficient, assuming uncorrelated noise in the Fisher information matrix, is given by where σ 2 ℜ(g 1 ) (z, τ m ) represents the variance of the real part of g 1 (z, τ m ), which is half of the variance of the complex g 1 (z, τ m ).It can be obtained with simulations or measurements and can also be derived analytically using where M is the total number of temporal sampling points (time series length).The first term represents the variance in the real part of the field autocorrelation and is equal to one half, the second term is a correction for the explained variance in the Pearson correlation coefficient [24][25][26][27][28], and the third term is a long-lag correction for the effect of the correlation magnitude on the variance [29].As expected, σ 2 ℜ(g 1 ) (z, τ m ) = 0 when ℜ Even though g 1 (z, τ m ) is complex in nature, we only consider its real part when fitting the diffusion coefficient.According to Eq. ( 1), in non-flowing suspensions, the imaginary part ℑ(g 1 (z, τ m )) is pure noise and contains no information about the particle diffusion.Therefore, in phase-stable systems it is always beneficial to use the real part ℜ(g 1 (z, τ m )) instead of the absolute value |g 1 (z, τ m )|.
A similar analysis can be performed with g 2 (z, τ m ).In this case we use Eq. ( 2) for fitting with model parameters D and A 2 (z).The Cramer-Rao lower bound on the variance of the diffusion coefficient, when using the second-order normalized autocovariance function, is with σ 2 g 2 (z, τ m ) being the variance of g 2 (z, τ m ) given by Here the higher order (3 rd and 4 th ) intensity autocorrelation functions were calculated using theory derived by Lemieux et al. [30].This analysis is limited to the Gaussian scattering process.Equation ( 14) incorporates the Siegert relation and relies on the fact that the diffusive g 1 (z, τ m ) is real-valued.

Precision of velocity estimation
The velocity of a non-dilute flowing particle suspension can be obtained by fitting Eq. ( 4) at every depth z to the measured g 2 (z, τ m ).In this case v(z) and A 2 (z) are the model parameters, while D, w 0 , w z , and θ are assumed to be exactly known a-priori (calibrated).The minimum achievable variance of the flow speed v is given by where σ 2 g 2 (z, τ m ) can also be computed analytically similar to Eq. ( 14) using [30].However, in this case phase terms of g 1 (z, τ m ) need to be taken into account.
For a dilute suspension the velocity can be determined by fitting Eq. ( 5) to the measured g 2 (z, τ m ).In this case v(z) and A 3 (z) are the model parameters, while D, w 0 , w(z), w z , N s (z) and θ are assumed to be known in advance.The minimum achievable variance of the flow speed is given by and where σ 2 g 2 (z, τ m ) can no more be calculated using an equation similar to Eq. ( 14) using [30] due to a non-Gaussian scattering process.
For comparison, the minimum achievable variance of the flow speed when using 3) is independent of the particle concentration and given by (21)

Bias of diffusion and flow estimation
Bias in DLS-OCT is determined by systematic errors in the correlation function and is a measure of the accuracy of the method.The bias is not random but leads to a consistent over or under estimation of the fit parameter.The systematic errors arise when the model parameters in DLS-OCT are determined by fitting the models from Sec. 2 to the measured g 1 (z, τ m ) and g 2 (z, τ m ).These error sources are: • Estimation bias in g 2 (z, τ m ) due to subtraction of the sample mean instead of the true mean from the OCT signal intensity.This introduces an almost constant offset to the normalized autocovariance function given by [31] g where g 2 biased (z, τ m ) and g 2 true (z, τ m ) are biased and unbiased correlation coefficients, respectively.Calculation of the bias requires a-priori knowledge of a true correlation coefficient, which is not possible in practice.The estimation bias term in Eq. ( 22) is independent of the correlation coefficient and depends only on the correlation decay rate, SNR, and M. It can be reduced by increasing the intensity time series length with respect to the characteristic decay time of g 2 (z, τ m ).As reported in literature [13], the estimation bias (also referred to as the statistical bias) also affects the amplitude and decay rate of g 2 (z, τ m ) and even contains a random component which is not included in Eq. ( 22).This randomness can be reduced by averaging multiple g 2 (z, τ m ).Theoretically both g 1 (z, τ m ) and g 2 (z, τ m ) suffer from the additional estimation bias when subtracting a DC term from the OCT interference spectra before inverse Fourier transformation, if the DC term is estimated from the interference time series itself.In our analysis this effect is neglected.
• Sampling distribution bias in g 1 (z, τ m ) and g 2 (z, τ m ) is caused by the slight skewness of the correlation coefficient distribution [28] leading to a bias of the sample mean with respect to the true value.This bias depends on the correlation coefficient and decreases with increasing time series length as the distribution better approaches a normal distribution.A simple correction factor for this bias is given by [25,26,28] which is negligible for long acquisitions.For example, the sampling distribution bias for our diffusion measurements in Sec.5.1 is of the order of 0.01% of the correlation coefficient and can be disregarded without any corrections.
• Curve fitting bias in g 1 (z, τ m ) and g 2 (z, τ m ) is caused by the use of incorrect fit models.For example, the exponential fit curves from Sec. 2 cannot be negative at any time delay, whereas measured g 1 (z, τ m ) and g 2 (z, τ m ) can be, due to their probability distributions.So, fits to g 1 (z, τ m ) and g 2 (z, τ m ) with asymmetric noise will be slightly biased, which can happen if random errors are correlated along τ m .This bias can be lowered by increasing the time series length and reducing the variance of g 1 (z, τ m ) and g 2 (z, τ m ).Additional fit bias arises if a wrong model is used in the analysis, i.e., a non-dilute model for very dilute suspensions, static fit models for flowing samples, etc.Furthermore, the models are based on assumptions that may be invalid.For example, we assume that the mean scattering wavenumber q is exactly known and that the scattering process is stationary.
The polydispersity in particle size and refractive index can also add to the bias.Estimation bias also adds to the fit bias since it introduces a constant offset that cannot be incorporated in the models.
The bias for a sufficiently large M is predominantly represented by estimation and curve fitting errors.Since it affects the normalized autocovariance functions and not directly the model parameters, its influence on the fitted diffusion coefficient and the flow speed cannot be evaluated analytically.Predicting the bias is also impossible experimentally due to a lack of a-priori knowledge and is only feasible using simulations.In addition, systematic errors make our estimators from Sec. 2 biased which affect the random errors and prevents us from reaching the maximum precision of our model parameters.

OCT signal simulations for flow and diffusion
In simulations we generate a complex OCT signal scattered from an ensemble of particles with random (diffusion) and directional (flow) motion.The signal intensity and the normalized autocovariance functions are subsequently calculated.This is repeated N b times which allows us to compute the autocorrelation variances at every time delay.

Simulation of diffusion
Noiseless time-dependent complex field scattered from N p diffusing particles is simulated using where z j (t) is the j th particle axial position at time t due to the Brownian motion generated using normally distributed steps with σ = √ 2D∆t.The initial particle positions z j (0) are uniformly distributed in space.In Eq. ( 25) we have neglected the effect of the PSF on the scattered field fluctuations [3], considered motion solely in the depth direction, and assumed that all particles have an identical scattering cross section.The effects of Brownian motion in the lateral direction can be disregarded since they occur on a timescale much longer than that assessed with DLS-OCT.Inclusion of the noise in the scattered field then leads to [32,33] where φ(t) is a time-dependent random phase angle uniformly distributed between 0 to 2π.Note that this is only a good approximation when SNR ≫ 1 because we assume that field magnitude fluctuations are negligible and the noise originates purely from the random phase angle.

Simulation of flow
The scattered field from suspensions flowing perpendicular to the beam optical axis is simulated by replacing E 0 (t) in Eq. ( 26) with where v(z) is the depth-dependent transverse bulk velocity, w(z) is the local beam waist, R(z) the radius of curvature of the Gaussian beam, and x j is the uniformly distributed initial random transverse position of j th particle.The additional factor 2 in the radial PSF function is due to coupling efficiency of the scattered light in a confocal setup [2,3,34].
As we solely model the transverse flow, leading to identical decorrelation regardless of its direction within the transverse plane, we consider only one lateral dimension, denoted as x.This allows us to reduce the computational complexity and the number of required particles to be modeled.The particles are simulated with random, uniformly distributed x-positions between −L(z) and +L(z) given by where N s (z) is the number of particles in the scattering volume and N p is the total number of simulated particles.Equation (28) guarantees that the number of particles in the scattering volume is always N s (z) and does not depend on N p .The scattering volume (length) for 1D flow simulations is an integral of the intensity PSF over all space [7,19] and equals to √ πw 0 /2.Since we confine our modeling to the x-dimension spanning from −L(z) to +L(z), as the particles move with the bulk velocity v(z) and exit the simulated space, they are reintroduced based on the periodic boundary conditions.Simulations using Eq. ( 27) are valid for arbitrary particle concentrations.For non-dilute particle suspensions with N s ≫ 1, the normalized flow autocovariance functions from Sec. 2.2 depend only on the beam waist in focus and are independent of w(z) and R(z).In this case Eq. ( 27) simplifies into with w 0 being the beam waist in focus.Therefore, the simulations of non-dilute flowing suspensions do not require a-priori knowledge of any other beam shape parameter except w 0 .
Even though the models in this section are given for the transverse flow, they can be readily extended to incorporate the flow component along the beam optical axis.

OCT system
The experiments were performed using a Thorlabs GANYMEDE II HR series spectral domain OCT system, which has been described in detail in our previous works [6,7].The acquisition rate was 5.5 kHz for diffusion and dilute flow measurements (low-speed), and 36 kHz for non-dilute flow measurements (high-speed).The acquired signal spectrum was measured with a spectrometer with 2048 pixels.The maximum imaging depth in air is 1.87 µm.After acquisition, the measured spectrum was first resampled to a linearly-sampled wavenumber domain and then apodized using a Gaussian filter.After the apodization, the measured coherence function waist in sample was w z = 2.11 µm.We have neglected the effect of a gradient of the axial velocity on the autocovariance function for two reasons [20].First, the Doppler angle in this work is essentially zero (θ ≈ 0).Second, our optical resolution is high both in axial and transverse directions compared to the flow channel dimensions.Hence, the flow velocity within PSF can be assumed to be constant.
The OCT system is operated with a scan lens (LSM04-BB, Thorlabs) in a confocal setup with a manufacturer provided focal spot size of w 0 = 6 µm in air which was validated by axial confocal response measurements [6,7].The system has NA = 0.05.Depending on the angle of incidence, refractive index contrast and Gaussian beam parameters, w 0 and w(z) vary somewhat because of the passage of the beam through the various interfaces [35].Therefore, for flow measurements w(z) and w 0 were calibrated using the procedures described in [6,7].Since for the given OCT setup the coherence length is small and the NA is very low, it can be assumed that the scattering angle is 180 • and the scattering wavenumber q in the correlation analysis is constant at q = 2nk 0 .

Averaging strategies
We acquire N b separate OCT interference time traces, S(k 0 , t m ), each having a length of M, to compute N b normalized autocovariance functions.Measurements conducted under identical conditions can be averaged.The processing and averaging steps of DLS-OCT are illustrated in Fig. 1.Non-linear curve fitting is utilized, truncating the autocorrelation functions when they reach zero and are fully decorrelated.From N b measurements, three distinct approaches are used to derive the average model parameters: • This first approach (method 1), represented by the orange arrow, involves fitting our models to each g 1 (z, τ m ) and g 2 (z, τ m ) and subsequently averaging the resulting N b parameters.We denote this method as D and v for diffusion and flow measurements, respectively, and it is also referred to as the standard method by us.
• The second approach (method 2), indicated by the red arrow, is to average N b normalized autocovariance functions and perform the curve fitting procedure on the averaged data.This approach is denoted as D g and v g for diffusion and flow measurements, respectively.
• The third approach (method 3), denoted by the blue arrow, is a mixing method developed by us.It begins by creating an N b × M matrix containing N b statistically independent autocorrelation functions along each row.Subsequently, the j th column is circularly shifted by an integer number j − 1, with the shifting being performed in the same direction for every column.This process generates a resultant 2D matrix where each row represents a mixed autocorrelation function containing only one correlation coefficient from any single time trace.Finally, our models are fitted to the N b mixed autocorrelation functions, and the obtained parameters are averaged.This method is denoted as mixed D and mixed v for diffusion and flow measurements, respectively.The mixing method is illustrated in Fig. 2 using simple autocorrelation functions with N b = M = 3.The superscript indicates the statistically independent measurements.
For a single time trace measurement the precision is given by the standard deviation (σ D and σ v ), but when averaging multiple measurements the precision is given by the standard error of the mean (σ average D and σ average v ).Averaging N b uncorrelated measurements lowers the variance by a factor of N b .For positively correlated measurements, the decrease is slower.
When averaging the fitted model parameters, σ 2 D and σ 2 v are reduced upon averaging, whereas when averaging the autocorrelation functions, σ 2 g 1 and σ 2 g 2 are scaled at all time delays and the noise is reduced.As Eq. (11,13,15,17) show, multiplying σ 2 g 1 or σ 2 g 2 by a constant N −1 b at every τ m results in the identical scaling of σ 2 D and σ 2 v .This holds true for both uncorrelated and correlated random variables.As the averaging sequence does not influence the noise correlation, for unbiased or even slightly biased estimators identical overall precision is expected for the first  For a single time trace measurement the precision is given by the standard deviation (  and   ),

319
but when averaging multiple measurements the precision is given by the standard error of the 320 mean ( average  and  average  ).Averaging   uncorrelated measurements lowers the variance 321 by a factor of   .For positively correlated measurements, the decrease is slower.

322
When averaging the fitted model parameters,   and second approaches.In the third approach, if the autocorrelation noise at different time delays is correlated, we anticipate that the mixing process will eliminate the correlation in the noise.
Bias cannot be improved by averaging.However, one of the sources of the bias is the estimation bias which, according to Sec. 2.4, exhibits some random nature.Hence, mixing or averaging the correlation functions can reduce this variability.Therefore, we study the bias for all three averaging schemes.The bias is generally very small compared to the measurement uncertainty, and due to a lack of a-priori information it can only be quantified using simulations.

Diffusion measurements
For diffusion measurements the acquisition rate was 5.5 kHz which is the highest sensitivity setting of our OCT system.All measurements are performed using monodisperse 50 nm radius silica particles with a volume fraction f v ≈ 1%, provided by CWK (Chemiewerk Bad Köstritz GmbH, Germany).In total 1100 depth-resolved measurements were performed with a time series length of 4096 points covering T = 0.74 s.From this data 1100 correlation functions g 1 (z, τ m ) and g 2 (z, τ m ) were calculated at every depth.For diffusion analysis only the real part of g 1 (z, τ m ) was used.Since the OCT sensitivity decreases in depth [36], a usable depth-range of 0.92 mm was chosen where the fitted diffusion coefficients were constant and where the single scattering regime holds.The reference diffusion coefficient was calculated by fitting Eq. ( 1) to the mean ℜ (︁ g 1 (z, τ m ) )︁ and averaging the resulting diffusion coefficients over the chosen depth range, resulting in D 0 = 4.02 ± 0.02 µm 2 /s which corresponded to a mean particle radius of 53 nm.The SNR at every depth was calculated based on Eq. (1) using the fitted A 1 (z).Simulations were performed using the same parameters as input.In total 10000 simulations were performed for each SNR.Simulated variances were used for calculating the Cramer-Rao lower bounds using Eq. ( 11) and Eq. ( 13).Precision as a function of SNR for both diffusion measurements and simulations was determined by computing the standard deviation and the standard error of the mean in the fitted diffusion coefficient both when using ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ).The bias for different averaging schemes was assessed using simulations.

Flow measurements
The flow was generated using a syringe pump with a variable discharge rate (Fusion 100, Chemyx, Inc.).The flow geometry was aligned so that θ ≈ 0. This minimized the velocity gradient effect on the correlation [20] and simplified the beam waist calibration.Other than that, the effect of a nonzero Doppler angle on DLS-OCT is minimal because of similar optical resolutions in axial and transverse directions.

Non-dilute flow measurements
For non-dilute flow measurements we used 20 mL syringe (Terumo Europe NV) and OCT acquisition rate of 36 kHz.Diluted Intralipid solutions were used as a sample with a particle volume fraction f v ≈ 0.6%.The flow passes through a quartz rectangular flow cell with internal dimensions of 0.2 mm depth and 10 mm width (type 45-F, Starna Scientific).M-scan 1D measurements were performed at the center width of the flow cell as a function of depth with discharge rates of 1.5 and 2 mL/min.The beam focus was moved away from the center depth of the sample to increase the number of particles in the scattering volume and reduce the effect of number fluctuations [18,19].The particle diffusion coefficient and the beam waist in focus were calibrated using a g 2 (z, τ m ) based on lateral scanning over the static sample [6] and were found to be D = 1.40 ± 0.09 µm 2 /s and w 0 = 8.08 ± 0.31 µm, respectively.The flow speed was determined by fitting Eq. (3,4) to the measured or simulated first and second-order autocovariance functions.The depth-dependent SNR was calculated from the fitted autocovariance magnitude.In total 1000 measurements and 10000 simulations were performed with a time series length of 4096.The average depth-dependent flow velocities and SNR values from measurements were used as the input for simulations.Simulated variances were used for calculating the Cramer-Rao lower bounds both for |g 1 (z, τ m )| and g 2 (z, τ m ).The bias for different averaging schemes was also determined using simulations.

Dilute flow measurements
Dilute measurements were performed with the second-order normalized autocovariance function using the number fluctuations analysis.The OCT acquisition rate was 5.5 kHz and a 5 mL syringe (BD Plastipak) was used.Experiments were performed using a monodisperse suspension of polystyrene particles with volume fraction f v ≈ 0.006% and expected particle radius of 230 − 250 nm, provided by InProcess-LSP.The flow passes through a quartz rectangular flow cell with internal dimensions of 0.5 mm depth and 10 mm width (type 45-F, Starna Scientific).M-scan 1D measurements were performed at the center width of the flow cell as a function of depth at discharge rates of 0.15 and 0.225 mL/min.In contrast to our previous work, where only the number fluctuation term was used [7], here we use the full model from Eq. ( 5).
The particle diffusion coefficient was determined (calibrated) in the stationary suspension using ℜ (︁ g 1 (z, τ m ) )︁ and Eq. ( 1).The obtained diffusion coefficient of D 0 = 1.00 ± 0.03 µm 2 /s corresponded to an average particle radius of 220 nm.The beam shape calibration procedure was identical to that in [7], but the analysis was slightly different because of using the full autocorrelation model.The SNR at every depth was calculated by fitting Eq. (3) to the measured |︁ |︁ g 1 (z, τ m ) |︁ |︁ .The kurtosis was determined using Eq. ( 8).Then A 3 (z) was found based on Eq. ( 7) using the known SNR and kurtosis.Finally, Eq. ( 6) was fitted to the measured g 2 (z, τ m ) with N s (z) and w(z) being the fit parameters and using A 3 (z), SNR, and the measured |︁ |︁ g 1 (z, τ m ) |︁ |︁ .The obtained beam shape w(z) was fitted using w(z) = w 0 √︂ 1 + (z − z 0 ) 2 /z 2 R , with w 0 , z 0 , and z R being the free parameters equaling 6.45 ± 0.02 µm, 187 ± 3 µm and 293 ± 3 µm, respectively.The particle volume fraction f v was computed using w(z), D 0 , w z and N s (z) according to [7].It is given in Fig. 3(a,b) along with the obtained w(z).
The obtained w z (z) and N s (z) were used in subsequent number fluctuation flow measurements.In total, 750 measurements were performed each with a time series length of 8192.The flow speed was determined by fitting Eq. ( 5) to the measured second-order autocovariance function with A 3 (z) and v(z) being the free parameters.Simulations for the dilute regime were not performed due to the adverse effects of boundary conditions and long computational times.Therefore, the measured variances were used for calculating the Cramer-Rao bounds.Even though number fluctuations are only present in g 2 (z, τ m ), for comparison the CRLB for |g 1 (z, τ m )| was also determined.

DLS-OCT diffusion measurement
Depth dependent diffusion measurements were performed on a static particle suspension.Examples of measured ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) for SNR = 100 are shown in Fig. 4(a).Observe here that the autocorrelation functions do not oscillate around the mean, but, instead, consistently remain above or below it for some time.From N b measured normalized autocovariance curves (signal ACF), the normalized autocovariances of errors (error ACF) in the measured and simulated and g 2 (z, τ m ) were calculated using where g(z, τ m ) is the measured normalized autocovariance function, g(z, τ m ) is its average from N b measurements, and σ 2 e (z) is the variance of g(z, τ m ) − g(z, τ m ) along the τ m axis.Figure 4(c,d) illustrates that in standard autocorrelation functions, errors at different time delays exhibit strong correlation.This means that the magnitude and direction of errors in or g 2 (z, τ m ) at small τ m significantly affect the errors at larger τ m , which can be observed from the fact that any single autocorrelation function is consistently above or below the mean autocorrelation function.As described in Sec.4.2, we implement a novel way of reducing these correlations or even completely eliminating them by mixing different autocorrelation functions at every τ m .The number of independent autocorrelation functions used for mixing is denoted by N mix .To fully get rid of the correlations N mix must be greater or equal to the number of sampling points it takes for ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) to go to zero.For g 2 (z, τ m ) it takes fewer realizations to eliminate correlations because it decays twice as fast.
Figure 4(b-d) show that for a single fully mixed ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) with sufficiently large N mix the errors at different τ m are completely uncorrelated.As a result, the measured mixed autocorrelation data oscillate approximately around the mean as a function of τ m .With standard (unmixed) measurements, as shown in Fig. 4(a), error ACF is nonzero and every measured and g 2 (z, τ m ) deviates from the mean autocorrelation function.Figure 4(c,d) show that once the autocorrelation functions are sufficiently mixed (e.g.N mix = 32), error ACF becomes a delta function.At intermediate mixing ratios (N mix = 4) error ACF becomes periodic as the correlations reappear at later time points.

Precision in diffusion estimation
Next, we look at the diffusion estimation precision from a single autocorrelation measurement.To estimate the precision, the autocorrelation variance at every time delay needs to be known.Figure 5(a,b) show measured, simulated, and analytical variances from Eq. ( 12), ( 14) for the diffusive ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) at different SNR values.All three match relatively well with each other.The variances are higher for larger SNR values, but, in this case the correlation coefficients are also larger.Even though the variances of ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) are comparable at larger time delays, σ 2 g 2 (z, τ m ) is significantly larger at small τ m .The different shape of σ 2 g 2 (z, τ m ) and its sharp increase at small τ m is due to the mean subtraction from the intensity time series.11) and ( 13, where the z-dependence is converted to SNR.The results are based on 10000 simulations and 1100 measurements.Here D 0 was calculated from the measured ℜ (︁ g 1 (z, τ m ) )︁ as described in Sec.4.3 and subsequently used as input for simulations.The obtained σ D strongly depends on whether the errors in our normalized autocovariance functions are correlated.The obtained σ D from the standard measurements is several factors larger than the Cramer-Rao lower bound.After mixing the resulting σ D matches very well with the Cramer-Rao bound and we achieve the most precise determination of the diffusion coefficient possible.The measured and simulated σ D overlap each other except for g 2 (z, τ m ) at very low SNR values.Here the simulations underestimates the error.The lowest spread in the fitted D, given by σ D , is obtained when mixing the normalized autocovariance functions for removing the error correlations.However, this requires many measurements, as shown in Fig. 6(a,b), and is therefore less relevant for practical applications with few measurements.Typically, the number of DLS-OCT measurements is on the order of 5-10.For conventional DLS-OCT measurements without mixing, σ D is constant as expected for independent measurements.When employing autocorrelation mixing we see that σ D decreases with increasing N b until reaching the CRLB and after which it remains constant with N b .For g 2 (z, τ m ) it takes fewer measurements for σ D to reach the CRLB.In DLS-OCT measurements are averaged to improve the precision.In this case the measurement precision is given by the standard error of the mean, σ average D . Figure 7(c for both averaging methods.However, when mixing the normalized autocovariance functions, we notice that as the fitted diffusion coefficients now become statistically interdependent.Here σ average D initially decreases very slowly with N b and then starts to reduce faster above a certain N b .This coincides with the N b at which the CRLB is reached in Fig. 7(a,b).In Fig. 7(d) there is a small discrepancy between measurements and simulations at large N b which caused by the lack of sufficient number of averaged measurements.

Bias in diffusion estimation
The bias in the average diffusion coefficient is determined by comparing the averaged D to a ground truth.We average 1100 and 10000 values of D in measurements and simulations, respectively.With this large set of averages we can assume that the random errors are negligible and only the systematic errors remain.Figure 8(a,b) shows the measured diffusion coefficients as a function of SNR for both ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) with different averaging schemes.Figure 8(c,d) show the same results but obtained using simulations.In this case the input D 0 is also displayed.Both in measurements and simulations we can observe that mixing different correlation functions slightly reduces the bias.However, the lowest bias is obtained for D ḡ when averaging the correlation functions.This is more evident for g 2 (z, τ m ).The differences for the three techniques are minimal for ℜ (︁ . Overall, the bias is much smaller and less affected by SNR when using ℜ (︁ The trends from both and simulations in Fig. 8(a,b) and Fig. 8(c,d) are similar, however for the measurements we lack a-priori knowledge of the ground truth D 0 .Therefore, it is not possible to compare the bias directly between measurements and simulations.This can be rectified on a relative basis by not comparing the absolute deviation from D 0 , but the relative bias with respect to the most accurate method (correlation averaging).Figure 9(a,b) displays the relative bias when averaging the fitted diffusion coefficients with respect to averaging ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) first and then fitting the diffusion coefficient.The obtained measurements and simulations are in perfect agreement except for g 2 (z, τ m ) at very low SNR values.This emphasizes the reliability of our simulations and that they can be used for determining the absolute bias.The bias of the mean diffusion coefficient from Fig. 8(c,d) is determined by averaging over a large number of simulations.Figure 10(a,b) show simulated systematic errors as a function of the number of averaged measurements the bias is very low and almost independent of N b , as expected.The bias is more significant for g 2 (z, τ m ). Figure 10(b) shows the estimation bias due to mean subtraction modeled using Eq. ( 22).This is the main source of bias for g 2 (z, τ m ) which is absent in ℜ . The estimation bias also contains a random component which averages out with increasing number of measurements.As a result, the bias in g 2 (z, τ) initially decreases with N b until reaching the minimum and then remains constant.Overall, the bias is minimal when averaging the autocorrelation functions.The second best technique is mixing the autocorrelation functions before the curve fitting and averaging the fitted diffusion coefficients.The worst approach is to the fit model parameters to the standard autocorrelation functions and then average them.

Non-dilute DLS-OCT flow measurement
Precision analysis of flow was performed for DLS-OCT |︁ |︁ g 1 (z, τ m ) |︁ |︁ and g 2 (z, τ m ) measurements of a laminar transverse flow corrected for diffusion.Figure 11(a,b) show measured and simulated variances at different flow speeds in the channel, good agreement between the two is observed.The variances are higher with lower flow speeds, but do not depend much on the SNR.The highest SNR value in the channel was above 100, and the lowest was 6.This is sufficiently high to neglect the SNR dependence of σ v for plotting purposes and display it only as a function of the flow speed.Velocity profiles, shown in Fig. 12(a), were obtained by averaging the measured correlation functions and were subsequently used as ground-truth input for simulations.Figure 12(c,d) show measured and simulated standard deviation in the fitted flow speed as a function of velocity.Results are obtained using both the standard method with correlated errors and the mixing technique with uncorrelated errors are given.Similar to diffusion, mixing decreases the spread in the obtained velocity.For comparison, the Cramer-Rao lower bounds for g 2 (z, τ m ) and |g 1 (z, τ)|, calculated using Eq. ( 15) and (21), are also displayed.The theoretical Cramer-Rao bound for |g 1 (z, τ)| is lower than for g 2 (z, τ m ).For both methods there is a good agreement between measurements and simulations.The obtained σ v from measurements and simulations is significantly improved when mixing removes the error correlations and matches well with the Cramer-Rao bound, especially for g 2 (z, τ m ).In general, the standard deviation in the flow speed decreases with lower velocities.However, it does not decrease to zero and has a certain minimum value.Below the threshold velocity σ v starts to increase because of diffusion.This is only clearly visible when the errors are not correlated.
Figure 12(b) displays the simulated bias for different averaging schemes both for |g 1 (z, τ)| and g 2 (z, τ m ).The bias was calculated with respect to the simulation input which included both depth-dependent velocity and SNR profiles.The simulated bias is lowest when averaging the complex g 1 (z, τ m ), taking the absolute value and then fitting v(z).This is denoted by the black curve.The second lowest bias is obtained when averaging g 2 (z, τ m ) and then fitting v(z) or when mixing g 2 (z, τ m ), fitting and averaging v(z).In this case the errors are identical and given by the red curve.This approach is only marginally less accurate than the former method.The third best method relies on using the standard g 2 (z, τ m ), fitting and then averaging v(z).This approach is less accurate than the second best method only at low flow speeds.The largest bias is obtained when using |g 1 (z, τ)| to fit v(z) and then average the obtained velocities.In this case the sample population is significantly skewed and the mixing method has virtually no effect.

Dilute DLS-OCT flow measurement
Similar to non-dilute, we analyse precision of dilute DLS-OCT flow measurements where the effect of diffusion is absent in g 2 (z, τ m ).Also for this case we obtained the variance of the correlation function for different flow speeds.The results look similar to the non-dilute case, except there is no overshoot at small τ m .Figure 13(a,b) show dilute flow measurements using the second-order normalized autocovariance function incorporating the number fluctuations.Velocity profiles are shown in Fig. 13(a) which were again obtained by averaging the measured correlation functions.The lowest SNR value in the channel was again around 6, and the highest was above 100.In Fig. 13 we have again neglected the SNR dependence for the same reasons as mentioned in Sec.5.2 and plot σ v only as a function of the velocity.17), and |g 1 (z, τ m )|, calculated with Eq. ( 21), are also displayed.Here we used the measured variances of σ 2 g 2 (z, τ m ) and σ 2 |g 1 | (z, τ m ) for calculating the Cramer-Rao lower bounds.The Cramer-Rao bound for g 2 (z, τ m ) is considerable lower than for |g 1 (z, τ m )| and agrees very well with our measurements.When using g 2 (z, τ m ) the velocity standard deviation decreases, with decreasing velocity, to zero.However, this is not the case with  3) is ultimately limited by particle diffusion and is independent of number fluctuations.The second-order normalized autocovariance function in Eq. ( 5), on the other hand, contains the number fluctuation term that is independent D. This increases the velocity precision for this method at extremely low flow speeds.For number fluctuation DLS-OCT no bias analysis is performed because of boundary condition limitations and increasing computational complexity of our simulations.

Discussion
In this work we review autocorrelation averaging and mixing techniques for reducing random and systematic errors.In the absence of an analytical model for error correlations, only using the mixing method suggested by us, it is possible to quantify the precision in DLS-OCT and verify whether the CRLB is reached.In diffusion measurements using the average real part of g 1 (z, τ m ) results in the highest accuracy and precision.In flow measurements it is more reliable and convenient to use the average g 2 (z, τ m ) because it is not affected by the sampling distribution bias, does not depend on phase and can be implemented in phase-unstable systems.It is important to note that mixing and/or averaging of the autocorrelation functions can only be performed using the autocorrelation functions with the same SNR (and identical optical or sample properties).For single scattering diffusion measurements where D is constant as a function of depth (SNR), we can also average the measurements with different signal-to-noise ratios.In this case the measured g 1 (z, τ m ) and g 2 (z, τ m ) first have to be noise-corrected [13] before any further processing.

DLS-OCT diffusion estimation
Random error correlations in ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) are the main reason for σ D not reaching the theoretical Cramer-Rao bound.Once these correlations are removed by mixing, the standard deviation in the fitted diffusion coefficient is reduced significantly and reaches the CRLB.The improvement is bigger for larger SNR values since ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) are lower in the presence of more noise and hence the random errors are less correlated.The improvement with mixing is less for g 2 (z, τ m ) because it decays faster and therefore correlations play less of a role in the precision.We can conclude that faster decorrelations and smaller autocovariance magnitudes, i.e., more noisy measurements, are less affected by mixing.We have not derived analytical expressions for ρ e , but our observations suggest that it depends on the autocorrelation decay rate and also experiences a delta function noise decorrelation at τ m = 0.
Even though the improvement in σ D due to mixing is significant, practical applications of the mixing technique for increasing the measurement precision are limited.In order to reduce some of the error correlations we need at least two measurements.Mixing just 2 measurements already provides us with a significant improvement in σ D .However, the precision of the average diffusion coefficient when using multiple measurements is given by the standard error of the mean, σ average D , which is identical for all averaging techniques with or without mixing the correlation functions.So, even though σ D is significantly lower when mixing the autocorrelation functions, it does not decrease as fast with averaging because the mixed autocorrelation functions are not statistically independent anymore.As a result, the fitted diffusion coefficients are also interdependent and because every fitted D is independent.With mixing we basically move the error correlations from ℜ (︁ g 1 (z, τ m ) )︁ and g 2 (z, τ m ) to the fitted D. Therefore, the actual measurement precision is identical for all averaging methods and cannot be improved any further.
We noticed that systematic errors are much larger when using g 2 (z, τ m ) compared to ℜ (︁ g 1 (z, τ m ) )︁ and are dominated by the estimation bias.Averaging N b measurements slightly reduces the bias in g 2 (z, τ m ) but only to a certain limit, beyond which the errors remain constant.This is probably caused by variability and randomness in the estimation bias which averages out and diminishes with increasing N b .We have mentioned that in phase-stable systems it is always preferable to utilize ℜ(g 1 (z, τ m )) instead of |g 1 (z, τ)|.Using |g 1 (z, τ)| adversely affects the bias because the autocorrelation coefficients are always positive.This leads to a higher sampling distribution bias which is otherwise negligible.Overall, random errors are larger than systematic errors unless we average thousands of measurements.
There is an asymptotic dependence of precision and bias on SNR.Therefore, the benefits of using extremely high signal-to-noise ratios are rather limited.As Fig. 6(b) and Fig. 9(b) show, there is a mismatch between simulations and measurements at very low SNR when using g 2 (z, τ m ).This is caused by the specific noise model used in our simulations [33] that is based on the assumption that the signal intensity is much larger than the noise intensity.However, this assumption is not valid at low SNR values.As a result, this noise model overestimates A 2 (z) at low SNR values, leading to the lower σ D .

DLS-OCT flow estimation
Both in dilute and non-dilute flow measurements, the standard deviation (spread) in the fitted velocity is significantly reduced when the error correlations are removed.In non-dilute flows, σ v decreases with decreasing velocity and then starts to increase at very low flow speeds.A minimum in the error occurs because the diffusive limit is reached.At this stage, any further reduction in the velocity does not affect autocorrelation functions or their variances.According to Eq. ( 15), the velocity is explicitly included in the denominator of σ 2 CRLB, v .As a result, σ v starts to increase when the diffusive limit is reached.This behaviour is overshadowed by error correlations and is only visible when employing correlation function mixing.
For non-dilute flow the Cramer-Rao bound for g 2 (z, τ) is slightly lower than the simulated σ v at larger velocities, which itself is marginally lower than the measured σ v .The difference between simulations and the Cramer-Rao bound is probably because of the Siegert approximation.In this case, the approximate number of particles in the scattering volume is around 49 both in simulations and measurements.This could be sufficiently low to violate the large number of particle assumption and have a small effect on g 2 (z, τ m ).Furthermore, the offset between measurements and simulations can be caused by galvo or pump instabilities.The Cramer-Rao bound of |g 1 (z, τ)| shows that on paper it is possible to achieve higher precision in non-dilute flow measurements when using g 1 (z, τ) instead of g 2 (z, τ).However, in practice this is more problematic.The mismatch among measurements, simulations and the Cramer-Rao bound at all velocities is considerably larger for |g 1 (z, τ)| compared to g 2 (z, τ).This is largely caused by the bias of the sampling distribution and the limited phase stability.
Theoretically, the lowest bias in non-dilute flow measurements can also be achieved when using |g 1 (z, τ)| by averaging the complex autocorrelation functions first, then taking the absolute value, and then fitting v(z).However, getting rid of the sampling distribution bias when taking the absolute value requires sufficient averaging of the complex g 1 (z, τ m ).In our work we averaged 1000 autocorrelation functions which in real-time flow measurements is impossible.Insufficient averaging will result in oscillations with |︁ |︁ g 1 (z, τ m ) |︁ |︁ >0.In this work we have not investigated the dependence of the sampling distribution bias on the number of measurements.However, it is clear that with fewer measurements it is preferred to use g 2 (z, τ) as it does not suffer from the population skewness.
Similar to diffusion, random velocity errors in flow measurements are also generally larger than systematic errors.However, at very low speeds both errors become comparable.Despite some advantages of |︁ |︁ g 1 (z, τ m ) |︁ |︁ , for practical reasons it is more convenient to use g 2 (z, τ m ) in non-dilute flows.In dilute low-speed flows it is always preferable to use g 2 (z, τ m ) because it is not limited by the particle diffusion.

Conclusion
We have investigated precision and bias of diffusion coefficient and flow speed measurements using DLS-OCT based on the first and second-order normalized autocovariance functions.We found that errors in the autocovariance functions are strongly correlated.This significantly reduces the precision and bias of fit parameters and prevents us from reaching the Cramer-Rao lower bound.We demonstrated that mixing different autocovariance functions at every time delay before the curve fitting procedure reduces the standard deviation in the fitted parameters and reaches the Cramer-Rao lower bound.This proves the validity of our DLS-OCT models.
When using the mixing technique the correlations are transferred to the fitted parameters and the standard error of the mean remains unchanged.We conclude that the precision in DLS-OCT is identical for all averaging techniques and ultimately limited by correlations between the random variables.The bias is lowest when averaging the measured normalized autocovariance functions before fitting the model parameters.

Fig. 1 .
Fig. 1.Data processing steps for obtaining D and v from OCT data with different averaging techniques.

Fig. 2 .
Fig. 2. Schematic overview of the calculation of standard and mixed autocorrelation functions.

Fig. 3 .
Fig. 3. Obtained (a) beam shape w(z) and (b) particle volume fraction f v as a function of depth.

Figure 6 (
Figure6(a,b) show measured and simulated standard deviations in the fitted diffusion coefficient from a single correlation function, as well as the Cramer-Rao lower bound calculated with Eq. (11) and (13, where the z-dependence is converted to SNR.The results are based on 10000 simulations and 1100 measurements.Here D 0 was calculated from the measured ℜ (︁ g 1 (z, τ m ) )︁ as described in Sec.4.3 and subsequently used as input for simulations.The obtained σ D strongly depends on whether the errors in our normalized autocovariance functions are correlated.The obtained σ D from the standard measurements is several factors larger than the Cramer-Rao lower bound.After mixing the resulting σ D matches very well with the Cramer-Rao bound and we achieve the most precise determination of the diffusion coefficient possible.The measured and simulated σ D overlap each other except for g 2 (z, τ m ) at very low SNR values.Here the simulations

Fig. 6 .
Fig. 6.Measured (points) and simulated (lines) standard deviation for estimating the diffusion coefficient from a single autocorrelation function and the corresponding Cramer-Rao lower bound for (a) g 1 (z, τ m ) and (b) g 2 (z, τ m ).

Figure 7 (
Figure 7(a,b) show the standard deviation (spread) in the diffusion coefficient obtained from fitting a single correlation function as a function of the number of measurements N b .For conventional DLS-OCT measurements without mixing, σ D is constant as expected for independent measurements.When employing autocorrelation mixing we see that σ D decreases with increasing N b until reaching the CRLB and after which it remains constant with N b .For g 2 (z, τ m ) it takes fewer measurements for σ D to reach the CRLB.

Fig. 7 .
Fig. 7. (a,b) Standard deviation and (c,d) standard error of the fitted diffusion coefficient as a function of the number of averaged correlation functions for SNR = 100.Points correspond to DLS-OCT measurements and lines denote simulations.
) shows simulated and measured standard error of the mean diffusion coefficient, σ average D , as a function of N b for different averaging techniques.The standard error of the mean is identical for all three averaging methods, even though σ D itself is different depending on error correlations.The decrease is proportional to N b −1/2 .Figure 7(d) shows how σ average D is related to σ D for the three different averaging techniques.For the standard analysis we see that σ average D = σ D √ N b

Fig. 8 .
Fig. 8. Diffusion coefficients as a function of SNR and averaging scheme determined using (a,b) measurements and (c,d) simulations.

Fig. 10 .
Fig. 10.Simulated bias of the fitted diffusion coefficient as a function of the number of used correlation functions at SNR = 100 for (a) ℜ (︁ g 1 (z, τ m ) )︁ and (b) g 2 (z, τ m ).

Figure 13 (
Figure13(b)  shows the measured standard deviation in the fitted flow speed as a function of velocity.Results using both the standard g 2 (z, τ m ) with correlated errors and the mixed g 2 (z, τ m ) with uncorrelated errors are given.Measurements using |︁ |︁ g 1 (z, τ m ) |︁ |︁ are not shown because they are too noisy and have a significantly larger σ v .For comparison, the Cramer-Rao lower bounds for g 2 (z, τ m ), calculated with Eq. (17), and |g 1 (z, τ m )|, calculated with Eq. (21), are also displayed.Here we used the measured variances of σ 2 g 2 (z, τ m ) and σ 2 |g 1 | (z, τ m ) for calculating the Cramer-Rao lower bounds.The Cramer-Rao bound for g 2 (z, τ m ) is considerable lower than for |g 1 (z, τ m )| and agrees very well with our measurements.When using g 2 (z, τ m ) the velocity standard deviation decreases, with decreasing velocity, to zero.However, this is not the case with|︁ |︁ g 1 (z, τ m ) |︁ |︁ .As the Cramer-Rao lower bound shows, σ v increases dramatically at lower flower speeds when using|g 1 (z, τ m )|.This happens because |︁ |︁ g 1 (z, τ m ) |︁ |︁ in Eq. (3) is ultimately limited by particle diffusion and is independent of number fluctuations.The second-order normalized autocovariance function in Eq. (5), on the other hand, contains the number fluctuation term that is independent D. This increases the velocity precision for this method at extremely low flow speeds.For number fluctuation DLS-OCT no bias analysis is performed because of boundary condition limitations and increasing computational complexity of our simulations.

|︁ |︁ g 1
Figure13(b)  shows the measured standard deviation in the fitted flow speed as a function of velocity.Results using both the standard g 2 (z, τ m ) with correlated errors and the mixed g 2 (z, τ m ) with uncorrelated errors are given.Measurements using |︁ |︁ g 1 (z, τ m ) |︁ |︁ are not shown because they are too noisy and have a significantly larger σ v .For comparison, the Cramer-Rao lower bounds for g 2 (z, τ m ), calculated with Eq. (17), and |g 1 (z, τ m )|, calculated with Eq. (21), are also displayed.Here we used the measured variances of σ 2 g 2 (z, τ m ) and σ 2 |g 1 | (z, τ m ) for calculating the Cramer-Rao lower bounds.The Cramer-Rao bound for g 2 (z, τ m ) is considerable lower than for |g 1 (z, τ m )| and agrees very well with our measurements.When using g 2 (z, τ m ) the velocity standard deviation decreases, with decreasing velocity, to zero.However, this is not the case with|︁ |︁ g 1 (z, τ m ) |︁ |︁ .As the Cramer-Rao lower bound shows, σ v increases dramatically at lower flower speeds when using|g 1 (z, τ m )|.This happens because |︁ |︁ g 1 (z, τ m ) |︁ |︁ in Eq. (3) is ultimately limited by particle diffusion and is independent of number fluctuations.The second-order normalized autocovariance function in Eq. (5), on the other hand, contains the number fluctuation term that is independent D. This increases the velocity precision for this method at extremely low flow speeds.For number fluctuation DLS-OCT no bias analysis is performed because of boundary condition limitations and increasing computational complexity of our simulations.
2  and  2  are reduced upon averaging, whereas 323 when averaging the autocorrelation functions,  2  1 and  2  2 are scaled at all time delays and the 324 noise is reduced.As Eq. (11,13,15,17) show, multiplying  2  1 or  2  2 by a constant  −1  at every 325   results in the identical scaling of  2  and  2  .This holds true for both uncorrelated and 326 correlated random variables.As the averaging sequence does not influence the noise correlation, 327 for unbiased or even slightly biased estimators identical overall precision is expected for the first 328