A.
Madanchi†
^{a},
Ji Woong
Yu†
^{b},
Won Bo
Lee
*^{b},
M. R.
Rahimi Tabar
*^{cd} and
S. H. E.
Rahbari
*^{be}
^{a}Department of Physics, McGill University, H3A2T8, Montreal, Canada
^{b}School of Chemical and Biological Engineering, Institute of Chemical Processes, Seoul National University, Seoul 08826, Korea
^{c}Department of Physics, Sharif University of Technology, Tehran, 11155-9161, Iran
^{d}Institute of Physics and ForWind, Carl von Ossietzky University, 26111, Oldenburg, Germany
^{e}School of Physics, Korea Institute for Advanced Study, Seoul 02455, Korea

Received
16th November 2020
, Accepted 19th March 2021

First published on 22nd March 2021

Owing to the local/heterogeneous structures in supercooled liquids, after several decades of research, it is now clear that supercooled liquids are structurally different from their conventional liquid counterparts. Accordingly, an approach based on a local probe should provide a better understanding about the local mechanical properties as well as heterogeneous structures. Recently, the superiority of active microrheology over global rheology has been demonstrated [Yu et al., Sci. Adv., 2020, 6, 8766]. Here, we elaborate this new avenue of research and provide more evidence for such superiority. We report on the results of an extensive molecular dynamics simulation of active microrheology of a model glass. We identify several time scales in time series of friction, and detect a transition in dynamical behavior of friction. We discuss the possible relation to structural heterogeneities—a subject of considerable interest in glass physics.

In a supercooled liquid, particles transiently trap in cages formed by their neighbors. The dynamics of a supercooled liquid is thus described by two main relaxation processes, which is mathematically described by a two-step relaxation function C(t), (1) a fast exponential relaxation, known as β process with C(t) ∼ e^{−(t/τβ)}, and (2) a slow stretched-exponential relaxation, known as α process with C(t) ∼ e^{−(t/τα)b} where b < 1. In the β process, a particle rattles within the cage created by its neighbors. This back and forth movement terminates when the particle eventually escapes the cage, and this results in the final α relaxation. Below the melting temperature, cages become permanent and the α process, within accessible laboratory time scales, becomes unattainable. This picture may suggest a homogeneous relaxation dynamics in a supercooled liquid. In contrast, it was revealed that in a supercooled liquid, particles displaced, let's say longer than the average particle displacement, form clusters; This is known as DH. As a result of DH, the dynamics in a supercooled liquid is heterogeneous.^{1,3} This is in contrast to liquids that are widely considered homogeneous. Even though its name suggests a dynamic scenario, it is not clear whether the origin of DH is static or dynamic. Unveiling the origin of DH is important, because it will help to understand whether glass transition is a phase transition (the static origin) or a dynamic phenomenon (the dynamic origin). There are conflicting reports regarding the origin of DHs: some researchers advocate a static origin,^{4–7} while others suggest a dynamic origin.^{8} For example, Kawasaki and Tanaka identified clusters of medium-range crystalline order (MCRO), and demonstrated that MCROs are slower than the rest of the particles. The size of these clusters grows upon approaching the glass transition temperature. This suggests a possible link between MCRO and DH, and the authors concluded that MCRO may be the origin of DHs. Tanaka and co-workers pioneered this line of research and identified many other static structures;^{2} These results demonstrate the existence of a hierarchical order in supercooled liquids. Most of these structures provide length scales that grow upon cooling. However, these length scales never become macroscopic. This is perhaps because of the the dramatic slowing down of the dynamics near the glass transition temperature.

Interestingly, the mechanical properties of a glass are very different from those of liquids, and share many similarities with conventional solids; though glasses have some special features such as the existence of excess vibrations, known as the boson peak. These similarities can be summarized into one main feature: both glasses and solids are rigid. These mechanical properties are investigated via global deformations, or shearing. However, a global deformation would only render an average response, and the local properties due to typical structural heterogeneities of supercooled liquids cannot be explored by a global deformation. A remedy for this problem is an approach based on local deformation: active microrheology. In this technique, local mechanical and transport properties can be obtained via the investigation of dynamics of a single-out probe particle driven into host medium.^{9} This technique, active microrheology, has been widely used for investigations on biological materials, complex fluids,^{10} and hard spheres as model systems for colloids.^{11–14}

Three of us have previously investigated active microrheology of model glasses in a recent paper.^{15} In this paper, a probe particle moves with a constant velocity, U, and the local response of the medium is obtained via friction force experienced by the probe particle. It was revealed that the time series of friction contains blueprints of the structural relaxations of the host medium. The study had three major outcomes:

First, the probe particle exhibited a second-order delocalization transition when the temperature was varied across the glass transition temperature. As a result, the friction, F, satisfied a scaling ansatz of the form:

(1) |

Second, auto-correlation of the friction showed a two-step relaxation, described by a compressed exponential C(t) ∼ e^{−(t/τα)b} where b > 1. Third, the study provided the first glance into the spatial distribution of micro-viscosity, η_{m}. The distribution was bell-shaped, and the skewness of the distribution function of η_{m} increased as the temperature was lowered. A possible connection to DH was discussed in ref. 15. We further elaborate this in the current manuscript and provide more clues about the dynamical behavior of the friction in this system.

We notice that the delocalization transition is related to a larger class of far-from equilibrium phase transitions of driven particles over random or periodic substrates; dubbed depinning transition. At low temperatures, where the thermal motion of the medium is seized, active microrheology resembles a probe particle driven over a random substrate. The depinning transition has far more ramifications in a much broader context of vortices in type-II superconductors, and skyrmions, and many other examples are given in a recent review at ref. 16.

In this work, we look at more clues in time series of friction. We will follow a usual practice in time series analysis to investigate fluctuations and the corresponding scaling behaviors across all characteristic time scales. We will ask whether time series of friction are stationary. By stationarity, we mean a time series in which two point correlation function decays as a function of the time lag, similar to the behavior of a correlation function in an equilibrium system. Another possibility is a reverse behavior in which the correlation increases as a function of time lag; this is similar to fractional Brownian Motion (fBM). This corresponds to a non-stationary time series, and results in the persistency of the underlying fluctuations. For this purpose, we will not directly use correlation functions owing to the inadequacy of the correlation functions for non-stationary processes. Instead, we will use a standard method, known as multifractal detrended fluctuation analysis,^{17} for the time series of friction. We will show a transition in time series from stationary to non-stationary in the dynamics of friction. This will help us to identify a length scale that grows by lowering the temperature.

(2) |

We used a state-of-the-art method, known as shear annealing, to stabilize the host medium.^{20,21} In this method, the system is periodically sheared, and below a threshold yielding amplitude, the system settles down into deeper basins in the energy landscape. The more the system is annealed, the more stable it becomes, and the fewer fluctuations happen. This step is crucial to reduce the fluctuations of the system when we perform active microrheology.

(3) |

For the given time series x(t_{i}) ≡ x(i) (here friction time series) for {I = 1,…,N}, we define the mean of x series as:

(4) |

_{i} = x_{i} − 〈x〉. | (5) |

(6) |

(7) |

C(s) ∼ s^{−γ}, | (8) |

2.3.1 Detrended and multifractal detrended fluctuation analyses.
It is known that the presence of trends in a time-series can result in fictitious correlations.^{17} Therefore, it is important to eliminate all possible trends in a time-series prior to any statistical computation. For this purpose, a technique, dubbed detrended fluctuation analysis (DFA), has been developed^{17} and has already been applied to various time series analyses.^{23–27} For a given time-series x(t_{i}) ≡ x_{i}, DFA is performed in these consecutive steps:

For simplicity from now on, we will refer to Y^{det.}_{s} without the superscript det. After the detrending step, we can proceed to the final step for the investigation of fluctuations. For this purpose, we will measure mean variance (fluctuation function) of data averaged over all segments.

The standard fluctuation function, F(s), is obtained by

The superiority of DFA becomes clear upon a re-inspection of eqn (11). Usually, a variance is defined as F_{s}^{2} = 〈(Y_{s}(i))^{2}〉 − 〈Y_{s}(i)〉^{2}, where the last term on the right-hand-side (rhs) of this equation can cause numerical instabilities during a scaling analysis. This hinders the detection of genuine scaling properties of the data. For clear reasons, in DFA that term is not required. To investigate the underlying scaling properties, we seek the following scaling relation

where the index h has the following relationship with exponent of the correlation function C(s) of eqn (8):

If h < 1, then γ > 0, and thus, the correlation function C(t) will decay as a function of the time lag, s. As a result, a times series with h < 1 is called stationary. However, if h > 1, then γ < 0 and C(s) ∼ s^{|γ|}: the correlation will grow as a function of time s, akin to fBM, and as a result a time series with a positive correlation is considered non-stationary:

For a stationary process, h < 1, the exponent h is identical to the well-known Hurst exponent H.^{28,29} For this reason, h is called generalized Hurst exponent. Depending on the value of H, the time series may have various dynamical states:

(i) Integration: compute integral of the original time series

(9) |

(ii) Division: the time series (with the length of N) is then divided into non-overlapping segments, where s is the box size. Because N/s has normally a non-zero remainder, especially when s is large, the remaining part of the time series is discarded. To prevent any bias as a result of division, the time-series is first divided from the beginning to end, and the remaining part is left, and second time the division is performed in the opposite direction. As a result 2N_{s} segments are produced. Now, one can perform the detrending step.

(iii) Detrending: for each segment s, a linear, or higher-order, regression is performed using least square algorithm. For a linear dependence, the resulting function is P_{s}(i) = a_{s} + b_{s}i, where a_{s} and b_{s} are constants obtained from the regression, and i is time. The detrending is performed by a simple subtraction from each data point within the interval:

Y^{det.}_{s}(i) = Y(i) − P_{s}(i). | (10) |

(iv) Scaling of fluctuations: variance of an interval of length s is defined as

(11) |

(12) |

F(s) ∼ s^{h}, | (13) |

γ = 2 − 2h. | (14) |

(15) |

(16) |

The implementation of DFA has been done by the MFDFA package in python provided by LRydin/MFDFA.

Fig. 1 Friction time series. The magnitude of the total force exerted by the medium to the probe particle, friction, is depicted versus the probe displacement, t × U, for different probe velocities. The medium is in the liquid state with temperature T = 0.45, which is above the delocalization transition temperature T_{mc} = 0.36875.^{15} The stochastic characteristics of the times series of friction depends on temperature and velocity. We will show that for U = 1 and 0.018, the corresponding generalized Hurst exponents are h < 1 and h > 1, respectively. However, for U = 0.223, some of simulated time series are stationary and some others are non-stationary. This represents a transition as a function of the control parameter U in the behavior of time series. |

In order to gain more insight about the dynamical nature of friction time series, similar to our previous work,^{15} we present auto-correlation of time series of friction, C(t) = 〈F(0)F(t)〉, in Fig. 2a and d, for T = 0.15 and 0.466, respectively. Each curve has a different color which corresponds to different probe velocity U whose values are given in the legend. One can see that C(t) faithfully represents the typical two-step structural relaxation of a supercooled liquid: (1) an exponential decay as a result of ballistic motion, followed by a plateau representing the cage effect akin to the β relaxation with time scale τ_{β}, and (2) final relaxation due to escape of a particle from the surrounding cage akin to the α relaxation with τ_{α}. These time scales can be obtained by a model-fitting function of the form:

C(t) = (1 − A)e^{−(t/τβ)} + Ae^{−(t/τα)b} | (17) |

We have previously shown that by rescaling C(t) versus x = t × U^{ν}, where ν = 1 − q = 0.95 being the force-thinning exponent, all curves collapse into a master curve given by a compressed exponential function e^{−(x/τα)b} with b > 1, resembling a time-force superposition similar to time-shear superposition principle predicted by mode-coupling theory^{31} and spin glass theory.^{32} These results established the fact that friction time series contains blueprint of structural relaxations in supercooled liquids.

An interesting point is the appearance of a second small plateau for low-temperature and glass phase and small enough probe velocities in panel-a, which might be related to the fBm nature of friction dynamics in those scales and propagation of waves in the system as a result of the probe motion. This secondary plateau does not exist in the liquid state in panel-d. Therefore, active microrheology reveals many more interesting time scales that are not visible for globally averaged quantities, such as the self part of the scattering function. These time scales cannot be obtained by a simple fitting procedure using eqn (17), and more elaborate techniques should be used. To distinguish these cross overs, we compute the magnitude of Fourier transform of time series, {F(t)}, depicted in Fig. 2 panel-b and -e. The existence of many cross overs at different time scales becomes much more clear in the frequency domain, f = 1/t. Below a frequency f_{×}, curves nearly level off at f → 0. For f > f_{×}, a scaling regime is observed. In panel-f in liquid phase, we observe a ripple mode on f_{x} line (zoom in required).^{33} This scaling region levels off at f_{α}, which is followed by another plateau, similar to the corresponding plateau in the time domain in C(t), then there is another scaling regime at the end of this region for f > f_{β}. We use these plots to identify cross over frequencies f_{×}, f_{α} and f_{β} for every curve individually. These cross overs are marked by dashed lines. The boundaries between these cross overs become blurred at large probe velocities for U > 2.33 × 10^{−1}.

We are now ready to perform DFA. This is done in Fig. 2 panel-c and -f. One can see that each curve crosses over to various scaling regimes. Generally, DFA provides more clean scaling behavior with respect to other methods. The corresponding cross overs can be found using both DFA or power spectrum analysis. Using the results of the Fourier transform analysis, we can determine cross overs in these curves as, s_{β} = 1/f_{β}, s_{α} = 1/f_{α} and s_{×} = 1/f_{×}. Using linear regressions for the ranges s_{β} < s < s_{α} and s_{α} < s < s_{x} in double logarithmic scale, we examine the following scaling behaviors

(18) |

More interesting behavior is found in panel-b for α region. Here, we discover that h_{α} crosses over from a region with h > 1 to h < 1. A transition boundary can be drawn at h = 1; this is marked by a dashed line. The area inside the cone represents the non-stationary regime. The non-stationary region expands as temperature is lowered. The area outside the cone represents the stationary region. To gain a better insight about the nature of this “stationary to non-stationary” transition, we look individually at each trajectory within our ensemble, and we calculate the ratio of cases where h < 1 to the total number of trajectories, or probability P(h_{α} < 1) in panel-c. It can be seen that deep inside the non-stationary region near the origin of the plot, almost all trajectories are non-stationary P(h_{α} < 1) ≃ 0; however, near the transition boundary almost half of trajectories are stationary. This may imply that the transition is second-order-type.

The transition boundary at h = 1 for a given T occurs at a critical probe velocity U = U_{c}. We define a critical length scale, l_{c}, given by l_{c} = (t_{x} − t_{α})U_{c}. We plot this length scale as a function of temperature in Fig. 4. One can see that l_{c} increases as T is lowered. The largest gap in l_{c} happens near the delocalization transition temperature at T_{mc} = 0.36875. In the liquids phase, l_{c} is almost zero. Similar to many other length scales that are proposed for supercooled liquids, this length scale remains microscopic across all temperatures. However, to gain a better understanding of this length scale, a finite size scaling is required to see how l_{c} grows with the system size. This is out of reach for the present work, but we will address it in near future.

Previous studies on the rheology of supercooled liquids by Berthier-Barrat,^{35} and Varnik^{36} demonstrated that, even though a supercooled liquid falls out of equilibrium near the melting temperature, however, the application of a shear flow into the system settles down the system into a dynamical stationary state akin to an equilibrium system. Here, active microrheology demonstrates that there exists a stationary state only above a critical threshold, U > U_{c}. Note that, we previously have shown that the probe velocity U is akin to a global flow rate via =U/D, where D is the probe diameter. It is possible that the onion structure in T–U plane is smeared out in macrorheology. This demonstrates once again the superiority of active microrheology with respect to macrorheology.

Our analysis allowed us to extract a length scale, l_{c}, from our data, which is defined based on the cross over probe velocity U_{c}, and the corresponding time scale, t_{x} − t_{α}. We found that this length scale is almost zero in the liquid phase. Near T_{mc}, the critical delocalization temperature, l_{c} has a sharp increase, and then grows monotonically as the temperature is lowered. Similar growth of a length scale has been also observed in dynamic heterogeneities and hierarchal structures.^{2} Our next goal is to find a relationship between this length scale and other known length scales in the literature. It will also be very interesting to see whether l_{c} will grow with the system size. This requires a systematic finite size scaling. We will address these questions in near future.

Microrheology is known to exhibit pronounced finite-size effects due to the finite size of system and the size of tracer. In its passive version, when the external drive on the tracer is switched off, results of simulations showed satisfactory quantitative agreement with continuum theories when the tracer size was about five times larger than that of the embedding medium particles.^{37,38} A recent study on active microrheology found good agreement between the hydrodynamic theory and simulations for the dependence of viscosity as a function of system size only at intermediate system sizes,^{39} and the the subject requires further investigation. Our scaling approach offers a sound framework for this problem, because the system size, L, can be cast as a relevant perturbation to the scaling ansatz proposed in our previous work.^{15} Because the delocalization transition in this system suggests a critical point at which the correlation length diverges, b = ξ ∝ δT^{−ν}, and the friction scales as a function of distance from the critical temperature, F ∝ δT^{y}, near the critical temperature δT → 0 a standard critical scaling can be written as

F = b^{−y/ν}0(Ub^{z},bL^{−1}), | (19) |

F = L^{−y/ν}_{1}(UL^{z}), | (20) |

- C. A. Angell, Science, 1995, 267, 1924–1935 CrossRef CAS PubMed.
- H. Tanaka, H. Tong, R. Shi and J. Russo, Nat. Rev. Phys., 2019, 1, 333–348 CrossRef.
- L. Janssen, Front. Phys., 2018, 6, 97 CrossRef.
- A. Widmer-Cooper and P. Harrowell, Phys. Rev. Lett., 2006, 96, 185701 CrossRef PubMed.
- H. Shintani and H. Tanaka, Nat. Phys., 2006, 2, 200 Search PubMed.
- T. Kawasaki, T. Araki and H. Tanaka, Phys. Rev. Lett., 2007, 99, 215701 CrossRef PubMed.
- T. Kawasaki and H. Tanaka, J. Phys.: Condens. Matter, 2010, 22, 232102 CrossRef PubMed.
- L. O. Hedges, R. L. Jack, J. P. Garrahan and D. Chandler, Science, 2009, 323, 1309–1313 CrossRef CAS PubMed.
- A. M. Puertas and T. Voigtmann, J. Phys.: Condens. Matter, 2014, 26, 243101 CrossRef CAS PubMed.
- T. A. Waigh, Rep. Prog. Phys., 2005, 68, 685 CrossRef.
- P. Habdas, D. Schaar, A. C. Levitt and E. R. Weeks, Europhys. Lett., 2004, 67, 477 CrossRef CAS.
- I. C. Carpen and J. F. Brady, J. Rheol., 2005, 49, 1483–1502 CrossRef CAS.
- A. S. Khair and J. F. Brady, J. Fluid Mech., 2006, 557, 73 CrossRef.
- L. G. Wilson, A. W. Harrison, A. B. Schofield, J. Arlt and W. C. K. Poon, J. Phys. Chem. B, 2009, 113, 3806–3812 CrossRef CAS PubMed.
- J. W. Yu, S. H. E. Rahbari, T. Kawasaki, H. Park and W. B. Lee, Sci. Adv., 2020, 6, eaba8766 CrossRef PubMed.
- C. Reichhardt and C. J. O. Reichhardt, Rep. Prog. Phys., 2016, 80, 026501 CrossRef PubMed.
- S. V. Buldyrev, A. L. Goldberger, S. Havlin, R. N. Mantegna, M. E. Matsa, C. K. Peng, M. Simons and H. E. Stanley, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 1995, 51, 5084 CrossRef CAS PubMed.
- W. Kob and H. C. Andersen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1993, 48, 4364 CrossRef CAS PubMed.
- W. Kob and H. C. Andersen, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 4626 CrossRef CAS PubMed.
- P. Das, A. D. S. Parmar and S. Sastry, 2018, arXiv preprint arXiv:1805.12476.
- D. J. Lacks and M. J. Osborne, Phys. Rev. Lett., 2004, 93, 255501 CrossRef PubMed.
- J. Feder, Fractals, Springer Science & Business Media, 2013 Search PubMed.
- P. C. Ivanov, M. G. Rosenblum, C. K. Peng, J. E. Mietus, S. Havlin, H. E. Stanley and A. L. Goldberger, Phys. A, 1998, 249, 587–593 CrossRef PubMed.
- A. Bunde, S. Havlin, J. W. Kantelhardt, T. Penzel, J. H. Peter and K. Voigt, Phys. Rev. Lett., 2000, 85, 3736 CrossRef CAS PubMed.
- M. S. Movahed, F. Ghasemi, S. Rahvar and M. R. R. Tabar, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 021103 CrossRef PubMed.
- A. Madanchi, M. Absalan, G. Lohmann, M. Anvari and M. R. R. Tabar, Sol. Energy, 2017, 144, 1–9 CrossRef.
- A. Madanchi, F. Taghavi-Shahri, S. M. Taghavi-Shahri and M. R. R. Tabar, Eur. Phys. J. B, 2020, 93, 1–8 CrossRef.
- P. C. Ivanov, L. A. N. Amaral, A. L. Goldberger, S. Havlin, M. G. Rosenblum, Z. R. Struzik and H. E. Stanley, Nature, 1999, 399, 461–465 CrossRef CAS PubMed.
- M. S. Movahed, G. R. Jafari, F. Ghasemi, S. Rahvar and M. R. R. Tabar, J. Stat. Mech. Theor. Exp., 2006, 2006, P02003 Search PubMed.
- S. Sengupta, S. Karmakar, C. Dasgupta and S. Sastry, J. Chem. Phys., 2013, 138, 12A548 CrossRef PubMed.
- W. Gotze and L. Sjogren, Rep. Prog. Phys., 1992, 55, 241 CrossRef.
- L. Berthier and G. Biroli, Glasses and Aging, A Statistical Mechanics Perspective on, 2009 Search PubMed.
- N. Abedpour, R. Asgari and M. R. R. Tabar, Phys. Rev. Lett., 2010, 104, 196804 CrossRef CAS PubMed.
- M. R. R. Tabar, M. Anvari, G. Lohmann, D. Heinemann, M. Wächter, P. Milan, E. Lorenz and J. Peinke, Eur. Phys. J.-Spec. Top., 2014, 223, 2637–2644 CrossRef.
- L. Berthier and J. L. Barrat, J. Chem. Phys., 2002, 116, 6228–6242 CrossRef CAS.
- F. Varnik, J. Chem. Phys., 2006, 125, 164514 CrossRef CAS PubMed.
- I. C. Yeh and G. Hummer, J. Phys. Chem. B, 2004, 108, 15873–15879 CrossRef CAS.
- R. O. Sokolovskii, M. Thachuk and G. N. Patey, J. Chem. Phys., 2006, 125, 204502 CrossRef CAS PubMed.
- F. Orts, G. Ortega, E. M. Garzón and A. M. Puertas, Comput. Phys. Commun., 2019, 236, 8–14 CrossRef CAS.

## Footnote |

† Madanchi and Yu have equal contribution. |

This journal is © The Royal Society of Chemistry 2021 |