Probing roto-translational diffusion of small anisotropic colloidal particles with a bright-field microscope

Abstract Soft and biological materials are often composed of elementary constituents exhibiting an incessant roto-translational motion at the microscopic scale. Tracking this motion with a bright-field microscope becomes increasingly challenging when the particle size becomes smaller than the microscope resolution, a case which is frequently encountered. Here we demonstrate squared-gradient differential dynamic microscopy (SG-DDM) as a tool to successfully use bright-field microscopy to extract the roto-translational dynamics of small anisotropic colloidal particles, whose rotational motion cannot be tracked accurately in direct space. We provide analytical justification and experimental demonstration of the method by successful application to an aqueous suspension of peanut-shaped particles. Graphic abstract


Introduction
Over a century ago, Jean Perrin conducted a series of experiments to mark the birth of quantitative microscopy [1,2]. Among his many intriguing results, the simultaneous measurement of the translational and rotational diffusion coefficients of micron-sized particles provided the first convincing experimental demonstration of the energy equipartition between the translational and rotational degrees of freedom, one of the most important predictions of the kinetic theory. Since then, the rotational and the roto-translational motion of colloidal objects has been repeatedly used as the tool of choice to investigate fundamental physical processes like mass transport or phase transitions [3,4], as a probe of material properties at the micro-scale [5,6], as a biosensor [7], to monitor the motility pattern of active swimmers and their interactions with passive particles [8,9], to characterize complex transport dynamics in cells [10], just to mention a few recent examples.
Nowadays, optical microscopy remains a key tool for probing roto-translational particle motion, especially in combination with specific (for example fluorescent) markers. For low particle density, single particle tracking (SPT) represents a rather straightforward way to extract a quantitative information on the translational dynamics of well-resolved particles by following their center of mass; by contrast, the accurate determination of their three-dimensional rotational dynamics with SPT is usually much more challenging and time cona e-mail: fabio.giavazzi@unimi.it (corresponding author) suming [11]. Moreover, most of the SPT algorithms are developed and optimized for spherical particles.
The study of the roto-translational dynamic of small anisotropic particles and/or of particles at high densities is thus the realm of ensemble-averaging techniques, such as depolarized dynamic light scattering (DDLS) [12,13], polarized fluorescence recovery after photobleaching (pFRAP) [14,15], and nuclear magnetic resonance (NMR). However, the main limitation of these approaches is that they cannot be used in complex, highly heterogeneous environments, which is a common scenario in biology and material science, and one for which microscopy is particularly utilized.
An intermediate approach is provided by a family of methods known as digital Fourier microscopy (DFM) [16]: As in microscopy, temporal image sequences of the sample are acquired in direct space but the sample intermediate scattering function f (q, τ ) is extracted via image correlation in the reciprocal (Fourier) space, as a function of the scattering wave-vector q and of the time delay τ . One of the most known implementations of DFM is differential dynamic microscopy (DDM), a technique that has found vast application with a variety of soft and biological matter systems [17,18]. DDM has been successfully used for studying the translational dynamics of shape-anisotropic or functionally anisotropic particles [19][20][21][22]. In combination with a suitable imaging mode, like polarized [23,24] or darkfield [25] microscopy, DDM enables the accurate measurement of the roto-translational dynamics of optically and/or shape-anisotropic particles.
While measuring the translational motion in DDM is based on the simple fact that a translating particle occupies different portions of an image, the key for DDM to be able to capture the particle rotational motion is to convert the latter into an intensity fluctuation of the particle image. For example, the image of optically anisotropic, uniaxial particle observed between cross-polarizers appears darker or brighter, according to the angle formed between its optical axis and the axes of the polarizers. In this condition, rotational diffusion, by promoting the random reorientation of the particle, leads to an intermittent "blinking" of the intensity associated to the particle image. If DDM analysis is now performed on an image sequence collected in these conditions, a q-independent decay is observed in f (q, τ ), which is particularly evident in the low-q regime. With proper fitting of f (q, τ ) the characteristic correlation time of this decay can be measured and used to determine the associated rotational diffusion coefficient [23]. A similar effect, although realized via a different optical mechanism, is exploited with dark-field microscopy to measure the characteristic reorientation time of micrometer-scale, shape-anisotropic particles [25].
Using the same strategy for small anisotropic objects imaged in bright-field or in non-polarized fluorescence microscopy remains much more challenging: An inplane rotation of the object simply leads to a rotated image, while out-of-plane rotations can introduce subtle shape changes, but with no effect on the overall associated intensity. In reciprocal space, this corresponds to the fact that the low-q portion of the image power spectrum is not affected by rotations. In theory, rotational diffusion contributions become important in scattering at large wave vectors q, where the object rotation can induce a fluctuation in the scattered intensity [12,26]. However, the relaxation of the intermediate scattering function for such large q is usually dominated by the translational term whose relaxation rate, in the case of a Brownian particle, rapidly grows as D T q 2 , where D T is the translational diffusion coefficient. To the best of our knowledge, the only study reporting a DDM-based measurement of the rotational dynamics of colloidal particles in bright-field is Ref. [27], in which the particles translational motion was suppressed by partially tethering them to a solid surface. A general recipe for using bright-field microscopy to simultaneously characterize the translational and rotational diffusive dynamics of Brownian particles remains thus unavailable.
In this work, we address this issue by combining the standard DDM analysis of bright-field movies with a simple image preprocessing step, whose function is to digitally reproduce the intensity change induced by particle rotation that is obtained by optical means in polarized or dark-field microscopy: For each image, we generate a local orientation map whose intensity level encodes the local particle orientation at those image pixels where a particle is present; we then perform a standard DDM analysis on these local orientation maps and obtain the accurate characterization of both translational and rotational diffusion. A detailed description of the implementation of the method is reported in Sect. 2.3, while a simple analytic justification is discussed in Appendix A.

Sample preparation and imaging
The sample is a highly diluted aqueous suspension of hematite-silica core-shell peanut-shaped particles, synthesized according to the protocol described in detail in Refs. [28,29]. The length and diameter of the lobes of these particles are (1723 ± 50) nm and (740 ± 50) nm, respectively. Before measurement, the sample is loaded in a glass capillary of thickness 100 µm, which is then sealed with vacuum grease. After about one hour at room temperature ((T = 22 ± 1) • C), the sample appears to be sedimented close to the bottom of the container as the sedimentation length l g ∼ 175 nm of these particles is smaller than their diameter [30]. The final number density is about 5 · 10 3 mm −2 , as obtained from particle counting (see Fig. 2).
The sample is observed with an inverted Nikon Ti-E bright-field microscope equipped with a Hamamatsu Orca Flash 4.0 camera (pixel size 6.5 µm). The condenser diaphragm (numerical aperture NA = 0.52) is kept completely open to achieve incoherent illumination of the sample. Two different objectives are used, a 40×, 0.60 NA objective and a 10× 0.25 NA objective. The size of the imaged region is 6.9 · 10 3 µm 2 and 2.8 · 10 4 µm 2 for the higher and the lower magnification objective, respectively. Image sequences of 10 4 frames each were recorded at 25 frames/s.

Differential dynamic microscopy
The translational dynamics of the particles is characterized by analyzing the sample movies with the standard DDM analysis scheme [17,31]. In short, we calculate the difference ΔI(x, t, τ) = I(x, t+ τ ) − I(x, t) between two images acquired at times t and t + τ . By averaging the spatial Fourier power spectrum of ΔI(x, t, τ) obtained for the same τ but different reference times t, we obtain the so called image structure function that captures the sample dynamics as a function of the two-dimensional scattering wavevector q and of the lag time τ . The symbol· indicates the two-dimensional digital Fourier transform, usually performed with a Fast Fourier Transform algorithm. In a last step, we take advantage of the circular symmetry of the problem to also perform an azimuthal average of D(q, τ), which leads to the one-dimensional function D(q, τ ) of the radial wavevector q = q 2 x + q 2 y . The obtained struc-ture function is connected to the intermediate scattering function (ISF) f (q, τ ) [12] by the relation where the term B(q) accounts for the camera noise and the amplitude A(q) depends on the scattering properties of the sample and the transfer function of the imaging system [16].

Squared-gradient differential dynamic microscopy
In order to probe also the rotational dynamics of the particles, we used a digital preprocessing step of the bright-field images to reveal information about the particles orientation. A simple implementation of this concept, which we call "squared-gradient" (SG), involves the calculation of the partial derivative of the image with respect to a given component μ, which is then squared to give While other implementations are also possible [32], in this work, the partial derivatives in Eq. 3 are calculated simply as the difference between the original image and its copy translated along one of the main axes by a single pixel. As it can be appreciated in Fig. 1, where a single elongated particle is considered, the obtained map has the required property, as the overall intensity changes according to the particle's orientation. Indeed, when the long axis of the particle is oriented along the μ axis (the vertical direction in Fig. 1), the contrast in the gradient map is lower compared to the case where the particle is oriented in the perpendicular direction. In the square-gradient map c μ (x, t), the particle is replaced by two smaller lobes. The most evident effect of the particle rotation on c μ (x, t) is a modulation of the intensity associated with the two lobes. As shown in Appendix A, a similar intensity modulation in c μ (x, t) is also produced by changes in the angle formed by the long axis of the particle with optical axis. Once a sequence of SG-maps is obtained from the original images, the standard DDM analysis is performed to calculate the corresponding structure function where we also average over the two orthogonal directions. As it is shown explicitly in Appendix A for the special case q = 0, D SG (q, τ ) can be expressed in terms of the normalized ISF via Eq. 2, and we employ the general expression Fig. 1 "Squared-gradient" maps for a single anisotropic particle whose long axis lays on the image plane. First row: simulated images of the same particle with bivariate Gaussian intensity profile and different orientations. Second row: "gradient" maps. Each map correspond to the partial derivative along the vertical direction of the image above it. Third row: "squared-gradient" maps. Each map corresponds to the square of the map above it. In each row, maps are shown with the same grayscale where f T is the translational ISF, f RT is the rototranslational one, and α(q) is the q-dependent relative amplitude of the roto-translational term. In Eq. 5, we neglect a term describing the coupling between orientational and translational diffusion, which is expected to be of little significance [3]. For monodisperse Brownian particles, the ISFs take the well-know expressions where D T and D R are the average translational and rotational diffusion coefficients, respectively [12]. It is worth stressing that the taking the square in Eq. 3 represents an essential step of the procedure, as it enables "translating" the orientation-dependent spatial modulation introduced by the gradient operation into a global intensity variation, which affects in a qindependent fashion all the Fourier components of the SG map.

Single particle tracking
To validate SG-DDM, we analyze the same images also with single particle tracking (SPT), by using the MATLAB code developed by the Kilfoil group at the University of Massachusetts [33] and available at https://github.com/dsseara/microrheology. We reconstruct individual trajectories of the center of mass of the particles, and we calculate their mean-squared displacement Δr 2 (τ ) as a function of the delay time τ . To track over time the orientation of the particles, we use a custom MATLAB code implementing a procedure similar to the one described in Refs. [34,35], which is described in detail in Appendix B. We applied the described procedure, which provides an independent estimate of the rotational diffusion coefficient, only to the image sequence recorded with the higher magnification (40×), as we were not able to successfully perform  panel (a). Each particle appears to be replaced by two bright spots whose intensity depends on the particle orientation the same orientational tracking-based analysis on the image sequences acquired with a 10× objective.

Results and discussion
We show in Fig. 3 representative ISFs obtained from both standard DDM (panel a) and SG-DDM (panel b) analyses of the same sequence of bright-field images (magnification 40×) of the sample. The ISFs obtained from standard DDM analysis exhibit a single exponential decay e −Γ (q)τ , whose relaxation rate Γ (q) displays a clean ∼ q 2 scaling (Fig. 4, gray triangles). A fit of Γ (q) with a quadratic model Γ (q) = D T q 2 provides an estimate for the translational diffusion coefficient D T = 0.24 ± 0.02 µm 2 /s of the particles, which is in excellent agreement with the value D T,P T = 0.23 ± 0.01 µm 2 /s obtained with SPT (Fig. 5a) via a linear fit Δr 2 (τ ) = 4D T,P T τ of the MSD of the particles.
Inspection of the ISFs obtained from SG-DDM (Fig.  3b) reveals instead the existence of a double decay, which is particularly evident for the lowest values of q. Fitting a double exponential model to the data allows to reliably estimate, at least for q 2.5 µm −1 , two well distinct decay rates: a fast decay with an almost q-indpendent rate Γ 1 , and a slow diffusive decay with rate Γ 2 ∼ q 2 that is fully compatible with the relaxation rate obtained from standard DDM (Fig. 4). For q > 2.5 µm −1 , the two relaxation rates are too close to be reliably separated and the fitting procedure provides a single relaxation rate.
According the discussion in Sec.2.3, we attribute the fast decay to the roto-traslational diffusion of the particles, Γ 1 (q) = 6D R + D T q 2 , and we obtain the estimate D R = 0.30 ± 0.05 s −1 . This result agrees with those obtained with SPT according to the procedures described in Section 2.4, in particular with the estimate D R,P T1 = 0.28 ± 0.05 s −1 obtained from a linear fit of the angular MSD Δφ 2 (τ ) = 2D R,P T1 τ (Fig. 5b) from an exponential fit of the effective aspect ratio autocorrelation function C (τ ) = e −6DR,P T 2 τ (Fig. 5c).
We perform the same analyses on an image sequence of the same sample, acquired with a 10X objective (NA = 0.25). In this imaging condition, the particle size is very close to the diffraction limit (d ∼ λ/(2NA) 1.2 µm, assuming λ 0.6 µm) and the effect of the convolution with the point-spread-function is to markedly reduce the apparent anisotropy in the particle image, that is now barely appreciable (inset of Fig. 4). As a consequence, the uncertainty in the frame-by-frame determination of single particle orientation is so large that the angular tracking procedure described in Sect. 2.4 no longer provides reliable results. On the contrary, the results of both DDM and SG-DDM are in excellent agreement with the one obtained with the 40X objective (NA = 0.60), as shown in Fig. 4.
As a further validity check for the proposed approach, we compare the obtained diffusion coefficients with the ones predicted by the theory for a rigid rod, which represents a fair model for our peanut-shaped particles. If we consider a cylinder of length 1723 nm and diameter 740 nm, in a Newtonian fluid of viscosity η = 1.0 mPa at temperature T = 22 • C, by using the theoretical expressions reported in Ref. [36] (Eqs. 1-4), we obtain the values D R,th = 0.44 s −1 , D T,th = 0.25· µm 2 /s for the rotational and the translational diffusion coefficient, respectively. We note that, while D T,th is in very good agreement with the corresponding values obtained with both SG-DDM and SPT, the calculated diffusion coefficient D R,th is off by about 30% with respect the experimental ones. This moderate discrepancy could be explained as a systematic effect of the adopted cylindrical approximation. However, it can be also a genuine effect due the hydrodynamic interactions between particles and the bottom plate of the container, which can slow down the particles' diffusivity compared to the bulk [37].

Conclusions
In this work, we introduced SG-DDM, a simple procedure enabling the tracking-free determination of the rotational and translational dynamics of anisotropic particles in microscopy images. SG-DDM was here demonstrated for bright-field microscopy experiments, but we trust that it could be applied without modifications to a broad range of imaging modes, including wide-field or confocal fluorescence microscopy.
Compared to SPT-based approaches, we highlight a number of significant advantages. First of all, at a variance with most SPT methods, the proposed approach does not require user-dependent parameters and does not rely on a specific model for the particle shape. Moreover, by exploiting the intrinsic ensemble-averaging capability of the method, we proved that SG-DDM is able to provide a quantitative information on the rotational dynamics even when the size of the particles is very close to the diffraction limit, and a low signal-to-noise ratio does not allow a reliable frameby-frame determination of single-particle orientation in real space. We thus expect the method to be compatible with a wide range of shape anisotropic particles, such as for example metallic nanowires [38], as well as with spherical, but optically anisotropic, objects, like Janus particles [27]. Computationally, SG-DDM is as effective as DDM, whose use becomes advantageous in particular when a large number of particles are considered. In fact, the computation time does not depend on the number of imaged particles, so that thousands of particles can be analyzed in parallel in a single movie.
Being based on a nonlinear preprocessing step, a condition under which the method is expected to provide the correct translational dynamics of the sample is the absence of a significant overlap between images of the different particles [25,31]. While this condition is almost automatically fulfilled for rigid objects in two dimen- and c-d show results obtained from the analysis of image sequences of the same sample, acquired with a 40×, 0.60 NA objective and a 10×, 0.25 NA objective, respectively. a Gray triangles: relaxation rate Γ (q) obtained from the single exponential fit of the structure functions obtained from standard DDM analysis. Orange circles and blue squares: relaxation rates Γ1(q) and Γ2(q) obtained from a double exponential fit of the structure functions obtained from SG-DDM analysis on the same image sequence. The continuous line Γ = DT q 2 is a best fitting curve to the blue squares. The analogous curve obtained from the gray triangle is indistinguishable. The dashed curve is a best fitting curve Γ1 = 6DR + DT q 2 to the orange circles. The vertical dashed line marks the q-value above which the two relaxation rates Γ1 and Γ2 cannot be resolved any more. The small inset shows a representative ROI of size 10.4 µm. b Orange circles: relative amplitude α(q) of the roto-translational contribution to the ISF (see Eq.5), obtained from a double exponential fit of the SG-DDM structure functions. The description of panels c and d is identical to the one of panels a and b, respectively c Symbols: normalized temporal correlation function of the apparent aspect ratio (t) = b 2 (t)/a 2 (t) of each particle, as obtained from a bivariate Gaussian fit of the associated intensity distribution (see main text for details). The continuous line is a best fitting curve to the data with an exponential decay C (τ ) = exp(−γτ ), with γ = (1.75 ± 0.2) s −1 sions, it could represent a more serious limitation on the particle concentration in the three dimensional case.
In addition, if the particles deviate significantly from a simple uniaxial symmetry, such as the one investigated here, other contributions (corresponding to highorder terms in the spherical harmonics expansion of the particles' shape) become relevant in the ISF. In those cases, a more refined system-dependent model for the ISF, including multiple exponential decays, should be used to properly interpret the SG-DDM signal.
Finally, also the presence of a strongly inhomogeneous optical background can affect the quality of the results, as the nonlinear processing can introduce a coupling between moving objects of interest and static features. A mitigation of this problem could be provided by a background subtraction before calculating the SGmap. For example, if the observation time window is long enough to enable a complete decorrelation in the particles positions, a good estimate of the background intensity could be obtained as the time average of the whole sequence.
In all the above-mentioned cases, a quite compelling internal quality test is provided by the comparison between the results obtained from standard DDM and SG-DDM on the same sequence, as the exact same translational dynamics is expected to be measured in DDM and SG-DDM. Any lack of correspondence between the translational ISFs obtained with the two methods is a proxy of a potential artifact introduced by the SG mapping, offering a simple and sensitive check of the validity of the mapping procedure.
Finally, beyond offering a novel way to determine the rotational dynamics of small anisotropic entities, SG-DDM represents also the first example of nonlinear preprocessing of microscope images prior to DDM analysis, which we believe will expand the range of applications of DDM and other DFM techniques beyond the current frontier.

Author contribution statement
FG and RC designed research. FG performed experiments and analyzed the data. AP synthesized the particles and prepared the samples. All authors discussed the experimental results. FG and RC wrote the paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix A. Analytical calculation of D SG (q = 0, τ )
In this Appendix, we provide analytical justification of SG-DDM as a probe of the rotational diffusion of anisotropic particles. In particular, we show that, in the q → 0 limit, the structure function DSG(q, τ ) is dominated by a qindependent relaxation mirroring the rotational dynamics.
We assume a simple model where each "particle" described by a 3D Gaussian profile with variances (a 2 , a 2 , b 2 ) along its three main axes (a < b). Let (x0, y0, z0) be the coordinate of its center of mass, while the orientation of the long axis is determined by the polar angle θ (i.e., the angle formed with the z-axis, corresponding to the optical axis) and the azimuthal angle φ.
In what follows, we consider for simplicity the case of a single particle, the generalization to a collection of identical, non-overlapping, particles being straightforward. We assume that the intensity distribution I(x, y) on the image plane is given by the projection on the x − y plane of the object density distribution (good for example in the case of incoherent bright-field or wide-field fluorescence microscopy of a particle close to the object plane [31]). Under this hypothesis, I(x, y) is a 2D Gaussian distribution where c is a dimensionless constant. We initially consider the case of a moderately anisotropic particle (i.e., one whose eccentricity e ≡ 1 − a 2 /b 2 is small), for which we obtain the following expression, valid to the first order in e: i0(θ, φ) = c 1 + e 2 sin 2 θ (1 + (1/2) cos 2φ) , (6) where c = (2c/a 4 ).
Combining Eq. 7 with the identity [12] Y m l Y m we obtain the following expression for DSG(q = 0, τ) This simple expression is valid under the hypothesis of small shape anisotropy. Within the same approximation, the effective aspect ratio = σy/σx is given by = 1 + 1 2 e 2 sin 2 θ and the corresponding time autocorrelation function C (τ ) = (t + τ ) (t) − (t) 2 / (t) 2 − (t) 2 can be easily calculated as Fig. 6 Schematic representation of the rotational tracking algorithm. The coordinates (xm, ym) of the centroid of each particle are determined via a standard particle tracking algorithm [33]. The shape and orientation of the particle's image are determined via a bivariate Gaussian fit to the local intensity profile. This enables estimating the length σx and σy of the short and of the long axis, respectively, and the azimuthal angle φ formed by the particle's long axis with the horizontal direction. The effective aspect ratio = σy/σx depends of the polar angle θ between the particle's axis and the vertical direction, as described in the text If the eccentricity is large, higher-order spherical harmonics are expected to appear in Eq. 7, leading to multiple exponential decay of the correlation functions [12].

Appendix B. Rotational tracking algorithm
In this Appendix, we present a detailed description of the algorithm used to track the rotational motion of the particles in real space. As described in Par. 2.4, we first apply a standard particle tracking algorithm [33] to determine, frame by frame, the position of the centroid of each particle. For each tracked particle i and for each time t, a square ROI of size 16x16 pixel (corresponding to area of about 27 µm 2 when the 40X objective is used) is obtained from the original image, centered on the particle's centroid; a bivariate Gaussian distribution with arbitrary orientation of the main axes f (x, , is fitted to the ROI intensity distribution, which provides an estimate for the variances σ 2 x (t) and σ 2 y (t) along the main axes, and for the azimuthal angle φi(t) identifying the direction of projection of the long axis on the x, y plane. The values obtained for σx(t) and σy(t) are quite consistent along each trajectory and display time-correlated fluctuation around their mean values, which we attribute to random variations in the polar angle θ formed by the long axis of the particle with the vertical direction (see Fig. 6).
We then calculate the autocorrelation function C (τ ) = (t + τ ) (t) − (t) 2 / (t) 2 − (t) 2 of the effective aspect ratio (t) = σy(t)/σx(t). For an elongated particle performing rotational Brownian motion, C (τ ) is expected to display a simple exponential relaxation C (τ ) = e −γτ , with γ = 6DR, as shown in Appendix A. This provides a (a) (b) Fig. 7 a Blue (orange) curve: time-dependent length of the major (minor) axis σy(t) (σx(t)) of a single particle, as obtained from a bivariate Gaussian fit of its intensity distribution. The black continuous line corresponds to the mode σy of the frequency distribution of σy(t), while the black dashed line corresponds to threshold value 0.95σy (see main text for details). b Frequency distribution of σy(t).
The vertical line corresponds to the mode σy = (0.90±0.02) µm of the distribution robust means to measure the rotational dynamics of the particles, namely by fitting an exponential function to the obtained C (τ ). In order to extract a similar information from the azimuthal degree of freedom, which is encoded in φ, the simultaneous knowledge of the polar angle θ is also needed. This is due to the fact that the same change in the azimuthal angle δφ can correspond to actual displacements of very different amplitudes according to the value of θ [11]. To account for this, a different variable ψ is often considered, which is defined by the relation dψ = dφ sin θ [6]. For a particle performing rotational diffusion, the mean squared value of ψ follows the simple relation Δψ 2 (τ ) = 2DRτ .
In our experiments, although we could in principle calculate θ(t) from σx(t), σy(t), we found that this direct inversion was too noisy to provide reliable results. As a consequence, we adopt a simplified yet more robust approach. For each particle, we consider the histogram of the values assumed by the long axis length σy(t), and we identify the value σy corresponding to the right peak of the distribution (see Fig.7). According to the model described in the Appendix, σy corresponds to the length of the long axis when it lies on the x − y plane, i.e., when θ = π/2. For each particle, we identify the portions of the trajectory such that (σy − σy(t))/σy < 0.05.