# Apparent finite-size effects in the dynamics of supercooled liquids

## Abstract

Molecular dynamics simulations are performed for a supercooled simple liquid with changing the system size from $`N=108`$ to $`10^4`$ to examine possible finite-size effects. Although almost no systematic deviation is detected in the static pair correlation functions, it is demonstrated that the structural $`\alpha `$ relaxation in a small system becomes considerably slower than that in larger systems for temperatures below $`T_c`$ at which the size of the cooperative particle motions becomes comparable to the unit cell length of the small system. The discrepancy increases with decreasing temperature.

As liquids are cooled toward the glass transition temperature $`T_g`$, a drastic slowing-down occurs in dynamical properties, such as the structural relaxation time, the diffusion constant, and the viscosity , while only small changes are detected in static properties. The goal of theoretical investigations on the glass transition is to understand the universal mechanism which gives rise to the drastic slowing-down. To this end, a great number of molecular dynamics (MD) simulations has been carried out for supercooled liquids . Several large scale simulations have been performed very recently and reviled that the dynamics in supercooled liquids are spatially heterogeneous ; rearrangements of particle configurations in glassy states occur cooperatively involving many molecules. We have examined bond breakage processes among adjacent particle pairs and found that the spatial distribution of broken bonds in an appropriate time interval ($`\tau _\alpha 0.1\tau _b`$ where $`\tau _\alpha `$ is the structural $`\alpha `$ relaxation time and $`\tau _b`$ is the average bond life time) is very analogous to the critical fluctuation in Ising spin systems. The structure factor is excellently fitted to the Ornstein-Zernike form , and the correlation length $`\xi `$ thus obtained grows rapidly with decreasing temperature. Furthermore, we demonstrated that $`\xi `$ is related to $`\tau _\alpha `$ trough the dynamical scaling law, $`\tau _\alpha \xi ^z`$ with $`z4`$ in 2D and $`z2`$ in 3D. The heterogeneity structure in our bond breakage is essentially the same as that in local diffusivity , which leads to a systematic violation of the Stokes-Einstein low in supercooled states.

To investigate long-time behavior of glassy materials by MD simulation, rather small systems typically composed of $`N=10^210^3`$ particles have been used with the periodic boundary condition (PBC). Such small systems have generally been considered to be large enough to avoid finite-size effects in the case of amorphas materials in which no long-range order exists. In fact, static properties such as the radial distribution function $`g(r)`$ or the static structure factor $`S(q)`$ of glassy materials are not seriously affected by the system size as long as reasonably large systems ($`N10^2`$) are used. However, this is not always the case for dynamical properties. For example, it is known that the use of a small system with PBC gives an manifest effect in relatively short-time behavior of the density-time correlation function. There appears an artifact in time scale of order $`L/c`$, where $`L`$ is the size of the simulation cell and $`c`$ is the sound velocity . As we already mentioned, the dynamical correlation length $`\xi `$ in supercooled liquids grows rapidly with lowering the temperature. It is thus possible that some kinds of finite-size effects may appear in the dynamics of supercooled liquids when $`\xi `$ becomes comparable to $`L`$ even if no such effect is detected in the static correlation functions. The main purpose of this paper is to examine carefully this point for a simple soft sphere mixture.

Our model mixture is composed of two soft sphere components $`1`$ and $`2`$ having the size ratio $`\sigma _1/\sigma _2=1/1.2`$ and the mass ratio $`m_1/m_2=1/2`$ while $`ϵ_1=ϵ_2=ϵ`$. The units of length, time, and temperature are $`\sigma _1`$, $`\tau _0=(m_1\sigma _1^2/ϵ)^{1/2}`$, and $`ϵ/k_B`$ in this paper. Details of simulations are given in our earlier paper . We presently performed MD simulations only in three-dimensional space with the systems composed of $`N=N_1+N_2=108`$, $`10^3`$, and $`10^4`$ particles while the density $`\rho =N/L^3=0.8`$ and the composition $`N_1/N=0.5`$ are fixed. The corresponding system linear dimensions are $`L^{N=108}=5.13`$, $`L^{N=10^3}=10.8`$, and $`L^{N=10^4}=23.2`$. Simulations were carried out at $`T=0.772`$, $`0.473`$, $`0.352`$, $`0.306`$ and $`0.267`$ with the time step $`\mathrm{\Delta }t=0.005`$. The PBC was used in all cases. At each temperature, the systems were carefully equilibrated in the canonical condition so that no appreciable aging effect takes place. Data are then taken in the microcanonical condition.

We first calculate the partial static structure factor,

$$S_{ab}(q)=\frac{1}{N}\underset{j=1}{\overset{N_a}{}}\underset{k=1}{\overset{N_b}{}}\mathrm{exp}(i𝒒(𝒓_j^a𝒓_k^b)),$$
(1)   
to investigate whether finite-size effects are detectable in static particle configurations. Here $`𝒓_j^a`$ and $`𝒓_k^b`$ are the positions of the $`j`$-th and $`k`$-th particles in the $`a`$ and $`b`$ components ($`a,b1,2`$) and $`\mathrm{}`$ indicates the ensemble average over different initial configurations. The dimensionless wave number $`q`$ is in units of $`\sigma _1^1`$. In Fig. 1, we plotted $`S_{11}(q)`$ for $`N=108`$, $`10^3`$, and $`10^4`$ at $`T=0.473`$ (a) and $`0.267`$ (b). One can find that $`S_{11}(q)`$ for all cases excellently agrees with each other both in (a) and (b); no systematic size dependence can be detected among them. We examined also $`S_{12}`$ and $`S_{22}`$ and confirmed the same tendency. Our results indicate that finite-size effects are very small or almost negligible in static pair correlations for $`N10^2`$, as that is generally believed.

Let us next consider finite-size effects in dynamical properties. The structure relaxation in glassy materials can be measured by calculating the coherent or the incoherent intermediate scattering functions, $`F(q,t)`$ or $`F_s(q,t)`$. The decay profiles of those two functions tend to coincide at the first peak wave number $`q_m`$ in $`F(q,0)`$. This has been confirmed for the present soft-sphere mixture  and also for a Lennard-Jones binary mixture . Since $`F_s(q,t)`$ can be more accurately determined via MD simulation, we here calculate the incoherent scattering function for the component 1,

$$F_s(q,t)=\frac{1}{N_1}\underset{j=1}{\overset{N_1}{}}\mathrm{exp}[i𝒒\mathrm{\Delta }𝒓_j^1(t)],$$
(2)   
where $`\mathrm{\Delta }𝒓_j^1(t)=𝒓_j^1(t)𝒓_j^1(0)`$ is the displacement vector. Although $`F_s(q,t)`$ decays monotonically in normal liquid states, it exhibits multi-step relaxations in highly supercooled states. This is due to the fact that at lower temperatures the particle motions are highly jammed and thus trapped considerably in effective cages formed by their neighbors. We then defined the $`\alpha `$ relaxation time $`\tau _\alpha `$, which corresponds to a characteristic life time of the effective cage, by $`F_s(q,\tau _\alpha )=e^1`$ at $`q=2\pi `$ for several temperatures. Fig. 2 shows the decay profiles of $`F_s(q=2\pi ,t)`$ obtained for $`N=108`$, $`10^3`$, and $`10^4`$ at $`T=0.473`$ and $`0.267`$. At $`T=0.473`$, we see that the two curves from $`N=10^3`$ and $`10^4`$ entirely coincide, and one from $`N=108`$ is also close to them. The relaxation times thus obtained are $`\tau _\alpha ^{N=108}3.5`$, $`\tau _\alpha ^{N=10^3}=\tau _\alpha ^{N=10^4}2.3`$. However, the situation is different at $`T=0.267`$ where $`F_s(q,t)`$ exhibits two-step relaxations. The faster and the slower parts of the decay are called the fast-$`\beta `$ (thermal) and the $`\alpha `$ relaxations, respectively. We note that the decay profiles for the three systems differ significantly in the $`\alpha `$ regime ($`t10^2`$), whereas they agree well in the fast-$`\beta `$ regime ($`t1`$). Here we determined $`\tau _\alpha ^{N=108}11000`$, $`\tau _\alpha ^{N=10^3}6500`$, and $`\tau _\alpha ^{N=10^4}2000`$ at $`T=0.267`$. Fig. 3 shows the temperature dependence of $`\tau _\alpha ^{N=108}`$, $`\tau _\alpha ^{N=10^3}`$ and $`\tau _\alpha ^{N=10^4}`$. At the highest temperature $`T=0.772`$, $`\tau _\alpha ^{N=108}`$, $`\tau _\alpha ^{N=10^3}`$ and $`\tau _\alpha ^{N=10^4}`$ are exactly equal. However, $`\tau _\alpha ^{N=108}`$ begins to deviate from the others around $`T=0.473`$ at which $`\xi 5`$  is comparable to $`L^{N=108}=5.13`$. The deviation increases with further decreasing temperature, and $`\tau _\alpha ^{N=108}`$ becomes almost one order larger than $`\tau _\alpha ^{N=10^4}`$ for $`T0.352`$. Furthermore, $`\tau _\alpha ^{N=10^3}`$ begins to deviate from $`\tau _\alpha ^{N=10^4}`$ around $`T=0.306`$, at which $`L^{N=10^3}<\xi <L^{N=10^4}`$. We thus suppose that the present finite-size effects are attribute to suppressions of cooperative particle motions due to insufficient system size. The structural relaxation time of smaller systems thus tend to show a stronger (super-Arrhenius) temperature dependence as a result of the finite-size effects. Remembering the fact that the static structure factors are almost identical among those three systems at all temperatures, the origin of this effect may be purely kinetic, or higher order correlations in particle configurations may relevant to this. In their recent paper , Horbach et al. found similar finite-size effects in a model silica glass which is known as a typical strong glass former, in addition to the present soft sphere mixture which is usually classified in fragile glass former.

To understand what happens in microscopic scale, we next visualize individual particle motions in $`N=10^4`$ system at $`T=0.267`$. First, we pick up mobile particles by the condition $`|\mathrm{\Delta }𝒓_j^a(t)|>l_c^a`$ in a time interval $`[t_0,t_0+t]`$, where $`t=0.125\tau _\alpha =250`$, and $`l_c^a`$ is defined separately for the component $`a1,2`$ such that the sum of $`\mathrm{\Delta }𝒓_j^a(t)^2`$ over the mobile particles covers $`66\%`$ of the total sum $`_j^{N_a}\mathrm{\Delta }𝒓_j^a(t)^2`$. Then we define clusters of the mobile particles by connecting $`i`$ and $`j`$ if $`|𝒓_i(t)𝒓_j(0)|<0.3(\sigma _i+\sigma _j)`$ or $`|𝒓_i(0)𝒓_j(t)|<0.3(\sigma _i+\sigma _j)`$ similar to Donati et al. . In Fig. 4, we show spatial distribution of the clusters having the size $`n5`$; those are all chain-like  and have large scale correlations. Although only $`5`$% of the total particles are shown in Fig. 4, the sum of $`\mathrm{\Delta }𝒓_j(t)^2`$ covers approximately $`40`$% of the total $`_{i=1}^N\mathrm{\Delta }𝒓_j(t)^2`$. This clearly indicates that the cooperative motions become dominant in glassy states. To investigate finite-size effects in cooperative motions quantitatively, we here introduce the distribution function,

$$P(n)=\underset{i=1}{\overset{N}{}}{}_{}{}^{}\mathrm{\Delta }𝒓_j(t)^2\delta (nn_i)/\underset{i=1}{\overset{N}{}}{}_{}{}^{}\mathrm{\Delta }𝒓_j(t)^2,$$
(3)   
where the sum runs over mobile particles only. $`n_i`$ is the size of the cluster in which the mobile particle $`i`$ belongs, and thus $`\delta (nn_i)`$ is $`1`$ if $`i`$ is a member of the cluster having the size $`n`$ and $`0`$ if not. The physical meaning of $`P(n)`$ is as follows; clusters having the size $`n`$ contribute $`P(n)`$ to the total squared displacements of the mobile particles. In Fig. 5, we show $`P(n)`$ for $`N=108`$, $`10^3`$, and $`10^4`$ at $`T=0.267`$ at which $`\xi 40`$ obtained for $`N=10^4`$ is even larger than $`L^{N=10^4}=23.2`$. We found that the cooperative motions in $`N=108`$ system are strongly suppressed. By comparing $`N=10^3`$ and $`10^4`$, it is found that larger scale cooperative motions ($`n>10`$) are considerably suppressed also in $`N=10^3`$ system. Infinitely large clusters which percolate the system though the PBC have never been found in all cases. The characteristic cluster size $`\overline{n}=_{i=1}^{\mathrm{}}nP(n)`$ thus obtained is 3.38, 6.11, and 7.73 for $`N=108`$, $`10^3`$, and $`10^4`$, respectively.

In the framework of conventional liquid theories , changes in static particle configurations upon cooling lead to a drastic slowing-down toward the grass transition. The mode coupling theory (MCT)  is the most successful self-consistent approach within this framework. MCT describes onset of glassy slowing down (or slow structural relaxations) in the density-time correlation functions. In the original MCT, however, a sharp Ergodic/Non-ergodic transition is predicted at a temperature $`T_0`$ which is considerably above $`T_g`$. Although such a tendency has been found in colloidal systems in which thermal activation processes are negligible, it has never been observed in real glassy materials. It is thus believed that the MCT has some difficulties for describing the true dynamics of supercooled liquids apparently below $`T_0`$. The main problem is the fact that the original MCT do not take into account the hopping motions of particles, which must be cooperative and thus long ranged as is seen in recent MD simulations . Unfortunately the problem has not yet been overcome in fully self-consistent way because efforts for including thermal activations make the theory more or less ad hoc. It is a interesting fact that the behavior of structural relaxations in our smallest system, in which cooperative hopping motions are highly suppressed, becomes somehow closer to the original MCT prediction.

It is worth mentioning several experimental attempts to find the evidence of the dynamical heterogeneity in glassy materials. One of the most interesting and useful approaches are the recent experiments on glass-forming thin films. The thickness $`d`$ dependence of film properties is the main interest in these studies . The motivations of those studies are quite similar to the present study; they aimed to control the size $`\xi `$ of the cooperative motions by changing $`d`$ and found that the relaxation time and $`T_g`$ considerably depend on $`d`$. They considered that the cooperative particle motions in the direction normal to the film may be truncated near the interface, and this effect become dominant when $`\xi d`$. Thus the system size restriction can enhance the particle motions. Note that this seems to contradict to the finite-size effect in the present MD simulations in which the size restriction suppresses the cooperative particle motions . The mechanism of this discrepancy is still an open question; we naively speculate that the situation is much more complicated in polymer films than in MD simulations. The system size restriction occurs only in one direction normal to the film and other two in-film directions are free in thin films, whereas all directions are equally restricted in the present simulations. We believe that investigating the microscopic relaxation mechanisms in glass-forming (both simple and polymeric) thin films are definitely important.

In summary, we have examined system-size effects in the dynamics of supercooled liquids by MD simulation. We found significant finite-size effects in the structural relaxation at lower temperatures, whereas no such effect is detectable in static pair correlation functions. The cooperative particle motions, which leads to the $`\alpha `$ relaxation in glassy states, are strongly suppressed in smaller systems for temperatures lower than $`T_c`$ at which $`\xi `$ becomes comparable to the system size. The present finite-size effects are regarded as a natural consequence of the dynamical heterogeneity appearing in supercooled liquids. We point out that finite-size effects are significant in the dynamics of highly supercooled liquids, and thus special attention must be payed to investigate true relaxation dynamics in computer simulations particularly in the $`\alpha `$ regime.

We thank Profs. T. Munakata and A. Onuki for helpful discussions. This work is supported by a Grant in Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture of Japan. Calculations have been carried out at the Supercomputer Laboratory, Institute for Chemical Research, Kyoto University and Human Genome Center, Institute of Medical Science, University of Tokyo.