Temperature dependence of the transition packing fraction of thermal jamming in a harmonic soft sphere system

We consider a soft sphere system in three dimensions where the particles repel each other with a force that is proportional to the overlap, i.e. the pair potential for the spheres i and j  is given by

$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:scaling} V(r_{ij})=\left\{\begin{array}{@{}c@{}} \epsilon\left(1-\frac{r_{ij}}{\sigma_{ij}}\right)^2/2, r_{ij}<\sigma_{ij}, \nonumber \\ 0, r\geqslant\sigma_{ij},\end{array}\right. \nonumber \end{align} \tag{ 1 }$

where r_ij is the distance between the spheres, $\sigma_{ij}$ is their mean diameter, and $\newcommand{\e}{{\rm e}} \epsilon$ is a positive constant that denotes the strength of the repulsion.

The temperature T is given in units of $\newcommand{\e}{{\rm e}} \epsilon/k_{\rm B}$ in the data from [11] that we analyze in section 3.2. Originally, the temperature T in [9, 10] is given in units of $\newcommand{\e}{{\rm e}} 2\epsilon/k_{\rm B}$, which we change to $\newcommand{\e}{{\rm e}} \epsilon/k_{\rm B}$ in our discussion in section 3.3.

In this article we reanalyze simulation data from [11] and from [9, 10]. In the following, we shortly comment on the methods employed in those works as well as on the quantities that are considered in the respective data sets. Further details are given in the original articles [9–11].

In [11] we employed molecular dynamics simulations of a monodisperse system at constant temperature T and constant pressure p . The timescale that describes rearrangements is given by the time $\tau$ where the mean square displacement reaches the squared diameter, i.e. $\left\langle r^2 (\tau)\right\rangle=\sigma^2$. For the mapping onto the hard sphere dynamics the dimensionless time $\tau^*=\tau\sqrt{p\sigma/m}$ is used.

In [9, 10] Berthier and Witten present molecular dynamics simulations at constant temperature and volume. A bidisperse system with a 50:50 mixtures with a ratio of diameters $\sigma_1/\sigma_2=1.4$ is used. Relaxation times $\tau_{\alpha}$ are measured as $\alpha$-relaxation time by considering the decay of the self-part of the intermediate scattering function and usually the dimensionless timescale $\newcommand{\e}{{\rm e}} \tau_\alpha^*=\sqrt{T\epsilon/m}\tau_\alpha/\sigma_2$ is employed. The data that we reanalyze in section 3.3 has been extracted from figure 2(a) of [9].

Note that different ways to measure the timescales are used in [11, 16] and [9, 10]. However, within each set of data we stick to one definition, namely the definition employed in the original publication. We assume that for the considered systems the scaling works for all definitions if it is based on the divergence that can be obtained by extrapolation.

In [16] we introduced a new approach to explore the energy landscape of a glassy soft sphere system. Motivated by the method employed to determine the athermal jamming transition [17, 18], we start with a random initial configuration and then minimize the energy without crossing energy barriers, e.g. by employing a conjugate gradient minimization or a steepest decent approach. Thermal fluctuations around the paths of these minimization methods are not considered. However, in order to denote the possibility that an energy barrier can be crossed, we introduce additional steps that occur with a small probability p  and that displace a particle in a random direction such that a barrier can be crossed. Note that we have checked that the details on how these additional steps are implemented do not matter [16] as long as they enable the rare random crossing of energy barriers. In our simulations with this new method we analyzed whether the ground state can be reached or not. In case the ground state cannot be reached within the time of the simulation the system is effectively non-ergodic. If the ground state is reached, the system is considered to be unjammed [17, 18]. Another approach to explore and characterize the energy landscape of glasses were employed in [19], where the inherent structures for equilibrated systems are determined and a sharp increase of the obtained jamming packing fraction is observed when the equilibrated systems enter the glass regime [19]. Note that our approach in [16] does not rely on determining inherent structures.

We have confirmed that within our approach the transition packing fraction as a function of the probability p  (see figure 1(a)) in the limit of $p\rightarrow 0$ does not approach the packing fraction $\phi_J\approx 0.64$ [17, 18] of athermal jamming that one obtains for p   =  0 [20]. Instead the glass transition line $\phi_g(\,p)$ for small but non-zero p  saturates at a lower value $\phi_{g,0}\approx 0.55$. The transition packing fraction $\phi_g(\,p)$ decreases with increasing p , i.e. if energy barriers are crossed more often, the relaxation of the system becomes more difficult. In [16] we showed that these properties can be explained by the number of particles that are affected by a random barrier crossing. Close to the packing fraction $\phi_g(\,p)$ a crossing of a barrier can affect almost all particles of the system due to a spatial percolation of the system that also has been observed in [21]. Therefore, above $\phi_g(\,p)$ almost all particles in the system have to restart their relaxation process after a barrier crossing event and as a consequence their relaxation is significantly delayed. Note that a similar mechanism of correlated relaxation processes is also observed in [22] and percolation transitions have been connected to the glass transition in other systems as well [23, 24].

Zoom In Zoom Out Reset image size

Note that in [16] we also showed that in case of non-instantaneous quenches or in smaller systems glass transition packing fractions above 0.55 are usually observed though the critical behavior remains unchanged. Therefore, in this article we will concentrate on how the glass transition is shifted due to changes of the temperature but not on the absolute values of the dynamical glass transition because the latter obviously are protocol dependent (see also [4, 6, 16]).

For the case that fluctuations within the valleys of the energy landscape are neglected, we determined how the glass transition packing fraction $\phi_g$ depends on the probability p  that determines how often energy barriers are crossed [16] (see also figure 1(a)). Obviously the probability p  is related to the temperature T of a soft sphere system. In principle, the probability can be approximated by Kramer's rate [25]. However this requires an extensive study of the energy barriers heights as well as the curvatures of the non-trivial energy landscape. Here we employ an alternative approach in order to relate the probability p  to the temperature T based on a comparison of properties of the pair distribution functions $g(r)$ close to the glass transition.

For the pair distribution functions $g(r)$ close to the glass transition we consider the peak heights g_m. The pair distribution functions $g(r)$ have been determined in [16] for various probabilities p  (see also figure 1(b)). As a result, the peak heights g_m(p ) are available as a function of p . By inverting the the glass transition line $\phi_g(\,p)$ shown in figure 1(a) we can express g_m(p ) as a function of the packing fraction $\phi_g$ (see figure 2). We find that $g_{m}(\phi_g)$ as a function of how much the glass transition packing fraction deviates from the $p\rightarrow 0$ result, i.e. as a function of $\phi_{g,0}-\phi_g(\,p)$, can be described by a power law $g_{m}\propto (\phi_{g,0}-\phi_g(\,p)){}^{-1/(2\beta)}$ with $\beta=0.202\pm 0.005$ (see figure 2). In the following, we will employ the approximation $\beta\approx 0.2$.

Zoom In Zoom Out Reset image size

According to [7], the peak heights g_m for soft harmonic spheres depend on the temperature as $g_{m}\propto T^{-1/2}$. As a consequence, we expect $\phi_{g,0}-\phi_g(\,p)\propto g_{m}^{-2\beta}\propto T^{\beta}$. Therefore, we are able to express the scaling behavior of the glass transition curve $\phi_g(\,p)$ not only in terms of the probability p , but also in terms of the temperature T, namely we expect $\phi_{g,0}-\phi_g(T)\propto T^{\beta}$ with $\beta\approx 0.2$. Note that $\phi_g(T)$ is smaller than $\phi_{g,0}$ for T  >  0.

While fluctuations around a valley of the energy landscape lead to an effective packing fraction that is smaller than the real packing fraction and as a consequence increases the distance from the glass transition for increasing temperature, the crossing of energy barriers can prevent the relaxation and thus leads to a decrease of the glass transition packing fraction as a function of temperature. If the effective packing fraction $\phi_{\rm ef\,\!f}(T)$ scales like $\phi-\phi_{\rm ef\,\!f}\propto T^{1/2}$ for harmonic interactions and $\phi_{g,0}-\phi_g(\,p)\propto T^{\beta}$ with $\beta\approx 0.2$ the distance to the glass line behaves like $\left(\phi_g(T)-\phi_{g,0}\right)-\left(\phi-\phi_{\rm ef\,\!f}(T)\right)\propto T^{0.5}-bT^{0.2}$ (corresponding to equation (1)) where b is a positive constant that weights the two scaling contributions. Note that $\phi_{\rm ef\,\!f}$ can also be determined by calculating effective diameters, e.g. by employing the Andersen–Weeks–Chandler method [12]. Furthermore, our choices to denote the fluctuations around the valley by an decreasing effective diameter instead of an increasing glass transition and to consider the barrier crossing by a deceasing glass transition density instead of an increasing effective packing fraction are arbitrary. The only physical relevant quantity is the difference of the effective packing fraction from the glass line at the respective temperature. This difference is determined such that the soft sphere system behaves similar as a hard sphere system or another soft sphere system with the same difference.

In this subsection we reanalyze the simulation data from [11]. Our main goal is not not cover a large region of packing fractions and especially we do not aim to get even close to the glass transition packing fraction. For such a task better simulation methods and data are available, especially since Monte Carlo simulations with particle swaps have been developed [6]. Our focus here is on the temperature dependence and specifically on the systematic deviations from a hard sphere-like dynamics that have been observed in [11] and that are obvious especially for larger temperatures, i.e. temperatures $T\geqslant 10^{-3}$. In order to quantify the dependence of these deviations on the temperature we employ fits to Vogel–Fulcher–Tammann functions. We want to emphasis that due to the limited range concerning the packing fraction the obtained values of the fit parameters might be somehow arbitrary. However, we are only interested on how these parameters change as a function of temperate. Over the large temperature range of the data from [11] these changes are systematic.

In [11] the dynamics of soft spheres with an diameter $\sigma$ were compared to the dynamics of hard spheres with a diameter $\sigma_{\rm ef\,\!f}$, where $\sigma_{\rm ef\,\!f}$ is the effective diameter as determined by using the Andersen–Weeks–Chandler method [12]. Originally, this method was developed to map the structure of soft sphere systems to the structure of hard spheres with the effective diameter. Concerning the dynamics, it turned out that for small overlaps the soft spheres dynamics can be well-mapped onto the hard sphere dynamics. In figure 3(a) the results for small temperatures leading to small overlaps are shown by black and blue circles and can in principle be nicely described by a single curve corresponding to the hard sphere behavior (not shown). However, significant deviations occur for larger overlaps, especially if $\sigma_{\rm ef\,\!f}/\sigma<0.9$ [11]. To be specific, the relaxation times of soft spheres at larger temperatures systematically deviate from the relaxation times of the corresponding hard spheres or data in the $T\rightarrow 0$ limit as can be seen in figure 3(a). For small effective packing fractions the relaxation times for larger temperatures (green and red data in figure 3(a)) lies below the curve expected from relaxation times that are obtained in the $T\rightarrow 0$ limit (black and blue data in figure 3(a)). Note that for larger packing fractions where no obvious trend is visible there might be larger uncertainties concerning the data points (for more details see discussion in [11]).

Zoom In Zoom Out Reset image size

While the effective diameter takes into account the fluctuations around the equilibrium-like average distance between soft spheres, the dynamical contribution of the crossing of energy barriers as studied in [16] are not considered. Note that the crossing of barriers slows down the relaxation process. In the following we show that the systematic deviations described in [11] can be explained in terms of temperature-dependence of the glass transition packing fraction that we have discussed in section 3.1.

The data of [11] is obtained for a monodisperse system at packing fractions where within the time of simulations no crystallization can be observed. As a consequence, the considered packing fractions cannot be close to the glass transition packing fraction. However, the glass transition packing fraction can be estimated, e.g. by fitting Vogel–Fulcher–Tammann functions

$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber \tau(\phi,T)(\,p\sigma/m)^{1/2}&=A(T)\cdot\nonumber \\ &\cdot\exp\left[-\frac{B(T)}{\left(\phi_{g}(T)-\phi_{\rm ef\,\!f}(T)\right)^{\delta}}\right] \label{eq:VFT} \nonumber \end{align} \tag{ 2 }$

to the relaxation time as function of the packing fraction for given temperatures leading to the glass transition packing fraction $\phi_g(T)$ at that temperature. Note that there are different ways how an extrapolated glass transition packing fraction can be determined. Concerning the Vogel–Fulcher–Tammann functions different exponents $\delta$ can be employed. While we cannot decide on which $\delta$ is the best choice, we show in the following, that the obtained glass transition packing fractions in all considered cases agree to the scaling proposed in section 3.1.

First we fit Vogel–Fulcher–Tammann functions as in equation (2) with $\delta=1$ to the data for temperatures within certain intervals (i.e. $\newcommand{\e}{{\rm e}} T/\epsilon\in [10.0^{m},10.0^{m+1}]$ for $m=-9,-8,$ $-7,-6,-5,-4,-3,-2$) in order to determine how $A(T)$, $B(T)$, and $\phi_{g}(T)$ depend on the temperature T. We find that within the considered temperature range $A(T)$ and $B(T)$ can be well fitted by linear functions of $\newcommand{\e}{{\rm e}} \log(T/\epsilon)$, i.e. $\newcommand{\e}{{\rm e}} A(T)\approx 0.16\log(T/\epsilon)+0.57$ and $\newcommand{\e}{{\rm e}} B(T)\,\approx$ $\newcommand{\e}{{\rm e}} -0.071\log(T/\epsilon)-0.039$. Furthermore, $\phi_{g}(T)$ can be fitted by a power law: $\phi_{g}(0)-\phi_{g}(T)\propto T^\beta$ with $\beta=0.177 \pm 0.030$ and $\phi_{g}(0)=0.72\pm 0.01$ (see figure 3(b)). Note that we fitted all data from within a whole temperature interval which might lead to mistakes. Therefore, in the following we extend our analysis by performing more extensive fits.

We now employ fits according to

$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \nonumber \tau(\phi,T)(\,p\sigma/m&)^{1/2}=\left(A_1\log(T/\epsilon)+A_2\right)\cdot\nonumber \\ &\cdot \exp\left[-\frac{\left(B_1\log(T/\epsilon)+B_2\right)}{\left(\phi_{g}(0)-cT^\beta-\phi_{\rm ef\,\!f}(T)\right)^{\delta}}\right] \label{eq:VFT-fit} \nonumber \end{align} \tag{ 3 }$

with fitting constants A₁, A₂, B₁, B₂, c, $\phi_{g}(0)$, and $\beta$. For $\delta$ we consider either $\delta=1$ (see figure 3(c)), $\delta=2$ (see figure 3(d)), or we also use $\delta$ as an fit parameter (see figure 3(e)). Note that due to the large number of fitting parameters the fitting process is not stable for all combinations of starting values. Therefore, for A₁, A₂, B₁, B₂ we employ the values that we have obtained by the method described in the previous paragraph as starting values. Our fits lead to the following results: for $\delta=1$: $\phi_{g}(0) = 0.663\pm 0.007$ and $\beta= 0.215\pm 0.009$, for $\delta=2$: $\phi_{g}(0) =0.690\pm 0.008$ and $\beta=0.208\pm 0.009$, for fitted $\delta$: $\delta=1.51\pm 0.22$, $\phi_{g}(0) = 0.677\pm 0.010$ and $\beta=0.206\pm 0.009$.

Therefore, in all cases the additional correction of $\phi_{g}(T)$ seems to be well-described by a power law with an exponent $\beta\approx 0.2$ as expected from our analysis presented in section 3.1. Note that while the temperature-dependence of our fitting curves relies on the large temperature-range of the data, due to the limited range in packing fraction we cannot deduce from our fits which exponent $\delta$ or which other details of the Vogel–Fulcher–Tammann curves as functions of the packing fraction would be best to describe glassy dynamics.

In order to test our scaling approach for simulation data close to the glass transition, we have a closer look on the data of of Berthier and Witten [9, 10], who demonstrated that the relaxation times of soft spheres can be quite well collapsed by suitable rescaling of the packing fraction and the temperature onto one universal function below the glass transition density and one function above the glass transition density. Here we show that the collapse of the data can be improved by employing the scaling relation that we propose in this article. To be specific, we show that with our approach it is not necessary to choose an exponent $\mu$ that differs from the theoretically expected value $\mu=1$.

As outlined in the introduction, according to Berthier and Witten [9, 10] the best collapse is obtained by assuming an effective packing fraction $\phi_{\rm ef\,\!f}=\phi-a T^{\mu/2}$ with $\mu=1.3$, which is intermediate between $\mu=1$ as expected for distinct overlaps between harmonic spheres (see also, e.g. [7]) and $\mu=3/2$ as claimed for ballistic, dilute particles [9]. Furthermore, the exponent of the Vogel–Fulcher–Tammann function is set to $\delta=2.2\pm 0.2$ and the best fit is claimed to occur for a glass transition at packing fraction $\phi_{g,0}=0.635\pm 0.005\pm$ [9, 10]. The collapse is obtained by plotting $\left|\phi_{g,0}-\phi\right|^{\delta}\log{\left(\sqrt{T}\tau_{\alpha}\right)}$ as a function of $\left(\phi_{g,0}-\phi\right){}^{2/\mu}/T$ [9, 10]. The collapse works well for packing fractions close to $\phi_{g,0}$ and for small temperature, but essential deviations occur for larger $\left|\phi_{g,0}-\phi\right|$ or larger T. We will show in the following that better collapses can be obtained (even with less free fitting parameters). Note that the limits of the scaling approach of Berthier and Witten in [9] can be easily seen by considering the curvature in log–log-plots of the curves of $\log{\left(\sqrt{T}\tau_{\alpha}\right)}$ as function of 1/T with constant packing fraction: according to the approach of Berthier and Witten all curves below $\phi_{g,0}$ should be concave in a log–log-plot while all curves above $\phi_{g,0}$ should be convex such that they can be collapsed onto the universal scaling functions. This is because the rescaling according to Berthier and Witten corresponds to affine transformations of the curves in the log–log-plot such that concave curves are always mapped on concave curves and convex curves always on concave ones. From looking at $\log{\left(\sqrt{T}\tau_{\alpha}\right)}$ as function of 1/T in a log–log-plot one therefore would expect the glass transition to occur at a packing fraction larger than 0.662 where there is a transition from concave curves to maybe convex curves (see red line in figure 4(a) where we replot the data extracted from figure 2(a) of [9] in the described log–log-representation). To be specific, while for packing fractions larger than 0.662 it is hard to tell whether the curves are convex or concave, below or at 0.662 all curves are concave, i.e. they cannot be shifted such that they suit the convex scaling function of the upper branch as claimed in [9, 10]. For comparison, the transition reported in [9, 10] is indicated by a green line in figure 4(a). All curves to the left of that line in [9, 10] are displaced such that they somehow arrange on a convex curve though some of the individual curves clearly are concave.

Zoom In Zoom Out Reset image size

Note that the rescaling of Berthier and Witten can also be achieved by plotting $\left|\phi_{g,0}-\phi\right|^{\delta}\log{\left(\sqrt{T}\tau\right)}$ as a function of $\left(\phi_{g,0}-\phi\right)/T^{\mu}$. In figure 4(b) we use this representation for the rescaling proposed by Berthier and Witten [9, 10]. Unlike Berthier and Witten in [9, 10] we plot all data points and not just those that collapse onto the scaling functions. We want to point out that similar collapses can be achieved for various choices of parameters, e.g. $\delta$ can also be smaller if $\mu$ is increased (see figure 4(c)). As a consequence, the parameters $\mu$, $\delta$, and $\phi_{g,0}$ are not that well determined as Berthier and Witten suggest. Concerning $\phi_{g,0}$ we show in figure 4(d) that a similar collapse can be achieved for a glass transition packing fraction of 0.662 that is expected by considering the transition between concave and maybe convex curves as explained in the previous paragraph.

Now we want to demonstrate that the data of Berthier and Witten agrees with the scaling that we propose in this article. Instead of rescaling the x-axis by $T^{-\mu/2}$ with $\mu=1.3$ as chosen in [9, 10] we choose the scaling of section 3.1, i.e. we rescale it by $(T^{\mu/2}-bT^{\beta}){}^{-1}$ with $\mu=1$, $\beta=0.2$, and a parameter b that we chose such to obtain the best fit. In principle other values for $\mu$ or $\beta$ could be chosen. We employ $\mu=1$ such that the effective packing fraction scales as in [7] and $\beta=0.2$ in agreement to [16] as explained in section 3.1. For the glass transition packing fraction we use $\phi_{g,0}=0.635$ as in [9, 10]. The y -axis is left unchanged in order to resemble the Vogel–Fulcher–Tammann behavior. In figure 4(e) we employ an exponent $\delta=2.2$ as in [9, 10], while in figure 4(f) we use $\delta=2$. In both cases the collapse of the data, especially concerning the upper branch, is superior to the collapse of Berthier and Witten and by further varying $\delta$ (not shown) we find that $\delta=2$ seems to be the best choice of $\delta$. Note that we can also use other values for $\phi_{g,0}$, but the best collapses seem to be reached for $\phi_{g,0}$ between 0.62 and 0.64.

Note that in contrast to the scaling proposed by Berthier and Witten, our way of rescaling is not an affine transformation and therefore one can map a concave curve onto a convex curve such that the collapse of the upper branch can be composed of only convex curves. The collapse of the lower branch can even be further improved if the logarithmic correction due to the prefactor of the assumed Vogel–Fulcher–Tammann extrapolation function is not neglected as in [9, 10] such that there is an additional fit parameter 1/c within the logarithm on the y -axis. To be specific $\left|\phi_{g,0}-\phi\right|^{\delta}\log{\left(\sqrt{T}\tau_{\alpha}/c\right)}$ can be used on the y -axis (see, e.g. figure 4(g)). In this case the best results are obtained for transition packing fractions $\phi_{g,0}$ between 0.64 and 0.65. As a consequence, a different transition packing fraction might be obtained if the collapse is further improved. Note that the data in [9, 10] was probably obtained in a way that simulations at larger packing fractions might have used the results of simulations at lower packing fractions as an input. As a consequence, in principle the ageing times for different data points might be different. Since the transition packing fraction is expected to increase with increasing ageing time [6, 16] this might slightly influence the analysis. Note that today much better equilibration methods are available [6]. However, we have no indication that the data in [9, 10] is not equilibrated sufficiently for our analysis.