Cross-correlation analysis of X-ray photon correlation spectroscopy to extract rotational diffusion coefficients.

Coefficients for translational and rotational diffusion characterize the Brownian motion of particles. Emerging X-ray photon correlation spectroscopy (XPCS) experiments probe a broad range of length scales and time scales and are well-suited for investigation of Brownian motion. While methods for estimating the translational diffusion coefficients from XPCS are well-developed, there are no algorithms for measuring the rotational diffusion coefficients based on XPCS, even though the required raw data are accessible from such experiments. In this paper, we propose angular-temporal cross-correlation analysis of XPCS data and show that this information can be used to design a numerical algorithm (Multi-Tiered Estimation for Correlation Spectroscopy [MTECS]) for predicting the rotational diffusion coefficient utilizing the cross-correlation: This approach is applicable to other wavelengths beyond this regime. We verify the accuracy of this algorithmic approach across a range of simulated data.

The analysis of Brownian motion of different types of particles is a classic problem in research fields such as molecular biology and materials science. For a dilute suspension, collisions from the solvent particles lead to random reposition and reorientation that can be decomposed into translational and rotational diffusion. These two types of diffusion can be characterized by coefficients, namely, the translational diffusion coefficient and rotational diffusion coefficient in two dimensions and the corresponding tensor in three dimensions. Knowledge of these parameters provides insight into the structure and dynamic properties of the particles, opening the door for understanding functions and transport process of proteins (1), synthesis and stability of materials (2), biomolecular reactions (3), etc.

As an emerging X-ray scattering technique, X-ray photon correlation spectroscopy (XPCS) is able to probe length scale down to nanometers and time scale from below microseconds to hours and is well-suited for studying dynamics of disordered systems, of which analyzing diffusion is a critical informative characteristic. In XPCS experiments, samples are illuminated by partially coherent X-ray beams, and time series of scattering images are collected by detectors. (See Fig. 1 for an illustration of the XPCS experiments.) The inhomogeneity of the electron density leads to spatial variation of the brightness of the images, which is called a ”speckle pattern,” and the fluctuation of the electron density leads to the temporal variation of the speckle patterns that contains valuable information of the sample dynamics, which can be revealed by analyzing these collected, time-dependent images.

Fig. 1.

Schematic illustration of the XPCS experiments. The translation and rotation of the particles within the scattering volume leads to variation of the speckle patterns shown on the right. (While the grainy, noise-like texture makes these images appear visually similar, the proposed algorithm is able to detect and analyze the contained variations.)

One of the most well-known tools for analyzing the images is the temporal autocorrelation function , which depends on the second-order degree of coherence (4); see, for example, refs. 5–7. While translational diffusion coefficients can be determined through , to the best of our knowledge, there is no current algorithm to extract from XPCS data.

The objective of this paper is to present a methodology for calculating the rotational diffusion coefficient in the two-dimensional (2D) case using XPCS data. By exploiting the more detailed angular-temporal cross-correlation function, we are able to discover more information than the autocorrelation function . Information about this cross-correlation function can be related to the static electron density of individual particles and the rotational diffusion coefficient . To estimate from the proposed cross-correlation, we first introduce a mathematical model and related set of equations, whose solution corresponds to the rotational diffusion. We then design a numerical algorithm, “Multi-Tiered Estimation for Correlation Spectroscopy (MTECS),” based on an approach introduced here, which solves the relevant equations by following a multitiered iterative projection (M-TIP) philosophy, first introduced in ref. 8. In particular, we construct several operators and apply them in an iteration to efficiently map the cross-correlation data into a form that can then be processed by the algorithm to converge to the solution of the underlying equations.

While we mainly discuss XPCS in this paper, the developed methodology does not rely upon the specific range of the wavelength essentially, so that it can be generalized to other experiment techniques, such as dynamic light scattering. Though dynamic light scattering has been applied to measure the rotational diffusion (9, 10), these methods are limited by their requirements of the type and structure of the particles. The MTECS algorithm requires almost no prior assumption about particle structure, but can take advantage of such information if a priori known. MTECS is robust against noisy data, as it includes additional filtering methods applied to the input cross-correlation data.

Background

Assume a monodisperse dilute sample of particles undergoing free translational and rotational diffusion. Here, we focus on the 2D case, where the rotational diffusion is around a single axis parallel to the incident beam. The angular perturbation of particles can be modeled by Gaussian variables (11),where is the angle indicating the orientation at time , and is the desired rotational diffusion coefficient.

Let and be the Cartesian coordinates in real and Fourier space, respectively, and in polar coordinates, they can be expressed as and , where are the radial coordinates, and are the angular coordinates. Let be the electron density of an individual particle whose center of mass is at origin and orientation is associated with angle zero. can be expressed in circular harmonic expansion asand the circular harmonic coefficients can be calculated by the circular harmonic transformThe form factor is the Fourier transform of :and the intensity is the squared modulus of the form factor :Similarly to Eq. , and can both be represented in circular harmonic expansion with coefficients , , respectively. Since the electron density is real, application of Friedel’s law means that for odd .

Let denote the scattering data, which is the squared modulus of the Fourier transform of the electron density within the scattering volume, recorded by the detectors at scattering vector and time . Throughout this paper, we consider only nonzero scattering vectors. Assuming a flat Ewald sphere (12), the scattering vector can be reduced to its component orthogonal to the incident beam. Then, we represent the images in polar coordinate system as . The temporal autocorrelation function is defined as follows:where the angle bracket indicates the ensemble average that is equivalent to the time average taken over for systems under thermal equilibrium. For isotropic systems, does not depend on in the limit of infinite number of images (11). In this case, the averages are also taken over , and the dependency of on can be dropped, thus becoming .

For spherical particles undergoing free diffusion (13),* and, hence, estimation of can be obtained according to this equation. For nonspherical particles, the particle structure plays an important role. Assuming anistropic structure of particles and decoupling of translation diffusion and rotational diffusion, we are able to show () thatin the limit of infinite number of images. Other than and , the above equation involves also the unknowns, . However, there may be insufficient information in Eq. to extract all of these quantities, which leads us to the analysis introduced in the next section.

Angular-Temporal Cross-Correlation

In the cases that the X-ray exposure time is below the rotational diffusion relaxation time, which can be achieved by slowing down the motion or with ultrafast X-ray light sources, the angular fluctuation on images can be recorded (14, 15), providing key information that can be used to estimate the rotational diffusion coefficient . However, this information is missing from the autocorrelation function due to the fact that it is an average over the system orientations for an isotropic system. More specifically, since correlates points in the same spatial position from different images, it is not able to capture information from photons that are scattered along different scattering vectors and correlated by rotational dynamics or symmetric structures of particles.

To address this limitation, we propose exploiting the angular-temporal cross-correlation function, defined aswhere the average is performed over and . This function can be viewed as a natural generalization of the autocorrelation . Taking and , the cross-correlation becomes the autocorrelation without the normalization, which implies that data given by are a subset of data given by . Moreover, by exploiting correlation between points at different spatial positions, is able to uncover information contained in the angular fluctuation of images to produce orders of magnitude more information than . As the arguments vary, generates only a data matrix, while the more involved generates a four-way data tensor. We note that in a non-XPCS context, reduced forms of Eq. with consideration of only angular effect have been applied to discovery of local symmetries (16, 17), ab initio structure determination (18, 19), and investigation of orientational distribution (20–22).

In practice, the cross-correlation is evaluated aswhere is the time difference between consecutive images, and is the number of images collected. Assuming or , in the limit of infinite number of images, i.e., , while keeping constant, we can derive the following limit:up to a constant. When or and , the cross-correlation manifests as sharp peaks (see Fig. 5). The above limit and the peaks are derived in . The second term on the left-hand side of Eq. is the angular average of the cross-correlation. The subtraction of this term can be replaced by the techniques outlined in ref. 23 to reduce the effect of background noise and detector malfunction.

Fig. 5.

An example of the cross-correlation data at with against angular coordinates ranging from . The data on and the data on are identical. The vertical black bars indicate the boundaries of the masking that masks out the peaks at and .

Intuitively, on the assumption that the motions of different particles are uncorrelated, the interparticle effects vanish in the ensemble averages. The main contribution arises from the effect of each individual particles. The rotational diffusion results in the decaying terms , which eventually lead to the decorrelation of the angular fluctuation.

In comparison to Eq. , the rates of the multicomponent exponential decay in Eq. depend solely on . Moreover, Eq. provides considerably more data to help handle the unknowns related to the structure of the particle.

MTECS Algorithm

Our focus is on anisotropic particles (otherwise, the right-hand side of Eq. is always zero). The rotational diffusion coefficient can be estimated by solving Eq. , which we do by building an algorithm, called MTECS, following the M-TIP approach of decomposing complex optimization into subparts that can be efficiently inverted/pseudoinverted or solved exactly and combining the solutions from the different subparts in an iteration (8, 24). This approach allows us to avoid solving a high-dimensional nonconvex problem directly and, instead, decompose the original problem into a few efficiently solvable subparts, narrowing the nonconvex part to a much lower-dimensional space.

Let be an upper bound of the diameter of the minimal circle enclosing the particle. This is the only prior knowledge required by MTECS. We will discuss this parameter in later subsections. Define as the right-hand side of Eq. truncated to a finite sum,Assume that the above function is measured at , , and , where is the number of radial coordinates, is the number of the angular coordinates, and is the number of lag time. Denote the maximum of as . These values give rise to a four-way data tensor , i.e., , whose values form the input to the algorithm. Here, we use the similar notations of tensors as ref. 25. The element of the four-way tensor is denoted by . The subarrays of are specified by using colons in subscripts; e.g., indicates a vector obtained by varying the third subscript, and represents a matrix obtained by varying the first two subscripts. Notice that the parameter should be chosen so that . A good heuristic for is approximately , which is computed by determining the number of Shannon elements that lie within the domain of the polar grid (26).

The flowchart in Fig. 2 outlines the procedure of the MTECS algorithm. Now, we turn to showing how to construct the subparts, each of which enforces a constraint induced by the mathematical description of the cross-correlation .

Fig. 2.

High-level flowchart of the MTECS algorithm. See the subsections Correlation Noise Projector, Band-Limiting Projector, Rank-One Tensor Decomposition, Exponential Fitting, and MTECS for more detailed description.

Correlation Noise Projector.

Let be a four-way tensor, whose entries areWe first design a correlation noise projector to seek the minimum perturbation of the current guess of to make the updated cross-correlation within a range from the input data.

One may notice that Eq. is a circular harmonic expansion. However, the coefficients cannot be attained by the circular harmonic transform Eq. due to the fact that we usually don’t have access to the full angular range in practice. This is caused by either detector artifacts or the peaks of the cross-correlation at or , .

Instead, for each , Eq. can be rewritten as a matrix-vector multiplication:whereSince and are known and , the above problem is an overdetermined linear system.

To solve Eq. , we formulate as a Tikhonov regularization, i.e., an -norm constrained linear regression:where for each , solveswhere are scaling parameters, and is a fitting parameter of the constraint. These parameters can be chosen based on how well the cross-correlation data fit the model Eq. (see for more details).

Band-Limiting Projector.

Starting with the Fourier–Bessel expansion of the density–density autocorrelation function (), we derive an basis expansion of the circular harmonic coefficients of the intensity ,where are Fourier–Bessel series coefficients of , andwhere are Bessel functions of the first kind, and are the -th zeros of . We don’t need to be precisely determined: It is only required to be larger than the true value of the diameter to ensure the convergence of Eq. . A special property of the set of basis functions is that its mass concentrates around , which implies that it is practicable to take truncation of Eq. based on the number of terms whose majority of mass is within the measured range. Here, we assume that the series can be accurately truncated to the first terms.

Combining Eqs. and , we havewhereSince the matrices as well as their pseudoinverses can be computed explicitly, we define the band-limiting projector to solve Eq. for as

Rank-One Tensor Decomposition.

Next, we obtain the coefficients and the term . According to Eq. , for each , the tensor can be viewed as an outer product:whereThe above equation implies that is at most rank one. Relying on this, and can be retrieved by using a partially Hermitian bound-constrained rank-one tensor decomposition:where are weights, and is the Frobenius norm. The purpose of the weights is to balance the exponentially decaying signal and the noise, whose magnitude remains the same with respect to different . If or is quite large, the unweighted Eq. tends to pick up the noise instead of true signal. A reasonable choice of is based on the decaying terms , where is the current estimation. This subpart can be expressed briefly as

Exponential Fitting.

Once the matrix of decaying terms is computed, we are able to estimate the rotational diffusion coefficient by solving the following minization:and then lettingIf the magnitudes of are too close to zero for some , which happens when the particles possess unknown radial symmetry, then could be too noisy for these . To deal with this situation, the algorithm allows assigning the weights according to the estimated magnitude of , e.g., . This subpart is denoted by

MTECS.

We now assemble the subparts described in the previous subsections to give an iterative algorithm, MTECS (see Fig. 3 for the flowchart).

Fig. 3.

Low-level flowchart of the MTECS algorithm. The blue arrows represent the subparts that can be computed exactly, and the red arrows stand for the subparts whose solutions can only be approximated iteratively.

With initialization and , the algorithm consists of the following steps:

Given the cross-correlation , apply the correlation noise projector to compute ;

For , compute by performing the band-limiting projector for ;

Compute by the rank-one tensor decomposition for ;

Make a current estimation of the rotational diffusion coefficient as ;

Form an updated matrix of the decay terms whose entries are for and ;

For , obtain the tensors ;

Update the tensor by calculating for and ;

Terminate if the difference between and is smaller than a preset tolerance; otherwise, set and go to 1.

The stopping criterion is that the relative change between and is smaller than a tolerance rather than hinged on the desired coefficient , as the former choice is more robust. (See for more details of selection of parameters.)

The upper bound of the diameter of the minimal bounding circle of the particle required by step 2 is the only necessary prior knowledge. When available, MTECS can be modified to capitalize on additional prior information of particle structure. Some additional constraints on the intensity and its circular harmonic coefficients can be derived accordingly. For example, if the particles have -fold symmetry, then for that is not a multiple of . Then, Eqs. and can be adapted accordingly.

Analogous to the methods introduced in ref. 24, a filtered cross-correlation tensor, , is acquired through Eq. using the last . Against the input , the signal-to-noise ratio of is significantly enhanced to benefit the aforementioned angular cross-correlation analysis for other X-ray diffraction techniques, e.g., structure reconstruction (8, 27).

Results

To illustrate the capabilities of this approach, we applied MTECS to synthesized images collected from simulated systems. All the particles of interest in the systems are assumed to be identical. Fig. 4 shows the electron density . Each particle is constructed by attaching two spheres to ends of a cylinder in three-dimensional (3D) spaces, then projected to 2D space by integrating along the axis. (We are currently exploring extending these techniques to systems with polydispersity.) The analytical expressions of and are provided in . We calculate the circular harmonic coefficients by numerical integration.

Fig. 4.

(A) The electron density of an individual particle, . (B) An example of generated XPCS images. The colors are mapped according to the logarithm of the collected intensity data. The shaded area enclosed by the two circles describes the domain of the polar grid on which the intensity is used. (C) Logarithm of angular average of the image shown in B versus the radial coordinate .

The systems are simulated on a 2,047 × 2,047 Cartesian grid equally spaced in domain , with 500 particles in each system with uniform distribution in initial position and orientation. The position replacements and orientation perturbations are treated as Gaussian random variables. The translation diffusion coefficient is set to be D = 1/2,0472, while the rotational diffusion coefficients are selected differently for different systems. Overlapping and collisions are ignored. The scattering fields are calculated via the 2D fast Fourier transform. Instead of the whole images on a Cartesian grid, we interpolate them to a polar grid with radial coordinates spaced equally on and 3,210 angular coordinates spaced equally on . (See Fig. 4 for an example of the images whose angular average against the radial coordinate is plotted in Fig. 4.) The time difference between consecutive images is 1.0.

For each system, we let , and calculated a cross-correlation data tensor . (See for more details.) The calculation of a such tensor from images took 4 h using 32 single-socket Intel Xeon Phi 7250 processors with 68 cores at 1.4 GHz. A curve of cross-correlation data against is shown in Fig. 5.

There are two sharp peaks at and . Additionally, the data on are identical to the data on due to symmetry of the intensity . Thus, we kept the data on , ending up with a data tensor .

We used in the truncation of the cross-correlation Eq. and the upper bound of the diameter of the minimal circle enclosing the particle that is 10% larger than the actual diameter. Such precision can be easily achieved in practice. The first subpart, Tikhonov regularization , was solved by finding the Tikhonov parameter (28). We solved the tensor decomposition through a modified alternating least-square algorithm (29). The exponential fitting was handled by limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (30). Please refer to for more details about solving the subparts. The MTECS algorithm terminates if the convergence is reached or after 200 iterations. It took at most 20 min to complete the algorithms on a single core of a 2.3-GHz 16-core Intel Xeon Processor E5-2698 v3.

For each of the rotational diffusion coefficient values ranging from 0.01 to 0.1, we collected images. Table 1 gives the ground truth, the estimated given by the algorithm and the relative error between them. In Fig. 6, we present two curves of the filtered cross-correlation data against the angular coordinate . The input data and ground truth calculated from the circular harmonic coefficients via Eq. are also shown for comparison.

Table 1.

Ground truth of , estimation of , and the relative errors

True D_r	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.10
Estimated D_r	0.010062	0.020086	0.030490	0.040521	0.050919	0.061311	0.071812	0.080488	0.092633	0.100093
Relative error, %	0.62	0.43	1.63	1.30	1.84	2.19	2.59	0.61	2.93	0.09

Fig. 6.

Filtered cross-correlation data at against angular coordinates. The cyan dots are the cross-correlation data computed from images, which are the input to MTECS. The black solid lines are the cross-correlation filtered by the algorithm. Ground truth calculated from is indicated by the red dashed lines.

As shown in Table 1, all the rotational diffusion coefficients are estimated with relative errors smaller than . Fig. 6 shows that the curves of filtered cross-correlation data almost completely overlap with the ground truth, indicating the ability of the MTECS algorithm to filter the cross-correlation data and, hence, boost the signal-to-noise ratio.

To understand how the computed accuracy depends on the number of images, for , we gradually reduced the number of images from to . As less numbers of the images are used, there will be more noise caused by the random motion of the particles. Due to the stochastic nature of the experiments, which is also true in practice, the outputs exhibit some randomness.

To characterize the results better, for each value of and each number of images, we repeated the experiments 10 times and, hence, generated 10 different estimations, whose statistics are shown in Fig. 7. We note that the MTECS algorithm is still able to give reasonable estimation while using less images. Despite the fluctuation, the trend shows that as the number of images increases, the quality of the estimations improves.

Fig. 7.

Quantiles of 10 different estimated values against number of images. Horizontal bars from top to bottom are maxima, third quartiles, medians, first quartiles, and minima. White boxes are averages. Horizontal lines behind boxes are ground truth.

We also conducted experiments in which we added Poisson noise to the images. Rather than the polar grid introduced before, we used a coarsened polar grid with only half amount of the radial coordinates and half amount of the angular coordinates. The domain of this grid is the same as before. is selected as 0.01, and the other parameters remain unchanged. A total of images were collected. We multiplied the intensity by a scaling factor before added Poisson noise. The level of Poisson noise can be controlled by tuning this factor, and, thus, the Poisson noise has a more severe effect when the magnitude of data are small. Results are shown in Table 2: The relative errors are all smaller than . Such results are quite close to the estimation in Table 1.

Table 2.

Estimation of and the relative errors against the factor controlling the level of Poisson noises

Factor	0.25	0.5	1.0	2.0	4.0
Estimated D_r	0.010088	0.010081	0.010093	0.010093	0.010091
Relative error, %	0.88	0.81	0.93	0.93	0.91

The smaller the factor, the more severe the noise. Ground truth is

Conclusions and Discussion

We have shown how to estimate rotational diffusion coefficients from XPCS. The cross-correlation provides significant information, which cannot be acquired through the autocorrelation function . In fact, viewed as a natural generalization of , the cross-correlation function has the potential of being applied to various other problems, despite the fact that it takes substantially more effort and resources to process and interpret the cross-correlation data, owing to the size of data.

The derivation of Eq. is the basis of application of the algorithm. Consequently, the outlined approach of the combination of cross-correlation functions and MTECS applies to other wavelengths and, hence, can be carried directly to other experimental regimes.

Moreover, the approach developed here provides a reliable methodology to measure the rotational dynamics beyond rotational diffusion. For noninteracting particles, one can derive a limit similar to Eq. for other dynamics. The variations of the circular harmonic coefficients of cross-correlation curves can be computed through the first three subparts of the MTECS algorithm and be studied by other operators that replace the last subpart, the exponential fitting . For interacting particles, e.g., in concentrated systems, the ensemble average of the cross-correlation contains terms that are related to the pair distribution of the particles. These terms can still be expressed in terms of product of the circular harmonic coefficients of or . Additional constraints should be explored to design operators handling these extra terms.

The current MTECS algorithm is designed for estimating the rotational diffusion coefficient in 2D rotational diffusion. A natural next step is to estimate the rotational diffusion tensor in 3D rotational diffusion. We are currently working on this extension: A limit similar to Eq. can be derived, and then the subparts can be modified accordingly. The overall procedure remains basically the same.