4.3: Viscoelastic Model

Introduction

The Generalized Langevin Equation and Mode-Coupling theory are subsets of molecular hydrodynamics, the theory that was developed to bridge the gap between hydrodynamics and molecular dynamics. Hydrodynamics, which we discussed in chapter 3, describes the macroscopic, long time behavior of systems in the limit as \(t \rightarrow \infty\) and \(k \rightarrow 0\). It uses the transport coefficients \(D, \lambda\), \(\eta\), and \(\eta_{B}\) to predict long time fluctuations. Molecular dynamics, which we discussed in section I of chapter 4 , describes the microscopic, short time behavior of systems in the limit as \(t \rightarrow 0\) and \(k \rightarrow \infty\). In this limit, systems behave as static liquid structures, and their dynamics are largely determined by the pairwise interaction potential.

In this section, we will use the GLE to derive the viscoelastic model for transverse current. By taking the appropriate limits, we can show that the results of the viscoelastic model are consistent with those of hydrodynamics and molecular dynamics, and that this model provides a successful bridge between the two limits.

Phenomenological Viscosity

Consider a constant shear force applied to a viscous liquid. At long times, the shear stress \(\sigma_{x z}\) in the liquid is related to the rate of strain \(\partial_{z} \vec{v}_{x}+\partial_{x} \vec{v}_{z}\) by

\[\sigma_{x z}=-\eta\left(\partial_{z} \vec{v}_{x}+\partial_{x} \vec{v}_{z}\right)\]

Liquids behaving in this fashion do not support shear waves. However, if the force is applied instantaneously, the system does not have the time to relax like a liquid. Instead, it behaves like an elastic solid. The stress is now proportional to the strain rather than the rate of strain. The short term response is

\[\sigma_{x z}=-G\left(\partial_{z} \vec{x}+\partial_{x} \vec{z}\right)\]

where \(G\) is the modulus of rigidity. When the liquid is behaving like a solid, it supports shear waves propagating at a speed of \(v_{s}=\sqrt{\frac{G}{\rho m}}\).

To determine the time scale on which the liquid behaves like an elastic solid, define the constant

\[\tau_{M}=\frac{\eta}{G}\]

This is the Maxwell relaxation time. For the timescales \(t\) when

\[\frac{t}{\tau_{M}} \ll 1\]

the system behaves like an elastic solid. For the timescales when

\[\frac{t}{\tau_{M}} \gg 1\]

the system behaves like a viscous liquid.

Viscoelastic Approximation

To interpolate between the two extremes, we can write

\[\left(\frac{1}{\eta}+\frac{1}{G} \frac{\partial}{\partial t}\right) \sigma_{x z}=-\left(\frac{\partial}{\partial z} \vec{v}_{x}+\frac{\partial}{\partial_{x}} \vec{v}_{z}\right)\]

The Laplace transform of this equation yields

\[\hat{\eta}(s)=\frac{G}{s+\frac{G}{\eta}}=\frac{G}{s+\frac{1}{\tau_{M}}}\]

In the steady-state limit, as \(s \rightarrow 0\)

\[\lim _{s \rightarrow 0} \hat{\eta}=\eta\]

and in the high-frequency limit, as \(s \rightarrow \infty\)

\[\lim _{s \rightarrow \infty} \hat{\eta}=\frac{G}{s}\]

Transverse Current Correlation Function

We will use the transverse current correlation function to demonstrate the viscoelastic approximation. In Section I, we defined the transverse current as

\[J_{T}(\vec{k}, t)=\sum_{j=1}^{N} \overrightarrow{\dot{x}}_{j} e^{i \vec{k} \vec{z}_{j}}\]

and the transverse current correlation function as

\[C_{T}(\vec{k}, t)=\frac{1}{N}\left\langle J_{T}(\vec{k}, t) \mid J_{T}(\vec{k}, 0)\right\rangle\]

We have studied the transverse current in both the hydrodynamic limit \((k \rightarrow 0, \omega \rightarrow 0)\) and the short-time expansion limit \((k \rightarrow \infty, \omega \rightarrow \infty)\). In chapter 3, we used the Navier Stokes equation to find an equation of motion for the transverse correlation function in the hydrodynamic limit

\[\dot{C}_{T}(\vec{k}, t)=-\nu k^{2} C_{T}(\vec{k}, t)\]

This has the solution

\[C_{T}(\vec{k}, t)=C_{T}(\vec{k}, 0) e^{-\nu k^{2} t}=v_{o}^{2} e^{-\nu k^{2} t}\]

where \(\nu=\frac{\eta}{m \rho}\) is the shear viscosity. Therefore, in the hydrodynamic limit, transverse current fluctuations decay exponentially with a rate determined by the shear viscosity \(\nu\).

In section I of this chapter, we used the short-time expansion approximation to show that in the \(k \rightarrow \infty, \omega \rightarrow \infty\) limit, the transverse current correlation function can be written as

\[C_{T}(\vec{k}, t)=v_{o}^{2}\left(1-\omega_{T}^{2} \frac{t^{2}}{2}\right)+\cdots\]

where the transverse frequency \(\omega_{T}\) is related to the transverse speed of sound \(c_{T}(k)\) by

\[\omega_{T}^{2}=k^{2} c_{T}^{2}(k)\]

And the transverse speed of sound is given by

\[c_{T}^{2}(k)=v_{o}^{2}+\frac{\rho}{m} \int g(\vec{r}) \partial_{x}^{2} U(\vec{r}) \frac{\left[1-e^{i k z}\right]}{k^{2}} d \vec{r}\]

where \(g(\vec{r})\) is the pairwise correlation function and \(U(\vec{r})\) is the pairwise interaction potential. This frequency term can also be written as [3]

\[\omega_{T}^{2}=\frac{\left(k v_{o}\right)^{2}}{n M} G_{\infty}(k)\]

where \(G_{\infty}(k)\) is the shear modulus. This indicates that at short times and wavelengths, the dissipation effects are diminished and transverse current fluctuations can propagate through the material with speed \(c_{T}(k)\).

Using the Generalized Langevin equation, we can generate a model for transverse current fluctuations that replicates the results of hydrodynamics and the short-time expansion when the appropriate limits are taken. Begin by writing the GLE for transverse current. Since the frequency matrix is zero, the GLE is written

\[\dot{J}_{T}(k, t)=-\int_{0}^{t} \kappa(t-\tau) J_{T}(k, \tau) d \tau+f(t)\]

where \(\kappa(t)\) is the memory function and \(f(t)\) is the noise term. Multiplying through by \(J_{T}(k, 0)\) and taking the average gives us the equation of motion for the transverse current correlation function

\[\dot{C}_{T}(k, t)=-\int_{0}^{t} \kappa(t-\tau) C_{T}(k, \tau) d \tau\]

Take a closer look at the memory kernel

\[\kappa(t)=\frac{\langle R(t) \mid R(0)\rangle}{\langle A \mid A\rangle}=\frac{\left\langle e^{i \mathcal{Q} \mathcal{L} t} i \mathcal{Q} \mathcal{L} J_{T}(k, 0) \mid i \mathcal{Q L} J_{T}(k, 0)\right\rangle}{\left\langle J_{T}(k, 0) \mid J_{T}(k, 0)\right\rangle}\]

The normalization factor is simply \(\frac{\beta m}{\rho}\). By writing the projection operator \(\mathcal{Q}\) as \(1-\mathcal{P}\) and eliminating terms, this can be written

\[\begin{aligned} \kappa(t) &=\frac{\beta m}{\rho}\left\langle e^{i(1-\mathcal{P}) \mathcal{L} t} i \mathcal{L} J_{T}(k, 0) \mid i \mathcal{L} J_{T}(k, 0)\right\rangle \\ &=\frac{\beta m}{\rho}\left\langle\dot{J}_{T}(k, t)\left|e^{i(1-\mathcal{P}) \mathcal{L t}}\right| \dot{J}_{T}(k, t)\right\rangle \end{aligned}\]

The equation of motion for the transverse current can be written as [5]

\[\dot{J}_{T}(\vec{k}, t)=i \frac{k}{m} \sigma^{z x}(\vec{k}, t)\]

where \(\sigma^{z x}(\vec{k}, t)\) is the zx-component of the microscopic stress tensor

\[\sigma^{z x}(\vec{k})=\sum_{i}\left\{m v_{i, z} v_{i, x}-\frac{1}{2} \sum_{j \neq i} \frac{z_{i j} x_{i j}}{r^{2}} \mathcal{P}_{k}\left(r_{i j}\right)\right\} e^{i k z_{i}}\]

Then the memory kernel becomes

\[\kappa(t)=k^{2} \frac{\beta}{\rho m}\left\langle\sigma^{z x}(\vec{k}, t)\left|e^{i(1-\mathcal{P}) \mathcal{L} t}\right| \sigma^{z x}(\vec{k}, t)\right\rangle=k^{2} \nu(k, t)\]

where \(\nu(k, t)\) is defined as

\[\nu(k, t)=\frac{\beta}{\rho m}\left\langle\sigma^{z x}(\vec{k}, t) e^{i(1-\mathcal{P}) \mathcal{L t}} \sigma^{z x}(\vec{k}, t)\right\rangle\]

This demonstrates that the memory kernel is proportional to \(k^{2}\). Then the transverse current correlation function can be written

\[\dot{C}_{T}(k, t)=-k^{2} \int_{0}^{t} \nu(k, t-\tau) C_{T}(k, \tau) d \tau\]

The memory kernel is the key element that links the two limits. In general, the presence of the propagator \(e^{i \mathcal{Q} L t}\) makes it very difficult to evaluate \(\nu(k, t)\) explicitly. However, the presence of \(\mathcal{Q}\) indicates that we can separate out fast and slow motions and use this to construct a form for \(\nu(k, t)\) that will bridge the short and long time limits. To find this form, the viscoelastic model starts my assuming that the memory kernel has an exponential form:

\[\nu(k, t)=\nu(k, 0) \exp \left[\frac{-t}{\tau_{M}(k)}\right]\]

where \(\tau_{M}(k)\) is the Maxwell relaxation time, discussed above. Before using this function, it is necessary to specify the values of the two parameters, \(\nu(k, 0)\) and \(\tau_{M}(k)\). These can be found by taking the short and long time limits of the GLE and comparing them to the short-time expansion and hydrodynamic results, respectively.

The Short Time Limit

The value of \(\nu(k, t)\) at short times can be obtained by comparing the GLE at time \(t=0\) to the short time expansion of the transverse correlation function. To find the GLE at time \(t=0\), take its time derivative

\[\begin{aligned} \ddot{C}_{T}(k, t) &=-k^{2} \frac{d}{d t} \int_{0}^{t} \nu(k, t-\tau) C_{T}(k, \tau) d \tau \\ &=-k^{2} \nu(k, 0) C_{T}(k, t) \\[4pt] &=-k^{2} \nu(k, 0) C_{T}(k, 0) \\[4pt] &=-k^{2} \nu(k, 0) v_{o}^{2} \end{aligned}\]

The first two terms of the short time expansion of the correlation function are

\[C_{T}(\vec{k}, t)=v_{o}^{2}\left(1-\omega_{T}^{2} \frac{t^{2}}{2}\right)+\cdots\]

The second derivative of this expansion gives

\[\ddot{C}_{T}(k, 0)=-k^{2} c_{t}^{2}(k) v_{o}^{2}\]

Comparison of equations (4.38) and (4.35) shows that

\[\nu(k, 0)=c_{t}^{2}(k)=\frac{G(k)}{\rho m}\]

Further, we see that in this limit the material supports propagating waves. The form of the waves can be found by solving the differential equation Eq.(4.38), and is given by

\[C_{T}(k, t)=v_{o}^{2} \cos \left(\omega_{t} t-k x\right)\]

where \(\omega=k c_{t}(k)\) and the speed of the waves are \(c_{t}=\sqrt{\frac{G(k)}{\rho m}}\).

The Hydrodynamic Limit

The value of \(\nu(k, t)\) at long times can be obtained by comparing the hydrodynamic equation to the long time limit of the GLE for \(C_{T}(k, t)\) :

\[\dot{C}_{T}(k, t)=-k^{2} \int_{0}^{t} \nu(k, t-\tau) C_{T}(k, \tau) d \tau\]

To take the long time limit of this equation, note that the memory function will generally be characterized by some relaxation time \(\tau_{\kappa}\). When the time \(t\) is much greater than this relaxation time, the major contribution to the integral will come when \(t \sim \tau\). Therefore, we can approximate \(C_{T}(k, \tau) \sim C_{T}(k, t)\). With this approximation, the correlation function can be taken out of the integral in the GLE:

\[\dot{C}_{T}(k, t)=-k^{2} C_{T}(k, t) \int_{0}^{\infty} \nu(k, t-\tau) d \tau\]

where the integration limit has been extended to \(\infty\) to indicate that we are taking the long time limit.

This result should be identical to the hydrodynamic solution in the long time and long wavelength limit. By taking the long wavelength limit \((k \rightarrow 0)\) and comparing to the hydrodynamic result (Eq.(4.35)), we see that this only holds when:

\[\int_{0}^{\infty} \nu(k, t-\tau) d \tau=\nu=\frac{\eta}{\rho m}\]

The Viscoelastic Solution

We now have the information we need to construct the explicit form of the viscoelastic memory kernel.

\[\nu(k, t)=\nu(k, 0)\left[\frac{-t}{\tau_{M}(k)}\right]\]

From the short time limit, we found that \(\nu(k, 0)=\frac{G(k)}{\rho m}\), which allows us to write

\[\nu(k, t)=\frac{G(k)}{\rho m} \exp \left[\frac{-t}{\tau_{M}(k)}\right]\]

From the long time limit, we know that

\[\int_{0}^{\infty} \nu(k, t-\tau) d \tau=\frac{\eta}{\rho m}\]

Now, plug in the exponential memory kernel for \(k=0\)

\[\int_{0}^{\infty} \frac{G(0)}{\rho m} \exp \left[\frac{-t}{\tau_{M}(0)} d \tau=\frac{\eta}{\rho m}\right]\]

The elastic modulus has no time dependence, so it can be taken out of the integral

\[\frac{G(0)}{\rho m} \int_{0}^{\infty} \exp \left[\frac{-t}{\tau_{M}(0)} d \tau=\frac{\eta}{\rho m}\right]\]

Finally, evaluate the integral to find the Maxwell relaxation time at \(k=0\). It is reasonable to assume that the Maxwell relaxation time remains constant over all \(k\) values. Therefore, the Maxwell relaxation time can be written as the ratio of the shear viscosity coefficient of the liquid to the modulus of rigidity of the elastic solid at \(k=0\).

\[\tau_{M}=\frac{\eta}{G(0)}\]

When \(\tau_{M}\) is small compared to the time \(t\), the viscosity term dominates and the system will behave as a viscous liquid. However, when \(\tau_{M}\) is large compared to the time \(t\), the system does not have time to respond to a stimulus as a viscous liquid. The modulus of rigidity dominates, and the material will behave as an elastic solid, supporting propagating shear waves.

Finally, we can use the Maxwell relaxation time to write the explicit form of the viscoelastic memory kernel.

\[\nu(k, t)=\frac{G(k)}{\rho m} \exp \left[\frac{-\eta t}{G(0)}\right]\]

With this memory kernel in hand, we can now go on to find an explicit solution to the transverse current correlation function.

To find the equation for the viscoelastic wave, we first find the Laplace transform of the transverse current correlation function

\[\hat{C}_{T}(k, s)=\frac{C(k, t=0)}{s+k^{2} \hat{\nu}(k, s)}=\frac{v_{o}^{2}}{s+k^{2} \hat{\nu}(k, s)}\]

Now, solve this equation using the exponential memory kernel \(\nu(k, t)=\frac{G(k)}{\rho m} \exp \frac{-t}{\tau_{M}}\). The Laplace transform of an exponential function is well defined

\[\mathcal{L}[\exp [-\alpha t u(t)]]=\frac{1}{s+\alpha}\]

Therefore, the Laplace transform of the viscoelastic memory kernel is

\[\hat{\nu}(k, s)=\frac{\nu(k, 0)}{s+\tau_{M}^{-1}}\]

Plug this into the the Laplace transform of the transverse current correlation function

\[\begin{aligned} \hat{C}_{T}(k, s)=\frac{v_{o}^{2}}{s+k^{2}\left(\frac{\nu(k, 0)}{s+\tau_{M}^{-1}}\right)} \\ \hat{C}_{T}(k, s)=\frac{v_{o}^{2}\left(s+\tau_{M}^{-1}\right)}{s\left(s+\tau_{M}^{-1}\right)+k^{2} \nu(k, 0)} \end{aligned}\]

Since function is quadratic, it is relatively easy to find the reverse Laplace transform, using the same method as that presented in section 4.2.C.2, or reference [1].

\[C_{T}(k, t)=v_{o}^{2} \frac{1}{\lambda_{+}+\lambda_{-}}\left(e^{-\lambda+t}-e^{-\lambda_{-} t}\right)\]

where the eigenvalues \(\lambda_{\pm}\)are given by the solutions to the quadratic equation \(s^{2}+\tau_{M}^{-1} s+k^{2} \nu(k, 0)\) :

\[\lambda_{\pm}=-\frac{-\tau_{M}^{-1}}{2} \pm \sqrt{\frac{\tau_{M}^{-2}-4 k^{2} \nu(k, 0)}{4}}\]

Complex eigensolutions exist if

\[\frac{1}{\tau_{M}^{2}}<4 k^{2} \nu(k, 0)\]

Recall that \(\nu(k, 0)=c_{t}^{2}(k)\). Then we can rewrite the above inequality in terms of the wavenumber

\[k>\frac{1}{2 \tau_{M} c_{t}(k)}\]

Define the critical wavenumber, \(k_{c}=\left(2 \tau_{M} c_{t}(k)\right)^{-1}\). For more information on the viscoelastic approximation and its application to transverse current, please see Chapter 6 of Molecular Hydrodynamics by Jean-Pierre Boon and Sidney Yip [3] and chapter 3 and chapter 6 of Dynamics of the Liquid State by Umberto Balucani [5].