Signature of Dynamical Heterogeneity in Spatial Correlations of Particle Displacement and its Temporal Evolution in Supercooled Liquids

The existence of heterogeneity in the dynamics of supercooled liquids is believed to be one of the hallmarks of the glass transition. Intense research has been carried out in the past to understand the origin of this heterogeneity in dynamics and a possible length scale associated with it. We have done extensive molecular dynamics simulations of few model glass-forming liquids in three dimensions to understand the temporal evolution of the dynamic heterogeneity and the heterogeneity length scale. We find that although the strength of the dynamic heterogeneity is maximum at a timescale close to characteristic $\alpha$-relaxation time of the system, dynamic heterogeneity itself is well-developed at timescale as short as $\beta$-relaxation time and survives up to a timescale as long as few tens of $\alpha$-relaxation time. Moreover, we discovered that temperature dependence of heterogeneity length remains the same in the whole time window although its absolute value changes over time in a non-monotonic manner.

The existence of heterogeneity in the dynamics of supercooled liquids is believed to be one of the hallmarks of the glass transition. Intense research has been carried out in the past to understand the origin of this heterogeneity in dynamics and a possible length scale associated with it. We have done extensive molecular dynamics simulations of few model glass-forming liquids in three dimensions to understand the temporal evolution of the dynamic heterogeneity and the heterogeneity length scale. We find that although the strength of the dynamic heterogeneity is maximum at a timescale close to characteristic α-relaxation time of the system, dynamic heterogeneity itself is well-developed at timescale as short as β-relaxation time and survives up to a timescale as long as few tens of α-relaxation time. Moreover, we discovered that temperature dependence of heterogeneity length remains the same in the whole time window although its absolute value changes over time in a non-monotonic manner.
Dynamic heterogeneity (DH) is ubiquitous in a vast variety of natural processes spanning from molecular systems to biological cells and tissues. Existence, characterization and its role in different dynamical processes particularly in the dynamics of glass forming liquid approaching glass transition, is an active field of research. The dynamics of the supercooled liquid become increasingly heterogeneous [1] as the system approaches the putative glass transition point. Different regions of the system relax at time scales that differ from each other by several orders of magnitude. It is argued that the slowing down of dynamics arises mainly due to the cooperative motion of the particles. The typical size of these cooperatively rearranging regions (CRR) is believed to be of the order of few particles diameter within which motions of particles are spatially and temporarily correlated. It is already shown both in experiments [2][3][4][5] and in computer simulations [6][7][8][9] that the mobile particles are non uniformly distributed in the system and tend to form a cluster. This inclination of clustering, as well as the size of the cluster increases as glass transition, is approached.
Extensive studies [1,3,10,11] have been performed in past to understand the behaviour of DH at the characteristic long relaxation time scale or the α-relaxation time scale, τ α (defined later) [12] and only a handful studies are done at shorter β-relaxation time scale [13]. In [13], it has been shown that β-relaxation time, τ β (defined later) has a strong finite size effect, which can be rationalized if one assumes the existence of a growing correlation length. It was surprisingly found that the temperature dependence of this growing correlation length at β-relaxation time is the same as that of the heterogeneity length scale obtained at α-relaxation time. This observation is very surprising as these time scales can differ by many orders of magnitude, especially at low temperatures.
The main goals of this work are two-fold. The first goal is to find signatures of DH in the dynamics at the short time scale of the order of τ β . Then we would like to understand the subsequent growth and temporal evolution of DH at timescales ranging from τ β to an order of magnitude larger than τ α . As most of the research works have focused on the characterization of DH in the α-relaxation timescale, it is very important to comprehend the time evolution of DH in the intermediate as well as long-time scale compare to τ α , in order to understand the role of DH in glass transition. The main results of this work are the observation of signature of DH in the displacement fields of particles over β time scale and the survival of DH at timescales that are larger than τ α by a factor of 30 − 40 in the studied temperature range. We have also discovered that temperature dependence of the heterogeneity length scale remains the same throughout the studied time window, but the region of heterogeneity or the spatial extent of heterogeneity changes with time in a non-monotonic way with its maximum appearing at or near τ α .
Although in [13], it was shown that DH seems to be already developed in the system at timescale close to τ β , it was not immediately clear how particle motions at this short time scale get affected due to the presence of the heterogeneity, in other words, whether particle motions at τ β are correlated over the length scale of dynamic heterogeneity, is not immediately clear. Following [13], we define τ β to be the time at which logarithmic derivative of MSD develops a minimum which corresponds to an inflection point in the log − log plot of MSD (see SI for further details).
To measure the spatial correlation and to extract the associated length scale in the displacements of particles at τ β , we have implemented the procedure given in Ref. [14][15][16]. Note that this measure of the spatially correlated motion in super-cooled liquids do not depend on arbitrary cutoff parameters as already conclusively shown in Refs. [14][15][16] for DH at τ α . The spatial correlation of the particle displacements g uu (r, ∆t) is defined as (1) where u i (t, ∆t) = r i (t + ∆t) − r i (t) is the vector displacement of the particle between time t and t+∆t. u(∆t) = 1 N N i=1 |u i (t, ∆t)| , is the average displacement of particles within time interval ∆t. r ij (t) = r j (t) − r i (t), is the distance between i th and j th particles. Note that our definition of displacement-displacement correlation is slightly different from the definition given in Ref. [14]. In Ref. [14], scalar displacement of the particles are considered, whereas we have considered the vector displacement of the particles [17]. It captures the orientational as well as translational correlation in the particle displacements within the time of observation.
We have performed extensive computer simulation of four well-studied model glass formers in three dimensions with different inter particle potential over a wide range of temperatures. The model systems studied are the following (i) 3dKA [18], (ii) 3dR10 [19], (iii) 3dIPL [20] and (iv) 3dHP [21]. The details of the models and simulations are given in the Supplemental Information (SI). We find that the growth of DH identified using displacementdisplacement correlation function show strong system size dependence at least for at τ β . This was not the case in [22,23] when the correlation function was computed at τ α . Thus we demonstrate the results of g uu (r, ∆t) for the above models for two different system sizes N = 10000 and N = 108000.
In Fig. 1 (left panel) we show the g uu (r, ∆t) for N = 108000. It is observed that g uu (r, ∆t) exhibits damped oscillation which is in agreement with previous numerical [16,24] as well as experimental studies [25]. The correlation function decays to zero exponentially as a function of distance r and with decreasing temperature, the dynamics of the liquid becomes more heterogeneous as the correlation between the particles' displacements in space extend up to a larger distance as shown in the left panel of Fig. 1. It physically means that particles in the liquids are moving in a co-operative fashion with a monotonically increasing size of the co-operative region as the temperature is lowered. We find a strong system size dependence in g uu (r, ∆t) as shown in the inset of the left panel of Fig. 1. We have computed the correlation for N = 108000 and N = 8000 for the 3dKA model as shown in the inset of Fig. 1 which shows that at low temperature, relative correlation increases with increase in system size. For the robustness of our results, we have calculated the correlation for the other two model systems. We found the results are qualitatively similar (see SI for further details).
In the right panel of Fig. 1, we show a plot of heterogeneity length scale as a function of temperature for the 3dKA model. As expected we observe strong finite size effects on the temperature dependence of DH length obtained from displacement-displacement correlation function. Heterogeneity length obtained from very large system size (N = 108000) grow very similarly with dynamical length scale obtained from finite size scale (FSS) of τ β [13]. For N = 8000 system size one observes that the heterogeneity length grows mildly, thus studies on smaller system size would have lead to a conclusion that DH is not present at β-relaxation time. After conclusively establishing the existence of DH at short time scale, we now focus on the temporal evolution of DH and the associated length scale across the whole range of timescales that can be accessed in simulation.
In [23], although DH was studied over different time scales but a systematic study on the temperature depen-  dence of the DH length scale was not done. Equipped with the method of block analysis, introduced in Ref. [26], a systematic study of the temperature dependence of the dynamical length scale across different time scale for different model glass-forming liquids became computationally feasible. Usually, four-point correlation function, g 4 (r, t) and the corresponding susceptibilities, χ 4 (t) [27] are used to study DH.χ 4 (t) is related to fluctuations in two-point function, Q(t). The Fourier transform of g 4 (r, t) is known as four-point structure factor S 4 (q, t) and related to χ 4 (t) as lim q→0 S 4 (q, t) ≡ χ 0 (t). τ α is defined as Q(t = τ α ) = 1/e, where . . . denotes ensemble averages. (See SI for further details and definitions). To perform the block analysis, we equilibrate a large system of N = 108000 particles and measure various quantities of interest by coarse-graining over different block sizes, L B . We then obtain the dynamic length scale ξ d by FSS analysis of χ 4 (L B , t).
In previous studies [26,28], dependence of χ P 4 (L B ), the maximum intensity (peak value) of dynamical susceptibility on block size was studied and the dynamic heterogeneity correlation length ξ d has been estimated by FSS analysis using the following scaling form where χ P 4 (∞, T ) is the L B → ∞ value of dynamical susceptibility at a temperature T . In this work we have done similar scaling analysis for the block size dependence of the intensity of dynamical heterogeneity at few particular time scales t = τ α /3, τ α /2, 3τ α to understand how the characteristic length scale and the exponent associ-ated with dynamical susceptibility and the heterogeneity length scale change with time at different temperatures.
In Fig. 2 we have plotted the results for the 3dKA model system. In top panels, we reported the block size dependence of χ 4 (L B , T ) (inset) and the scaling collapse of χ 4 (L B , T ) at τ α /3 (left panel) and τ α /2 (right panels) time scales respectively. A similar analysis also is shown in the bottom panels at time scales τ α (left panel) and 3τ α (right panel). The scaling observed in all these four cases are indeed very good and the calculated length scales from FSS analysis of the block method is found to be in good agreement with the ξ d obtained from the wave vector dependence of S 4 (q, t) [29] (discussed in subsequent paragraph). In FSS, ξ d (T ) is known up to a multiplicative factor which is the same for all temperatures. In order to fix this uncertainty, ξ d obtain from FSS is scaled to match with ξ d obtained from S 4 (q, t) at one temperature.
By fitting the q dependence of S 4 (q, t) for small q values to the Ornstein-Zernike (OZ) form S 4 (q, t) χ 0 (t)/[1 + (qξ) 2 ], one can also obtain ξ d . It has already been shown that heterogeneity length scale obtained from FSS of block method is in good agreement with the same obtained from S 4 (q, t). In top panels of Fig. 3 we plot the wave vector dependence of the inverse of four point structure factor S 4 (q, t) for the 3dKA model for two different times τ α /3 (left), 3τ α ) (right). One can clearly see that OZ form fits the data very well, thus the extracted length scale will be quite accurate. In the bottom panels of the same figure, the temperature dependence of the length scales computed by different methods are compared for t = τ α /3 (left) and t = 3τ α (right). The legend "ξ d from FSS" refers to the length scale obtained from finite size scaling at t = τ β and taken from [28]. Note that ξ d from FSS are scaled at T = 0.80. It is worth highlighting that these results suggests that temperature dependence of ξ d is same across timescales starting from τ β to at least 3τ α . To check the robustness of our results, length scales for other two models (3dR10 and 3dIPL) are also computed and the temperature dependence of the length scale are found to be same over the time interval (τ β , 3τ α ) (see SI for further details). Although variation of heterogeneity with decreasing temperature remains same over studied timescales across model systems, the spatial extent of the heterogeneity is observed to increases up to some timescale and then starts to decrease. We are not able to estimate the dynamical length from FSS of block methods at very early time (t < τ α /3) and very long time (t > 3τ α ) as the intensity of χ 4 (t), itself becomes extremely small to do any computation reliably.
Next we examine the power law relation between χ 4 (T ) and ξ d (T ) across these studied timescales. According to Inhomogeneous mode coupling theory (IMCT) [30][31][32] there exists a power relation between χ 3 (τ α ) (a three point correlator which is believed to be similar to χ 4 at least in the scaling behaviour) and ξ d (τ α ) as χ 3 (τ α ) = where the theory predicts exponent 2 − η = 4 in the α regime and 2 for the β regime [33]. Following Ref. [26] we have done the scaling analysis of χ t * 4 (T ) to obtain the exponent 2 − η at different times, t * . In the large system size limit, L B >> ξ d , the L B dependence should disappear in the scaling relation in Eq.2 and it should approach a constant value for x >> 1. On the other hand for ξ d → ∞ and L B remains finite the dependence of χ 4 on ξ d should go away. This implies that scaling function f (x) should be proportional to x 2−η at x → 0 and χ t 4 (L B , T ) should grow as L 2−η B . We show the results of such an analysis for the 3dKA model in Fig. 4. The exponent value is found to be 2 − η 2 for both t = τ α /3 as well as at 3τ α . Note this is completely different from the exponent (2 − η 4) predicted by IMCT in the α regime but in very good agreement with the prediction at β regime. This observation can be rationalized if one assumes that there will be less activated relaxation at a short time scale as compared to α relaxation time scale and MCT approximation should then be reasonable. Thus we can expect to have a reasonable agreement with MCT predictions at short timescales but not so good at longer timescales. Our results very nicely highlights this agreement with good quality data. We have done similar analysis for other models systems (3dIPL and 3dR10) and found that the exponent 2 − η is very close to 2 (see SI for details).
Next, we look at the time evolution of the DH length scale. In [29], it was shown that time dependence of ξ d (t) is same as χ 4 (t), which is in contradiction with the results reported in [34]. In [34], ξ d is found to increase monotonically with time, in partial agreement with the results reported in [35] for hard sphere systems. Moreover in [35] it is found that ξ d saturates to a plateau at a later time. On the other hand, IMCT predicts ξ d to remain constant in between τ β and τ α . We then look at the In top panel of Fig. 5, the time dependence of ξ d (t) for the 3dR10 model is shown. It is clear that ξ d (t) grows up to τ α and decreases at later time, in agreement with Ref. [29]. We then looked at the results obtained for 3dHP model (see SI) [21], which is one of the paradigmatic models in the context of jamming physics. For the 3dHP model, ξ d (t) shows a peak at timescale close to 4τ α (bottom panel of Fig. 5), which implies that ξ d (t) increases even though the overall strength of the heterogeneity decreases after τ α . This suggests that a hard sphere like model are probably different from those models where the particles are treated as point particles. These observations although corroborate the previous observations, the reason for it is not immediately clear. A different choice of the steepness of the potential does change the qualitative results (see SI for detailed discussion). Although for 3dHP model ξ d (t) ∼ log(t) in agreement with Ref. [35] for hard sphere systems, one does not see similar dependence for 3dR10 model. The mutual dependence between ξ d (t) and χ 4 (t), as χ 4 (t) ∼ ξ d (t) 2−η , seems to have two different regime -up to t ∼ τ α it is power law like with exponent 2 − η ∼ 4 for all the models [4.18 (3dKA), 3.50 (3dR10) and 3.39 for 3dHP model] but a extremely different one with very large exponent for t > τ α as shown in the insets of Fig. 5. Thus it suggests that one will not be able to extract even approximately the dynamic heterogeneity length ξ d from the measurement of χ 4 (t) alone as the exponent 2 − η varies across model systems as well as over the time window of calculation.
Finally, to conclude, we have shown the presence of dynamic heterogeneity in the displacement field of particles at τ β timescale and also highlighted the strong system size effect. Then we showed that although the absolute value of the dynamic heterogeneity length changes with time, the temperature dependence of this length scale across timescale spanning from τ β to couple of times τ α remain same. This is really surprising and may become an important test for validation of different theories of glass transition. We also found that absolute value of ξ d reaches its maximum value at t ∼ τ α for 3dKA, 3dIPL and 3dR10 models but it does so at ∼ 4τ α for 3dHP models. These observations are indeed extremely intriguing and warrants further research as to why heterogeneity length scale continues to increase at timescales where χ 4 already decays to small values for soft sphere and hard sphere model systems.

Signature of Dynamical Heterogeneity in Spatial Correlations of Particle Displacement and its Temporal Evolution in Supercooled Liquids -Supplementary Information
Indrajit Tah

MODELS AND SIMULATION DETAILS
We have studied four different model glass forming liquids in three dimensions. The model details are given below: 3dKA Model: This is the well known 80:20 binary mixture of Lennard-Jones particles [1]. This model was first introduced by Kob-Anderson to simulate N i 80 P 20. The interaction potential is given by where α, β ∈ {A, B} and AA = 1.0, AB = 1.5, BB = 0.5; σ AA = 1.0, σ AB = 0.80, σ BB = 0.88. the number density is ρ = 1.20. The interaction potential was cut off at 2.50σ αβ . We use a quadratic polynomial to make the potential and its first two derivatives smooth at the cutoff distance. Length, energy and time scale are measured in units of σ AA , AA and 3dR10 Model: This is a 50:50 binary mixture [2] with repulsive interaction potential, defined as Here αβ = 1.0, σ AA = 1.0, σ AB = 1.18, σ BB = 1.40. The number density is ρ = 0.81. The interaction potential is cut-off at 1.38σ αβ .
3dIPL Model: The model was first studied in [3]. The interaction potential is given by All the parameters and interaction cut-offs are the same as 3dKA model. Though this model has a purely repulsive interaction but the interaction domain is much larger than that in the 3dR10 model. The interaction range plays an major role in determining both dynamic and mechanical properties of the system [4][5][6].
3dHP Model: We have simulated a model system that interpolates between finite-temperature glasses and hard-sphere glasses and has been studied extensively in the context of jamming physics. This is a 50:50 binary mixture with diameter ratio of 1.4 [7]. The interaction potential is given by . We simulated this model at a constant density ρ = 0.82 and temperature as a control parameter. The value of is chosen to be 2.0. We also did simulation by varying the α value. The α value we have used are the following, (α = 2 for harmonic sphere, α = 3 and α = 6). Many studies of this model have been done at finite temperatures [8][9][10] and it finds experimental realizations in soft colloids, emulsions and grains [11]. For this model we have done simulation of a very large system with N = 108000 particles.
We have done NVT simulations for all the model systems using velocity-Verlet integration scheme. For all the simulations we have first equilibrated our systems at least for 200τ α (defined later) and stored data for similar simulation time. The four-point dynamical susceptibility χ 4 (t) [12] is defined in terms of the fluctuations of the two point overlap correlation function Q(t) as The function Q(t) is defined as This function measures the overlap of a configuration of particles at a given initial time (t = 0) with the configuration at a later time t. w(x) is a window function defined as is w(x) = 1 for x < a and 0 otherwise. The window function w(x) is introduced to remove the apparent de-correlation that might happen due to the vibrational motion of particles inside the cages formed by their neighbours. To calculate τ β , we follow the procedure of Ref. [13,14], where τ β is calculated from the inflection point in the loglog plot of mean squared displacement (MSD) r 2 (t) . Mean square displacement ( r 2 (t) ) is defined as In Fig. 1, we plot the MSD in log-log as a function of time and log-derivative of MSD with time, d log |∆r(t)| 2 /d log t. From the minimum in the derivative, one can estimate the τ β time scale without much uncertainty.

DISPLACEMENT-DISPLACEMENT CORRELATION
To validate our results we measure the spatial correlation of particle displacement and associated length scale in the displacements correlation of particles at τ β time scale for 3dIPL model system. We find strong system size dependence of the displacement-displacement correlation as shown for two different system sizes (N = 108000 and N = 10000) in middle panel of Fig. 2. In right panel of Fig. 2, we plot the heterogeneity length scale as a function of temperature for 3dIPL model. Like 3dKA model, we also observe strong system size effects on the temperature dependence of DH length obtained from displacementdisplacement correlation function. Heterogeneity length obtained from very large system size (N = 108000) grow very similarly with dynamical length scale obtained from finite size scale (FSS) of τ β [14].

BLOCK ANALYSIS TO CALCULATE THE DYNAMICAL LENGTH SCALE
To calculate the dynamical length scale, we perform the finite size scaling (FSS) analysis of dynamical susceptibility χ 4 (L B , t) following the procedure of Ref. [15]. χ 4 (L B , t) is computed from the fluctuations of Q(L B , t) as Q(L B , t) is the two-point density-density correlation function computed for the particles that are residing in the block of size L B at the chosen time origin t = 0. (5) where N B is the number of blocks of size L B , n i is the number of particles in the i th block at time t = 0. The dynamic heterogeneity correlation length ξ d has been estimated by FSS analysis using the following scaling form where χ t 4 (∞, T ) is the L B → ∞ value of dynamical susceptibility at temperature T and time t.

TEMPERATURE DEPENDENCE OF THE DYNAMICS HETEROGENEITY LENGTH SCALE OVER DIFFERENT TIME FOR 3DR10 AND 3DIPL MODEL SYSTEM
In Fig. 3 we plot the block size dependence of χ 4 (L B , T ) and its scaling collapse at time scale τ α /3 and 3τ α for 3dR10 and 3dIPL model systems for different temperatures. For all these models observed scaling collapse are reasonably good and the calculated length scales are very good agreement with the previous reported length scales [15]. We found that temperature dependence of dynamical length scale ξ d remains same in between time scale τ α /3 to 3τ α as shown in Fig. 4. This results confirm the validity of the results over different model glass forming liquids.

SCALING FUNCTION AND THE SCALING EXPONENT FOR 3DR10 AND 3DIPL MODEL
We now examine the power law dependence of dynamical susceptibility χ 4 (T ) and the correlation length ξ(T ) at different time interval (e.g. t = τ α /3 and 3τ α ). The finite size scaling form used as following Here χ t 4 denotes the susceptibility at time t. χ 0 (t) = lim L B →∞ χ t 4 (L B , T ). As shown in [15], we do the similar scaling analysis of χ t 4 (T ) to obtain the exponent (2−η) at different time interval as detailed in the main article for 3dR10 and 3dIPL model systems. The exponent value is found to be 2 − η 2. This analysis is shown in Fig. 5. Growth of dynamic length scale as a function of temperature for 3dR10 model at time interval τα/3 and 3τα obtained from block analysis method and compared with the corresponding quantities which is taken from Ref. [15]. Bottom Panel: Similar comparison for 3dIPL model system.

CALCULATION OF DYNAMICAL LENGTH SCALE FROM FOUR POINT STRUCTURE
FACTOR S4(q, t) Most of the important information about dynamical heterogeneity length scale has been obtained from studies of multi-point correlation function and its associated susceptibilities [12]. To extract the dynamical length scale many numerical studies consider the four point time dependent correlation function defined as: g 4 (r, t) = δρ(0, 0)δρ(0, t)δρ(r, 0)δρ(r, t) − δρ(0, 0)δρ(0, t) δρ(r, 0)δρ(r, t) (8) δρ(r, t) is the deviation of local density ρ(r, t) at position r and time t from its average value ρ 0 (r, t), ... represents the thermal or time average. The fourier transform of g 4 (r, t) is known as four point time dependent structure factor S 4 (q, t) and its associated susceptibilities defines as lim q→0 S 4 (q, t) ≡ χ 0 (t).
To compute the dynamical heterogeneity length scale we compute four point structure factor S 4 (q, t) [16], which is defined as where where w(x) is the same window function defined before. T . This behaviour is similar of that observed in static structure factor S(q) near liquid-gas critical point where two point density fluctuation diverge at small q. Near the glass transition temperature or density, usual two-point density fluctuation does not show any diverging peak. In Fig. 6, we schematically represents the wave vector (q) dependence of S 4 (q, t) over different time interval. The intensity of heterogeneity increases at low q limit. The four-point dynamic correlation length is extracted by fitting S 4 (q, t) at small q limit with the Ornstein-Zernike (OZ) form S 4 (q, t) = S 4 (q → 0, t) 1 + (qξ) 2 In left panel of Fig. 6 we have plotted the wave vector (q) dependence of four point structure factor S 4 (q, t) over different time interval in the deep supercooled regime. To obtain the length scale we fit S 4 (q, t) to the Ornstein-Zernike (OZ) form in the range qξ ≤ 1.0. Now by using the value of ξ obtained from the OZ fit we did scaling plot S 4 (q, t)/S 4 (q → 0, t) vs qξ in right panel of Fig. 6 and we find a very good scaling collapse for all the model systems.

EXPONENT: η
We have calculated the exponent associated with χ 4 (t) and ξ d (t), as χ 4 (t) ∼ ξ 2−η d for all the models. Results for the 3dR10 and the 3dHP are given in the main article. Here we show the results for the 3dKA model system in Fig. 7. We obtained the exponent 2 − η = 4.175 up to t ∼ τ α and a very different behavior for t > τ α in agreement with the results obtained for 3dR10 model. For the 3dHP model we have shown that ξ d (t) has a peak at time t 4τ α . On the other hand in Ref. [17] it was shown that ξ d (t) saturates to a plateau for hard sphere systems. To test whether one will get similar results if one systematically tunes the interaction potential such that one asymptotically approach hard sphere like behaviour by changing the parameters of the potential, we did a similar analysis for the same 3dHP model but this time by changing the α parameter of the potential. It is clear that if one takes limit α → ∞, as lim α→∞ α 1 − r σ αβ α , then the soft sphere potential will tend towards the hard sphere case. Thus we did a similar analysis for two different choices of α, α = 3 and 6 to see if one obtains hard sphere like results. In the bottom panel of Fig.7, we have shown the data once more to highlight that within the studied range of values of α, the results do not change qualitatively. It is quite possible that if one takes much larger values of α, then probably one will get results similar to the hard sphere, but this is beyond the scope of the present work and we intend to do this analysis in future.