{'fraction_non_alphanumeric': 0.03408647382800642, 'fraction_numerical': 0.017129970944099918, 'mean_word_length': 4.318265682656826, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical Heterogeneities Below the Glass Transition

We present molecular dynamics simulations of a binary Lennard-Jones mixture at temperatures below the kinetic glass transition. The ``mobility'' of a particle is characterized by the amplitude of its fluctuation around its average position. The 5% particles with the largest/smallest mean amplitude are thus defined as the relatively most mobile/immobile particles. We investigate for these 5% particles their spatial distribution and find them to be distributed very heterogeneously in that mobile as well as immobile particles form clusters. The reason for this dynamic heterogeneity is traced back to the fact that mobile/immobile particles are surrounded by fewer/more neighbors which form an effectively wider/narrower cage. The dependence of our results on the length of the simulation run indicates that individual particles have a characteristic mobility time scale, which can be approximated via the non-Gaussian parameter.


I. INTRODUCTION
Although glasses have already been studied for a long time, their complete understanding is still an open problem due to the complex behavior of their static and dynamic quantities [1]. We focus here on their dynamics which has been found to relax non-exponentially in the supercooled liquid and shows strong history dependence below the glass transition. The question arises whether this behavior is due to spatially homogeneous non-exponential dynamics or spatially heterogeneous dynamics (for review articles see [2][3][4]). Since a glass is an amorphous solid and therefore not all atoms are structurally equivalent, one should expect that also their dynamics differs, i.e. that the system has a heterogeneous dynamics. Recently the answer to the question of dynamic heterogeneity has been addressed both by means of experiments [5][6][7] and computer simulations of two-dimensional [8][9][10] and three-dimensional systems [11][12][13][14][15][16]. Here we study a binary Lennard-Jones mixture in three dimensions which has been investigated before extensively [15] and which shows clear dynamic heterogeneity above the glass transition. Similar dynamics has been found experimentally with confocal microscopy of a colloidal suspension in the supercooled fluid and in the glass [6,7]. Whereas most simulations have been done at relatively high temperatures above the calorimetric glass transition, and the experiments of atomic systems were performed near the glass transition, we simulate, in this paper, below the glass transition (which has so far only been done experimentally by Weeks et al. [7] and in simulations by Oligschleger et al. [16]). We find here dynamic heterogeneity via simulations of the same binary Lennard-Jones system as [15] but below the glass transition. In contrast to the previous simulations we have the picture of a solid in mind, instead of coming from the liquid. We use the "localization length" of the work of reference [17] to define the mobility of a particle as the mean fluctuation around its average position. To address the question of what allows or inhibits a particle to be mobile, we study the surrounding of these particles. Using different lengths of simulation runs we also learn about the time scale over which mobile and immobile particles sustain their character.
We review in section II the model used and give details of the simulation. In section III we present the mean square displacement and the mean fluctuations of a particle around its average position and define what we mean by mobile and immobile particles. We then study their spatial distribution (sec. IV), surrounding (sec. V) and time scale (sec. VI) and conclude with section VII.

II. SIMULATION DETAILS
We study a binary Lennard-Jones (LJ) mixture of 800 A and 200 B particles. Both A and B particles have the same mass. The interaction between two particles α and where ǫ AA = 1.0, ǫ AB = 1.5, ǫ BB = 0.5, σ AA = 1.0, σ AB = 0.8 and σ BB = 0.88. We truncate and shift the potential at r = 2.5σ αβ . From previous investigations [18] it is known that this system is not prone to crystallization and demixing. In the following we will use reduced units where the unit of length is σ AA , the unit of energy is ǫ AA and the unit of time is mσ 2 AA /(48ǫ AA ). We carry out molecular dynamics (MD) simulations using the velocity Verlet algorithm with a time step of 0.02. The volume is kept constant at V = 9.4 3 = 831 and we use periodic boundary conditions. We are interested in the dynamics of the system below the glass transition. Since recent simulations [18] showed that for present day computer simulations the system falls out of equilibrium around T ≈ 0.44, we run (NVE)-simulations at temperatures T = 0.15/0.2/0.25/0.3/0.35/0.38/0.4/0.41/0.42 and 0.43. To do so we start with a well equilibrated configuration at T = 0.466. After an instantaneous quench to T = 0.15 we first run a (NVT)-simulation [19] for 10 5 MD steps to let the system anneal, and then run the production run with a (NVE) simulation for also 10 5 MD steps. We then increase the temperature to T = 0.2 and then again run a (NVT)-simulation followed by a production run each with 10 5 MD steps; then increase the temperature to T = 0.25 and so forth. In this paper we refer to the so obtained production runs as "short runs," of which some preliminary results have been published elsewhere [20]. We present here mainly results of the so-called "long runs" for which we use the configurations at the end of the equilibration period of the short runs but the production runs are for 5 · 10 6 MD steps. To improve the statistics we run 10 independent configurations for both long and short runs and for each temperature.
As it has been demonstrated in earlier work [21], the structural properties of glasses studied by molecular dynamics simulation do depend on the preparation history quite distinctly. Since we study the system out of equilibrium and at finite temperature, the resulting configurations show some "aging phenomena" during the time intervals used for the production of the present data. However, for those temperatures where even the mean square displacements of the 5% fastest particles are still small in comparison to typical interparticle distances over the whole time span of the averaging, this change of the glass structure due to aging should have a relatively small effect on our data. As will be seen below, this is the case for T ≤ 0.35, while for somewhat higher temperatures aging effects can be expected to become important. The reason for this is that τ , the typical relaxation time of the system in equilibrium, is at T = 0.446 around 800, 000 time units (= 4 ·10 7 MD steps) [22], thus only a factor 10 longer than the long runs in the present work. Since in equilibrium the α−relaxation time corresponds to the typical time scale on which a substantial fraction of the particles have moved (40-70%), it can be expected that quite a few particles will show relaxation even on time scales significantly shorter than τ . Similar relaxation processes are also expected in the out of equilibrium situation, i.e. in the glass and thus we do indeed expect aging effects at temperatures slightly below the (kinetic) glass transition. We shall comment on this problem that there is some aging of the glass structure occurring when appropriate.

III. MOBILE AND IMMOBILE PARTICLES
Similar to previous work on dynamic heterogeneities [6,7,14,15], we study the dynamics of the system by observing the fastest (mobile) and the slowest (immobile) particles. Since the focus of this work is, however, on the dynamics of the glass below the glass transition, our definition of mobile and immobile particles is different. We have in mind the picture of a harmonic solid for which the vibrational amplitude carries essential information about the local dynamics. We therefore characterize the mobility of each particle i by where the bar denotes an average over a certain time interval. We call the 5% A particles and separately the 5% B particles with the largest/smallest d 2 i the mobile/immobile particles. With "mobile" we intend to indicate that these particles are relative to all other particles more mobile, they are however in the results, presented here, in most cases still bound to their site. All results presented below are qualitatively the same if the 5% are replaced by 10% particles and are therefore independent of the specific percentage used. The results depend however on the length of the time average, as will be discussed in detail in section VI. We use in this paper the non-Gaussian parameter [9,15,23] to determine the length of the time average as follows. The non-Gaussian parameter is defined as where · corresponds not to the canonical average since the system is out of equilibrium. Instead we mean by · here and in the following an average over particles and initial configurations, i.e.
where the sum goes over either all A particles (to obtain α 2A ) or all B particles (to obtain α 2B ). This parameter vanishes if the van Hove correlation function is equal to a Gaussian [24] G s (r, t) = 3 2π r 2 (t) We expect Eq. (6), and therefore α 2 = 0, to be a good approximation for t → 0 (because this corresponds to the ballistic regime where r(t) ∝ v · t which is Maxwell distributed [23]) as well as for t → ∞ (diffusive behavior). For intermediate times α 2 (t) = 0 (see Fig. 1). Although the non monotonous temperature dependence of α 2 (t) in Fig. 1 indicates that our statistics is not very good, we use the time t max where α 2 (t) reaches its maximum as an estimate for the characteristic time scale of mobility. (Note that for T ≤ 0.25 we use t max = 10 5 .) In all following results which involve time averages, we use t max , if not otherwise stated, as the time length over which we average. Note that t max is much larger than the microscopic oscillation time which is of the order of 1.0. We obtain thus for each temperature and for either all A or all B particles the distribution of d 2 i defined in Eq.(2). Fig. 2 shows for the A particles that with increasing temperature the distribution shifts to the right and develops a longer tail. (Similar results are obtained for the B particles.) The curves are zero for d 2 i < 0.0025, which reflects the fact that all particles are oscillating somewhat. The tail of P (d 2 i ) extends for high temperatures to values of d 2 i that are twice as large as the d 2 i at the peak position, which shows that the dynamics is rather heterogeneous. We also mention that the P (d 2 i ) for different T can be collapsed onto a single curve by rescaling the distribution to [25]). Fig. 3   more than linearly. The discrepancies from a straight line through the origin for A and B particles at temperatures T > ∼ 0.25 show that anharmonic effects become important already at small temperatures. As suggested in [26], this onset of anharmonicity may be related to the calorimetric glass transition. The log-log plot of the mean square displacement where we average over all A particles (see Fig. 4), shows that the slopes at large times are smaller than one. Therefore most particles never reach the diffusive region (m = 1) and are trapped in their cages during the whole simulation run at least for T ≤ 0.35. We find the same for the B particles where m < ∼ 0.93. If we average over only the 5% particles with the largest r 2 (t end = 10 5 ) , the late time slopes are at T = 0.43 m ≈ 1.4 for the A and m ≈ 1.2 for the B particles. This transient behavior of m > 1 might be due to jump processes. Jumps are clearly visible for the B-particles in Fig. 5, which shows the 5% particles with the smallest r 2 (t end ) . At low temperatures the slowest A particles are trapped at their site as can be seen in Fig. 5, since r 2 (t) < 10 −2 over the whole simulation run. Note also that the slowest B particles are faster than the average A particles. Fast B particles at T = 0.43 are reaching values of even r 2 (t end ) > 10 (see Fig. 6).

IV. SPATIAL DISTRIBUTION
With the definition of mobile and immobile particles as given in the last section, we study now how they are spatially distributed. Fig. 7 and 8 show the spatial distribution of the mobile (light spheres) and immobile (dark spheres) at temperatures T = 0.15 and T = 0.43, and at the beginning of the time interval for which their mobilities are determined. For clarity all other 900 particles are not shown. In these snapshots, and similarly for all other temperatures and times, the particles are clearly distributed in a heterogeneous way. We therefore find dynamic heterogeneity for all investigated states in the glass phase. The number of particles in the largest cluster [27] is at all investigated temperatures for mobile particles about 30 particles and for immobile particles about 22 particles. These clusters are smaller than the clusters of Weeks et al. [7]. The likely reason for this discrepancy is that we study a smaller system and with less good statistics. To quantify the spatial heterogeneity we plot similar to [15] the ratio g mAmA /g AA between the radial pair distribution [28,29] of solely mobile particles and that of all particles ( Fig. 9) with In the case of randomly selected 5% particles from all A particles, this ratio would be one. We find however that this ratio is not a constant with respect to r (and similarly for the corresponding ratios of AB and BB), which confirms the dynamic heterogeneity. Since g mAmA /g AA > 1.0 for distances r < ∼ 3.2, mobile particles tend to be near each other. We can conclude from the position of the first peak of g mAmA /g AA (see inset of Fig. 9) that separation distances which are for average particles very unlikely, as at the left wing of g AA , occur for mobile particles more often.
We can draw similar conclusions for the immobile particles. Fig. 10 shows that the ratio of the radial pair distribution of immobile particles to that of all particles, g iAiA /g AA , is also larger than one for small distances. The inset, which includes g AA for comparison, reflects that also for immobile particles very small separation distances are more likely.

V. SURROUNDING
In the last section we found that the mobile/immobile particles form clusters and that the typical distances between the particles are different from those in the bulk. We now address the question of the reason for mobility by comparing the surrounding of mobile/immobile particles to the one of average particles.

A. Coordination Numbers
To probe the immediate neighborhood of the mobile and immobile particles we count their number of nearest neighbors (coordination number z) where a particle j is defined to be a neighbor of particle i if their distance | r ij | = | r i − r j | is smaller than the position of the first minimum r min of the corresponding (average) radial pair distribution function (r min = 1.4 for AA, 1.2 for AB and 1.07 for BB, independent of temperatures). Fig. 11 shows that a mobile A particle is on average surrounded by fewer particles, z mAA+B , than an average A particle, z AA+B . This sug- gests that one of the reasons for a particle to be mobile is that it is on average caged by fewer particles. The same is true for mobile B particles (see Fig. 12). However, we cannot make the stronger statement that any A particle with z mAA+B less than, say, 13.4 is mobile, because the distribution of coordination numbers is quite broad with a standard deviation of σ ≈ 2.0. Mobile A particles are furthermore surrounded by a lower than average percentage of B neighbors (see Fig. 13), because the latter trap A particles both energetically (ǫ AB > ǫ AA ) as well as geometrically (σ AB < σ AA ).
Similarly immobile particles have the property to have more neighbors than average particles (see Fig. 14) and of a higher percentage of B particles than usual (see Fig. 15 and 16). Notice that the latter is true both for A and B particles (see Fig. 15 and 16) due to the tighter packing with the smaller B particles.
We also find that the percentage of mobile/immobile neighbors of a mobile/immobile particle is significantly larger than 5%, i.e., z mAmA+mB > 0.05 · z mAA+B and z iAiA+iB > 0.05 · z iAA+B , which reflects once more the spatial heterogeneity discussed in Fig. 9 and 10.

B. Radial Pair Distribution Functions
Next we use the radial pair distribution function to study the environment of the mobile and immobile particles beyond the nearest neighbor shell. Fig. 17 shows the radial pair distribution function of a mobile A particle with any B particle, g mAB , in comparison with g(r) of any A and B particles, g AB , at T = 0.15. We find that g mAB has smaller maxima and broader peaks than g AB which corresponds, specifically for the first neighbor    Radial pair distribution function of a mobile A particle and any B particle (solid line) and for comparison the corresponding radial distribution function of any A and B particles (dashed line). The inset shows the similar g(r) for the immobile A particle with any B particle. Both graphs are for T = 0.15.
shell, to an effectively wider cage around the mobile particles [30]. The wider cage allows larger distances and thus larger d 2 i which corresponds by definition to mobile particles. We see the same effect for g mAA (for g mBB the statistics is not sufficient) and for all other temperatures.
Immobile particles, in contrast, are surrounded by an effectively narrower cage, as can be concluded from the more pronounced peak of the first maximum (see inset of Fig. 17). This is observed for all temperatures and all radial pair distribution functions characterizing the surrounding of an immobile particle. Notice that the change of the neighborhood is larger around an immobile particle than around a mobile particle, as the comparison of Fig. 13 and Fig. 15 and the comparison of Fig. 17 and its inset show. This is probably due to our definition of mobility: since the distribution of d 2 i (see Fig. 2) is very asymmetric, 5% particles with the smallest d 2 i cover a much smaller range of d 2 i than 5% particles with the largest d 2 i . Immobile particles are thus more distinct than mobile particles.

VI. TIME SCALE OF MOBILE AND IMMOBILE PARTICLES
In this section we get back to Eq. (2), which is essential for the definition of mobile and immobile particles. We vary the time length over which we average. Specifically, we average over the simulation time for the long and short runs (see Sec. II), rather than using α 2 to determine a temperature-dependent time for calculating the average. We now investigate the influence of this averaging time on the results presented in the previous three sections.   18 shows a comparison of the total number of nearest neighbors of a mobile A particle for the short runs, z s mAA+B , and the long runs, z l mAA+B . All coordination numbers have been averaged over the 10 independent initial configurations, and no data from later times are included. The difference in z s mAA+B and z l mAA+B is solely due to the different definition of mobile particles. As before (see Fig. 11) we find that the mobile particles are surrounded by fewer particles. For the long runs, however, this effect vanishes at higher temperatures. The likely reason for this decreasing difference is that we used in Eq. (2) a time average over the entire long simulation run, while a fast particle might be fast only over some fraction of the simulation run. While the particle is fast, its environment is different than that of a regular particle. However, when the time average includes also times when the particle is not fast, then we dilute the average with environments that are not special. This mixing with the environment of average particles happens more readily at temperatures T > ∼ 0.4 when we approach the glass transition because the typical time scale over which a particle is fast is shorter than at lower temperatures.
Other coordination numbers show the same behavior with the exception of immobile A particles (see Fig. 19). The latter distinguish themselves from the average particle even at high temperature and thus are immobile over the whole simulation run (consistent with Fig. 5).
Another quantity which is also strongly dependent on the time length of the run is the distribution of 1/d 2 i (see Fig. 20 for long runs and its inset for short runs). In the case of long runs a double peak structure develops for increasing temperature. We tentatively associate the particles in the peak with large 1/d 2 i with localized particles and those in the peak with small 1/d 2 i with mobile parti- cles. Further work is required to check this hypothesis.
The peak at small 1/d 2 i starts to dominate with increasing temperatures, which reflects that at high temperatures most particles have become fast at some time during the simulation run.
As the figures of Secs. IV and V demonstrate, these effects of averaging over fast and regular particles can be avoided via appropriate choice of averaging time t max which we chose to be when α 2 (t) reaches its maximum. Therefore t max gives us a rough estimate about the lifetime of fast particles.

VII. SUMMARY AND OUTLOOK
We investigate the mobile and immobile particles of a glass. Note that our definition of mobile and immobile is different than the definition given in [15] since we have the picture of a solid in mind (similar to G. Johnson et al. who study solid-like particle clusters [10]).
We find below the glass transition temperature a clear dynamic heterogeneity which has previously been seen in simulations above T g and experiments below T g .To address the question why certain particles are more/less mobile than others, we study their surrounding. As one might have expected, the mobile/immobile particles are surrounded by fewer/more neighbors, forming a cage which is effectively wider/narrower than the one of a regular particle. In addition mobile/immobile particles are trapped by fewer/more B particles, which are smaller than the A particles and therefore allow closer packing [31]. We expect that, similarly, a surrounding specific to the mobility of the central particle might be found in the future in experiments. Both the dynamic heterogeneity as well as the particular surrounding of mobile and immobile particles are consistent with collective behavior as it has been found above the glass transition temperature T g . A more detailed analysis below T g is left for future work.
The characteristics of mobility show a time dependence which is well estimated with the time t max of the maximum of α 2 . We conclude that mobile particles are "fast" only for a certain time window of the simulation run whereas the immobile A particles seem to stay mostly immobile over the range of our simulation runs. This raises the question of a more precise criterion for the time scale of fast and slow processes, which we leave for future work.

VIII. ACKNOWLEDGMENTS
KVL gratefully acknowledges financial support from the SFB 262.