Viscoelasticity and dynamical gaps: rigidity in crystallization and glass-forming liquids

TL;DR

This paper explores how viscoelasticity and dynamical gaps relate to rigidity in glass-forming liquids and crystals, proposing a new method to assess rigidity through dynamical structure factors and simulations.

Contribution

It introduces a novel approach linking viscoelasticity and dynamical gaps to rigidity, including defining a dynamical average effective coordination number.

Findings

Dynamical gaps in transverse vibrational modes relate to rigidity.

Numerical simulations support the link between dynamical gaps and relaxation times.

Viscoelasticity can serve as an alternative measure of glass rigidity.

Abstract

Rigidity plays an important role on the relaxation properties of glass forming melts, yet it is usually determined from the average co-ordination number through the chemical composition. A discussion is presented on how viscoelasticity can be used as an alternative way to determine glass rigidity and to give clues about the relaxation processes. In particular, it is shown that the transverse current dynamical structure factor of dense glass and crystal forming fluids contain rich information about rigidity that can be related with the presence of a dynamical-gap for transversal vibrational modes. As a result, the number of floppy modes can be related with the size of the dynamical gap and with the liquid relaxation time. From this, a dynamical average effective coordination number can be defined. As an example, numerical simulations for hard-disks in a dense fluid phase are provided. Finally, a discussion is presented on the need to improve glass viscoelasticity models to describe consistently non-exponential stress and strain relaxation.

Equations

$$S(k,\omega )={\int }_0^{\infty }dt{e}^{-ı\omega t}C(k,t).$$

$$C(k,t)=⟨J_T^{\ast }(k,t)J_T(k,0)⟩,$$

$$J_T(k,t)=\frac{1}{\sqrt{2N}k}\sum _{i=1}^N\mathbf{k}\times {\mathbf{v}}_i(t)\exp \left(ı\mathbf{k}\cdot {\mathbf{r}}_i(t)\right).$$

$$N_f(k_c)\approx 2{\int }_0^{k_c}4\pi k^{2}dk=\frac{8\pi }{3}{k}_c^3$$

$$f\approx \frac{2}{3}{\left(\frac{k_c}{k_D}\right)}^3$$

$$f\approx \frac{2}{3{\omega }_D^3}\left(\frac{1}{{\tau }^3(T)}\right)$$

$$<r>=\frac{12}{5}\left(1-\frac{1}{3}{\left(\frac{k_c}{k_D}\right)}^3\right)$$

$$<r>=6\left(1-\frac{2}{3}{\left(\frac{k_c}{k_D}\right)}^3\right)$$

Full text

One of the most important problems in glass formation is the understanding of structural relaxation mechanisms near glass transition [1, 2, 3, 4, 5, 6, 7], as well as how supercooled liquid relaxation wins over crystal nucleation [8]. Certainly a huge body of research has been focused on the subject (see [9, 10, 11, 12, 13, 14, 15, 16] and in particular [17] and references therein), yet there is not a definitive consensus. As is well known, experiments and simulations still have many feats to achieve [18].

Relaxation is related with one of the key features of glass formation: the minimal speed required in order to make a glass, a property known as glass forming ability [3]. Phillips and Thorpe’s rigidity theory gives good insights on how this problem is related to network topology [19, 20]. These ideas can be extended to include non-directional potentials [21, 22, 23]. Eventually, the pioneering work of Gupta and Mauro [1] led ridigidty theory to produce a new and highly accurate viscosity model, known in the literature as the Mauro-Yue-Ellison-Gupta-Allison Model (MYEGA) [24]. This allows us to understand the chemical composition and temperature effects on the viscosity of glass-former melts [24]. As a result, we are closer than ever to an age of glasses obtained by design [25, 26] . This goes together with the advances made by Micoulaut and Bauchy who had extensively studied how to define rigidity for realistic potentials (see for instance [27]). Previous efforts were made in simple models by Huerta et. al.[21, 22, 28, 29]. Stochastic models also provided a different pathway to include chemical composition effects [30, 31]. From an experimental point of view, Boolchand and coworkers have extensively studied the optical, mechanical and thermodynamical properties in terms of rigidity [32, 33]. Theoretical models allows an understanding of some general properties of thermodynamics in terms of rigidity [34, 35] and there are suggestions of a connection with the boson peak [36, 37, 38].

As a matter of fact, any symmetry-breaking thermodynamic phase transition involves the development of some kind of generalized rigidity by the system [39]. This allows the given system to preserve the phase order against thermal fluctuations [39]. In spite of this fundamental character, it is surprising to find that in general such observation is not emphasized when phase transitions are studied. A fluid is different from a solid precisely due to its rigidity, and thus a simple first-order fluid-solid phase transition must also contain a rigidity transition as its main signature. Moreover, the lack of rigidity is the defining property of a Newtonian fluid, i.e., the absence of elastic behavior against shear stress. This leads to the absence of transversal waves in a fluid.

The main aim of this work is to emphasize the need to decode how the rigidity of glass forming melts depends upon the time and spatial scales in which the system is probed or perturbed. Moreover, to accomplish this feat we need to understand rigidity transitions not only in glasses, but for crystallization and in cluster nucleation at the kinetic spinodal temperature [40].

A lot of knowledge on these aspects could be obtained by looking at the similarities and differences between rigidity in organic and inorganic glasses, all of them above $T_{g}$ and close to the mechanical vitrification point [41]. For organic glasses there is a vast literature concerning flexible and rigid polymer models in which relaxation has been tested [41, 42, 43]. On the experimental side, modulated differential scanning calorimetry (MDSC), dielectric relaxation and rehology measurements have been very succesful for the understanding of relaxation processes [44]. For inorganic glasses, a series of different experiments such as MDSC and Raman scattering among others, as well as computational simulations reveal interesting aspects of the rigidity transition [32, 45, 46, 47]. However, results on chalcogenide glasses frequency-dependent rheology are recent [48, 49, 50, 51].

One may wonder what is the fundamental difference between rigidity in organic and inorganic glasses above $T_{g}$. Many years ago the answer to this question was not clear. Above $T_{g}$ and due to their polymeric nature, organic glasses display transitions from the folded to the stretched chain forms and thus present viscoelasticity [44]. On the other hand, inorganic glass melts were thought to be purely Newtonian fluids [41]. Yet, inorganic glasses, as polymeric systems, were expected to display viscoelasticity. This apparent paradox was solved by the observation made by G. M. Bartenev, who started by adscribing the prominent differences between the $T_{g}$ of inorganic and organic glasses to the much higher rotation flexibility of the C-C bonds [41]. Therefore, it was concluded that viscoelasticity was also possible for inorganic glasses, although happens to be smaller than in its organic counterparts due to their somewhat limited angular bond excursions. This sole fact explained why the viscoelastic response in inorganic glasses remained for a long time unnoticed [41].

For organic glasses, the key to understand the relationship between relaxation and rigidity is given by measuring the viscoelasticity using rheological experiments [44]. In viscoelasticity, the relationship between the stress $\sigma(\omega)$ and strain $\epsilon(\omega)$ is measured as a function of the frequency $\omega$. For $\epsilon(\omega)=\cos\omega t$, we have $\sigma(\omega)=G^{\prime}\cos\omega t+G^{\prime\prime}\sin\omega t$, and thus a complex modulus $G(\omega)$ is obtained [49]. The real part of $G(\omega)$, denoted by $G^{\prime}(\omega)$, is the storage modulus while the imaginary part $G^{\prime\prime}(\omega)$ gives the loss modulus. The phase lag between strain and stress is given by $\tan\theta(\omega)=G^{\prime\prime}(\omega)/G^{\prime}(\omega)$, while a frequency-dependent viscosity is obtained from $\eta(\omega)=\sqrt{G^{\prime\prime}(\omega)^{2}+G^{\prime}(\omega)^{2}}/\omega$.

Above the glass transition and for low frequencies[48, 49, 50, 51], the viscosity $\eta(\omega)$ is strongly frequency-dependent and $G^{\prime}(\omega)<<G^{\prime\prime}(\omega)$. At these frequencies, the system behaves as a viscous fluid as $\eta(\omega)\approx G^{\prime\prime}(\omega)/\omega$. In the limit of high frequencies, denoted by $\omega=\infty$, we have that $G(\infty)\approx G^{\prime}(\infty)>G^{\prime\prime}(\infty)$ and mainly a purely elastic response is obtained. The lesson taken from these rheology experiments is that at high-frequencies, the system is rigid-like, while at low frequencies is non-rigid. Thus, rigidity in melts involves the time-scale in which the system is probed. Moreover, this aspect means that there must also be a transition concerning the propagation of transversal waves. As the dispersion relationship of waves involves $\omega$ as a function of the wavevector $\bm{k}$, is clear that rigidity involves time and space density-density fluctuations. A striking demonstration of this phenomena is the report of transversal-wave branchs in the dynamical structure factor [52, 53, 54]. The transversal part of the dynamical structure factor is defined as [52],

[TABLE]

where $C(k,t)$ is the transversal current density correlation function,

[TABLE]

and the brackets $\langle...\rangle$ represent an ensemble average. The function $J_{T}(k,t)$ is the transversal density current averaged over the different directions of $\bm{k}$ given the wavenumber $k=|\bm{k}|$,

[TABLE]

Here, $\bm{v}_{i}(t)$ and $\bm{r}_{i}(t)$ are the velocity and position of the $i^{th}$ particle of a given system at time $t$. The $1/\sqrt{2}$ factor takes into account the two transverse currents in three-dimensional systems, and is replaced by one in two dimensions.

As an example, in Figure 1 we present the transversal part of the dynamical structure factor $S(k,\omega)$ for the simplest imaginable system: hard-disks. This result was obtained from a molecular dynamical simulation of $2500$ hard-disks. Once the simulation was thermalized, we ran the simulation 2000 times for different velocities and positions. The transversal current density correlation function, Eq. (2) was averaged over these 2000 simulation samples in order to reduce the noise. We used the event driven molecular dynamics simulation called DynamO [55]. It is important to remark that the presented results in Figure 1 were obtained for a packing fraction $\phi=0.68$, where the system is in a very dense fluid phase, close to the freezing point which is known to be at $\phi_{C}\approx 0.72$. [56, 57]

In the upper panel of Fig. 1, we show the resulting contour plot of $S(k,\omega)$. We can see that for small wavenumbers, shear waves do not propagate as expected for the fluid phase. However, Fig. 1 reveals a threshold $k_{c}$. Whenever $k>k_{c}$, shear waves indeed propagate in the system. In the lower panel of Fig. 1, we show the transversal part of the dynamical structure factor vs $\omega$ for different wavenumbers $k$ given in terms of the lowest wavenumber $k_{min}=4\sqrt{\eta\pi/N}$. Notice how as $k$ increases, the peaks in $S(k,\omega)$ shift to larger values of $\omega$. Furthermore, there is a gap between the peaks for $k\leq 3k_{min}$ and $k=4k_{min}$. From the upper panel in Fig. 1 we can see that $\omega(k)\approx\sqrt{k^{2}-k_{c}^{2}}$ in agreement with a recent theoretical solid-state approach to liquids [53, 58, 59].

Fig. 1 shows another viewpoint to look at viscoelasticity, but here the change from a fluid-like to a solid-like behavior is revealed by the presence of a dynamical gap [52, 53, 60]. Transversal wave propagation is only possible for modes with $k>k_{c}$. For $k<k_{c}$, in Fig. 1 we observe that $S(\bm{k},\omega)\approx\delta(\omega)$, where $\delta(\omega)$ is the Dirac delta function. As for $k<k_{c}$ we have $\omega=0$, we can consider these states in terms of rigidity as floppy, i.e., the system is flexible.

In general we can estimate a relationship between $k_{c}$ and the number of floppy modes as follows. Since the fluid is isotropic, the number floppy modes in three dimensions is,

[TABLE]

The fraction of floppy modes ($f$) with respect to the total number of modes is then,

[TABLE]

The normalization factor $k_{D}\;(\gg k_{c})$ is the Debye wavevector [60]. We thus arrive to the conclusion that floppy modes are related with a dynamical gap. Moreover, as $k_{c}=1/c\tau(T)$, where $c$ is the transverse sound speed and $\tau(T)$ is the average time at temperature $T$ it takes for a molecule to diffuse a distance equal to the inter-atomic separation [60], we can further relate floppy modes with this characteristic time,

[TABLE]

where $\omega_{D}=ck_{D}$. Although floppy modes in principle reduce the internal energy [61, 60], this will not happen in all cases, as entropy has two sources, vibrational and configurational [34].

As a matter of fact, floppy regions favor the maximization of vibrational entropy [34] and thus under certain conditions domains of floppy regions appear [35]. This in turn has huge consequences for relaxation [62, 63, 64, 65] and it becomes difficult to characterize rigidity using a mean-field approach above glass transition. Nevertheless, following the spirit of a mean field , we can define a mean coordination number $<r>$ of an effective topological lattice [20]. The fraction of floppy modes is $f=(3N-c)/3N$, where $c$ is the number of constrains. When angular and radial forces are present, this results in $f=2-5<r>/6$, while $f=1-<r>/6$ for radial forces. By using Eq. (5) we arrive to a possible and alternative definition for a “dynamical” mean coordination number in the melt when angular forces are present,

[TABLE]

and for pure radial forces,

[TABLE]

We remark that here $k_{c}=0$ implies $<r>=2.4$ whenever angular forces are present. In a similar way, $k_{c}=0$ implies $<r>=6$ for pure radial forces. These are the magical coordinations for rigidity transitions [19] and thus contain and highlight what we expect for a transition from a liquid to a solid. When there is a hierarchy of forces, these coordination numbers are not intended to necessarily caracterize glasses below $T_{g}$ as the solidified network can be already classified as floppy, isostatic or rigid. This task requires a more involved treatment, yet the present ideas suggest a path to be followed.

Let us discuss these dynamical results in the context of the usual invoked arguments relating relaxation time ($\tau$) and Newtonian viscosity ($\eta_{0}$) above $T_{g}$. This characteristic time at glass transition is estimated by setting $\eta_{0}\approx 10^{12}\ Pa\ s$ in the Maxwell relationship $\tau\approx\eta_{0}/G^{\prime}(\infty)$. This comes from the simplest model of viscoelasticity: a spring with a dashpot connected in series. However, the Maxwell model automatically implies exponential stress relaxation [66]. Glasses and glass-forming melts are known to have non-exponential relaxation [66], as for example, streteched exponential relaxation $\sigma(t)=\sigma_{0}\exp[-t/\tau]^{\beta}$ where $\beta$ depends upon the range of the interaction [6], as happens for relaxation in other topologically connected lattices [67, 68]. To be consistent, is paramount to search beyond the Maxwell picture. The task can be performed by using an extensive collection of models [41, 66]. Several paths are envisioned which include the use of fractional derivatives and generalized Maxwell-Voigt-Zener models with many spring dash-pots circuits to accurately reproduce all frequency decades [69]. This is in agreement with the use of Prony series to represent many relaxations times in order to obtain an accurate SER [70]. However, even for organic glasses is difficult to obtain models able to reproduce all kind of possible protocols for elasticity measurements [69]. Moreover, for inorganic glasses the relationship between rigidity and elasticity protocols is still a work in progress. For chalcogenide glasses, recent works add to viscoelasticity a plastic response [50] or a delayed elasticity [48] to account for the results on specific protocols. Any advance in this area is essential, as elastic stresses are related to thermodynamic driving forces for crystallization [71]. In this regard, Grassia et al. have made significant progress [72, 73, 74, 75, 76]. By linking viscoelasticity and the phenomenological KAHR model for structural relaxation, developed by Kovacs, Aklonis, Hutchinson, and Ramos [77, 78], they were able to characterize amorphous polymers and, in particular, predict the isobaric and isothermal glass transition for polystyrene [79].

Finally, we conclude by observing that for the system presented in Fig. 1, the dynamical gap goes to zero ($k_{c}\rightarrow 0$) as the hexatic to solid second order phase transition is approached. It remains to determine how the transition to rigidity occurs in glass forming melts, for example, by considering polydisperse disks. Also, we need to perform simulations on realistic Hamiltonians with angular dependent potentials. For organic glasses, it is known that such contributions increase relaxation times by steric shielding [43]. Cuts of the polymer chains and therefore, chain length, is an important parameter for relaxation in organic glasses [41, 42], yet is a factor that still needs to be addressed in time-dependent constraint theory for inorganic glasses.

Acknowledgments

This work was partially supported by DGAPA-UNAM project IN102717. J.Q.T.M. acknowledges a doctoral fellowship from CONACyT.