## Microscopic Origin of Shear Relaxation in Strongly Coupled Yukawa Liquids

∗

Ashwin J. and Abhijit Sen Institute for Plasma Research, Bhat, Gandhinagar-382428, India. (Dated: August 4, 2014)

We report accurate molecular dynamics calculations of the shear stress relaxation in a twodimensional strongly coupled Yukawa liquid over a wide range of the Coulomb coupling strength Γ and the Debye screening parameter κ . Our data on the relaxation times of the ideal- , excess- and total shear stress auto-correlation ( τ id M , τ ex M , τ M respectively) along with the lifetime of local atomic connectivity τ LC leads us to the following important observation. Below a certain crossover Γ ( c κ ), τ LC → τ ex M , directly implying that here τ LC is the microscopic origin of the relaxation of excess shear stress unlike the case for ordinary liquids where it is the origin of the relaxation of the total shear stress. At Γ &gt;&gt; Γ ( c κ ) i.e. in the potential energy dominated regime, τ ex M → τ M meaning that τ ex M can fully account for the elastic or 'solid like' behavior.

PACS numbers: 52.27.Lw,51.35.+a,52.65.Yy

The Yukawa liquid is routinely used to model a wide variety of strongly coupled systems such as laboratory and astrophysical dusty plasma systems, charged colloids and liquid metals [1, 2]. The particle interaction potential used in this liquid has the form φ r ( ) = Q 2 4 πglyph[epsilon1] 0 r e -r λ D . Here Q and λ D refer to the particle charge and Debye shielding distance respectively. The thermodynamic state point of the Yukawa liquid is completely characterized by two dimensionless quantities namely, the screening parameter κ = a λ D and the coupling strength Γ = Q 2 4 πglyph[epsilon1] 0 ak B T . Here a is the Wigner-Seitz radius such that πa n 2 = 1 with n being the areal number density. In a strongly coupled Yukawa liquid (SCYL), the average Coulomb interaction energy can exceed the average kinetic energy per particle thus leading to Γ &gt; 1. A direct consequence of this is the emergence of solid-like features such as sustaining low frequency shear modes [3] originally predicted in theoretical works [4, 5] and later realized in laboratory experiments [6]. This makes SCYL a good model system to study a range of collective phenomena in dusty plasma systems. SCYL's have also been shown to be excellent test-beds for modeling hydrodynamic flows [79], self-organization phenomena such as clustering [10] and lane formation [11] in complex plasmas. A significant amount of numerical work has also been done on the viscosity measurements in SCYL using both equilibrium [12-14] and non-equilibrium molecular dynamics simulations [15, 16] . These studies have confirmed the existence of a viscosity minimum at some crossover value Γ ( c κ ). The minimum arises due to the competition between ideal and the excess part of the stress tensor. In later works, visco-elasticity was also quantified using both experiments [17, 18] and numerical simulations [19] where the crossover frequency for the real and imaginary parts of the complex viscosity was shown to be empirically related to the inverse Maxwell time. It was shown in Ref [19] that this crossover frequency develops a maximum at some Γ ( c κ ) implying a minimum in the Maxwell time τ M at the same Γ ( c κ ). However, a systematic study of the microscopic origin of shear relaxation and Maxwell time in these liquids at various (Γ , κ ) values is still lacking. The aim of this letter is to fill this gap and provide from first principles, an atomistic study that would provide answers to the following two fundamentally important questions: (a) What is the connection between the microscopic world and the macroscopic shear relaxation in these SCYL? (b) Is τ M the duration of dominant elastic response for these liquids especially at Γ glyph[lessmuch] Γ ( c κ )? To address the first question, we have calculated the lifetime of local atomic connectivity τ LC [20] and found that it converges to the relaxation of the excess part of shear stress auto-correlation τ ex M at Γ &lt; Γ ( c κ ). Since τ LC corresponds to the time duration where the topology of nearest neighbors remains intact, the fact that τ ex M → τ LC for Γ &lt; Γ ( c κ ) directly indicates that τ ex M is the duration of dominant elastic response at these temperatures. This is different from the results reported in Ref [19] where τ LC was shown to converge to τ M below the cross-over temperature for a range of ordinary liquids. It is important to note that none of those liquids showed a nonmonotonous behavior in τ M with temperature. We also show that at Γ &gt; Γ ( c k ), τ LC deviates from τ ex M directly indicating a crossover from a kinetic regime to potential energy dominated regime. Our study further shows that the infinite frequency shear modulus G ∞ does not have a minimum at any temperature as opposed to the viscosity data which has a well known minimum at Γ = Γ ( c κ ). In the following we provide the details of our numerical work and also explain the procedure used to extract quality data.

Numerical simulations: We have performed molecular dynamics (MD) simulations on a two-dimensional (2D) Yukawa liquid in a canonical ensemble under periodic boundary conditions. All distances are normalized to Wigner-Seitz radius a , energies are normalized to Q 2 4 πglyph[epsilon1] 0 a and times are normalized to ω -1 pd . Here ω pd is the 2D nominal plasma frequency given by ω pd = √ Q 2 2 πglyph[epsilon1] 0 ma 3 .

The simulation box contains 5016 particles at a reduced number density n = π -1 . The dimensions of rectangular box were chosen to be 125.705291 × 125.358517 which allows for the formation of a perfect triangular lattice below the freezing transition. The interaction potential is truncated smoothly to zero along with its first two derivatives by employing a fifth order polynomial as a switching function in the range ( r m &lt; r &lt; r c ) where r m and r c are the inner and the outer cutoff respectively. We chose r m and r c subject to the criteria: φ r ( m ) ≈ 2 27 . × 10 -6 and φ r ( c ) ≈ 1 49 . × 10 -7 thus ensuring negligible perturbation to the bare Yukawa potential. A Nose-Hoover thermostat [21] with a time constant of 1 √ 2 is employed to maintain the temperature at a desired Γ. To improve statistics, we have averaged our data over an ensemble of 3200 independent realizations. This was necessary to reduce the fluctuations present in the long time tail of the stress relaxation function originating from the long range nature of the interaction potential. We now proceed to explain the computation of various dynamical quantities reported in this paper.

The stress relaxation function used in our work is the auto-correlation of the shear stress tensor and is given as

$$G ( t ) = \frac { 1 } { A k \cap T } \langle \sigma _ { x y } ( t ) \sigma _ { x y } ( 0 ) \rangle$$

$$\L ^ { ( t ) } = \frac { 1 } { A k _ { B } T } \langle \sigma _ { x y } ( t ) \sigma _ { x y } ( 0 ) \rangle \quad \quad ( 1 )$$

with initial value of this auto-correlation giving the infinite frequency shear modulus G ∞ = G (0). The angular brackets here denote the average over the entire ensemble, with A being the area and σ xy ( ) t the microscopic stress tensor being defined as

$$\sigma _ { x y } ( t ) = \sum _ { i = 1 } ^ { N } m _ { i } v _ { i } ^ { x } ( t ) v _ { i } ^ { y } ( t ) - \sum _ { i = 1 } ^ { N } \sum _ { j > i } ^ { N } x _ { i j } ( t ) y _ { i j } ( t ) \frac { \phi ^ { \prime } ( r _ { i j }, t ) } { r _ { i j } ( t ) } \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
T h o \ f o w t \text{ to the thead} \text{ in and thead} \text{ in to}$$

The first term on the right hand side has purely kinetic origins and is the dominating term at high temperatures whereas the second term is the excess part and has its origins in particle interactions. The second term is named excess as it would be absent in an ideal gas (gas of 'hard spheres' ) and arises only in a real gas. Figure 1 shows the normalized stress relaxation function G t ( ) of a 2D Yukawa liquid at various coupling strengths Γ &gt; Γ ( c κ ). It is clear that at short times, G t ( ) has a zero slope meaning that liquid response is dominantly elastic. This is followed by a region of fast decay where both the elastic and viscous effects are comparable. At large times, G t ( ) has become much smaller (within statistical noise) indicating a regime dominated by viscous response and negligible elastic effects. It is also seen from the figure that the relaxation time increases with Γ implying that elastic response will dominate for longer times as Γ increases. This is mainly due to growing structural order as shown in the inset of Figure 1. Next we turn our attention to the hydrodynamic (or zero frequency) shear viscosity η which is an important dynamical property responsible for viscous dissipation in liquids.

FIG. 1. (color online). Normalized stress relaxation function G t ( ) for a Yukawa liquid at various Γ. The dashed line shows the location of 1 /e falling time. It is clear from the figure that the relaxation time increases with Γ for Γ &gt; Γ ( c κ ). Note that the relaxation time (or e-folding time) increases with Γ. [Inset: shows the growing hexagonal order in the liquid with Γ.]

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0141v1-figure-750d7bfa-e2e9-45b9-8f02-a059ad5f6a31.jpg)

All liquids tend to gradually resist the deformation imposed by an external shear stress. A measure of this resistance is the hydrodynamic shear viscosity η (bulk viscosity in case of a compressive external force). For a liquid at thermal equilibrium, this shear viscosity can be obtained by the long time integral of the stress relaxation function -a procedure well known as the Green-Kubo formula [22].

$$\eta = \int _ { 0 } ^ { \infty } G ( t ) d t$$

In calculating η from Eq. (3), the upper limit of the integral is set to the first zero crossing time of the ensemble averaged shear relaxation function G t ( ). In Figure 2, we have shown our data on the shear viscosity measurement at three different values of κ . Our results for η are both qualitatively and quantitatively similar to the existing report on 2D equilibrium simulations [14]. Next we show our data on the G ∞ calculation in Figure 3 at various κ . It is interesting to see that unlike shear viscosity, G ∞ does not have a minimum at any temperature. The contributions coming from both the ideal and excess part of G ∞ are shown in the inset of Figure 3. At high temperatures (low Γ), kinetic effects are dominant and the ideal part G id ∞ is the dominating term and becomes a major contributor to the overall G ∞ . At low temperatures (high Γ) the excess term G ex ∞ begins to dominate and the ideal term becomes very small. It is interesting to note that the excess part G ex ∞ saturates as Γ increases.

FIG. 2. (color online). Hydrodynamic viscosity (normalized to η 0 = a nmω 2 pd ) computed from Eq. (3) at various κ . The minimum at Γ = Γ ( c κ ) is due to the competition of the ideal part and the excess part of the stress tensor and is well known in the literature [12, 14].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0141v1-figure-dfa1a914-7735-4029-99bb-6b8fce19968e.jpg)

FIG. 3. (color online). Infinite frequency shear modulus G ∞ vs. Γ at various κ . Unlike viscosity which has a minimum, G ∞ does not have a minimum at any temperature. [Inset: The total, ideal and excess contributions for the case κ = 0 6 . only. It is clear that at low Γ, the entire contribution comes from the ideal term whereas at high Γ, the excess term is responsible for the entire G ∞ .]

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0141v1-figure-ad543ff2-0cb3-48ca-a7e9-4d05806a5e16.jpg)

In the following, we will show the effect of quantities addressed so far in the determining the Maxwell shear stress relaxation time τ M .

It is generally accepted that all liquids are visco-elastic in nature meaning that their mechanical response to an external force will be both viscous and elastic at the same time. This can be explained through the concept of the Maxwell relaxation time scale τ M such that at times t &lt;&lt; τ M , the response of the liquid will be dominantly elastic ( reversible ) and at t &gt;&gt; τ M the response will be dominantly viscous ( irreversible ). At intermediate time

FIG. 4. (color online). Maxwell relaxation time τ M as calculated from Eq. (4). The data shows a clear minimum around a crossover Γ ( c κ ) as previously shown in Ref [19]. We explain this non-monotonous behavior by measuring independently the relaxation times of the ideal and excess parts of the shear stress auto-correlation function [See text and also Fig. 5 for details].

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0141v1-figure-122302c3-367d-4a10-b15b-1fde8b5c1021.jpg)

scales, both elastic and viscous features will be comparable. The ratio of the hydrodynamic shear viscosity η to the infinite frequency shear modulus G ∞ defines an average shear relaxation time or the Maxwell time [23],

$$\tau _ { M } = \frac { \eta } { G _ { \infty } } = \frac { \int _ { 0 } ^ { \infty } G ( t ) d t } { G _ { \infty } } \, \quad \quad \quad ( 4 )$$

FIG. 5. (color online). The relaxation times of the stress relaxation function G t ( ) along with its ideal G id ( ) and excess t parts G ex ( ) for the case of t κ = 0 6. . It is interesting to note that while the excess term τ ex M continues to rise all the way to the freezing transition, the ideal part saturates to a constant value around Γ ( c κ ) indicating the onset of potential energy influenced regime. [Inset: Data plotted on a log-linear axis. A dashed exponential is drawn to aid the eye of the reader.]

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0141v1-figure-9a8fb4ad-4652-46a9-a816-ea0374bf66aa.jpg)

FIG. 6. (color online). Comparison of τ LC with time scales τ M , τ ex M and τ id M for the cases (a) κ = 0 6, (b) . κ = 1 0, (c) . κ = 2 0. . At Γ &lt; Γ ( c κ ), τ ex M becomes equal to the lifetime of atomic connectivity τ LC . At Γ &gt; Γ ( c κ ) these two time scales deviate as the liquid enters a landscape dominated regime where the interaction between local networks becomes significant and leads to cancellation of long range elastic fields [20]

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0141v1-figure-e94083c7-2448-4800-b3ce-0dc93670b06c.jpg)

.

The normalized integral given above provides the duration of the auto-correlation and is hence a good measure of the time duration of G t ( ). Figure 4 shows a plot of τ M vs. coupling strength Γ for three values of κ . Our data shows the minimum in τ M at Γ = Γ ( c κ ) similar to the results reported in Ref [19]. This means that if one considers τ M (the relaxation of the total stress auto-correlation) as a measure of elastic response then it would imply that elasticity will persist for longer times as the weakly coupled limit or Γ → 0 is approached. We thus propose that τ M would be an incorrect measure of elasticity in these liquids as Γ → 0. To seek for an appropriate measure of elasticity we turn our attention to the components of the total stress auto-correlation G t ( ) namely the ideal part G id ( ) and the excess part t G ex ( ). t There is a cross term as well but its value is negligible at all values of ( κ, Γ). In figure 5 we show the relaxation times of the total total stress auto-correlation and its two components over the range of Γ that spans from the very weakly coupled fluid limit Γ = 0 25 to the highest Γ just below freez-. ing transition. As the liquid nears freezing (Γ → Γ m ), τ ex M → τ M and the shear relaxation is dominantly due to the particle interactions. We find that the relaxation time of the G ex ( ) denoted here by t τ ex M serves as a good measure of elasticity as it goes to zero monotonically as Γ → 0. On the high Γ side, it rises exponentially with Γ (see inset: Figure 5) as the liquid approaches the freezing transition. On the other hand, as Γ → 0, relaxation time of G id ( ) denoted here as t τ id M becomes dominant and is entirely responsible for the origin of τ M . This means that in the weakly coupled fluid limit, stress relaxation is purely kinetic and has very little contribution coming from particle interactions. It is very important to note here that our assertion on τ ex M being the right indicator for elasticity is not just based on the right asymptotic behavior with Γ but also due to its convergence to the lifetime of local atomic connectivity τ LC as explained in the following section.

In solids, phonons are weakly scattered (long lived) and hence they are the microscopic origin for vibrations. In liquids however, they exhibit highly marginalized behavior as they are strongly scattered due to lack of any underlying long range structural order. As a result phonons cannot be used to explain the microscopic origins of viscosity and the infinite frequency shear modulus. One is thus led to the following fundamentally important question: What is the microscopic origin for shear relaxation in SCYL? To answer this question, we measure the lifetime of the local atomic connectivity τ LC [20] in SCYL. To calculate τ LC one first assigns a set of bonds between a central atom and its nearest neighbors at some reference time t 0 . The nearest neighbors here correspond to atoms situated at a distance lesser than the location of the first minima in the radial distribution function. Once a reference state is assigned we track the change in local atomic connectivity i.e. the central atoms losing or gaining some neighbors. As the simulation proceeds the average coordination number will reduce and τ LC is extracted as the time duration (elapsed since t 0 ) in which the average coordination falls by 1. To improve statistics, we averaged our data over several independent realizations. In this sense, τ LC is the average time during which any atom will lose one neighbor. It was shown in Ref [20] that above a certain crossover temperature, τ LC

converges to the Maxwell time τ M and is thus the microscopic origin of shear relaxation in a range of liquids. It should be noted that none of the liquids used in that work showed a non-monotonic behavior in τ M . In Figure 6 we show our data on τ LC calculated at various values (Γ , κ ) and make the following interesting observation. At Γ &lt; Γ ( c κ ), τ LC → τ ex M making τ LC the microscopic origin of excess part of stress at these temperatures. Since τ LC corresponds to the time duration where the topology of nearest neighbors remains intact it is also the duration of the dominant elastic or 'solid like' response. These observations directly imply that the relevant shear relaxation time scale which quantifies the elastic response in these liquids is τ ex M which becomes equal to τ LC at Γ &lt; Γ ( c κ ) and approaches the value of τ M at Γ &gt; Γ ( c κ ). The reader should also note that at Γ &gt; Γ ( c k ), τ LC deviates from τ ex M directly indicating a crossover from a kinetic regime to potential energy dominated regime. This is a regime where interaction between local networks becomes significant and leads to cancellation of long range elastic fields [20].

Summary: We have calculated the infinite frequency shear modulus G ∞ , Maxwell relaxation times ( τ M , τ id M and τ ex M ) and the lifetime of the local atomic connectivity τ LC of a 2D SCYL using accurate MD simulations. The simulations were done under a canonical ensemble and quality data is obtained by averaging the runs over 3200 independent realizations. Our major finding is the microscopic origin of shear relaxation in SCYL. We find that the lifetime of topology of nearest neighbors τ LC is responsible for the relaxation of the excess part of the shear stress relaxation τ ex M . Wealso provide a solution to the riddle concerning the minimum in τ M which renders it unsuitable to quantify the elastic properties of a SCYL approaching the weakly coupled fluid limit Γ → 0. We resolve this issue by proposing that the correct object that quantifies elastic response of such strongly coupled liquids is the relaxation time of the excess part of the shear stress auto-correlation τ ex M which goes to zero as Γ → 0 and becomes exponential as the liquid nears the freezing transition. Our assertion is backed by careful measurements of the lifetime of the local atomic connectivity τ LC which clearly shows that τ ex M → τ LC below the crossover Γ ( c κ ) directly indicating that τ ex M is the duration of dominant elastic response. The relaxation of the ideal part of shear stress auto-correlation τ id saturates with Γ &gt; Γ ( c κ ) implying a crossover from the kinetic regime to the potential energy regime. We also find that the infinite frequency shear modulus G ∞ does not have a minimum at any temperature as opposed to the viscosity data which has a well known minimum around Γ ( c κ ).

We believe the present work is fundamentally important as it provides the microscopic origin of shear relaxation time scales in a SCYL. The results presented here may also apply to other models of visco-elastic liquids exhibiting a non-monotonic behavior of Maxwell time with temperature. It will be very interesting to extend the present study to the binary super-cooled Yukawa liquids.

Discussions with P.K. Kaw and R. Ganesh are gratefully acknowledged. This work was supported under the INSPIRE faculty program, Department of Science and Technology, Ministry of Science and Technology, Government of India.

- ∗ ashwin@ipr.res.in
- [1] G. E. Morfill and A. V. Ivlev, Rev. Mod. Phys. 81 , 1353 (2009).
- [2] V. Fortov, A. Ivlev, S. Khrapak, A. Khrapak, and G. Morfill, Physics Reports 421 , 1 (2005).
- [3] H. Ohta and S. Hamaguchi, Phys. Rev. Lett. 84 , 6026 (2000).
- [4] P. K. Kaw and A. Sen, Physics of Plasmas 5 , 3552 (1998).
- [5] G. Kalman, M. Rosenberg, and H. E. DeWitt, Phys. Rev. Lett. 84 , 6030 (2000).
- [6] J. Pramanik, G. Prasad, A. Sen, and P. K. Kaw, Phys. Rev. Lett. 88 , 175001 (2002).
- [7] Ashwin J. and R. Ganesh, Phys. Rev. Lett. 104 , 215003 (2010).
- [8] Ashwin J. and R. Ganesh, Phys. Rev. Lett. 106 , 135001 (2011).
- [9] Ashwin J. and R. Ganesh, Physics of Fluids (1994present) 24 , 092002 (2012).
- [10] H. Totsuji, C. Totsuji, and K. Tsuruta, Phys. Rev. E 64 , 066402 (2001).
- [11] A. Wysocki and H. L¨ owen, Phys. Rev. E 79 , 041408 (2009).
- [12] T. Saigo and S. Hamaguchi, Physics of Plasmas 9 , 1210 (2002).
- [13] G. Salin and J.-M. Caillol, Phys. Rev. Lett. 88 , 065002 (2002).
- [14] B. Liu and J. Goree, Phys. Rev. Lett. 94 , 185002 (2005).
- [15] K. Y. Sanbonmatsu and M. S. Murillo, Phys. Rev. Lett. 86 , 1215 (2001).
- [16] Z. Donk´ o, J. Goree, P. Hartmann, and K. Kutasi, Phys. Rev. Lett. 96 , 145003 (2006).
- [17] Y. Feng, J. Goree, and B. Liu, Phys. Rev. Lett. 105 , 025002 (2010).
- [18] Y. Feng, J. Goree, and B. Liu, Phys. Rev. E 85 , 066402 (2012).
- [19] J. Goree, Z. Donk´, o and P. Hartmann, Phys. Rev. E 85 , 066401 (2012).
- [20] T. Iwashita, D. M. Nicholson, and T. Egami, Phys. Rev. Lett. 110 , 205504 (2013).
- [21] W. G. Hoover, Phys. Rev. A 31 , 1695 (1985).
- [22] J. P. Hansen and I. McDonald, Theory of Simple Liquids: with applications to soft matter , 4th ed. (Academic, Oxford, 2013).
- [23] R. D. Mountain and R. Zwanzig, The Journal of Chemical Physics 44 (1966).
Ashwin J. Abhijit Sen
Microscopic Origin of Shear Relaxation in Strongly Coupled Yukawa Liquids
We report accurate molecular dynamics calculations of the shear stress
relaxation in a two-dimensional strongly coupled Yukawa liquid over a wide
range of the Coulomb coupling strength $\Gamma$ and the Debye screening
parameter $\kappa$. Our data on the relaxation times of the ideal- , excess-
and total shear stress auto-correlation ($\tau^{id}_M, \tau^{ex}_M, \tau_M$
respectively) along with the lifetime of local atomic connectivity $\tau_{LC}$
leads us to the following important observation. Below a certain crossover
$\Gamma_c(\kappa)$, $\tau_{LC} \rightarrow \tau^{ex}_M$, directly implying that
here $\tau_{LC}$ is the microscopic origin of the relaxation of excess shear
stress unlike the case for ordinary liquids where it is the origin of the
relaxation of the total shear stress. At $\Gamma >> \Gamma_c(\kappa)$ i.e. in
the potential energy dominated regime, $\tau^{ex}_M\rightarrow \tau_M$ meaning
that $\tau^{ex}_M$ can fully account for the elastic or "solid like" behavior.