Elucidating fluctuating diffusivity in center-of-mass motion of polymer models with time-averaged mean-square-displacement tensor

TL;DR

This paper introduces a tensor-based method to analyze fluctuating diffusivity in polymer models using trajectory data, revealing distinct fluctuation patterns across different models and applicable to experimental SPT data.

Contribution

A novel tensor-based approach to extract and analyze fluctuating diffusivity from trajectory data, capturing both magnitude and orientation fluctuations.

Findings

Different polymer models show unique fluctuation patterns.

The method distinguishes between magnitude and orientational diffusivity fluctuations.

Applicable to single-particle-tracking experimental data.

Abstract

There have been increasing reports that the diffusion coefficient of macromolecules depends on time and fluctuates randomly. Here, a novel method to elucidate the fluctuating diffusivity from trajectory data is developed. The time-averaged mean square displacement (MSD), a common tool in single-particle-tracking (SPT) experiments, is generalized into a second order tensor, with which information on both magnitude fluctuations and orientational fluctuations of the diffusivity can be clearly detected. This new method is utilized to analyze the center-of-mass motion of four polymer models, the Rouse model, the Zimm model, a reptation model, and a rigid rodlike polymer, and it is found that these models exhibit distinctly different types of magnitude and orientational fluctuations of the diffusivity. This method of the time-averaged MSD tensor can be applied also to the trajectory data obtained in SPT experiments.

Equations

$$\overline{\delta {\mathbf{r}}^2}(\Delta ,t):=\frac{1}{t-\Delta }{\int }_0^{t-\Delta }\delta {\mathbf{r}}^2(\Delta ,t^{\prime })d{t}^{\prime },$$

$$\frac{d\mathbf{r}(t)}{dt}=\sqrt{2}\mathbf{B}(t)\cdot \mathbf{\xi }(t),$$

$$⟨\mathbf{\xi }(t)\mathbf{\xi }(t^{\prime })⟩=\mathbf{I}\delta (t-t^{\prime }),$$

$$\mathbf{D}(t)=\mathbf{B}(t)\cdot {\mathbf{B}}^T(t),$$

$$\overline{\delta \mathbf{r}\delta \mathbf{r}}(\Delta ,t):=\frac{1}{t-\Delta }{\int }_0^{t-\Delta }\delta \mathbf{r}(\Delta ,t^{\prime })\delta \mathbf{r}(\Delta ,t^{\prime })d{t}^{\prime },$$

$${\left[\overline{\delta \mathbf{r}\delta \mathbf{r}}(\Delta ,t)\right]}_{ij}=\frac{1}{t-\Delta }{\int }_0^{t-\Delta }\delta r_i(\Delta ,t^{\prime })\delta r_j(\Delta ,t^{\prime })d{t}^{\prime }.$$

$$⟨\overline{\delta \mathbf{r}\delta \mathbf{r}}(\Delta ,t)⟩$$

$$=2{\int }_0^{\Delta }d{t}_1{\int }_0^{\Delta }d{t}_2⟨\mathbf{B}(t_1)\cdot \mathbf{\xi }(t_1)\mathbf{B}(t_2)\cdot \mathbf{\xi }(t_2)⟩$$

$$=2⟨\mathbf{D}⟩\Delta ,$$

$$⟨B_{ik}(t_1){\xi }_k(t_1)B_{jl}(t_2){\xi }_l(t_2)⟩$$

$$=⟨B_{ik}(t_1)B_{jl}(t_2)⟩{\delta }_{kl}\delta (t_1-t_2),$$

$$⟨\overline{\delta {\mathbf{r}}^2}(\Delta ,t)⟩=2\mathrm{tr}⟨\mathbf{D}⟩\Delta .$$

$$\mathbf{\Phi }(\Delta ,t):=$$

$$=$$

$$\mathbf{\Phi }(\Delta ,t)\approx {\mathbf{\Phi }}^{\mathrm{id}}(\Delta ,t)+{\mathbf{\Phi }}^{\mathrm{ex}}(\Delta ,t),$$

$${\Phi }_{ikmp}^{\mathrm{id}}(\Delta ,t)=\frac{2{\Delta }^3}{3t}\left(4-\frac{\Delta }{t}\right)\left(⟨D_{im}{D}_{kp}⟩+⟨D_{ip}{D}_{km}⟩\right),$$

$$\displaystyle\bm{\Phi}^{\mathrm{ex}}(\Delta,t)=\frac{8\Delta^{2}}{t^{2}}\!\int_{0}^{t}\!d\tau(t-\tau)\Bigl{[}\left\langle\bm{D}(\tau)\bm{D}(0)\right\rangle_{s}-\left\langle\bm{D}\right\rangle\!\left\langle\bm{D}\right\rangle\Bigr{]}.$$

$${⟨\mathbf{D}(\tau )\mathbf{D}(0)⟩}_s:=\frac{⟨\mathbf{D}(\tau )\mathbf{D}(0)⟩+⟨\mathbf{D}(0)\mathbf{D}(\tau )⟩}{2}.$$

$${\mathbf{\Psi }}^1(\Delta ,t)$$

$${\mathbf{\Psi }}^2(\Delta ,t)$$

$${\Psi }_{ikmp}^1(\Delta ,t)$$

$$+$$

$${\mathbf{\Psi }}^2(\Delta ,t)$$

$${\Phi }_{ikmp}(\Delta ,t)$$

$$+\frac{8{\Delta }^2}{t^2}$$

$${\Phi }_1(\Delta ,t)=$$

$$\approx$$

$${\Phi }_1^{\mathrm{id}}(\Delta ,t):=$$

$${\Phi }_1^{\mathrm{ex}}(\Delta ,t):=$$

$${\hat{\Phi }}_1(\Delta ,t)\approx {\hat{\Phi }}_1^{\mathrm{id}}(\Delta ,t)+{\hat{\Phi }}_1^{\mathrm{ex}}(t).$$

$${\hat{\Phi }}_1^{\mathrm{id}}(\Delta ,t)$$

$${\hat{\Phi }}_1^{\mathrm{ex}}(t)$$

$${\varphi }_1(\tau )$$

$${\varphi }_2(\tau )$$

$$C:=n\frac{\mathrm{tr}\left(⟨\mathbf{D}⟩\cdot ⟨\mathbf{D}⟩\right)}{{\left(\mathrm{tr}⟨\mathbf{D}⟩\right)}^2}.$$

$${\Phi }_2(\Delta ,t)=$$

Full text

I Introduction

Macromolecular diffusion in cytoplasm and cell membranes has received much attention in recent years, because it controls chemical kinetics and information processing in cells Berg et al. (1981). Single-particle-tracking (SPT) techniques have been used to study macromolecular diffusion in living systems, and remarkably complicated phenomena such as anomalous diffusion, weak ergodicity breaking (EB), and sample-to-sample fluctuations of the diffusion coefficient have been reported Golding and Cox (2006); Weigel et al. (2011); Parry et al. (2014); Jeon et al. (2011); Tabei et al. (2013). In such SPT experiments, a time average is commonly used to obtain the mean square displacement (MSD); the time-averaged MSD (TMSD) of a tagged particle is defined by Nordlund (1914); He et al. (2008); Metzler et al. (2014)

[TABLE]

where $\Delta$ is a lag time and $t$ is the total measurement time. In addition, a displacement vector $\delta\bm{r}(\Delta,t^{\prime})$ is defined as $\delta\bm{r}(\Delta,t^{\prime}):=\bm{r}(t^{\prime}+\Delta)-\bm{r}(t^{\prime}),$ where $\bm{r}(t^{\prime})$ is the position vector of the tagged particle at time $t^{\prime}$. Thus, the TMSD $\overline{\delta\bm{r}^{2}}(\Delta,t)$ can be obtained from a single trajectory $\bm{r}(t^{\prime})$.

In SPT experiments of macromolecules in living systems, sample-to-sample fluctuations of the diffusion coefficient have been observed frequently Golding and Cox (2006); Weigel et al. (2011); Parry et al. (2014); Jeon et al. (2011); Tabei et al. (2013). As stated above, the TMSD curve $\overline{\delta\bm{r}^{2}}(\Delta,t)$ (as a function of $\Delta$) is obtained from a single trajectory $\bm{r}(t^{\prime})$, and then, from this TMSD curve, the diffusion coefficient for that trajectory can be estimated. The values of this diffusion coefficient vary from trajectory to trajectory, but, for long trajectories (namely, at $t\to\infty$), they converge to a single value if the system is ergodic. In some SPT experiments, however, the values of the diffusion coefficient are scattered even for long trajectories, and this phenomenon cannot be explained by the ordinary Brownian motion Golding and Cox (2006); Jeon et al. (2011); Weigel et al. (2011); Tabei et al. (2013); Parry et al. (2014).

To explain such sample-to-sample fluctuation in the diffusivity, much effort has been devoted to investigating simple theoretical models such as the continuous-time random walk (CTRW) He et al. (2008); Lubelski et al. (2008); Neusius et al. (2009); Miyaguchi and Akimoto (2011, 2013); Thiel and Sokolov (2014a), fractional Brownian motion Deng and Barkai (2009); Thiel and Sokolov (2014b), and the random walk on fractals Miyaguchi and Akimoto (2015); Meroz et al. (2010). In these studies, the variance of the TMSD, which is commonly referred to as the EB parameter, has been used to characterize the fluctuation in the diffusivity. In particular, it was shown that the EB parameter for the CTRW converges to a non-vanishing value as $t\to\infty$. In other words, the TMSD behaves as a random variable even for long measurement times. Therefore, CTRW-like dynamics have been considered to be a factor in the sample-to-sample fluctuation of the diffusivity observed in SPT experiments He et al. (2008); Lubelski et al. (2008); Neusius et al. (2009).

However, fluctuations in diffusivity originate also from correlated dynamics of inner degrees of freedom. In Ref. Uneyama et al. (2015), the authors studied a reptation model (a tagged polymer model in entangled polymer solutions) and showed that the EB parameter of the center-of-mass (COM) motion is non-vanishing for quite a long measurement time. In other words, the system exhibits sample-to-sample fluctuations in diffusivity, that originate from non-Markovian dynamics of the end-to-end vector. Another important finding of Ref. Uneyama et al. (2015) is that the EB parameter is related to a correlation function of magnitude of diffusivity. Unfortunately, it was also found that much of the information contained in the trajectory data $\bm{r}(t)$ is lost in the EB parameter. Therefore, it is necessary to develop an efficient method to extract more information from the trajectory data.

In this paper, a novel method is developed for elucidating the fluctuating diffusivity of macromolecules from trajectory data $\bm{r}(t)$. More precisely, a TMSD tensor, a generalization of the TMSD [Eq. (1)], is proposed, and it is shown that correlation functions of this TMSD tensor contain plenty of information including a magnitude correlation and an orientation correlation of the fluctuating diffusivity. Moreover, by using this tensor analysis, four fundamental polymer models are investigated: the Rouse and Zimm models (polymer models in dilute solutions), a reptation model (a polymer model in concentrated solutions), and a rigid rod-like polymer (an extreme case of non-flexible polymers). It is shown that the COM motion of these polymer models exhibits distinctly different types of the fluctuating diffusivity. For example, it is shown that the COM motion of the Zimm and reptation models exhibits both magnitude and orientation fluctuations of the diffusivity, whereas that of the rigid rod-like polymer exhibits only orientation fluctuations. The tensor analysis presented in this article could be used to analyze the trajectory data obtained in SPT experiments.

This paper is organized as follows. In Sec. II, a Langevin equation with fluctuating diffusivity (LEFD) is defined. In Sec. III, the TMSD tensor is defined and its correlation functions are studied for the LEFD. It is also shown here that these correlation functions are related to a non-Gaussian parameter. In Secs. IV– VII, the COM motion of each of the aforementioned polymer models is studied with the TMSD tensor. Finally, Sec. VIII is devoted to a discussion. In the Appendices, we summarize some technical matters, including the simulation details.

II Langevin equation with fluctuating Diffusivity

As shown in subsequent sections, the COM of polymer models such as the Zimm and reptation models can be described by the following Langevin equation with time-dependent and fluctuating diffusivity Stein and Stein (1991); Łuczka et al. (1995); Rozenfeld et al. (1998); Łuczka et al. (2000); Chubynsky and Slater (2014); Massignan et al. (2014); Manzo et al. (2015); Cherstvy and Metzler (2016); Uneyama et al. (2015); Miyaguchi et al. (2016); Chechkin et al. (2017):

[TABLE]

where $\bm{r}(t)$ is an $n$-dimensional position vector of a tagged particle at time $t$, and the $n\times n$ matrix $\bm{B}(t)$ is a stochastic process. Moreover, $\bm{\xi}(t)$ is white Gaussian noise that satisfies

[TABLE]

where $\bm{I}$ is the identity matrix. Equation (2) is referred to as the LEFD.

In this study, it is assumed that $\bm{\xi}(t)$ and $\bm{B}(t)$ are mutually independent stochastic processes. Consequently, the diffusion coefficient tensor $\bm{D}(t)$ is given by

[TABLE]

where $\bm{B}^{T}$ is the transpose matrix of $\bm{B}$. It follows that $\bm{D}(t)$ is a symmetric tensor: $\bm{D}(t)=\bm{D}^{T}(t)$. In addition, $\bm{D}(t)$ is assumed to be a stationary process.

III TMSD tensor

In this section, the TMSD tensor is defined and its general properties are presented. In particular, it is shown that the TMSD tensor of the LEFD exhibits only normal diffusion, even though the density profile is non-Gaussian. Moreover, to extract information on the fluctuating diffusivity, correlation functions of the TMSD tensor are studied. In particular, a novel method to extract magnitude and orientation correlations of the diffusivity is presented.

III.1 TMSD tensor exhibits normal diffusion

As a generalization of the TMSD [Eq. (1)], a TMSD tensor (a second-order tensor) is defined as

[TABLE]

where the integral is taken for each element of the tensor in the integrand as

[TABLE]

Here $\delta r_{i}(\Delta,t^{\prime})$ is an element of $\delta\bm{r}(\Delta,t^{\prime})$, and $\left[\bm{H}\right]_{ij}$ represents an element of a second-order tensor $\bm{H}$: $\left[\bm{H}\right]_{ij}:=H_{ij}$.

Note that the TMSD tensor $\overline{\delta\bm{r}\delta\bm{r}}(\Delta,t)$ is the time-averaged counterpart of the ensemble-averaged MSD tensor Dhont (1996). Taking the trace of Eq. (5), we obtain the TMSD given in Eq. (1), and thus it is possible to extract more information with the TMSD tensor than with the TMSD. Moreover, taking the ensemble average in Eq. (5) and using Eqs. (2)–(4), we have

[TABLE]

where $\left\langle\dots\right\rangle$ is the ensemble average. For the first equality in Eq. (7), we used the stationarity of the system, and for the final equality, we used the fact that $\bm{B}(t)$ and $\bm{\xi}(t)$ are independent in the sense that

[TABLE]

where we have employed the Einstein summation convention. In particular, if the system is statistically isotropic, we have $\left\langle\bm{D}\right\rangle=D\,\bm{I}$. Taking the trace in Eq. (7), we obtain the TMSD again Uneyama et al. (2015)

[TABLE]

Surprisingly, all the elements of the ensemble-averaged TMSD tensor in Eq. (7) exhibit only normal diffusion (i.e., proportional to the lag time $\Delta$), even though the diffusion coefficient fluctuates. In other words, it is impossible to detect the fluctuating diffusivity with the first moment of the TMSD tensor [Eq. (7)], and so higher-order moments of the TMSD tensor are studied in the following subsections.

III.2 Correlation function of TMSD tensor

To extract information about the fluctuating diffusivity from trajectories $\bm{r}(t)$, we study a correlation function of the TMSD tensor

[TABLE]

where $\bm{\Phi}(\Delta,t)$ is a fourth-order tensor. Note that, in time-series analysis, Eq. (10) should be used instead of Eq. (11) to reduce numerical errors. In fact, Eq. (10) was used in all of the numerical simulations reported here.

If we assume that $\Delta$ is much shorter than a characteristic time scale $\tau_{D}$ of the fluctuating diffusivity, we can decompose $\bm{\Phi}(\Delta,t)$ into two parts (see below for a derivation) as

[TABLE]

where the fourth-order tensors $\bm{\Phi}^{\mathrm{id}}(\Delta,t)$ and $\bm{\Phi}^{\mathrm{ex}}(\Delta,t)$ are defined respectively as

[TABLE]

Here, $\left\langle\dots\right\rangle_{s}$ is a symmetrization given by

[TABLE]

Equation (12) can be derived as follows. First, $\bm{\Phi}(\Delta,t)$ is expressed as $\bm{\Phi}(\Delta,t)=\bm{\Psi}^{1}(\Delta,t)-\bm{\Psi}^{2}(\Delta,t)$, where $\bm{\Psi}^{1}(\Delta,t)$ and $\bm{\Psi}^{2}(\Delta,t)$ are fourth-order tensors defined [see Eq. (11)] as

[TABLE]

After a lengthy calculation, the elements of $\bm{\Psi}^{1}(\Delta,t)$ can be expressed (see Appendix A for detail) as

[TABLE]

where $\Phi_{ikmp}^{\mathrm{id}}(\Delta,t)$ is the ideal part defined in Eq. (13), and $\left[\bm{H}\right]_{ikmp}$ represents an element of a fourth-order tensor $\bm{H}$, i.e., $\left[\bm{H}\right]_{ikmp}:=H_{ikmp}$. On the other hand, from Eqs. (7) and (17), we have

[TABLE]

By subtracting Eq. (19) from Eq. (III.2), the elements of the fourth-order tensor $\bm{\Phi}(\Delta,t)$ are obtained as

[TABLE]

The second term in the right-hand side is equivalent to $\Phi_{ikmp}^{\mathrm{ex}}(\Delta,t)$ [see Eq. (14)], and hence Eq. (III.2) coincides with Eq. (12).

As can be seen from Eq. (14), the tensor $\bm{\Phi}^{\mathrm{ex}}(\Delta,t)$ is related to the autocorrelation function of the diffusivity tensor $\bm{D}(t)$. Thus, in contrast to the first moment of the TMSD tensor given in Eq. (7), the second moment $\bm{\Phi}(\Delta,t)$ can be used to characterize the fluctuating diffusivity. In particular, if $\bm{D}(t)$ does not fluctuate, then $\bm{\Phi}^{\mathrm{ex}}(\Delta,t)\equiv 0$; therefore, $\bm{\Phi}^{\mathrm{ex}}(\Delta,t)$ is hereinafter referred to as an excess part. In contrast, the qualitative features of $\bm{\Phi}^{\mathrm{id}}(\Delta,t)$ in Eq. (13) are independent of the fluctuating diffusivity, and therefore this part is referred to as an ideal part.

An important point is that the TMSD tensor $\overline{\delta\bm{r}\delta\bm{r}}(\Delta,t)$ and its correlation function $\bm{\Phi}(\Delta,t)$ can be calculated from the trajectory data $\bm{r}(t)$ alone, and there is no need to measure $\bm{D}(t)$. Since the trajectory data $\bm{r}(t)$ is available in many single-particle-tracking experiments, the TMSD tensor and its correlation function are useful tools for elucidating the fluctuating diffusivity. Note however that in the derivation of Eq. (12), it is assumed that $\Delta$ is shorter than a characteristic time scale $\tau_{D}$ of the fluctuating diffusivity. This means that the observation interval should be much shorter than $\tau_{D}$.

III.3 Correlation functions of diffusion coefficient

To obtain more specific information of the fluctuating diffusivity, two scalar functions $\Phi_{1}(\Delta,t)$ and $\Phi_{2}(\Delta,t)$ are derived from $\bm{\Phi}(\Delta,t)$. It is shown that these are related to a magnitude and orientation correlations, respectively, of the fluctuating diffusivity $\bm{D}(t)$.

III.3.1 Magnitude correlation of diffusion coefficient

Firstly, $\Phi_{1}(\Delta,t)$ is defined as a scalar quantity obtained by taking contractions in Eqs. (11) or (12) between the first and second indices, and also between the third and fourth indices. It follows that $\Phi_{1}(\Delta,t)$ is given by

[TABLE]

where the two scalar functions $\Phi_{1}^{\mathrm{id}}(\Delta,t)$ and $\Phi_{1}^{\mathrm{ex}}(\Delta,t)$ are defined by

[TABLE]

As can be seen from Eq. (III.3.1), $\Phi_{1}(\Delta,t)$ is the variance of the TMSD [Eq. (1)].

Furthermore, Eq. (III.3.1) can be made dimensionless by dividing it by $\langle\overline{\delta\bm{r}^{2}}(\Delta,t)\rangle^{2}=4\Delta^{2}\langle\mathrm{tr}\bm{D}\rangle^{2}$; this dimensionless quantity is denoted as $\hat{\Phi}_{1}(\Delta,t)$ and is given by

[TABLE]

Note that $\hat{\Phi}_{1}(\Delta,t)$ is the relative variance of the TMSD, which is equivalent to the EB parameter He et al. (2008); Deng and Barkai (2009); Uneyama et al. (2015); Miyaguchi et al. (2016). The two scalar functions $\hat{\Phi}_{1}^{\mathrm{id}}(\Delta,t)$ and $\hat{\Phi}_{1}^{\mathrm{ex}}(t)$ are defined respectively as

[TABLE]

Here, $n$ is the space dimension, $\phi_{1}(\tau)$ and $\phi_{2}(\tau)$ are magnitude and orientation correlation functions, respectively, of the diffusivity $\bm{D}(t)$:

[TABLE]

and $C$ is a constant defined by

[TABLE]

If the system is statistically isotropic, then we have $\left\langle\bm{D}\right\rangle=D\bm{I}$ and hence $C=1$.

As seen from Eq. (26), $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)$ is related to the magnitude correlation function $\phi_{1}(\tau)$ of the diffusivity. For example, if the magnitude of the diffusivity is constant [i.e., $\mathrm{tr}\bm{D}(t)\equiv\mathrm{const.}$] and only its direction fluctuates, we have $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)\equiv 0$ from Eqs. (26) and (27); thus, no information about the fluctuating diffusivity can be detected with $\hat{\Phi}_{1}(\Delta,t)$. This is actually the case for the COM motion of the rigid rod-like polymer (Sec. VII), and it is necessary to study a different quantity to elucidate the orientation fluctuation.

III.3.2 Orientation correlation of diffusion coefficient

To extract information about the orientation fluctuation, another scalar function $\Phi_{2}(\Delta,t)$ is defined by taking contractions in Eqs. (11) or (12) both between the second and third indices, and also between the first and fourth indices. Consequently, $\Phi_{2}(\Delta,t)$ is given by

[TABLE]

where a double dot product ”$:$” is defined by $\bm{A}\!:\!\bm{B}=\sum_{ij}A_{ij}B_{ji}$, and $\Phi_{2}^{\mathrm{id}}(\Delta,t)$ and $\Phi_{2}^{\mathrm{ex}}(\Delta,t)$ are scalar functions defined respectively as

[TABLE]

Again, let us make Eq. (31) dimensionless by dividing it by $\langle\overline{\delta\bm{r}\delta\bm{r}}(\Delta,t)\,\rangle:\left\langle\overline{\delta\bm{r}\delta\bm{r}}(\Delta,t)\right\rangle=4\Delta^{2}\mathrm{tr}\left(\langle\bm{D}\rangle\cdot\langle\bm{D}\rangle\right)$; we denote this dimensionless quantity as $\hat{\Phi}_{2}(\Delta,t)$, which is given by

[TABLE]

where the two scalar functions $\hat{\Phi}_{2}^{\mathrm{id}}(\Delta,t)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(\Delta,t)$ are defined as

[TABLE]

The function $\phi_{2}(\tau)$, which is defined in Eq. (28), represents an orientation correlation of the diffusivity, and hence information about the orientation correlation can be extracted by using $\hat{\Phi}_{2}(\Delta,t)$. Note however that, for the case in which the diffusivity tensor $\bm{D}(t)$ is given by a scalar function $D(t)$ as $\bm{D}(t)=D(t)\bm{I}$, the two functions $\hat{\Phi}_{1}^{\mathrm{ex}}(t)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ are equivalent: $\hat{\Phi}_{1}^{\mathrm{ex}}(t)=\hat{\Phi}_{2}^{\mathrm{ex}}(t)$. In this sense, $\phi_{2}(\tau)$ includes information about the magnitude correlation of the diffusivity as well as its orientation correlation; therefore $\tilde{\phi}_{2}(\tau):=\phi_{2}(\tau)-\phi_{1}(\tau)/C$ may be more suitable as an orientation correlation. In what follows, however, $\phi_{2}(\tau)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ are referred to as orientation correlation functions for simplicity. The special case in which $\bm{D}(t)=D(t)\bm{I}$ was studied extensively in Ref.Miyaguchi et al. (2016).

III.4 Non-Gaussian parameter

A non-Gaussian parameter of the displacement vector $\delta\bm{r}(t)=\bm{r}(t)-\bm{r}(0)$ is defined as Kegel and van Blaaderen (2000); Arbe et al. (2002); Ernst et al. (2014); Cherstvy and Metzler (2014)

[TABLE]

In Ref. Uneyama et al. (2015), it was shown that the non-Gaussian parameter $A(t)$ for the LEFD [Eq. (2)] is given by

[TABLE]

For isotropic systems, we have $C=1$; and hence the first term vanishes. Equation (38) shows that the non-Gaussian parameter $A(t)$ can be decomposed into two parts; one originates from the magnitude correlation of the diffusivity, and the other from its orientation correlation. Although Eq. (38) was derived previously in Ref. Uneyama et al. (2015), it was not known then how to calculate $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ from the trajectory data $\bm{r}(t)$. Therefore, the method for obtaining $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ as presented in the previous subsection is one of the main results of this article.

III.5 Isotropic case

If the system is statistically isotropic, $\bm{\Phi}(\Delta,t)$ is a fourth-order isotropic tensor. Moreover, from its definition [Eq. (10)], $\bm{\Phi}(\Delta,t)$ has the following symmetry properties: $\Phi_{ijkl}=\Phi_{jikl}$, $\Phi_{ijkl}=\Phi_{ijlk}$, and $\Phi_{ijkl}=\Phi_{klij}$. It follows that $\bm{\Phi}(\Delta,t)$ can be expressed as

[TABLE]

where $\lambda(\Delta,t)$ and $\mu(\Delta,t)$ are scalar functions (these functions are analogous to the Lamé coefficients in the theory of elasticity for isotropic bodies Landau and Lifshitz (1986)). Thus, in the isotropic case, the fourth-order tensor $\bm{\Phi}(\Delta,t)$ is characterized completely by $\lambda(\Delta,t)$ and $\mu(\Delta,t)$. Taking contractions in Eq. (39) between the first and second indices (i.e., $i$ and $j$) and between the third and fourth indices (i.e., $k$ and $l$), we have

[TABLE]

Similarly, taking contractions between the first and fourth indices (i.e., $i$ and $l$) and between the second and third indices (i.e., $j$ and $k$), we have

[TABLE]

Thus, we reach a significant conclusion that the two scalar functions $\Phi_{1}(\Delta,t)$ and $\Phi_{2}(\Delta,t)$ determine $\bm{\Phi}(\Delta,t)$ entirely for an isotropic system. For anisotropic systems, however, $\Phi_{1}(\Delta,t)$ and $\Phi_{2}(\Delta,t)$ may represent a small part of the information contained in $\bm{\Phi}(\Delta,t)$. For example, if the spatial dimension $n$ is $3$, as many as 21 elements of $\bm{\Phi}(\Delta,t)$ are independent.

III.6 Crossover

As seen from Eqs. (27) and (28), the correlation function $\phi_{i}(\tau)\,\,(i=1,2)$ satisfies $\lim_{\tau\to\infty}\phi_{i}(\tau)=0$. If $\phi_{i}(\tau)$ has a characteristic time scale $\tau_{i}$, then, from Eqs. (26) and (36), we have

[TABLE]

Thus, at the characteristic time scale $\tau_{i}$, $\hat{\Phi}_{i}^{\mathrm{ex}}(t)$ shows a crossover. For the polymer motion studied here, this crossover time $\tau_{i}$ corresponds roughly to the longest relaxation time of each polymer model as shown in the subsequent sections.

IV Rouse model

In this and the following three sections, the method of the TMSD tensor developed in the previous section is applied to the four polymer models stated in the Introduction. Here, the Rouse model is studied as the first example; although this is a very simple model of a flexible polymer chain in dilute solutions, it is the basis of many mathematical models of biopolymers Weber et al. (2010); Gong and van der Maarel (2014); Shinkai et al. (2016).

The Rouse model is composed of $N$ equivalent beads, the dynamics of which are subject neither to the excluded-volume nor hydrodynamic interaction Rouse (1953); Doi and Edwards (1986):

[TABLE]

where $\bm{R}_{n}(t)$ is the position of bead $n$, $k$ is the spring constant, and $\zeta$ is the friction coefficient. The spring constant $k$ is related to the mean bond length $b$ as $k=3k_{B}T/b^{2}$. The random force $\bm{f}_{n}(t)$ satisfies $\left\langle\bm{f}_{n}(t)\right\rangle=0$ and the fluctuation-dissipation relation $\left\langle\bm{f}_{n}(t)\bm{f}_{m}(t^{\prime})\right\rangle=2\zeta k_{B}T\delta_{nm}\delta(t-t^{\prime})\bm{I}$.

The equation of motion for the COM $\bm{R}_{G}(t):=\sum_{n=1}^{N}\bm{R}_{n}(t)/N$ is given by

[TABLE]

where $D_{G}=k_{B}T/N\zeta$ is the diffusion coefficient of the COM. Comparing with Eq. (2), we have $\bm{B}(t)=\sqrt{D_{G}}\bm{I}$. Because the diffusion coefficient $D_{G}$ is independent of time $t$, we have $\phi_{1}(t)\equiv\phi_{2}(t)\equiv 0$ from Eqs. (27) and (28). Consequently, the excess parts also vanish, namely $\hat{\Phi}_{1}^{\mathrm{ex}}(t)\equiv\hat{\Phi}_{2}^{\mathrm{ex}}(t)\equiv 0$, and, from Eqs. (25) and (35), the ideal parts are given by

[TABLE]

Note that the ideal parts decay simply as $1/t$ and do not exhibit crossover.

In Fig. 1, these formulas [Eqs. (45) and (46)] are displayed by the solid lines, and results of the numerical simulations by the circles and the triangles; the theoretical curves are in excellent agreement with the simulation results. These numerical results were obtained from trajectory data $\bm{R}_{G}(t)$ that were generated through Brownian dynamics simulations of the Rouse model [Eq. (43)].

V Zimm model

In this section, we study the Zimm model without the excluded volume interaction (i.e., the Zimm model in the $\Theta$ condition). Some scaling properties of the Rouse model are known to be inconsistent with experiments Doi and Edwards (1986), which is because the hydrodynamic interaction is disregarded entirely in the Rouse model. In contrast, the hydrodynamic interaction is taken into account in the Zimm model, which is another model of a flexible polymer chain in dilute solutions.

V.1 Model definition

As in the case of the Rouse model, the Zimm model consists of $N$ equivalent beads, and the equation of motion for bead $n$ is given by Zimm (1956); Doi and Edwards (1986); Ermak and McCammon (1978)

[TABLE]

where the hydrodynamic interaction is represented in terms of the mobility matrix $\bm{H}_{nm}$ defined by

[TABLE]

Here, $\eta$ is the viscosity of the solvent and $a$ is the radius of each bead. Moreover, $\bm{r}_{nm}(t)$ and $r_{nm}(t)$ are defined as $\bm{r}_{nm}(t):=\bm{R}_{n}(t)-\bm{R}_{m}(t)$ and $r_{nm}(t):=|\bm{r}_{nm}(t)|$, respectively. The thermal noise $\bm{f}_{m}(t)$ satisfies the fluctuation-dissipation relation

[TABLE]

The non-diagonal elements $\bm{H}_{nm}$ [Eq. (49)] are known collectively as the Oseen tensor, the nonlinearity of which makes theoretical analysis of the Zimm model considerably difficult.

A simple approximation that is commonly adopted is a pre-averaging approximation Doi and Edwards (1986) in which $\bm{H}_{nm}$ is replaced with its equilibrium average $\left\langle\bm{H}_{nm}\right\rangle=:h(n-m)\bm{I}$. In this approximation, the equation of motion for bead $n$ is expressed as

[TABLE]

Although this approximation works well for predicting the MSD of the COM motion Rey et al. (1989), it is impossible to use it to elucidate the fluctuating diffusivity. This is because the fluctuating diffusivity is disregarded entirely when replacing $\bm{H}_{nm}(t)$ in Eq. (50) with $h(n-m)\bm{I}$ [see Eq. (52)].

V.2 Equation of COM motion

To elucidate the effect of the fluctuating diffusivity, the pre-averaging approximation is applied to the internal modes only, whereas the COM motion is treated without pre-averaging.

The normal mode $\bm{X}_{p}(t)$ ($p=0,1,\dots$) of $\bm{R}_{n}(t)$ is defined by Doi and Edwards (1986)

[TABLE]

Note here that $\bm{X}_{0}(t)$ is equivalent to the COM position: $\bm{X}_{0}(t)=\bm{R}_{G}(t)$. Under the pre-averaging approximation [Eqs. (51) and (52)], the equations of motion for the normal modes are given by

[TABLE]

where $\hat{\bm{f}}_{p}(t)$ are random forces defined by

[TABLE]

with

[TABLE]

For $p\neq 0$, $\hat{h}_{pq}$ can be approximated further as $\hat{h}_{pq}\approx{\delta_{pq}}/{\zeta_{p}}$, where $\zeta_{p}:=(12\pi^{3}pN)^{1/2}\eta b$ Doi and Edwards (1986). Consequently, the Langevin equations for the internal modes, Eq. (56), are mutually independent because

[TABLE]

Moreover, in Eq. (56), $\tau_{p}$ is the relaxation time of the $p$-th mode, and given by

[TABLE]

where $\tau_{1}$ is the longest relaxation time.

Here, the COM equation of motion in Eq. (55) is rewritten as

[TABLE]

where $\bm{\xi}(t)$ is the white Gaussian noise given by Eq. (3). By comparing Eq. (61) with Eqs. (55), (57), and (58), $\bm{B}(t)$ is given by

[TABLE]

where we restored the time dependence of the diffusivity by formally replacing $h(n-m)\bm{I}$ with $\bm{H}_{nm}(t)$. In the following analysis, Eqs. (56) and (59) are used for the internal modes, whereas Eqs. (61) and (62) are used for the COM motion. Thus, the diffusion coefficient of the Zimm model, in contrast to that of the Rouse model, depends on time $t$ and fluctuates because of the hydrodynamic interaction.

From Eqs. (49) and (62), we have the ensemble average of the diffusion coefficient tensor as

[TABLE]

where $c:={k_{B}T}/(2\pi\eta N^{2})$ is a constant and we used the mutual independence of the magnitude $r_{nm}$ and direction $\bm{r}_{nm}/r_{nm}$ as follows Doi and Edwards (1986):

[TABLE]

From Eq. (63), we have

[TABLE]

The validity of Eq. (65) has been studied intensively Liu and Dünweg (2003), and it is shown that Eq. (65) is equivalent to the short-time diffusion coefficient of the COM and that it is also a good approximation to the long-time diffusion coefficient. In the next subsection, however, we have to study the second moment of the diffusion coefficient $\bm{D}(t)$.

Here, $\bm{r}_{nm}(t)$ follows three-dimensional Gaussian distribution with a covariant matrix $\bm{\Sigma}_{3}=\bm{I}|n-m|b^{2}/3$,

[TABLE]

Thus, $\left\langle 1/r_{nm}(t)\right\rangle$ is obtained by integrating over $f_{3}(\bm{r})$ in spherical coordinates as Doi and Edwards (1986)

[TABLE]

From Eqs. (63) and (68), we have an explicit expression of the ensemble-averaged diffusivity,

[TABLE]

It follows that

[TABLE]

V.3 Correlation functions of diffusion coefficient

In this subsection, we calculate the magnitude and orientation correlation functions $\phi_{1}(t)$ and $\phi_{2}(t)$, respectively, of the diffusivity [Eqs. (27) and (28)]. In the following derivation, we use crude approximations such as a single-mode approximation and a perturbation expansion of the Gaussian distribution. Nevertheless, the final results exhibit relatively good agreement with those of numerical simulations.

V.3.1 Magnitude correlation function of diffusion coefficient

We begin by deriving the magnitude correlation function $\phi_{1}(t)$ of the diffusivity. From Eqs. (49) and (62), we have

[TABLE]

To evaluate the ensemble average in the integrand, we define a six-dimensional vector $\bm{X}:=(x,x^{\prime},y,y^{\prime},z,z^{\prime})$, where $(x,y,z):=\bm{r}_{nm}(t)$ and $(x^{\prime},y^{\prime},z^{\prime}):=\bm{r}_{n^{\prime}m^{\prime}}(0)$. It can be shown that $\bm{X}$ follows six-dimensional Gaussian distribution (see Appendix B for a derivation), namely

[TABLE]

Here, $\bm{\Sigma}_{6}$ is a $6\times 6$ covariant matrix defined by

[TABLE]

where $\bm{0}$ is the $2\times 2$ zero matrix; $\alpha,\alpha^{\prime}$ and $\beta$ are defined by (see Appendix B)

[TABLE]

Hereinafter, we take only the longest relaxation mode ($p=1$) into account and ignore all the other modes (i.e., a single-mode approximation):

[TABLE]

Consequently, the determinant of the covariant matrix $\bm{\Sigma}_{6}$ is given by

[TABLE]

where $\epsilon:=\beta^{2}/(\alpha\alpha^{\prime})$ and $\tilde{\alpha}:=\alpha(1-\epsilon)$.

Using these quantities in Eq. (73), we have

[TABLE]

where $\bm{r}:=(x,y,z)$, and $\bm{r}^{\prime}:=(x^{\prime},y^{\prime},z^{\prime})$. For $t\gg\tau_{1}$, we have $\epsilon\ll 1$ and the above integrand can be approximated further as

[TABLE]

Integrating Eq. (79) in spherical coordinates, we have a perturbation expansion upto order $\epsilon^{1}$ as

[TABLE]

where we used Eqs. (68) and (75). Inserting this equation into Eq. (V.3.1) and taking Eq. (65) into account, we obtain

[TABLE]

where $K$ is a constant defined by

[TABLE]

Finally, from Eqs. (70) and (V.3.1), we have the magnitude correlation function $\phi_{1}(t)$ of the diffusivity [Eq. (27)] as

[TABLE]

V.3.2 Orientation correlation function of diffusion coefficient

We move on to a derivation of the orientation correlation function $\phi_{2}(t)$ of the diffusivity [Eq. (28)]. From Eqs. (49) and (62), we have

[TABLE]

where $\hat{\bm{r}}$ is the unit vector in the direction of $\bm{r}$. The ensemble average in Eq. (85) can be carried out in a way similar to the calculation of Eq. (79). In fact, under the approximation in Eq. (80), we obtain

[TABLE]

Inserting Eq. (86) into Eq. (85) and taking Eqs. (65) and (66) into account, we have

[TABLE]

Finally, from Eqs. (71) and (87), we have the orientation correlation function $\phi_{2}(t)$ of the diffusivity [Eq. (28)] as

[TABLE]

V.4 Correlation functions of TMSD tensor

Here, we derive the correlation functions $\hat{\Phi}_{1}(\Delta,t)$ and $\hat{\Phi}_{2}(\Delta,t)$ of the TMSD tensor. From Eqs. (25), (26), (84), and (88), we have

[TABLE]

where we used $C=1$ because the system is statistically isotropic. Similarly, from Eqs. (35), (36), (84), and (88), we obtain

[TABLE]

In contrast to the Rouse model, these correlation functions $\hat{\Phi}_{i}^{\mathrm{ex}}(t)$ $(i=1,2)$ for the Zimm model show crossovers. For example, from Eq. (42), $\hat{\Phi}_{1}^{\mathrm{ex}}(t)$ behaves as

[TABLE]

From Eq. (93), the crossover time $\tau_{c}$ can be estimated as $\tau_{c}=\tau_{1}$, i.e., the crossover time is equivalent to the longest relaxation time $\tau_{1}$. Also, $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ shows a crossover at $t=\tau_{1}$, because of Eq. (92).

As can be seen in Fig. 2, the theoretical predictions (the solid lines) [Eqs. (89)–(92)] are in good agreement with the results of the numerical simulations (the symbols). The slight deviations are due to the approximations used in the theoretical analysis. For example, in the simulations, the Rotne-Prager-Yamakawa tensor [Eq. (132)] was utilized as the mobility matrix instead of the Oseen tensor [Eq. (49)] to regularize the singularity in the Oseen tensor at $r_{nm}=0$. Moreover, we also applied the pre-averaging approximation to the inner degrees of freedom, and used the perturbation expansion in Eq. (80). However, incorporating a higher order term ($\epsilon^{2}$) in Eq. (80) improves the theoretical predictions only slightly (its contribution is less than 15 % of the leading term; data not shown).

VI Discrete reptation model

In this section, the focus is on the discrete reptation model, which describes tagged polymer motion in entangled polymer solutions Doi and Edwards (1978, 1986). Because of the entanglement, the tagged polymer chain of the reptation model is temporarily trapped in a virtual tube comprised of surrounding chains, and moves only in the longitudinal direction of the tube. Such reptation dynamics are an essential ingredient in modeling DNA molecules at high concentration Gong and van der Maarel (2014).

In the reptation model, the centerline of the tube, which is called a primitive chain, is considered instead of the real chain of the tagged polymer. The primitive chain is assumed to consist of $N$ tube segments $\bm{R}_{1}(t),\dots,\bm{R}_{N}(t)$ connected by bonds of constant length $b$. The primitive chain is allowed to move only in the longitudinal direction of the tube as a result of the entanglement. A single step of the primitive-chain dynamics is given as follows; one of the two end segments, $\bm{R}_{1}(t)$ or $\bm{R}_{N}(t)$, is chosen with equal probability; the chosen end segment hops with step length $b$ in a random direction; and each of the other segments slides to one of the positions of its neighboring segments [i.e., if $\bm{R}_{1}(t)$ is chosen, $\bm{R}_{n}(t)$ slides to $\bm{R}_{n-1}(t)$ ($n=2,\dots,N$); if $\bm{R}_{N}(t)$ is chosen, $\bm{R}_{n}(t)$ slides to $\bm{R}_{n+1}(t)$ ($n=1,\dots,N-1$)].

The COM $\bm{R}_{G}(t)$ of this primitive chain follows the LEFD [Eq. (2)] with $\bm{B}(t)$ given by Doi and Edwards (1978); Uneyama et al. (2015)

[TABLE]

where $D_{G}$ is the ensemble-averaged diffusion coefficient of the COM, and $\bm{p}(t)$ is the end-to-end vector of the primitive chain. It follows that the diffusion coefficient is obtained from Eq. (4) as

[TABLE]

Because the system is statistically isotropic, $\left\langle\bm{p}\bm{p}\right\rangle=A\bm{I}$ with a constant $A$. Taking the trace, we have $A=\left\langle\bm{p}^{2}\right\rangle/3$. It follows that $\left\langle\bm{D}\right\rangle=D_{G}\bm{I}$.

By using Eqs. (27), (28) and (95), the magnitude and the orientation correlation functions $\phi_{1}(\tau)$ and $\phi_{2}(\tau)$ of the diffusivity can be expressed as

[TABLE]

In Ref. Uneyama et al. (2015), $\phi_{1}(\tau)$ was obtained explicitly as

[TABLE]

where $\tau_{d}$ is the longest relaxation time of the reptation model, and $E_{m}(x)$ is the generalized exponential integral of order $m$ Olver et al. (2010). Furthermore, it is shown in Appendix D that

[TABLE]

From Eqs. (25) and (26), we have the correlation functions of the TMSD tensor as

[TABLE]

where we used $\phi_{2}(0)=6\phi_{1}(0)=4$. Similarly, from Eqs. (35), (36) and (99), we have

[TABLE]

As in the case of the Zimm model, both functions $\hat{\Phi}_{i}^{\mathrm{ex}}(t)$ $(i=1,2)$ show crossovers. For example, $\hat{\Phi}_{1}^{\mathrm{ex}}(t)$ behaves as Uneyama et al. (2015)

[TABLE]

Also, $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ shows a crossover at $t=\tau_{d}$, because of Eq. (103). From Eq. (104), this crossover time can be estimated as

[TABLE]

which is close to the longest relaxation time $\tau_{d}$.

In Fig. 3, results of the numerical simulations for the discrete reptation model are displayed; they exhibit remarkable agreement with the theoretical predictions [Eqs. (100)–(103)]. Moreover, $\hat{\Phi}_{2}^{\mathrm{ex}}(\Delta,t)$ far exceeds $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)$ in the reptation model, in contrast to the Zimm model for which the two functions are comparable. This means that the orientation fluctuation of the diffusivity is more prominent in the reptation model than in the Zimm model.

VII Rigid rod-like polymer

Finally, the rigid rod-like polymer in a dilute condition is investigated as an extreme example of non-flexible polymers Burgers (1995); Doi and Edwards (1986); Dhont (1996); Berne and Pecora (2000). In general, it is more difficult to observe rotational diffusion of an anisotropic particle than it is to observe its translational diffusion Han et al. (2006). With the TMSD tensor analysis, however, the rotational diffusion coefficient can be estimated by measuring translational motion of the COM.

Let us denote the COM of the rod as $\bm{R}_{G}(t)$, and assume that the rod is cylindrically symmetric along the long axis. Consequently, the COM position $\bm{R}_{G}(t)$ follows the LEFD [Eq. (2)] with $\bm{B}(t)$ given (see Appendix E) by

[TABLE]

where $\hat{\bm{u}}(t)$ is a unit vector in the direction of the rod’s long axis, and $D_{\parallel}$ and $D_{\perp}$ are the diffusion coefficients along and perpendicular to the long axis, respectively. Moreover, it is assumed that the rod is long and thin so that rotational motion around the long axis is disregarded. The time evolution of the rod’s direction $\hat{\bm{u}}(t)$ is given by Dhont (1996)

[TABLE]

where $\bm{\eta}(t)$ is white Gaussian noise, and $D_{r}$ is the rotational diffusion coefficient. The three diffusion coefficients $D_{\parallel}$, $D_{\perp}$ and $D_{r}$ can be expressed in terms of the length $L$ and diameter $b$ of the rod as Doi and Edwards (1986); Dhont (1996)

[TABLE]

These formulas are obtained through hydrodynamic calculations for a long thin rod, i.e., $L/b\gg 1$.

Firstly, we consider the magnitude correlation function $\hat{\Phi}_{1}(t)$. From Eqs. (4) and (106), we have the fluctuating diffusivity as

[TABLE]

Taking the trace, we obtain $\mathrm{tr}\bm{D}(t)=D_{\parallel}+2D_{\perp}$, i.e., the magnitude of the diffusivity is constant in time. It follows that the magnitude correlation of the diffusivity vanishes, i.e., $\phi_{1}(\tau)\equiv 0$; hence we have from Eq. (26) that

[TABLE]

Thus, for the rigid rod-like polymer, in contrast to the Zimm and reptation models, it is impossible to extract information about the fluctuating diffusivity by using $\hat{\Phi}_{1}(t)$ .

Therefore, to elucidate the fluctuating diffusivity of the rod, it is necessary to study $\hat{\Phi}_{2}(t)$. From Eqs. (28) and (111), we obtain the orientation correlation function as

[TABLE]

where we used $\left\langle\bm{D}\right\rangle=(D_{\parallel}+2D_{\perp})\bm{I}/3$. Because the rotational motion given by Eq. (107) is independent of the translational motion $\bm{R}_{G}(t)$, the correlation function $\langle[\hat{\bm{u}}(\tau)\cdot\hat{\bm{u}}(0)]^{2}\rangle$ can be calculated by employing the Smoluchowsky equation for the rotational motion as Doi and Edwards (1986); Berne and Pecora (2000)

[TABLE]

and hence we have

[TABLE]

From Eq. (36), the excess part $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ of the orientation correlation function is obtained as

[TABLE]

Moreover, by using Eqs. (25) and (35), the ideal parts are given by

[TABLE]

where $C=1$ is used because the system is statistically isotropic. In particular, from Eq (42), we have the following crossover:

[TABLE]

An estimate for the rotational relaxation time $1/D_{r}$ can be obtained from this crossover time, despite the fact that we observe only the translational motion of the rod. In fact, we have the crossover time $\tau_{c}$ from Eq. (119) as

[TABLE]

Thus, the crossover time $\tau_{c}$ gives an estimate of the rotational relaxation time $1/D_{r}$.

Results of the numerical simulations for $\hat{\Phi}_{1}(t)$ and $\hat{\Phi}_{2}(t)$ are presented in Fig. 4 (the circles and triangles). As predicted, $\hat{\Phi}_{1}(t)$ shows no crossover because the excess part is absent, whereas $\hat{\Phi}_{2}(t)$ exhibits a clear crossover. The numerical results are consistent with the theoretical predictions (the solid lines).

VIII Discussion

The sample-to-sample fluctuation of the diffusivity observed both in SPT experiments and theoretical models has been studied intensively for a decade. In such studies, the sample-to-sample fluctuation is usually characterized by the EB parameter He et al. (2008); Deng and Barkai (2009); Jeon et al. (2011); Tabei et al. (2013); Metzler et al. (2014); Miyaguchi and Akimoto (2011, 2013); Thiel and Sokolov (2014a, b); Uneyama et al. (2012, 2015); Miyaguchi and Akimoto (2015); Miyaguchi et al. (2016). However, when calculating the EB parameter from trajectory data $\bm{r}(t)$, much of the information originally contained in the data is lost. In this study, to obtain more information from the trajectory data, the EB parameter is generalized into the fourth-order tensor $\bm{\Phi}(\Delta,t)$, which is a correlation function of the TMSD tensor. Moreover, the two scalar functions $\hat{\Phi}_{1}(\Delta,t)$ and $\hat{\Phi}_{2}(\Delta,t)$ are derived from $\bm{\Phi}(\Delta,t)$; these functions are closely related to the magnitude and orientation correlation functions of the diffusivity, and in particular $\hat{\Phi}_{1}(\Delta,t)$ is equivalent to the EB parameter. It is also worth noting that a linear combination of the excess parts $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(\Delta,t)$ gives the non-Gaussian parameter $A(t)$ [Eq. (38)]. In other words, the non-Gaussianity can be decomposed into two parts: one originating from the magnitude fluctuation of the diffusivity, and the other from the orientation fluctuation.

Furthermore, by using the TMSD tensor analysis, it is shown that the four polymer models exhibit distinctly different types of fluctuating diffusivity in terms of the correlation functions $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(\Delta,t)$. For example, $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)\approx\hat{\Phi}_{2}^{\mathrm{ex}}(\Delta,t)$ in the Zimm model, $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)<\hat{\Phi}_{2}^{\mathrm{ex}}(\Delta,t)$ in the reptation model, and $\hat{\Phi}_{1}^{\mathrm{ex}}(\Delta,t)\ll\hat{\Phi}_{2}^{\mathrm{ex}}(\Delta,t)$ in the rigid rod-like polymer. This is in contrast to the non-Gaussian parameter $A(t)$, whose behavior is qualitatively similar for these three models; hence the polymer models are barely distinguishable with $A(t)$.

From these results, it seems that the fluctuating diffusivity might be ubiquitous in polymer motions from dilute to concentrated solutions and from flexible to non-flexible polymers. This is because the Zimm and the reptation models are flexible polymer models in dilute and concentrated solutions, respectively; in contrast, the rigid rod-like polymer is an extreme case of non-flexible polymers; each of these three models exhibits fluctuating diffusivity.

Moreover, the rotational relaxation time $\tau_{r}=1/D_{r}$ of the rigid rod can be obtained from the crossover time of $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ [Eq. (119)]. As a more direct approach, $\tau_{r}=1/D_{r}$ of an anisotropic particle was obtained in Ref. Han et al. (2006) by measuring the particle’s direction. Also, with the results of Refs. Brenner (1963); Cichocki et al. (2012), $\tau_{r}$ of the rigid rod can be estimated from the ensemble-averaged MSD of a reference point on the rod other than its COM. For both methods, however, it is necessary to measure at least one reference point other than the COM. In contrast, with the method proposed here, $\tau_{r}$ can be estimated by measuring only the translational motion of the COM.

Of course, the same information of $\hat{\Phi}_{1}^{\mathrm{ex}}(t)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ would be obtained from the ensemble-averaged quantities. In fact, the functions $\hat{\Phi}_{1}^{\mathrm{ex}}(t)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ are related to the non-Gaussian parameter $A(t)$ [Eq. (38)], which is defied by a fourth moment. Thus, essentially the same information as $\hat{\Phi}_{1}^{\mathrm{ex}}(t)$ and $\hat{\Phi}_{2}^{\mathrm{ex}}(t)$ might well be obtained from the translational correlation tensor of fourth order, which might be analyzed by the traditional approach with the Smoluchowski equation Rallison (1978); Cichocki et al. (2012, 2015). However, It should be noted that to calculate fourth moments such as $A(t)$ accurately, a large number of trajectories are necessary in general. In contrast, the present method, in which the time and ensemble averages are combined, works for a relatively small number of trajectories (typically, from tens to hundreds of trajectories), and therefore it would be useful in single-particle-tracking experiments, in which much effort is required to obtain a large number of trajectories.

Although the TMSD tensor analysis for the polymer models is based on the fact that the COM of these models can be described in terms of the LEFD [Eq. (2)], there are many phenomena that cannot be described with the LEFD. For example, the motion of a single bead in the Zimm and reptation models does not follow the LEFD because the bead shows anomalous subdiffusion, whereas the LEFD exhibits only normal diffusion as shown in Eq. (7). A candidate for describing such complex dynamics might be a generalized Langevin equation or fractional Brownian motion with fluctuating diffusivity, but the physical validity of such models should be clarified in future work.

Moreover, only two scalar functions, namely $\hat{\Phi}_{1}(\Delta,t)$ and $\hat{\Phi}_{2}(\Delta,t)$, were used here to analyze the isotropic polymer models. However, there must still be useful information in the fourth-order tensor $\bm{\Phi}(\Delta,t)$ for the case of anisotropic systems (see Sec. III.5). Future work should therefore include a full characterization of this tensor $\bm{\Phi}(\Delta,t)$.

Acknowledgements.

The author would like to thank T. Akimoto and T. Uneyama for fruitful discussions and comments. This work was supported by JSPS KAKENHI for Young Scientists (B) (Grant No. JP15K17590).

Appendix A Decomposition of fourth-order tensor $\bm{\Phi}(\Delta,t)$ into ideal and excess parts

In this appendix, the expression for $\bm{\Psi}^{1}(\Delta,t)$ given in Eq. (III.2) is derived. First, using Eqs. (2) and (5), we obtain

[TABLE]

where $t-\Delta$ is approximated as $t-\Delta\simeq t$, and $\bm{\Omega}(\Delta,t^{\prime},t^{\prime\prime})$ is another fourth-order tensor defined by

[TABLE]

By using the Heaviside step function $\Theta(t)$ and Wick’s theorem, namely

[TABLE]

the elements of $\bm{\Omega}(\Delta,t^{\prime},t^{\prime\prime})$ for $t^{\prime\prime}<t^{\prime}$ is obtained as

[TABLE]

where approximations such as $D_{ik}(s)\approx D_{ik}(t^{\prime})$ for $s\in[t^{\prime},t^{\prime}+\Delta)$ are applied; these approximations are justified by the assumption that $\Delta$ is much shorter than a characteristic time scale of the fluctuating diffusivity $\bm{D}(t)$. In addition, an expression similar to Eq. (124) can be obtained also for $t^{\prime}<t^{\prime\prime}$. By putting these equations into Eq. (121) and using the stationarity, the elements of $\bm{\Psi}^{1}(\Delta,t)$ can be expressed as Eq. (III.2).

Appendix B Derivation of six-dimensional covariant matrix $\bm{\Sigma}_{6}$ for Zimm model

Here, the covariant matrix $\bm{\Sigma}_{6}$ of the six-dimensional Gaussian distribution in Eq. (73) is derived. Firstly, let us denote a transition probability density function (PDF) for the normal mode $\bm{X}_{p}(t)$ as $P_{p}(\bm{X}_{p},t|\bm{X}^{\prime}_{p},0)$; more precisely, $P_{p}(\bm{X}_{p},t|\bm{X}^{\prime}_{p},0)d\bm{X}_{p}$ is the transition probability from $\bm{X}^{\prime}_{p}$ at time [math] to an interval $[\bm{X}_{p},\bm{X}_{p}+d\bm{X}_{p})$ at time $t$. From Eqs. (56) and (59), the PDF for $\bm{X}_{p}(t)\,(p=1,2,\dots)$ is given by

[TABLE]

where $\sigma_{p}^{2}(t)=(k_{B}T\tau_{p}/\zeta_{p})(1-e^{-2t/\tau_{p}})$ is the variance. For example, Eq. (125) can be derived by using Chandrasekhar’s theorem Dhont (1996). In particular, by taking the limit $t\to\infty$, the equilibrium PDF for $\bm{X}_{p}$ is obtained as

[TABLE]

where $\sigma_{p}^{2}=k_{B}T\tau_{p}/\zeta_{p}$.

By using these PDFs, the joint PDF of $\bm{r}_{nm}(t)$ and $\bm{r}_{n^{\prime}m^{\prime}}(0)$ is expressed as

[TABLE]

where $P(\bm{r},t;\bm{r}^{\prime},0)d\bm{r}d\bm{r}^{\prime}$ is the joint probability that $\bm{r}_{nm}(t)$ is in $[\bm{r},\bm{r}+d\bm{r})$ and $\bm{r}_{n^{\prime}m^{\prime}}(0)$ is in $[\bm{r}^{\prime},\bm{r}^{\prime}+d\bm{r}^{\prime})$. Note here that $\bm{r}_{nm}(t)$ and $\bm{r}_{n^{\prime}m^{\prime}}(0)$ are written as $\bm{r}_{nm}(t)=\sum_{p=1}^{\infty}c_{p}\bm{X}_{p}(t)$ and $\bm{r}_{n^{\prime}m^{\prime}}(0)=\sum_{p=1}^{\infty}c^{\prime}_{p}\bm{X}_{p}(0)$ with

[TABLE]

because of Eq. (54) and $\bm{r}_{nm}(t)=\bm{R}_{n}(t)-\bm{R}_{m}(t)$.

With Fourier transformation of Eq. (127) with respect to $\bm{r}$ and $\bm{r}^{\prime}$, we have a characteristic function

[TABLE]

where $\alpha$, $\alpha^{\prime}$, and $\beta$ are defined in Eqs. (75) and (76). Moreover, $\bm{k}$ and $\bm{k}^{\prime}$ are the Fourier variables conjugate to $\bm{r}$ and $\bm{r}^{\prime}$, respectively; their elements are defined as $\bm{k}:=(k_{x},k_{y},k_{z})$ and $\bm{k}^{\prime}:=(k^{\prime}_{x},k^{\prime}_{y},k^{\prime}_{z})$. To derive Eq. (130), we used a Fourier series $\sum_{p=1}^{\infty}\cos(px)/p^{2}=(x-\pi)^{2}/4-\pi^{2}/12$ for $x\in[0,2\pi]$. If we define a variable $\bm{K}$ as $\bm{K}:=(k_{x},k_{x}^{\prime},k_{y},k_{y}^{\prime},k_{z},k_{z}^{\prime})$, the right-hand side of Eq. (130) can be rewritten as

[TABLE]

This is a characteristic function of six-dimensional Gaussian distribution; consequently, Fourier inversion of $\hat{P}(\bm{K})$ gives Eq. (73).

Appendix C Rotne-Prager-Yamakawa tensor

To carry out numerical simulations of the Zimm model, it is necessary to regularize the singularity of the Oseen tensor at $r_{nm}=0$ [Eq. (49)]. A commonly employed regularization method is the Rotne–Prager–Yamakawa tensor $\tilde{\bm{H}}_{nm}\,(n\neq m)$ Rotne and Prager (1969); Yamakawa (1970):

[TABLE]

where $a$ is the bead radius. In our numerical simulations for the Zimm model, $\tilde{\bm{H}}_{nm}$ was used for the mobility matrix $\bm{H}_{nm}$ in Eq. (47). The Langevin equation [Eq. (47)] was solved numerically by using the Ermak–McCammon algorithm Ermak and McCammon (1978).

Appendix D Derivation of correlation functions for discrete reptation model

In this Appendix, the relation presented in Eq. (99) is derived for the discrete reptation model. The end-to-end vector $\bm{p}(t)$ of the reptation model can be expressed with a bond vector $\bm{u}(s,t)$ as

[TABLE]

where $s$ is the segment index and $N$ is the number of segments. The bond vector $\bm{u}(s,t)$ follows a Gaussian distribution with zero mean, and any two bond vectors $\bm{u}(s,t)$ and $\bm{u}(s^{\prime},t)$ are mutually independent. Thus, the first and second moments of $\bm{u}(s,t)$ in equilibrium are given by

[TABLE]

where $b$ is the bond length.

To derive an explicit formula for $\phi_{1}(\tau)$ and $\phi_{2}(\tau)$ [Eqs. (96) and (97)], we use the survival probability $\Psi(s;t)$ of segment $s$; more precisely, $\Psi(s;t)$ is the probability that segment $s$ at time [math] survives until time $t$ Doi and Edwards (1978). Also, we define a survival joint probability $\Psi(s,s^{\prime};t)$ of two segments $s$ and $s^{\prime}$ Uneyama et al. (2015). Namely, $\Psi(s,s^{\prime};t)$ is the probability that both segments $s$ and $s^{\prime}$ at time [math] survive until time $t$. In particular, $\Psi(s,s;t)=\Psi(s;t)$ is satisfied. Although an explicit expression for $\Psi(s,s^{\prime};t)$ was derived in Ref.Uneyama et al. (2015), it is not required here.

Correlation functions of the end-to-end vector $\bm{p}(t)$ can be expressed with $\bm{u}(s,t)$. For example, a fourth-order correlation function (tensor) of $\bm{p}(t)$ is written as

[TABLE]

The elements of the tensor in the integrand can be rewritten as

[TABLE]

where $\Psi(s^{\prime};t)-\Psi(s,s^{\prime};t)$ is the probability that only segment $s^{\prime}$ survives, and $1-\Psi(s;t)-\Psi(s^{\prime};t)+\Psi(s,s^{\prime};t)$ is the probability that neither of segments $s$ and $s^{\prime}$ survive. By using Eq. (134), the second and third terms on the right-hand side vanish. Meanwhile, the ensemble averages in the first and fourth terms can be rewritten as

[TABLE]

where we used Wick’s theorem Doi and Edwards (1986) and Eq. (134). Putting Eqs. (136), (137) and (138) into Eq. (135), we have

[TABLE]

where we used $\Psi(s,s;t)=\Psi(s;t)$.

Taking contractions in Eq. (139) between the first and second indices, and also between the third and fourth indices, we obtain

[TABLE]

where we used $\left\langle\bm{p}^{2}\right\rangle=b^{2}N$. Inserting Eq. (140) into Eq. (96), we obtain Uneyama et al. (2015)

[TABLE]

Similarly, taking contractions in Eq. (139) between the first and fourth indices, and also between the second and third indices, we obtain

[TABLE]

By inserting Eq. (142) into Eq. (97), $\phi_{2}(\tau)$ can be expressed as

[TABLE]

By comparing Eqs. (141) and (143), we obtain $\phi_{2}(\tau)=6\phi_{1}(\tau)$ [Eq. (99)].

Appendix E Langevin equation of COM motion for rigid rod-like polymer

Here, Eq. (106) for the rigid rod-like polymer is derived. The overdamped COM motion of the rod is described as follows Dhont (1996):

[TABLE]

where $\bm{\Gamma}_{f}^{-1}$ is the inverse of the friction matrix, namely

[TABLE]

and $\zeta_{\parallel}$ and $\zeta_{\perp}$ are the friction coefficients parallel and perpendicular to the rod’s long axis, respectively.

Note that the thermal noise $\bm{f}(t,\hat{\bm{u}}(t))$ depends on the direction $\hat{\bm{u}}(t)$ of the rod. This noise term can be decomposed as

[TABLE]

where $\bm{f}_{\parallel}(t,\hat{\bm{u}}(t))$ and $\bm{f}_{\perp}(t,\hat{\bm{u}}(t))$ represent equilibrium thermal noise in the parallel and perpendicular directions of the rod:

[TABLE]

Here, $\bm{\xi}(t)$ is the three-dimensional white Gaussian noise defined in Eq. (3). Note that $\bm{\xi}(t)$ is independent of $\hat{\bm{u}}(t)$, in contrast to $\bm{f}(t,\hat{\bm{u}}(t))$ in Eq. (144). Inserting Eqs. (145)–(148) into Eq. (144) and using the Einstein relations $D_{\parallel}=k_{B}T/\zeta_{\parallel}$ and $D_{\perp}=k_{B}T/\zeta_{\perp}$, we have Eq. (106).