Local versus global stretched mechanical response in a supercooled liquid near the glass transition

TL;DR

This study uses molecular dynamics simulations to analyze local and global mechanical responses in a supercooled liquid near the glass transition, revealing consistent Cole-Davidson spectra and size-dependent relaxation behaviors.

Contribution

It demonstrates that local and global spectra follow the Cole-Davidson formula with size-independent stretching, and explains size effects via the elastic shoving model.

Findings

Spectra fit Cole-Davidson formula across scales

Stretching exponent remains constant regardless of region size

Size dependence of relaxation explained by elastic shoving model

Abstract

Amorphous materials have a rich relaxation spectrum, which is usually described in terms of a hierarchy of relaxation mechanisms. In this work, we investigate the local dynamic modulus spectra in a model glass just above the glass transition temperature by performing a mechanical spectroscopy analysis with molecular dynamics simulations. We find that the spectra, at the local as well as on the global scale, can be well described by the Cole-Davidson formula in the frequency range explored with simulations. Surprisingly, the Cole-Davidson stretching exponent does not change with the size of the local region that is probed. The local relaxation time displays a broad distribution, as expected based on dynamic heterogeneity concepts, but the stretching is obtained independently of this distribution. We find that the size dependence of the local relaxation time and moduli can be well explained by the elastic shoving model.

Equations

$$I$$

$$J$$

$$G^{\prime }(\omega )=\frac{I\omega }{N\pi };G^{\prime \prime }(\omega )=\frac{J\omega }{N\pi }$$

$$G^{\ast }(\omega )\equiv G^{\prime }(\omega )+i{G}^{\prime \prime }(\omega )$$

$$G^{\ast }(\omega )=G_{\infty }-G_{\infty }(1+i\omega \tau {)}^{-\beta },$$

$$\frac{G^{\prime \prime }(\omega )}{G^{\prime }(\omega )}=\frac{\sin (\beta \theta )}{[1+{\omega }^2{\tau }^2{]}^{\frac{\beta }{2}}-\cos (\beta \theta )},$$

$$\frac{G^{\prime \prime }(\omega )}{G^{\prime }(\omega )}\approx \sin (\frac{\pi }{2}\beta ){\omega }^{-\beta }{\tau }^{-\beta }.$$

$$G_{\infty }(r_c)=G_{\infty }(\text{bulk})(1+\frac{A}{r_c}+\frac{B}{r_c^3})$$

$$\tau (r_c)=e^{\frac{G_{\infty }(r_c)}{k_{B}T\rho }}$$

$$⟨\ln (\tau (r_c))⟩=\frac{⟨G_{\infty }(r_c)⟩}{k_{B}T\rho }\equiv \mu$$

$$⟨\ln (\tau (r_c){)}^2⟩-⟨\ln (\tau (r_c)){⟩}^2=\frac{S{D}^2(G_{\infty }(r_c))}{k_{B}T\rho }\equiv {\sigma }^2$$

$$\ln ⟨\tau (r_c)⟩=\mu +{\sigma }^2/2$$

$$SD(\tau (r_c))=\sqrt{e^{{\sigma }^2}-1}⟨\tau (r_c)⟩$$

$$⟨\ln \tau (r_c{)}^{-\beta }⟩=-\beta \mu$$

$$SD(\ln \tau (r_c{)}^{-\beta })=\beta \sigma$$

$$\ln ⟨\tau (r_c{)}^{-\beta }⟩=-\beta \mu +{\beta }^2{\sigma }^2/2$$

$$SD(\tau (r_c{)}^{-\beta })=\sqrt{e^{{\beta }^2{\sigma }^2}-1}⟨\tau (r_c{)}^{-\beta }⟩$$

Full text

Nonexponential or stretched exponential relaxation is ubiquitous in amorphous materials, and is recognized as one of the key features in supercooled liquid and glassy states phillips1996stretched ; RevModPhys.78.953 . It appears in many relaxation processes at equilibrium or out of equilibrium, such as aging, stress relaxation and dielectric or mechanical relaxation spectra gotze1992relaxation ; Wang2012a . However, the origin of the stretching is still controversial cavagna2009supercool . Two hypotheses are typically put forward to explain the stretching: one identifies the stretched relaxation as resulting from dynamic heterogeneity in different regions of space, the other assumes that the relaxation in amorphous material is uniform, with stretched relaxation being a local feature ediger1996supercooled ; richert2002heterogeneous .

These different views can to some extent be reconciled within the now widely accepted concept of dynamical heterogeneity, which has been confirmed both in experiment and molecular simulation berthier2011overview . The supercooled liquid, for example, can be separated into fast regions of high mobility and slow regions with lower mobility, with a “slow” or “fast” character that persists over times comparable to the total $\alpha$ relaxation time. Mathematically, stretched exponential relaxation can be described as a superposition of simple exponential relaxation processes edholm2000stretched . It is then a natural hypothesis to assume that the slow and fast regions associated with dynamical heterogeneity each have a simple exponential relaxation, and that the global stretching results from the different relaxation times associated with different regions, which may be broadly distributed. In fact, this natural assumption was recently formalized in a series of works by Masurel et al. masurel2015role ; PhysRevLett.118.047801 ; masurel2017dynamical , who developed a mesoscale model to describe the viscoelastic spectrum in a polymer model near the glass transition temperature. In their model, every local region is described as a single Maxwell Voigt element, with a single relaxation time assigned randomly from a broad (log normal) distribution. Based on the idea that dynamic and elastic heterogeneity are related, SchirmacherPhysRevLett.115.015901 also uses a local Maxwell model to describe the relaxation spectra within a mean field theory.

However, this assumption that the stretched exponential relaxation arises from simple exponential relaxation in local regions has not, to our knowledge, been proven in direct investigations. Only indirect consequences, as in the work of Masurel et al, have been explored. In this work, we question directly the validity of this assumption for mechanical properties, using the flexibility offered by molecular dynamics simulations. We build on previous explorations of static properties such as local elastic constants PhysRevE.87.042306 or thermomechanical Baoshuang2018 properties and develop a methodology that allows us to obtain the dynamic modulus spectrum in a supercooled liquid near the glass transition at different length scales. We find that different dynamical spectra can be well fitted by a Cole-Davidson expression, with a distribution of relaxation times that evolves with the measurement scale. However, surprisingly, the stretching exponent does not change with increasing the spatial scale, and is nontrivial at the smallest scale investigated. Furthermore, We find a strong correlation between the local modulus and relaxation time, which can be rationalized within an elastic shoving model RevModPhys.78.953 at the local scale, and the size dependence of the average relaxation time and shear modulus can be well explained by confinement effects, which reflect the nature of elastic interactions in supercooled liquids.

Our study of local viscoelastic properties is based on molecular dynamics (MD) simulations of a classical 80:20 binary Lennard-Jones (LJ) glass model PhysRevE.52.4134 using the LAMMPS package plimpton1995fast . The interaction potential was truncated and force shifted toxvaerd2011communication . LJ units based on the mass $m$, the interaction diameter $\sigma$ and energy $\epsilon$ of the large particles are used throughout, with a time unit $\sqrt{m\sigma^{2}/\epsilon}$. Ten independent simulation samples containing 80,000 atoms each were generated to improve statistics. Fig. 1(a) illustrates the potential energy per particle during quenching of the sample at various rates and constant volume in periodic boundary conditions. The quench rate used in our study is $1.67\times 10^{-6}$. The temperature was controlled with a Berendsen thermostat berendsen1984molecular , the time step was set to 0.001, and the number density $\rho$ was fixed at 1.2.

To obtain the dynamic shear modulus of the whole sample (bulk) and of local regions for different frequencies, we performed a numerical analog of a mechanical spectroscopy experiment PhysRevB.90.144201 at a fix temperature $T=0.435$. To determine the local properties, we used a modified version of the frozen matrix method, which was previously used to measure local moduli PhysRevE.87.042306 or local yield stress PhysRevLett.117.045501 at zero temperature. As shown in Figure 1(b), we first choose a spherical region of radius $r_{c}$, then shear the whole sample with sinusoidal strain $\gamma=\gamma_{A}\sin(\omega t)$ and only let the selected region relax under NVT conditions, while the outside region was affinely deformed according to the imposed sinusoidal shear strain. As a result, the local shear stress $\Sigma(t)$ acquires an oscillatory component at the same frequency, which is extracted from the noise using the numerical analog of a lock-in amplification. As shown in Fig. 1(c), the storage shear modulus $G^{\prime}(\omega)$ and loss shear modulus $G^{\prime\prime}(\omega)$ were derived from the local stress signal $\Sigma(t)$ using:

[TABLE]

[TABLE]

[TABLE]

where $N$ is the number of strain cycles and $\gamma_{A}$ is the amplitude of strain. Here we choose $N=5$, $\gamma_{A}=0.02$; the influence of the choice of $N$ and $\gamma_{A}$ is discussed in Supplementary Material (SM)supp . In order to sample the space, the center points of the free regions were selected from a regular grid. For $r_{c}=5$ , the sample was meshed as $8\times 8\times 8$, for $r_{c}=6,7,10$, the grid was $6\times 6\times 6$ and for $r_{c}=15$, the grid was $4\times 4\times 4$. The bulk dynamic modulus was obtained by applying an oscillatory strain to the whole simulation box.

For a series of different frequencies, a probability distribution of the storage and loss moduli can be obtained from the statistics over different zones. As shown in Fig. 2(a) and (b), this probability distribution of the local moduli shows a distinct frequency dependence. As the frequency of the loading increases, the most probable value of the local storage modulus shifts to higher values and the one of the local loss modulus shifts to lower values. This trend is general, independent of the size of the local region (for different sizes see Fig. S2 for $r_{c}=5$, Fig. S3 for $r_{c}=6$, Fig. S4 for $r_{c}=10$ in the SM).

In Fig. 3, the average values of the local storage and loss modulus (obtained from the probability distribution function) are now plotted as a function of frequency and compared with the bulk values obtained from dynamical mechanical analysis on the whole sample. Both storage and loss moduli are notably influenced by size. However, all dynamic moduli frequency spectra, regardless of size, can be well fitted by a Cole-Davidson form davidson1951dielectric ,

[TABLE]

where $i=\sqrt{-1}$, $G_{\infty}$ is the high frequency shear modulus, $\tau$ is the relaxation time, and $\beta$ is Cole-Davidson stretching exponent. Since our simulations are performed in the supercooled liquid state, where we expect $G^{\prime}(\omega)|_{\omega\to 0}=0$, the Cole-Davidson formula contains only three independent parameters. The formula reduces to the usual Debye model with a single relaxation time for $\beta=1$. We will see below that this value of $\beta$ is clearly outside the uncertainty range on the fit parameters.

Figure 4 reports the parameters obtained from fitting the frequency dependent mechanical response, averaged over all sampled regions as shown in Fig. 3, as a function of the size of the region. Surprisingly, the stretching exponent $\beta$ does not change with size (see Fig. 4(a)). However, since the response shown in Fig. 3 is averaged over many different regions, this feature can still be explained by two different hypotheses: either the average value is the superposition of individual relaxations, and the each individual region still follows a Debye relaxation with $\beta=1$ in the Cole-Davidson formula, or the stretching is a feature of every individual region.

In order to distinguish between these two possibilities, it would in principle be adequate to perform an individual fit of eq. (5) to the frequency dependence of the mechanical response for each region. Unfortunately, such a procedure turns out to be difficult in view of the relatively large statistical uncertainty of individual spectra. Instead, we choose to investigate the consistency of the above hypotheses with the observed statistical properties of the locally measured $G^{\prime}(\omega)$ and $G^{\prime\prime}(\omega)$, as characterized by the probability distribution functions shown in Fig. 2.

To this end, we calculate for every zone and loading frequency the ratio of loss and storage modulus $G^{\prime\prime}(\omega)/G^{\prime}(\omega)$. Within a Cole-Davidson model, this ratio reads

[TABLE]

where $\theta\equiv\arctan(\omega\tau)$. Considering that the loading frequency is generally such that $\omega\tau\gg 1$, this formula can be simplified as

[TABLE]

Figure 2(c) shows the probability distribution of this ratio for different loading frequencies. As expected from eq. (7), the probability distribution function is sensitive to loading frequency, with a most probable value that decreases with increasing frequency. The width of the distribution also decreases with increasing $\omega$. However, if one now rescales the data to obtain the probability distribution of $G^{\prime\prime}(\omega)/G^{\prime}(\omega)\omega^{\beta}$, which according to the Cole Davidson model is $\tau^{-\beta}$, a very good collapse of the different distributions, independent of frequency, is obtained, as shown in Fig. 2(d). Here the value of $\beta=0.30$ was set to the one obtained from fitting the average response (see Fig. 4)(a)), which shows that the average response is relevant to describe the statistical properties of individual zones. This proves that every individual region actually follows the same stretched relaxation process as the bulk sample. Note that the data in Fig. 2 corresponds to a particular size of the local region, namely $r_{c}=7$. However, different sizes lead to similar conclusionssupp .

Following this analysis, one must conclude that the collapsed distributions shown in Fig. 2(d) represent the probability distribution of $\tau^{-\beta}$ (up to a factor $\sin(\frac{\pi}{2}\beta)$). From the present data, it follows that the distribution of $\tau^{-\beta}$ is relatively narrow and it would slightly increase with $r_{c}$ within our investigation regime (shown in Fig. S5 in SM), in contrast with the common assumption of a broadly distributed relaxation time, and the width of the local relaxation time would decrease with $r_{c}$supp . As expected from the dynamic heterogeneity picture, relaxation is heterogeneous, but the heterogeneity does not explain the stretching of the relaxation, which is present at the local scale, and rather homogeneous.

Another interesting conclusion from the analysis shown in figures 3 and 4 concerns the size dependence of the average relaxation times and high frequency moduli. In contrast to the result for the stretching exponent, the size effect is very pronounced for these quantities. We explain this dependence by assuming that the relaxation time $\tau$ is related to a local free energy barrier $\Delta F$ through an Arrhenius law, $\tau\sim e^{\frac{\Delta F}{k_{b}T}}$. Assuming further, following the general ideas of elastic “shoving” models RevModPhys.78.953 , that the free energy barrier is mainly due to elastic energy, then the high frequency modulus $G_{\infty}$ is directly correlated with the free energy barrier, as $G_{\infty}\sim\Delta F$. The consistency of these assumptions can be directly checked using a parametric plot of $\ln\tau$ versus $G_{\infty}$ with the size of the region as a parameter. Such a plot is shown in figure 4(d) and reveals an excellent correlation between these quantities.

In order to understand the origin of the size effect on the energy barrier, we may assume that this energy barrier is associated with shear transformations taking place within the zone, and compare the situation in which the outside region is affinely deformed with the one in which this region is allowed to relax. Two effects will contribute to increasing the elastic energy in the frozen matrix configuration. Firstly, the local shear modulus in the vicinity of the frozen boundary will have a smaller nonaffine contribution, as nonaffine relaxation is partially prevented by the boundary. It is well known that the nonaffine contributions decrease the shear modulus compared to the purely affine (Born) value PhysRevLett.93.175501 , so that a shell in the vicinity of the boundary will effectively have an increased modulus compared to the bulk. Such a surface contribution is expected to contribute as $1/r_{c}$ to the shear modulus of the zone. In addition, since the shear transformations in supercooled liquid shows a Eshelby inclusion-like pattern illing2016strain ; PhysRevLett.113.245702 ,the energy barrier of shear transformations can be represented by the Eshelby external field energyeshelby1957determination ; picard2004elastic . For the local region, the displacement field associated with a shear transformation has to vanish at the boundary $r_{c}$. This effectively increases the average shear strain, and an order of magnitude estimate leads to an increase in the elastic energy scaling as $1/r_{c}^{3}$ due to this effect li2007eshelby . As a result, we propose to fit the size dependence of the shear modulus (or equivalently of the energy barrier) using the form:

[TABLE]

The result of this fit yields $A\simeq 1.0$, $B\simeq 5.4$. In the absence of any other relevant length scale, these values are consistent with the particle size and volume, respectively. Together with the correlation between relaxation time and shear modulus shown in Fig. 4(d), this fits provide evidence for the elastic nature of the local energy barriers.

We have investigated the scale and frequency dependence of the viscoelastic moduli of a generic glass forming system near the glass transition. In the system we investigated, the stretching of the relaxation cannot be simply assigned to the superposition of dynamic heterogeneities, but already exists at a very local scale. The smallest size we investigated, $r_{c}=5$, contains a few hundreds atoms. It corresponds to the typical size below which the mechanical response of a local region no longer follows Hooke’s law PhysRevE.80.026112 , and can be considered as the possible starting point of a coarse graining approach. However, the complexity of this “elementary volume” is already such, that a naive coarse graining starting from a simple Maxwell description at the local scale is not possible. In fact, this result is consistent with established results concerning the potential energy landscape of systems with a small number of particles doliwa2003finite ; rehwald2012coupled , which already for a few tenths of particles has a complexity comparable to that of larger systems. It can also be understood by considering the actual length scales involved in the problem. Recent theories of glasses introduce the dynamical length $\xi_{d}$ and the static length $\xi_{s}$. $\xi_{d}$ quantifies the fact that, in a bulk system, two regions at points $r$ and $r^{\prime}$ have different mobilities at times $t$ and $t^{\prime}$. This length grows rather rapidly near the glass transition, and is comparable to $r_{c}=5$ for our system, according to earlier work karmakar2014growing . In contrast, $\xi_{s}$ characterizes the amorphous order, grows much more modestly in the temperature range we are studying karmakar2014growing , and is less than 3 particle diameters. As a result, the smallest value of $r_{c}\approx 5$ at which we can define a shear modulus already contains several amorphous regions in the spirit of random first order theory Wolynesbook . In view of this ordering of the relevant length scales, the elastic barrier defined on the scale $r_{c}$ is not able to probe directly the expected growth of $\xi_{d}$. One may however speculate, as pointed out by Bouchaud and Biroli in chapter 2 of Wolynesbook , that the relevant volume in the elastic barrier involves $\xi_{d}^{3}$ rather than the atomic volume.

While the complexity that determines the stretching exponent is essentially insensitive to the scale, a nontrivial scale dependence emerges due to the elastic nature of energy barriers that govern relaxation in supercooled liquids PhysRevLett.113.245702 ; PhysRevLett.78.2020 . As discussed above, this dominance of the elastic aspect appears as we are operating on a scale larger than the static correlation length PhysRevLett.119.195501 . A different regime may be observed for $r_{c}/\xi_{s}<1$, a regime that we are not able to probe. Understanding the global relaxation on the basis of a coarse graining approach between elastically interacting elements seems therefore a promising approach. It would, however, require to model each element by a complex behavior, an approach which, to our knowledge has not been attempted until now, and may be challenging within the framework of classical finite element codes.

This work is supported by (BSS and PGF) the NSF of China (Grant Nos.51601009,Nos.51571011), the MOST 973 Program (No.2015CB856800) and the NSAF joint program (No.U1530401). BSS and PFG acknowledge the computational support from the Beijing Computational Science Research Center (CSRC). J-L. Barrat is supported by Institut Universitaire de France.

I Supplementary Material

S1 Robustness of the calculation of local dynamic modulus

We have investigated the influence of the number and amplitudes of strain cycles. In particular, we checked that for the selected amplitude of the strain the response is still linear (see Fig. S1). We also checked that the mean value of the local dynamic modulus was not sensitive to the number of cycles. The variance in the probability distribution decreases with the number of cycles, but beyond 5 cycles becomes dominated by the intrinsic disorder rather than by the noise on the measurement.

S2 local dynamic modulus for regions of size $r_{c}=5,7,\mathrm{and}\,10$

S3 Local relaxation time distributions

First, let us recall that we find that the relaxation time is related to the local shear modulus through:

[TABLE]

If we assume (consistent with many previous findings) that the modulus follows a Gaussian, it follows that the distribution of $\tau(r_{c})$ is a log-normal distribution with the properties:

[TABLE]

[TABLE]

where $SD(x)\equiv\sqrt{\langle x^{2}\rangle-\langle x\rangle^{2}}$ indicates the standard deviation. From the properties of the log normal distribution it follows that

[TABLE]

[TABLE]

The analysis can be extended to the quantity whose distribution is shown in fig 2(d), which is actually $\tau(r_{c})^{-\beta}$. This quantity should also follow a log-normal distribution, shown in Fig. S5, with

[TABLE]

[TABLE]

Again from the properties of the distribution,

[TABLE]

[TABLE]

As a result, the standard deviation depends on $r_{c}$ through the average value $\langle\tau(r_{c})^{-\beta}\rangle$, which increases when $r_{c}$ increases, and the prefactor $\sqrt{e^{\beta^{2}\sigma^{2}}-1}$, which is expected to decrease when $r_{c}$ increases. In the range of values of $r_{c}$ we have investigated, the net effect is an increase of the standard deviation with $r_{c}$, as illustrated in Fig. S6. At larger $r_{c}$, we expect the variation of $\langle\tau(r_{c})^{-\beta}\rangle$ to saturate and the standard deviation to decrease again with $r_{c}$.