Effect of Inertia on Linear Viscoelasticity of Harmonic Dumbbell Model

TL;DR

This paper investigates how weak inertia influences the linear viscoelastic behavior of the harmonic dumbbell model, revealing kinetic coupling between momentum and bond relaxation modes and deriving an overdamped Langevin equation that accounts for inertia effects.

Contribution

It introduces a theoretical analysis showing inertia's impact on bond relaxation and derives an overdamped Langevin equation incorporating weak inertia effects.

Findings

Inertia can affect bond relaxation when momentum relaxation is slow.

Derived an overdamped Langevin equation with inertia effects.

Predicted bond relaxation dynamics accurately with weak inertia.

Abstract

The overdamped (inertialess) dumbbell model is widely utilized to study rheological properties of polymers or other soft matters. In most cases, the effect of inertia is merely neglected because the momentum relaxation is much faster than the bond relaxation. We theoretically analyze the effect of inertia on the linear viscoelasticity of the harmonic dumbbell model. We show that the momentum and bond relaxation modes are kinetically coupled and the inertia can affect the bond relaxation if the momentum relaxation is not sufficiently fast. We derive an overdamped Langevin equation for the dumbbell model, which incorporates the weak inertia effect. Our model predicts the bond relaxation dynamics with the weak inertia effect correctly. We discuss how the weak inertia affects the linear viscoelasticity of a simple harmonic dumbbell model and the Rouse model.

Equations

$$\mathcal{H}=\frac{1}{2M}({\mathbf{P}}_1^2+{\mathbf{P}}_2^2)+U({\mathbf{R}}_1-{\mathbf{R}}_2),$$

$$\frac{d{\mathbf{R}}_i(t)}{dt}$$

$$\frac{d{\mathbf{P}}_i(t)}{dt}$$

$$⟨{\mathbf{W}}_i(t)⟩=0,⟨{\mathbf{W}}_i(t){\mathbf{W}}_j(t^{\prime })⟩={\delta }_{ij}1\delta (t-t^{\prime }),$$

$$\mathcal{H}=\frac{1}{2\bar{M}}{\bar{\mathbf{P}}}^2+\frac{1}{2m}{\mathbf{p}}^2+U(\mathbf{r}),$$

$$\frac{d\bar{\mathbf{R}}(t)}{dt}$$

$$\frac{d\bar{\mathbf{P}}(t)}{dt}$$

$$\frac{d\mathbf{r}(t)}{dt}$$

$$\frac{d\mathbf{p}(t)}{dt}$$

$$\frac{d\bar{\mathbf{R}}(t)}{dt}$$

$$\frac{d\mathbf{r}(t)}{dt}$$

$$m\frac{d^2\mathbf{r}(t)}{d{t}^2}=-\frac{\partial U(\mathbf{r}(t))}{\partial \mathbf{r}(t)}-\zeta \frac{d\mathbf{r}(t)}{dt}+\sqrt{2\zeta k_{B}T}\mathbf{w}(t).$$

$$\frac{d\mathbf{r}(t)}{dt}={\left(1+{\tau }_m\frac{d}{dt}\right)}^{-1}\left[-\frac{1}{\zeta }\frac{\partial U(\mathbf{r}(t))}{\partial \mathbf{r}(t)}+\sqrt{\frac{2k_{B}T}{\zeta }}\mathbf{w}(t)\right],$$

$${\left(1+{\tau }_m\frac{d}{dt}\right)}^{-1}\left[\left(1+{\tau }_m\frac{d}{dt}\right)f(t)\right]=f(t).$$

$$\begin{array}{rl}\frac{d\mathbf{r}(t)}{dt}& \approx -\frac{1}{\zeta }\frac{\partial U(\mathbf{r}(t))}{\partial \mathbf{r}(t)}+\sqrt{\frac{2k_{B}T}{\zeta }}{\left(1+{\tau }_m\frac{d}{dt}\right)}^{-1}\mathbf{w}(t)\\ & +\frac{{\tau }_m}{\zeta }\left[\mathbf{K}(\mathbf{r}(t))\cdot \frac{d\mathbf{r}(t)}{dt}+\frac{k_{B}T}{\zeta }\frac{\partial }{\partial \mathbf{r}(t)}\cdot \mathbf{K}(\mathbf{r}(t))\right],\end{array}$$

$$\mathbf{K}(\mathbf{r})\equiv \frac{{\partial }^{2}U(\mathbf{r})}{\partial \mathbf{r}\partial \mathbf{r}},$$

$$\begin{array}{rl}\frac{d\mathbf{r}(t)}{dt}& \approx -\frac{\mathbf{\Lambda }(\mathbf{r}(t))}{\zeta }\cdot \frac{\partial U(\mathbf{r}(t))}{\partial \mathbf{r}(t)}+\frac{k_{B}T}{\zeta }\frac{\partial }{\partial \mathbf{r}(t)}\cdot \mathbf{\Lambda }(\mathbf{r}(t))\\ & +\sqrt{\frac{2k_{B}T}{\zeta }}\mathbf{\Lambda }(\mathbf{r}(t))\cdot {\left(1+{\tau }_m\frac{d}{dt}\right)}^{-1}\mathbf{w}(t),\end{array}$$

$$\mathbf{\Lambda }(\mathbf{r})\equiv 1+\frac{{\tau }_m}{\zeta }\mathbf{K}(\mathbf{r})=1+\frac{{\tau }_m}{\zeta }\frac{{\partial }^{2}U(\mathbf{r})}{\partial \mathbf{r}\partial \mathbf{r}}.$$

$$\frac{d\mathbf{r}(t)}{dt}=-\frac{\mathbf{\Lambda }(\mathbf{r}(t))}{\zeta }\cdot \frac{\partial U(\mathbf{r}(t))}{\partial \mathbf{r}(t)}+\frac{k_{B}T}{\zeta }\frac{\partial }{\partial \mathbf{r}(t)}\cdot \mathbf{\Lambda }(\mathbf{r}(t))+\sqrt{\frac{2k_{B}T}{\zeta }}\mathbf{B}(\mathbf{r}(t))\cdot \mathbf{w}(t),$$

$$G(t)=\nu k_{B}T{e}^{-2t/{\tau }_b},$$

$$\begin{array}{rl}\hat{\mathbf{\sigma }}(\{{\mathbf{R}}_i\},\{{\mathbf{P}}_i\})& \equiv k({\mathbf{R}}_2-{\mathbf{R}}_1)({\mathbf{R}}_2-{\mathbf{R}}_1)-\frac{1}{M}({\mathbf{P}}_1{\mathbf{P}}_1+{\mathbf{P}}_2{\mathbf{P}}_2)\\ & =k\mathbf{r}\mathbf{r}-\frac{1}{\bar{M}}\bar{\mathbf{P}}\bar{\mathbf{P}}-\frac{1}{m}\mathbf{p}\mathbf{p}.\end{array}$$

$$G(t)=\frac{\nu }{k_{B}T}⟨{\hat{\sigma }}_{xy}(t){\hat{\sigma }}_{xy}(0)⟩.$$

$$\begin{split}G(t)&=\frac{\nu}{k_{B}T}\bigg{[}\frac{1}{\bar{M}^{2}}\langle\bar{P}_{x}(t)\bar{P}_{x}(0)\rangle^{2}+\frac{1}{m^{2}}\langle p_{x}(t)p_{x}(0)\rangle^{2}+k^{2}\langle r_{x}(t)r_{x}(0)\rangle^{2}\\ &\qquad-\frac{k}{m}\left[\langle r_{x}(t)p_{x}(0)\rangle^{2}+\langle p_{x}(t)r_{x}(0)\rangle^{2}\right]\bigg{]}.\end{split}$$

$$⟨\bar{\mathbf{P}}(t)\bar{\mathbf{P}}(0)⟩=\bar{M}{k}_{B}T{e}^{-t/{\tau }_m}1,$$

$$\frac{d}{dt}\left[\begin{array}{c}\mathbf{r}(t)\\ \mathbf{p}(t)\end{array}\right]=-\left[\begin{array}{cc}0& -1/m\\ k& \zeta /m\end{array}\right]\cdot \left[\begin{array}{c}\mathbf{r}(t)\\ \mathbf{p}(t)\end{array}\right]+\left[\begin{array}{c}0\\ \sqrt{2\zeta k_{B}T}\mathbf{w}(t)\end{array}\right].$$

$$⟨\mathbf{r}(t)\mathbf{r}(0)⟩$$

$$⟨\mathbf{p}(t)\mathbf{p}(0)⟩$$

$$⟨\mathbf{r}(t)\mathbf{p}(0)⟩$$

$${\lambda }_{\pm }\equiv \frac{\zeta \pm \sqrt{{\zeta }^2-4km}}{2m}.$$

$${\lambda }_{\pm }=\frac{1\pm \sqrt{1-4{\tau }_m/{\tau }_b}}{2{\tau }_m},$$

$$G(t)=\nu k_{B}T\left[e^{-2t/{\tau }_m}+e^{-(1+\sqrt{1-4{\tau }_m/{\tau }_b})t/{\tau }_m}+e^{-(1-\sqrt{1-4{\tau }_m/{\tau }_b})t/{\tau }_m}\right].$$

$$G(t)\approx \nu k_{B}T{e}^{-2(1+{\tau }_m/{\tau }_b)t/{\tau }_b}.$$

$$\mathbf{K}(\mathbf{r})=k1,\mathbf{\Lambda }(\mathbf{r})=\left(1+\frac{{\tau }_m}{{\tau }_b}\right)1.$$

Full text

1 Introduction

To study molecular-level dynamics of polymers or other soft matters, simple Langevin equation models are useful. From the viewpoint of rheology, the dynamics of molecular-level microscopic degrees of freedom can be related to the macroscopic rheological functions such as the shear relaxation modulus. For example, simple models such as the Rouse model and the reptation model[1] can reproduce experimentally obtained shear relaxation modulus data, and provide the insights on the dynamics of polymers which cannot be directly observed by microscopes.

In many cases, the overdamped Langevin equations without the inertia terms are employed as the dynamic equations for polymers or soft matters. This treatment is reasonable because the dynamical modes related to the inertia (momentum relaxation modes) relax rapidly. If we are interested in relatively long-time scales, we can assume that the inertia effect is not important, and simply take the overdamped (or inertialess) limit. For example, the mean-square displacement (MSD) of a Brownian particle (which obeys the Ornstein-Uhlenbeck process[2]) shows ballistic behavior at the short-time scale. But except the short-time scale region, the MSD can be well approximated by a simple diffusion motion without inertia. However, the effect of the inertia on the dynamics is generally not so simple[3, 4]. The underdamped Langevin equations with the inertia terms are the second order differential equations whereas the overdamped Langevin equations are first order ones. The difference of the differential orders qualitatively affects some properties of the model. Therefore it is important to understand the effect of the inertia on the dynamics.

Here we limit ourselves to one of the simplest model, so-called the dumbbell model[5, 6]. Some researchers have investigated the effects of the inertia on rheological properties of the dumbbell model[7, 8, 9, 10, 11]. Roughly speaking, the overdamped Langevin equations are justified when the inertia effect is very weak. Fortunately, this condition is satisfied for many of polymeric solutions[8, 11], and thus the overdamped Langevin equations can be safely utilized to describe the dynamics. However, this does not mean that the overdamped Langevin equations are always justified. For example, if we model the relaxation process of relatively small objects such as segments, or if we consider the higher order Rouse modes (which essentially behave as statistically independent multiple dumbbell models), the inertia effect may not be simply neglected.

In this work, we investigate the effect of the inertia on the linear viscoelasticity of a simple dumbbell type model. We first derive the approximate Langevin equation for the dumbbell model which is valid for the weak (but non-zero) inertia case. Starting from the underdamped Langevin equation with the inertia term, we show that an overdamped Langevin equation with the lowest order correction by the inertia can be obtained. The dynamics is slightly accelerated by the inertia. Then we apply the obtained approximate dynamic equation to the harmonic dumbbell for which we can obtain analytic solutions. We calculate the shear relaxation modulus and the longest relaxation time, with and without approximations, and discuss how the inertia affects the linear viscoelasticity. We also discuss how the inertia effect accelerates the Rouse model, and a similar acceleration effect by the memory effect in the generalized Langevin equation.

2 Theory

2.1 Model

We start from a dumbbell model[5, 6] with an arbitrary bond potential and the inertial terms. We assume that a dumbbell consists of two beads with the same mass, $M$. We describe the positions and momentum of two beads as $\bm{R}_{1}$ and $\bm{R}_{2}$, and $\bm{P}_{1}$ and $\bm{P}_{2}$, respectively. The Hamiltonian is

[TABLE]

where $U(\bm{R}_{1}-\bm{R}_{2})$ is the bond potential (which is a harmonic function in the case of the harmonic dumbbell model). If we employ the Langevin thermostat to control the temperature of the system, the dynamic equations for the dumbbell is given as:

[TABLE]

for $i=1$ and $2$. Here, $Z$ is the friction coefficient for a bead, $k_{B}$ is the Boltzmann constant, $T$ is the temperature, and $\bm{W}_{i}(t)$ is a Gaussian white noise. The first and second moments of the noise is given as follows:

[TABLE]

where $\langle\dots\rangle$ represents the statistical average and $\bm{1}$ is the unit tensor.

We may rewrite the dynamic equation by introducing the center of mass $\bar{\bm{R}}\equiv(\bm{R}_{1}+\bm{R}_{2})/2$ and the bond vector $\bm{r}\equiv\bm{R}_{2}-\bm{R}_{1}$. We describe the momenta for the center of mass and the bond vector as $\bar{\bm{P}}\equiv\bm{P}_{1}+\bm{P}_{2}$ and $\bm{p}\equiv(\bm{P}_{2}-\bm{P}_{1})/2$. The Hamiltonian (1) can be rewritten as

[TABLE]

with the total mass $\bar{M}\equiv 2M$ and the reduced mass $m\equiv M/2$. Then we have statistically independent two set of dynamic equations:

[TABLE]

and

[TABLE]

where $\bar{Z}\equiv 2Z$, $\zeta=Z/2$, $\bar{\bm{W}}(t)\equiv(\bm{W}_{1}(t)+\bm{W}_{2}(t))/\sqrt{2}$ and $\bm{w}(t)\equiv(\bm{W}_{2}(t)-\bm{W}_{1}(t))/\sqrt{2}$. Here, it should be noticed that $\bar{\bm{W}}(t)$ and $\bm{w}(t)$ are statistically independent Gaussian white noises, and both of them have the zero mean and unit variance (in the same way as eq (4)).

At the overdamped limit, we set $d\bar{\bm{P}}/dt=0$ and $d\bm{p}/dt$ = 0 in eqs (7) and (9). Then eqs (6)-(9) reduce to the well-known dynamic equations for the fully overdamped dumbbell model:

[TABLE]

This approximation is usually justified due to the fact that the characteristic relaxation time of the inertia, $\tau_{m}\equiv m/\zeta=\bar{M}/\bar{Z}$, is much faster than that of the bond (by the competition between the restoring force and the friction force). The momentum relaxation time is common for the center of mass and the bond. Various properties including the linear viscoelasticity of eqs (10) and (11) have been studied, and can be found in literates[5, 6]. The center of mass motion does not contribute to the linear viscoelasticity, and the linear viscoelasticity is fully determined by eq (11). The properties of the overdamped dumbbell model with various bond potential models such as the finitely extensible nonlinear elasticity (FENE) model have been studied extensively[5, 6].

2.2 Approximation for Weak Inertia Case

To investigate the effect of the inertia, we consider an approximate dynamic equation which is valid for small $\tau_{m}=m/\zeta$. We rewrite eqs (8) and (9) as the second order stochastic differential equation for $\bm{r}$:

[TABLE]

We want a first order stochastic differential equation as an approximate form for eq (12). For this purpose, we further rewrite eq (12) as

[TABLE]

where $(1+\tau_{m}d/dt)^{-1}$ is the inverse operator of $(1+\tau_{m}d/dt)$, and is defined via the following relation for a given function $f(t)$:

[TABLE]

For small $\tau_{m}$, the operator can be expanded into the power series of $\tau_{m}$, and thus we have a simple approximate form, $(1+\tau_{m}d/dt)^{-1}\approx 1-\tau_{m}d/dt$. (We retain only the first order term in $\tau_{m}$.) Then we have the approximate form for the dynamic equation (13) as

[TABLE]

where $\bm{K}(\bm{r})$ is the potential curvature tensor defined as

[TABLE]

and we have utilized the Ito formula[12] to calculate $(d/dt)\partial U(\bm{r}(t))/\partial\bm{r}(t)$. In eq (15), we have terms which contain $d\bm{r}(t)/dt$ in both sides. Rearranging these terms and we have

[TABLE]

with the effective (dimensionless) mobility tensor defined as

[TABLE]

The third term in the right-hand side of eq (17) can be interpreted as the random noise. This noise is not well-defined (due to the fact that time derivative of the noise $\bm{w}(t)$ is not well-defined) and the fluctuation-dissipation relation seems not to be satisfied. It would be natural for us to require the thermodynamic equilibrium state can be correctly realized by the approximate dynamic equation. This requirement can be achieved by modifying the noise term in eq (17) to recover the fluctuation-dissipation relation. (Here, it should be noticed that such an ad-hoc modification cannot be fully justified in general. For the case of the harmonic dumbbell, however, the noise property does not directly affect the linear viscoelasticity. Thus this approximation is not bad for the analyses in what follows.) Finally, we have the following overdamped Langevin equation as an approximate dynamic equation:

[TABLE]

with $\bm{B}(\bm{r})$ being the noise coefficient matrix which satisfies $\bm{B}(\bm{r})\cdot\bm{B}^{\mathrm{T}}(\bm{r})=\bm{\Lambda}(\bm{r})$. Eq (19) is the approximate dynamic equation where the inertia effect is weak, but not absent. Even at the long-time scale, eq (19) does not reduce to eq (11). Only at the fully overdamped case where $\tau_{m}=0$ (or $m=0$), the mobility tensor simply becomes the unit tensor ($\bm{\Lambda}(\bm{r})=\bm{1}$) and eq (19) reduces to eq (11). The fact that the overdamped limit corresponds to the inertialess case ($m=0$) is consistent with the result obtained by Schieber and Öttinger[9].

2.3 Harmonic Dumbbell Model

Here we consider the linear viscoelasticity of the harmonic dumbbell model with the inertia effect and study how the inertia affects the linear viscoelasticity. We set $U(\bm{r})=k\bm{r}^{2}/2$ and assume that the spring constant $k$ satisfies $4km<\zeta^{2}$. (This condition is required to suppress the oscillatory behavior, which is out of the scope of the current work.) At the fully overdamped case where $m=0$, the shear relaxation modulus simply reduces to the single-mode Maxwell model:

[TABLE]

where $\nu$ is the number density of the dumbbell and $\tau_{b}\equiv\zeta/k$ is the bond relaxation time. In this subsection, we consider how the shear relaxation modulus deviates from eq (20) in the presence of the inertia effect.

We start from the exact solution for the case of the underdamped dumbbell model. The stress of the single harmonic dumbbell can be expressed as the Kramers form:

[TABLE]

Then, the shear relaxation modulus is given by the auto-correlation function of the stress tensor, by the Green-Kubo formula[13]:

[TABLE]

Here $\hat{\sigma}_{xy}(t)$ represents the $xy$ component of the stress tensor (21) at time $t$. From eqs (21) and (22), the shear relaxation modulus for the harmonic dumbbell can be expressed in terms of correlation functions. Because $\bar{\bm{P}}$ is statistically independent of $\bm{r}$ and $\bm{p}$, and the system is isotropic, the expression for the relaxation modulus can be simplified as

[TABLE]

Here we have utilized some relations for correlation functions such as $\langle r_{y}(t)r_{y}(0)\rangle=\langle r_{x}(t)r_{x}(0)\rangle$ and $\langle r_{x}(t)r_{y}(0)\rangle=0$.

The correlation function for $\bar{\bm{P}}$ can be easily obtained. From eq (7), we have

[TABLE]

with $\tau_{m}=\bar{M}/\bar{Z}$ being the momentum relaxation time. Also, the correlation functions for $\bm{r}$ and $\bm{p}$ can be calculated from eqs (8) and (9). For the harmonic potential, the restoring force $-k\bm{r}$ is linear in $\bm{r}$, and we can rewrite eqs (8) and (9) as

[TABLE]

Since eq (25) is linear in $[\bm{r},\bm{p}]^{\mathrm{T}}$, we can solve it by diagonalizing the coefficient matrix. To calculate the linear viscoelasticity, we need the two time correlation functions. The result is:

[TABLE]

with the eigenvalues of the coefficient matrix defined as

[TABLE]

See Appendix A for the detailed calculations. Eq (29) can be rewritten in terms of the momentum relaxation time $\tau_{m}=m/\zeta$ and the bond relaxation time $\tau_{b}=\zeta/k$:

[TABLE]

and if two relaxation times are well separated, $\tau_{b}\gg\tau_{m}$, we can approximate eigenvalues as $\lambda_{+}\approx 1/\tau_{m}$ and $\lambda_{-}\approx(1+\tau_{m}/\tau_{b})/\tau_{b}$.

By substituting eqs (24) and (26)-(28) into eq (23), we finally have the explicit expression for the relaxation modulus:

[TABLE]

The detailed calculations are shown in Appendix A. From eq (31), the relaxation modulus consists of three Maxwell models. At the long-time region, only the longest relaxation mode in eq (31) survives. In addition, if the momentum relaxation is sufficiently fast, $\tau_{b}\gg\tau_{m}$, the relaxation time can be approximated as a simple form. Thus we have the following approximate form:

[TABLE]

Eq (32) is similar to eq (20) but the relaxation time is slightly different. The relaxation time of eq (32) is $\bar{\tau}\equiv\tau_{b}/2(1+\tau_{m}/\tau_{b})$. This relaxation time $\bar{\tau}$ coincides to $\tau_{b}/2$ if $\tau_{m}/\tau_{b}=0$, but it is generally different from $\tau_{b}/2$. Therefore, we conclude that the overdamped limit cannot be justified unless $\tau_{m}/\tau_{b}$ is sufficiently small.

Now we apply the approximate dynamic equation (19) to the harmonic dumbbell case. By substituting $U(\bm{r})=k\bm{r}^{2}/2$ to eqs (16) and (18), we have

[TABLE]

From eq (33), we find that for the harmonic dumbbell model, the mobility tensor is independent of the bond vector $\bm{r}$. Therefore, eq (19) becomes the following overdamped Langevin equation with the additive noise:

[TABLE]

Eq (34) is a simple linear Langevin equation (or the Ornstein-Uhlenbeck process[2]), and can be solved easily. From eq (34), we have the correlation function as

[TABLE]

The single dumbbell stress tensor is given as the partial average over the momenta for eq (21):

[TABLE]

and the Green-Kubo formula (22) gives the following expression for the relaxation modulus:

[TABLE]

Finally, from eqs (35) and (37), we have the following explicit expression for the relaxation modulus:

[TABLE]

Eq (38) has the same form as eq (32). Therefore, from the approximate dynamic equation (19), we can successfully reproduce the relaxation time $\bar{\tau}=\tau_{b}/2(1+\tau_{m}/\tau_{b})$.

3 Discussions

3.1 Relaxation of Harmonic Dumbbell Model

The long-time dynamics of a Brownian particle with inertia can be well approximated by that without inertia. The diffusion coefficients calculated from the mean square displacement the long-time region, with and without inertia, are exactly the same. Thus for a diffusion process, we can reasonably approximate the underdamped Langevin equation by the overdamped one at the long-time region where $t\gtrsim\tau_{m}$. In contrast, the long-time dynamics of a dumbbell model with inertia does not reduce to that without inertia. The exact expression for the relaxation modulus with inertia (eq (31)) has three relaxation modes: The center of mass momentum relaxation model, the bond momentum relaxation mode, and the bond relaxation mode. The degrees of freedom of the center of mass ($\bar{\bm{R}}$ and $\bar{\bm{P}}$) give only a single mode, which is essentially the same as the relaxation mode of a single Brownian particle. On the other hand, the degrees of freedom of the bond ($\bm{r}$ and $\bm{p}$) have two relaxation modes, unlike the case of the center of mass. The relaxation times of these two modes depend both on the momentum relaxation time $\tau_{m}$ and the bond relaxation $\tau_{b}$. We may interpret that these two modes are kinetically coupled. Thus even at the long-time region, the relaxation time $\bar{\tau}$ depends on both $\tau_{b}$ and $\tau_{m}$.

We show the relaxation modulus of the harmonic dumbbell with the inertia effect ($\tau_{m}/\tau_{b}=0.1$) in Figure 1. At the short-time scale ($t\lesssim\tau_{m}$), we observe the momentum relaxation modes. The modulus at the overdamped limit (eq (20)) deviates from the exact one (eq (31)) even at the long-time region ($t\gtrsim\tau_{b}$). (The deviation seems not large at the logarithmic scale for $\tau_{m}/\tau_{b}=0.1$.) In contrast to the overdamped limit, the approximate form for weak inertia (eq (38)) agrees reasonably with the exact one. Thus we find that the approximate dynamic equation (19) derived in this work is a reasonable approximation where the inertia effect is small but not fully negligible.

To study the effect of the inertia, we consider how the zero shear viscosity $\eta_{0}$ and the longest relaxation time $\tau_{d}$ depend on the ratio $\tau_{m}/\tau_{b}$. From eq (31), the longest relaxation time is $\tau_{d}=\tau_{m}/(1-\sqrt{1-4\tau_{m}/\tau_{b}})$ (for $\tau_{m}/\tau_{b}<1/4$). The zero shear viscosity at the overdamped limit is simply given as $\nu k_{B}T\tau_{b}/2$. On the other hand, from eq (31), the zero shear viscosity with the inertia effet becomes

[TABLE]

Figure 2 shows the $\tau_{m}/\tau_{b}$ dependence of the zero shear viscosity $\eta_{0}$ and the longset relaxation time $\tau_{d}$. (For convenience, $\eta_{0}$ and $\tau_{d}$ are normalized by their overdamped limits.) As $\tau_{m}/\tau_{b}$ increases, the zero shear viscosity increases wheareas the longest relaxation time decreases. The increase of the viscosity is due to two relatively fast relaxation modes (the momentum relxation modes) which disappear at the overdamped limit.

We can understand the effect of the inertia term on the dynamics via the Langevin equation (19). In eq (19), inertia term modulates the the mobility tensor from $\bm{1}$ to $\bm{\Lambda}(\bm{r})$. Intuitively, this can be interpreted as the acceleration effect. The acceleration depends on the curvature tensor $\bm{K}(\bm{r})$ which generally depends on the bond vector $\bm{r}$. In the case of the harmonic potential, the curvature is isotropic and given as a positive constant $k$, and therefore all the dynamics is simply accelerated by the factor $(1+\tau_{m}/\tau_{b})$. This acceleration factor is consistent with the fact that the ratio of the relaxation time $\bar{\tau}$ to one at the overdamped limit $\tau_{b}/2$ is $\bar{\tau}/(\tau_{b}/2)=(1+\tau_{m}/\tau_{b})^{-1}$. If the potential is not harmonic, the acceleration depends on the bond vector. In such a case, the relaxation will be much more complex compared with the harmonic dumbbell model. However, any bond potential potentials will exhibit the acceleration effect on the bond dynamics, on average. Thus we expect that the acceleration of the viscoelastic relaxation time will be common and independent of the details of the bond potential.

3.2 Rouse Model with Inertia Effect

The acceleration effect by the inertia term will be non-negligible if we study relatively small scale relaxation modes. A simple yet nontrivial example is the Rouse model[1]. Usually, the Rouse model is defined as the the overdamped Langevin equation with the free energy for an ideal chain which consists of $N$ beads. For simplicity, we assume that $N$ is sufficiently large. We employ the underdamped Langevin equation instead of the overdamped one:

[TABLE]

where $\bm{R}_{k}$ represents the position of the $k$-th bead ($k=1,2,\dots,N$), $m$ is the mass of a bead, $\zeta$ is the friction coefficient for a bead, $\mathcal{F}(\{\bm{R}_{k}\})$ is the free energy for an ideal chain, and $\bm{w}_{k}(t)$ is the Gaussian white noise. The free energy is given as

[TABLE]

with $b$ being the segment size, and the noise satisfies the following relations:

[TABLE]

By using the Rouse mode defined as $\bm{X}_{p}\equiv\sqrt{2/N}\sum_{k}\cos(p\pi k/N)\bm{R}_{k}$ ($p=1,2,3,\dots$), we can rewrite the dynamic equation (40) as

[TABLE]

where $\tau_{R}\equiv N^{2}b^{2}\zeta/3\pi^{2}k_{B}T$ is the Rouse time and $p=1,2,3,\dots$ is the mode index (we ignore the zeroth mode $p=0$, because it does not contribute to the stress for the weak inertia case). $\bm{W}_{p}(t)$ is the Gaussian white noise which satisfies

[TABLE]

We apply eq (19) to the Rouse model, eq (43). Then, the approximate dynamic equation becomes

[TABLE]

where $\tau_{m}=m/\zeta$ is the momentum relaxation time. Eq (45) is valid only when the inertia effect is weak. Roughly, this condition is estimated as $p\lesssim p^{*}$, where $p^{*}$ is given as $p^{*}=\sqrt{\tau_{R}/\tau_{m}}=Nb\zeta/(\pi\sqrt{3mk_{B}T})$.

The discussions above implies that we should introduce the cutoff for the index $p$ unless the inertia effect is not fully negligible. In molecular dynamics simulations such as the Kremer-Grest model[14], we usually employ the unit where the mass, energy, and the bead size becomes unity. Other quantities such as the friction coefficient $\zeta$ and the temperature $k_{B}T$ are typically of the order of unity. Then the maximum index $p^{*}$ can be estimated as $p^{*}\approx N/(\pi\sqrt{3})\approx 0.18N$. Thus, we find that only less than $20\%$ of the Rouse modes work as naively expected. (Higher order modes will largely deviate from the prediction of the simple Rouse model.)

Even if the cutoff effect is minor and thus not important, the acceleration by the inertia effect depends on the index $p$. The shear relaxation modulus can deviate from the simple Rouse type behavior. From eq (45), the relaxation modulus becomes

[TABLE]

where $\nu$ is the chain density. Therefore, the longest relaxation time is given as $\bar{\tau}=\tau_{R}/2(1+\tau_{m}/\tau_{R})$. As the case of the harmonic dumbbell, the relaxation is accelerated by the factor of $1+\tau_{m}/\tau_{R}$. For the intermediate time region where $\tau_{m}\lesssim t\lesssim\tau_{R}$, we may replace the discrete sum over $p$ by a continuum integral (as the usual approximation for the Rouse model[1]), we have

[TABLE]

where $K_{1/4}(x)$ is the modified Bessel function of the second kind with the order $1/4$[15]. (The summation over $p$ for $p\geq 1$ in eq (46) is replaced by the integral for $p\geq 0$, not for $p\geq 1$, in eq (47). The integral over the range $0\leq p\leq 1$ is much smaller than that over the range $p\geq 1$, and thus we simply neglected it.) We can analytically calculate the complex modulus of the model based on eq (46) or eq (47). The complex moduli given by eqs (46) and (47) show reasonable agreement, and thus we find that eq (47) works as a good approximation. The detailed calculations are shown in Appendix B. At relatively long-time scale ($\tau_{m}\ll t\lesssim\tau_{R}$), the asymptotic expansion[15] can be utilized and eq (47) can be simplified as

[TABLE]

The first term in the parenthesis in the right-hand side of eq (48) corresponds to the usual (inertialess) Rouse type terms[1], and the second term can be interpreted as the correction. Eqs (47) and (48) imply that the shear relaxation modulus could slightly deviate from the well-known Rouse type behavior $G(t)\propto t^{-1/2}$, if the inertia effect is not sufficiently weak. We show the relaxation modulus by eq (47) (with $\tau_{m}/\tau_{R}=0.01$) together with the usual expression at the fully overdamped limit in Figure 3. For comparison, we show the relaxation modulus data directly calculated with the discrete sum expression, eq (46), with $\tau_{m}/\tau_{R}=0.01$ and [math]. With this parameter setting, the longest relaxation time $\bar{\tau}$ is almost the same as that at the overdamped limit $\tau_{R}/2$ ($\bar{\tau}/(\tau_{R}/2)\approx 0.99$ for $\tau_{m}/\tau_{R}=0.01$). The relaxation modulus by eq (47) looks almost the same as that at the overdamped limit at the relatively long-time region in Figure 3. However, we observe a clear deviation between eq (47) and the overdamped limit at the relatively short-time region. We observe a similar discrepancy for moduli calculated with the discrete sum expression. As the case of the dumbbell model, the inertia effect accelerates the relaxation.

Apparently, these results seem not to be consistent with the fact that the experimental and simulation data of shear relaxation modulus for unentangled polymer melts can be fitted well to the Rouse model[16, 17]. However, we should recall that the polymer chains in melts are strongly interacting with each other by the short-range repulsive interaction between segments, even if they are not entangled. From this point of view, it is not plausible for unentangled polymer melts to obey the simple Rouse model. Also, a recent study on the Langevin equation with a fluctuating diffusivity suggested that the Rouse type relaxation can emerge if the diffusivity (or the friction coefficient) is heterogeneous and fluctuates in time[18]. Therefore, even if the shear relaxation modulus exhibits the Rouse type behavior, we cannot conclude that the molecular level dynamics fully obeys the simple Rouse model. Careful and detailed analyses would be required to further investigate the molecular level dynamics of polymers.

3.3 Memory Effect

To model complex relaxation behavior, the generalized Langevin equation[3, 19, 20] is widely employed. In the generalized Langevin equation, the memory effect is expressed by using a memory kernel. Formally, the generalized Langevin equation is obtained by eliminating the microscopic fast degrees of freedom, and the projection operator formalism[19] gives the expression of the memory kernel. Although, in principle, the memory kernel for a specific system can be determined by the dynamics of the fast degrees of freedom in the system, we cannot calculate its explicit form in most of practical cases. Thus rather phenomoenolgical and simple memory functions such as a single exponential form are employed in some cases. In addition, if the relaxation time of fast degrees of freedom, $\tau_{f}$, is much faster than the characteristic time scale of the slow degrees of freedom, the memory effect may be neglected. (This relaxation time $\tau_{f}$ can be interpreted as the memory relaxation time.) Then the usual Langevin equation can be used as an approximate dynamic equation (the Markovian approximation).

In some aspects, the Markovian approximation explained above is similar to the overdamped limit for the Langevin equation with the inertia effect. Thus it would be informative to discuss how the memory kernel affects the relaxation time of the dumbbell model. Actually, the Langevin equation (12) can be rewritten as the generalized Langevin equation, and thus we can interpret our approximation method as an approximation for the generalized Langevin equation. We introduce the Green function $\Gamma(t)$ for the operator $(1+\tau_{m}d/dt)$:

[TABLE]

It is straightforward to see the following function $\Gamma(t)$ satisfies eq (49) and thus is the Green function:

[TABLE]

where $\Theta(t)$ is the Heaviside step function. The application of the operator $(1+\tau_{m}d/dt)^{-1}$ to a given function $f(t)$ can be interpreted as the convolution of $\Gamma(t)$ and $f(t)$:

[TABLE]

Then eq (12) (or eq (13)) can be exactly rewritten as

[TABLE]

where $\bm{\xi}(t)$ is a Gaussian colored noise which is defined as

[TABLE]

Clearly, the first moment of $\bm{\xi}(t)$ is zero, $\langle\bm{\xi}(t)\rangle=0$. The second moment reduces to a simple form:

[TABLE]

Eq (54) is nothing but the fluctuation-dissipation relation, and thus we can interpret eq (52) as the generalized Langevin equation with the memory kernel $\Gamma(t)$ given by eq (50).

As we discussed, the inertia effect accelerates the relaxation. This acceleration effect is also caused in the case of the generalized Langevin equation. Therefore, in general, the memory effect may not be fully neglected even if we consider relatively long-time scale dynamics. To reproduce the correct relaxation behavior, we should employ the effective mobility tensor. We expect that the analyses for the weak inertia cases can be used to analyze the dumbbell model with a weak memory effect, by replacing the momentum relaxation time $\tau_{m}$ by the memory relaxation time $\tau_{f}$.

4 Conclusions

In this work, we studied the effect of the inertia on the dynamics of the dumbbell model. We derived the approximate dynamic equation to describe the dynamics for the weak inertia case. The inertia effects modulate the mobility tensor, and thus the dynamics is slightly accelerated. As a result, the dynamical properties such as the shear relaxation modulus deviates from those at the overdamped limit.

As an example, we analyzed the harmonic dumbbell model which can be solved analytically. We calculated the exact expression for the shear relaxation modulus and showed that there are three relaxation modes. One mode is the momentum relaxation of the center of mass, and the other two modes are the momentum relaxation of the bond and the bond relaxation. The relaxation times of two modes for the bond depend both on the momentum relaxation time $\tau_{m}$ and the bond relaxation time $\tau_{b}$. The longest relaxation time is approximated as $\bar{\tau}\approx\tau_{b}/2(1+\tau_{m}/\tau_{b})$, which is different from that at the overdamped limit, $\tau_{b}/2$. On the other hand, the approximate dynamic equation (19) gives the correct long relaxation time. Although the inertia effect may not be important in many practical cases, it can be non-negligible for some limited cases. We consider that the analyses based on the approximate dynamic equation and the detailed and careful analyses of the experimental data will be required when the inertia effect becomes non-negligible.

Acknowledgment

The authors thank an anonymous reviewer to point that the sum in eq (46) can be evaluated analytically by utilizing the Laplace-Fourier transform. This work was supported by Grant-in-Aid (KAKENHI) for Scientific Research C 16K05513 and Grant-in-Aid (KAKENHI) for Scientific Research A 17H01152.

Appendix

Appendix A Detailed Calculations for Harmonic Dumbbell Model

In this appendix, we show detailed calculations for the harmonic dumbbell model with inertia effect, without approximations. The dynamic equation is given as eq (25) in the main text. We describe the coefficient matrix as

[TABLE]

$\bm{C}$ has the eigenvalues $\lambda_{\pm}$ in eq (29) and the corresponding eigenvectors are

[TABLE]

Then we can diagonalize eq (25) with the transformation matrix $\bm{V}\equiv[\bm{\psi}_{+}\,\bm{\psi}_{-}]$ and its inverse $\bm{V}^{-1}$. Multiplying $\bm{V}^{-1}$ from the left side to eq (25) gives

[TABLE]

After integrating eq (57) from [math] to $t$, we multiply $\bm{V}$ to both sides and have

[TABLE]

Now, $\bm{r}(0)$ and $\bm{p}(0)$ are statistically independent of the noise $\bm{w}(t)$ for $t>0$. In addition, in equilibrium, $\bm{r}(0)$ and $\bm{p}(0)$ obey the Boltzmann distribution with the Hamiltonian (5). Thus the averages over $\bm{r}(0),\bm{p}(0)$, and $\bm{w}(t)$ can be took separately. The two time correlation functions for $\bm{r}$ and $\bm{p}$ can be obtained as

[TABLE]

where we have utilized the relation for the eigenvalues, $\lambda_{+}\lambda_{-}=k/m$. Eq (59) gives eqs (26)-(28) in the main text.

The relaxation modulus $G(t)$ can be calculated by eqs (23), (24) and (26)-(28):

[TABLE]

where we have utilized the relation $k/m=(k/\zeta)(\zeta/m)=1/\tau_{b}\tau_{m}$. Some factors in eq (60) can be simplified as follows:

[TABLE]

[TABLE]

[TABLE]

Thus we have

[TABLE]

Eq (64) together with eq (29) gives eq (31) in the main text.

Appendix B Complex Modulus of Rouse Model

In this appendix, we calculate the complex modulus $G^{*}(\omega)$ of the Rouse model with inertia effect. The complex modulus is obtained as the Fourier transform of the relaxation modulus $G(t)$. From the view point of experiment, the complex modulus is much easier to measure and thus it would be informative to show some explicit expressions for the complex modulus.

Firstly we start from eq (46). The complex modulus $G^{*}(\omega)$ is calculated as

[TABLE]

where we have defined $\alpha_{\pm}$ as

[TABLE]

The sums over $p$ in eq (65) can be rewritten by utilizing the partial fraction expansion for $\coth z$[15]:

[TABLE]

By combining eqs (65) and (67), we have

[TABLE]

Eq (68) cannot be reduced to simpler form in general. At the intermediate angular frequency range where $\tau_{R}^{-1}\ll\omega\ll\tau_{m}^{-1}$, we can approximate $\alpha_{\pm}$ as

[TABLE]

$\cosh(z)$ can be approximated as $\cosh(z)\approx 1$ for $\mathrm{Re}(z)\gg 1$. Then the complex modulus can be approximated as follows:

[TABLE]

The first term in the parenthesis in the last line of eq (70) corresponds to the complex modulus of the fully overdamped Rouse model. Eq (70) means that the complex modulus is decreased (or the relaxation is accelerated) by the inertia effect, and this behavior is consistent with one for the relaxation modulus.

Although eq (68) is justified only for $\omega\tau_{m}\lesssim 1$, it would be informative to see how it behave at the high frequency region ($\omega\tau_{m}\gg 1$). $\alpha_{\pm}$ can be approximated as

[TABLE]

We can approximate eq (68) as

[TABLE]

Thus, for large $\omega$, the complex modulus depends on the angular frequency as $G^{*}(\omega)\propto\omega^{1/4}$. (However, this behavior cannot be observed, because we have momentum relaxation modes in such a region.)

Secondly, we use the approximate expression for the relaxation modulus, eq (47). The approximate expression for the complex modulus becomes

[TABLE]

where we have defined $u\equiv t/4\tau_{m}$ and $\beta\equiv 1-0-4i\omega\tau_{m}$ (we introduced $-0$ to satisfy the condition $\mathrm{Re}(\beta)<1$). The integral over $u$ in eq (65) can be simplyfied by utilizing the following formula[21] for $\mathrm{Re}(s)>-1$:

[TABLE]

From eqs (73) and (74), we have

[TABLE]

Eqs (68) and (75) look rather different, at least apparently. As before, we consider the case where $\tau_{R}^{-1}\ll\omega\ll\tau_{m}^{-1}$. For $\omega\tau_{m}\ll 1$, we have

[TABLE]

and we can approximate eq (75) as

[TABLE]

Eq (79) has the same form as eq (70), except an $\omega$-independent constant term. For the case of $\omega\tau_{m}\gg 1$, we can approximate

[TABLE]

[TABLE]

and thus

[TABLE]

Eq (82) has the same angular frequency dependence as eq (72) (although the absolute value is slightly different by the numerical factor $2^{-1/4}\approx 0.8409$.). Therefore, we conclude that the approximation employed to calculate eq (47) is reasonable.

Figure Captions

Figure 1: The shear relaxation modulus of a harmonic dumbbell with the inertia effect. The momentum relaxation time is $\tau_{m}/\tau_{b}=0.1$. The dashed black curve is the exact result by eq (31). The solid black curve and the dotted gray curve are the approximations for weak inertia (eq (38)) and the fully overdamped limit (eq (20)), respectively.

Figure 2: The zero shear viscosity $\eta_{0}$ and the longest relaxation time $\tau_{d}$ of a harmonic dumbbell with the inertia effect. Both $\eta_{0}$ and $\tau_{d}$ are normalized by their overdampled limits ($\nu k_{B}T\tau_{b}/2$ and $\tau_{b}/2$, respectively).

Figure 3: The shear relaxation modulus of the Rouse dumbbell with the inertia effect, at the intermediate time region $\tau_{m}\lesssim t\lesssim\tau_{R}$. The momentum relaxation time is $\tau_{m}/\tau_{R}=0.01$. The solid black curve is the approximation for weak inertia (eq (47)), and the dotted gray curve is modulus at the fully overdamped limit $G(t)=\nu k_{B}T\sqrt{\pi\tau_{R}/2t}$. Black circles are results of the direct numerical calculations of the discrete sum without the integral approximation (eq (46)). Gray diamonds are the results for the discrete sum at the fully overdamped limit (eq (46) with $\tau_{m}/\tau_{R}=0$).

Figures