Coarse-graining strategy for modeling effective, highly diffusive fluids with reduced polydispersity: A dynamical study

TL;DR

This paper introduces a coarse-graining method to simplify multi-component fluids into effective monodisperse systems, preserving key dynamic features and enabling easier experimental modeling of complex fluids.

Contribution

A novel coarse-graining strategy that reduces polydispersity in fluid models while maintaining dynamic and structural properties, facilitating experimental realization.

Findings

Effective monodisperse system mimics mixture dynamics

Similar glass transition temperature and critical exponents

Radial distribution function resembles original mixture

Abstract

We present a coarse-graining strategy for reducing the number of particle species in mixtures to achieve a simpler system with higher diffusion while preserving the total particle number and characteristic dynamic features. As a system of application, we chose the bidisperse Lennard-Jones-like mixture discovered by Kob and Andersen [Phys. Rev. Lett. 73, 1376 (1994)] possessing a slow dynamics due to the fluid's multi-component character with its apparently unconventional choice for the pair potential of the type-A-type-B arrangement. We further established in a so-formed coarse-grained and temperature-independent monodisperse system an equilibrium structure with a radial distribution function resembling its mixture counterpart. This one-component system further possesses similar dynamic features like glass transition temperature and critical exponents while subjected to Newtonian mechanics. This strategy may finally lead to the manufacturing of new nanoparticle/colloidal fluids by experimentally modeling only the outcoming effective pair potential(s) and no other macroscopic quantity.

Equations

$$\begin{gathered} U_{\alpha \beta }(R)= \\ =\{\begin{array}{ll}{U}_{\mathrm{LJ}}(R,{ϵ}_{\alpha \beta },{\sigma }_{\alpha \beta })-U_{\mathrm{LJ}}(2.5{\sigma }_{\alpha \beta },{ϵ}_{\alpha \beta },{\sigma }_{\alpha \beta })& ,R\le 2.5{\sigma }_{\alpha \beta }\\ 0& ,\mathrm{otherwise}\end{array} \end{gathered}$$

$${ϵ}_{\mathrm{AA}}=ϵ{ϵ}_{\mathrm{AB}}=1.5ϵ{ϵ}_{\mathrm{BB}}=0.5ϵ$$

$${\sigma }_{\mathrm{AA}}=\sigma {\sigma }_{\mathrm{AB}}=0.8\sigma {\sigma }_{\mathrm{BB}}=0.88\sigma \text{.}$$

$$g_{\mathrm{KA}}\stackrel{!}{=}{g}_{\mathrm{CG}}.$$

$$g(R)=\frac{2}{\rho \cdot (N-1)\cdot 4\pi \cdot R^2}⟨\sum _{i<j\in \{1,\dots ,N\}}\delta (\left|{\mathbf{r}}_j-{\mathbf{r}}_i\right|-R)⟩\text{,}$$

$$u_{i+1}(R)=u_i(R)-{\alpha }_i(R)\ln \left(\frac{g_i(R)}{g_{\mathrm{KA}}(R)}\right).$$

$$\begin{gathered} u_0(R)=u_{\mathrm{HNC}}(R)= \\ =-k_{\mathrm{B}}T\text{ln}(g_{\mathrm{KA}}(R))+k_{\mathrm{B}}T\cdot (g_{\mathrm{KA}}(R)-c_{\mathrm{KA}}(R)-1), \end{gathered}$$

$${\alpha }_i(R)={\alpha }_0\cdot k_{\mathrm{B}}T\cdot \exp (-R^2/(2{\kappa }^2))$$

$$MSD(t)=\frac{1}{N}⟨\sum _i{\left|{\mathbf{r}}_i(t)-{\mathbf{r}}_i(0)\right|}^2⟩$$

$$n(R)=4\pi \rho {\int }_0^{R}d{R}^{\prime }{R}^{\prime 2}\cdot g(R^{\prime })$$

$$g_{xy}(w)$$

$$g_{yz}(w)$$

$$g_{zx}(w)$$

$$g_{\perp }(w,\hat{\mathbf{n}})$$

$$\begin{gathered} F_s(k,t)=⟨{\hat{F}}_s(\mathbf{k},t)⟩\text{with}{\hat{F}}_s(\mathbf{k},t)= \\ =\frac{1}{N}\sum _i\cos (\mathbf{k}\circ \left({\mathbf{r}}_i(t)-{\mathbf{r}}_i(0)\right))\text{.} \end{gathered}$$

$$F_s({\tau }_{\alpha })=1/e.$$

$${\tau }_{\alpha }\propto \exp [(T^{\ast }{)}^{-p}]$$

$${\tau }_{\alpha }={\tau }_0{e}^{\frac{\Delta \mu \cdot z}{k_{\mathrm{B}}T}}$$

$${\tau }_{\alpha }\propto (T^{\ast }-T_c^{\ast }{)}^{-{\gamma }^{\tau }}$$

$$D=\frac{k_{\mathrm{B}}T}{6\pi a\eta }$$

$$\eta =\eta (T,{\tau }_{\alpha })=\frac{{\tau }_{\alpha }{k}_{\mathrm{B}}T}{4\pi a^3}.$$

$$D_{\mathrm{KA}}={\alpha }_{\mathrm{KA}}^D(T_{\mathrm{KA}}^{\ast }-T_{c,\mathrm{KA}}^{\ast }{)}^{{\gamma }_{\mathrm{KA}}^D}\propto D_{\mathrm{CG}}={\alpha }_{\mathrm{CG}}^D(T_{\mathrm{CG}}^{\ast }-T_{c,\mathrm{CG}}^{\ast }{)}^{{\gamma }_{\mathrm{CG}}^D}.$$

$$\eta =\eta ({\tau }_{\alpha })\propto {\tau }_{\alpha }$$

$${\chi }_4(k,t)$$

$${\gamma }^{\tau }$$

$$\lambda =\frac{\Gamma (1-a{)}^2}{\Gamma (1-2a)}=\frac{\Gamma (1+b{)}^2}{\Gamma (1+2b)}.$$

$$\lambda \approx S(k_0)\cdot \frac{k_0{A}^2}{32{\pi }^2\rho },$$

$${\gamma }_{\mathrm{KA}}^{\tau }\approx {\gamma }_{\mathrm{CG}}^{\tau }.$$

$$\frac{\partial }{\partial t}\rho (\vec{r},t)$$

$$\rho (\vec{r},t)={\rho }_0+(\rho \cdot {\chi }_T)\cdot p(\vec{r},t).$$

$$\rho \cdot {\chi }_T\cdot \frac{\partial p}{\partial t}(\vec{r},t)$$

$${\chi }_T=\frac{1}{\rho \cdot k_{\mathrm{B}}T}\left[1+\rho {\int }_0^{\infty }\mathrm{d}R4\pi R^2\left(g(R)-1\right)\right].$$

Full text

I Introduction

In this work we introduce a novel type of coarse-graining procedure for isotropic, polydisperse, Newtonian systems for achieving an effective system with fewer components possessing higher diffusion while the total particle number, the Newtonian mechanics, as well as relevant characteristic dynamic features are preserved. By coarse-graining (CG), we understand the systematic treatment of microscopic details and dynamics from a coarser perspective. Its purpose is to achieve a simpler system as well as having enlarged length or time scales in computer simulations by effectively considering microscopic details. Common CG methods in the literature include force matching schemes,Löwen and Kramposthuber (1993); Ercolessi and Adams (1994); Izvekov and Voth (2005a, b) the relative entropy method,Shell (2008) the conditional reversible work method,Brini, Marcon, and van der Vegt (2011) the inverse Monte Carlo method,Lyubartsev and Laaksonen (1995, 1997) the iterative Boltzmann inversion methodSoper (1996); Müller-Plathe (2002) or hybrid schemes.Rühle and Junghans (2011) These methods are primarily used to develop effective force fields among a small set of coordinates, called reaction coordinatesHeinemann et al. (2014); Heinemann and Klapp (2017) or collective variables, that form projections of microscopic degrees of freedom. A frequent choice of such coordinates is given by the molecules’ center-of-mass positions.Likos (2001)

Here, our presented coarse-graining strategy differs from the previously described standard strategy since we aim at effectively reducing the number of species while not reducing the particles, i.e. the number of collective variables equals the number of microscopic degrees of freedom. We recognized that top-down coarse-graining involving species reduction has been performed already.Ebrahimi, Whittle, and Pellenq (2016); Peters et al. (2018) We, on the contrary, perform coarse-graining in a bottom-up philosophy and thereby achieve a systematic reduction or even a neglect of multi-component character, which is regarded as microscopic detail. By eliminating such detail, one can achieve a higher diffusion and accordingly larger length and time scales in in-silico simulations as a consequence of a smoother energy landscape (see Ref. 17). A higher diffusion might lead to an improved wetting behavior as required in real-world applications like motor oil, hot glue, or inks. Besides better wetting, also the friction can be reduced with higher diffusion as required in machinery, e.g., via a single-component nanoparticle oil adhesion aiming to imitate a multi-component oil adhesion which, from experience, tends to be less tribological.Shahnazar, Bagheri, and Abd Hamid (2016) However, increasing diffusion might also lead to a change of relevant characteristic dynamic features like vitrification and glass transition. Preserving vitrification properties might be even considered as a necessity, e.g., in hot glue and ink-like applications. A common solution for controlling the dynamics, and thus preserving dynamic features, in an effective particle system is to modify the equations of motion with stochastic terms (e.g., Langevin dynamics, dissipative particle dynamics,Hoogerbrugge and Koelman (1992); Koelman and Hoogerbrugge (1993); Español and Warren (1995) …), i.e., density fluctuations stemming from a multi-component character could be emulated by density fluctuations from the reduced species set. We, on the contrary, stay with the Newtonian mechanics also in the CG system and do not intend to make the dynamics more complex, i.e. in our microscopic as well as our CG system we consider all particles to obey the Newtonian equations of motion (besides a weak temperature coupling). In that way, we aim to get all targeted dynamic aspects by having a proper choice for the effective pair potential and take possible structural deviations into account. From the experimental point of view, one can, thus, focus on how to imitate that pair potential.

With respect to the application point of view, there is also a direct demand for reducing polydispersity in melts towards one single component as long as the effective pair potential is not too complex for manufacturing. Potential examples are given by single-component polymer gels,Daniel et al. (2016) which do not require a solvent and thereby avoid unwanted species separations, or sealants,Chronister which do not have to be in a de-mixed state before being used at high temperatures.

To test our method, we chose the Kob-Andersen (KA) mixture,Kob and Andersen (1994, 1995) which is a bidisperse Lennard-Jones mixture with glass-forming ability whose Lennard-Jones potential range and strength parameters do not fulfill the common Lorentz-Berthelot mixing rule.Boda and Henderson (2008) The KA mixture, when supercooled but still being warmer than its glass transition temperature, has a very long-lasting isotropic regime before it crystallizes. Mixtures with similar dynamic features also exist for other models,Pedersen, Schrøder, and Dyre (2010); Kim and Saito (2013) even in different dimensions.Brüning et al. (2008); Flenner and Szamel (2015) By creating a coarse-grained or effective KA system consisting of a one-particle species, as schematically shown in Fig. 1, a higher diffusion and lower relaxation time may follow. These properties, however, act negatively on the glass-forming ability or influence its critical temperature due to less/no entropy of mixing.Pinal (2008) Therefore, the question to be investigated is: Can we yield a monodisperse glass-former resembling in structure but possessing a higher diffusion than the KA mixture? Prior model approaches for monodisperse glass-formers include particle shape anisotropy,Pfleiderer, Milinkovic, and Schilling (2008) pinned particles,Bhowmik, Das, and Karmakar (2016) ultrasoft interactions,Ikeda and Miyazaki (2011) an attractive pair potential with a shoulderDzugutov (1992) or the use of double wellHoang and Odagaki (2008) as well as undulating pair potentials.Elenius, Oppelstrup, and Dzugutov (2010) As we see later, our CG approach also leads to undulating pair potentials.

In order to further specify our polydispersity reduction, we require that the all-particle radial distribution function (RDF), which does not distinguish between the particle type, is set to remain unchanged. During the course of this study, we will slightly mitigate this requirement. As the principal RDF-matching method, we chose the iterative Boltzmann inversion (IBI) methodSoper (1996) which, when converged, leads to an effective model entirely build on an effective pair potential. In contrast to other CG approaches aiming at conserving dynamic features,Izvekov and Voth (2006); Davtyan et al. (2015) the IBI method conserves with the RDF an equilibrium quantity. Thus, it is interesting to scrutinize if, or under which modifications, the here presented CG method is still able to maintain the ability of the so-formed monodisperse system to appear in a glass phase. Within the dynamic context, the stability of the isotropic phase is another aspect we set out to scrutinize since monodisperse systems tend to have shorter relaxation times and crystallize at a fast paceNarumi and Tokuyama (2007); Tokuyama (2011) when supercooled. Even though the particles’ kinetic energies seem quite low for establishing a fast crystallization, a large diffusion, however, could also mean that crystallization barriers could be overcome easier. As we later point out in Sec. IV, we could also imagine ways of moderately reducing the high gain in diffusion in our effective model at the cost of a larger structural mismatch.

The remainder of this article is organized as follows. In Sec. II, we present the theory of the RDF-matching CG approach by describing the transformation from the bidisperse KA mixture towards an effective monodisperse fluid. Then in Sec. III, we propose temperature-independent versions of this approach by introducing the models M0.9 and M0.45. Later on, we present in Sec. IV an investigation of the dynamics covering the mixture and systems created through both temperature-independent models. Finally, we conclude our findings in Sec. V.

II RDF-matching coarse-graining approach

We motivate in the following the RDF-matching coarse-graining (CG) approach using the Kob-Andersen (KA) mixtureKob and Andersen (1994, 1995) in Subsec. II.1 and perform initial assessments on the outcoming effective monodisperse fluid in Subsec. II.2. In these assessments we point out the method’s limitations, which we set to overcome later in Sec. III.

II.1 Tranforming the KA-mixture into a monodisperse fluid

The KA mixture consists of Lennard-Jones particles of mass $m_{p}$ comprising two different species, labeled with A and B, whereas there are four times more A than B particles. It is these few B particles which prevent the other particles, that form the $80\%$ majority, from quickly crystallizing. The set of pair potentials as a function of the inter-particle distance $R$ in that mixture is defined through

[TABLE]

with $\alpha,\beta\in\{A,B\}$, and the Lennard-Jones parameters for the well depth

[TABLE]

and contact distance

[TABLE]

Corresponding potential curves are depicted in Fig. 2 as dashed lines. In accordance with the original KA system, we consider in the effective monodisperse system the same particle density $\rho$ being fixed at $\rho\sigma^{3}=1.2$.

The first trial within our coarse-graining approach consists in identifying at each considered temperature the all-particle RDF (denoted with $g$) of the bidisperse KA system with the RDF of the coarse-grained (CG) or effective monodisperse system, i.e.

[TABLE]

This is equivalent to identifying the mean force profiles.Kirkwood (1935) The definition of the all-particle RDF based on the $N$ particle positions $\{\mathbf{r}_{i}\}$, that is,

[TABLE]

is designed to not distinguish among particle species.

In the next step, we set out to find a monodisperse isotropic system with the same particle density, whose RDF is equivalent to the one of the KA system. Moreover, for computational reasons, it would be desired if the total potential energy in that effective system has only pair-wise contributions. The existence of such underlying pair potential, however, is not proven, but as soon as a system with only pair-wise energy contributions with a pair potential yielding the reference RDF exists, it is unique according to Henderson.Henderson (1974) We further want to point out that such a pair potential not only depends on temperature but also on the particle density. For the interested reader, we recommend the work of LouisLouis (2002) thematizing a discussion concerning density-dependent pair potentials.

In order to obtain the effective pair potential, we chose the iterative Boltzmann inversion (IBI) method.Soper (1996) The corresponding procedure requires that we simulate at each iteration a monodisperse system with an approximate pair potential $u_{i}$ and determine an improved pair potential $u_{i+1}$ by making use of the obtained RDF $g_{i}$ according to the following iterative formula

[TABLE]

As an initial pair potential, we choose the solution of the hypernetted chain approximation (HNC),van Leeuwen, Groeneveld, and de Boer (1959) that is

[TABLE]

with $c_{\rm KA}$ being the direct correlation function (implicitly defined through the Ornstein-Zernike equationOrnstein and Zernike (1914) but explicitly defined in k-spacevan Leeuwen, Groeneveld, and de Boer (1959)) of the KA system. If the functions $\alpha_{i}(R)$ in Eq. (4) are carefully chosen, we can reach convergence after a certain iteration index $j$, i.e. $g_{i\geq j}\approx g_{\rm KA}$. For this purpose, we implemented in accordance with the original literatureSoper (1996) the following mixing function

[TABLE]

with $k_{\rm B}$ being the Boltzmann constant and the remaining parameters were set to $\alpha_{0}=0.05$ and $\kappa=4\sigma$. The latter parameters are fixed for all $i$’s and represent the main mixing parameter as well as the effective range. The mixing parameter was chosen quite low, but its low-valued choice is necessary in contrast to standard IBI applications since we have to dampen condensation and crystallization effects appearing during the iteration in this supercooled fluid. The overall particle number in all systems was set to 35000, which is rather large compared to the original work of Kob and Andersen from 1994,Kob and Andersen (1994) but allows, on the one hand, to correctly track the generally long-range character of the RDFs $g_{i}$ and avoids system size effects on the other. With respect to the pair forces, a cut-off distance of $5\sigma$ turned out to be sufficient. All simulations for each iteration were performed with the Gromacs simulation package 4.6.7 covering $60000$ time steps per iteration, whereas a time step has a length of $0.005$ time units (t.u.) being $\sigma\cdot(m_{p}/\epsilon)^{1/2}$, i.e. each iteration covers a time interval of $300$ t.u. In order to obtain the reference RDF $g_{\rm KA}$ for all considered temperatures, we simulated each bidisperse KA system with even half the time step length and extracted the required positional data from the time interval $[7000~{}t.u.,15000~{}t.u.]$. As initial systems for the bidisperse as well as all initial monodisperse systems, we used random configurations that were quickly equilibrated using the steepest descent method, and in order not to disturb the systems at each iteration too much, we additionally took the last system snapshot of the $i$-th iteration as input for iteration ${i+1}$. Considering the temperature control, we used the Nosé-Hoover thermostat with a time constant of $1$ t.u.

After about 1200 iterations, convergence has been achieved in all our systems in the considered (dimensionless) temperature range covering $T^{*}=k_{\rm B}T/\epsilon=0.025,0.05,0.1,0.15,\dots,1$. At higher temperatures, the bidisperse system is above the freezing point and accordingly out of supercooling.Pedersen, Schrøder, and Dyre (2018) The total iteration number of 1200 seems to be quite large but its high value stems from the low choice of $\alpha_{0}$. The RDFs $g_{i}$ are obtained using the last third of the positional data of each iteration run. Though this only covers $100$ t.u., it leads to the convergence of the IBI scheme as exemplary shown in Fig. 3 for the $g_{i}$ and $u_{i}$ at temperatures being with (a) $T^{*}=0.6$ well above and with (b) $T^{*}=0.4$ below the critical (glass transition) temperature of $T_{c}^{*}\approx 0.435$ (see Ref. 24) within the framework of the mode-coupling theory (MCT).Bengtzelius, Gotze, and Sjolander (1984); Leutheusser (1984)

At both temperatures, the first iteration of the RDF, $g_{1}$ (being the RDF corresponding to $u_{\rm HNC}$) in Figs. 3(a-b) reveals a low-valued first peak (at $R\approx 0.85\sigma$) and a high-valued second peak (at $R\approx 1.05\sigma$) with respect to the specific reference RDF $g_{\rm KA}$ from the bidisperse KA system. In spite of this observed similar shape characteristic for both $g_{1}$ at these two temperatures, we detect during the first 400 iterations a significant difference in the convergence behavior that can be observed by analyzing the corresponding manifold of solutions marked in red. At $T^{*}=0.6$ (Fig. 3(a)) that manifold and the series of manifolds reveals a good-natured convergence, which is reflected by their location between $g_{1}$ and $g_{1200}$, whereas at $T^{*}=0.4$ (Fig. 3(b)) the manifold covering the first 400 iterations is quite broad and shows a strong tendency to over- and underestimate peak heights. Such an extreme overestimation of the first peak is observed since the system shows cavitation effects driven by a negative pressure (see snapshot in Fig. 3(b)). However, these effects become less pronounced the more steps we took into account, and as we later see, the converged RDFs all stem from isotropic phases in which the system cannot overcome critical cavitation within a reasonable simulation time of each iteration and no crystallization is present.

Another aspect we recognize during the IBI convergence is that in contrast to the initially fluctuating convergence behavior of the $g_{i}$ functions at low temperatures, the associated pair potentials $u_{i}$ show a rather well-behaved convergence at all temperatures as exemplarily shown in Figs. 3(c-d). This is a consequence of our low-valued mixing parameter $\alpha_{0}$. Even though the value of $\alpha_{0}$ appears small, it provides a feedback that is strong enough to prevent even supercooled systems close to the glass transition from quickly crystallizing. By scrutinizing the deviation between the converged pair potential $u_{1200}$, and the initial HNC pair potential $u_{\rm HNC}$ in Figs. 3(c-d), we detect a strong alignment at large inter-particle distances ($R\gg\sigma$) which is a known feature of the HNC approximation. However, the HNC approximation failed to correctly predict the effective pair potential $u$ at small inter-particle distances $R$. Therefore, the IBI scheme–or perhaps some promising alternative–turned out to be necessary, even though its implementation and convergence in a large simulation system is involved with respect to computational requirements. We further would like to point out that the values for $u_{i}$ below $R=0.8\sigma$ were extrapolated by cumulatively integrating the pair force along $R$, which we have assumed being equal to the force between A and B particles in that close-contact regime. Regarding the curves’ shape of the converged pair potentials $u_{1200}$ from the IBI scheme (Figs. 3(c-d)), we can identify an alternating behavior along the inter-particle distance $R$; especially a strong repulsive shoulder appearing at around $R=1.5\,\sigma$. To have a comparison of this effective pair potential with pair potentials of the original KA system, we have already included $u_{1200}$ from Figs. 3(c-d) in Fig. 2.

All in all, our IBI convergence behavior seems quite weak compared to the original work of Soper,Soper (1996) who used IBI to simplify water models with a mixing parameter of $\alpha_{0}=1$. Our reference system, however, is situated in the supercooled phase, which is rather unstable at temperatures slightly above the glass transition temperature. Even though the system below the glass transition temperature appears stable for a very long time (up to minor changes affecting the compressibility as pointed out in the appendix), approximate solutions during the IBI scheme might not (see Fig. 3(b)). Another problematic convergence issue in the IBI scheme is the high particle density which is increasing the likelihood of clumping due to a higher number of next neighbors.

II.2 Initial assessments of the RDF-matching approach

We start our approach’s assessments by investigating the pressure, which is in part also negative for our converged simulations as can be seen through the pressure progression along the temperature in Fig. 4. Corresponding ensemble averages (indicated by $\left<\dots\right>$) were ascertained within the framework dubbed as long sampling scheme comprising an evenly sampling in the time interval $[7000~{}t.u.,30000~{}t.u.]$ (forming the standard sampling interval within our investigation).

We recognize that the effective monodisperse system reveals a far lower pressure than the bidisperse system due to less packing fraction as a result of a small contact distance which is similar to the one of the A-B configuration (see Fig. 2). Effective particles thus can get very close before feeling core repulsion which would significantly increase the pressure in a close packed system. As a consequence of this volume effect, particles can cluster more easily and the phase characterized by the given density becomes meta- or unstable up to a point at which cavitation is unstoppable and the pressure turns negative.

Our next objective is the investigation of the influence of (over- and undercritical) temperature on the reference RDF $g_{\rm KA}$ and its associated converged pair potential $u_{1200}$ (see Fig. 5). As expected, these possess a slight progression towards the value $1$ (ideal gas limit) the higher the temperature. But as depicted in the inset, this progression seems to slow down at high temperatures since the system resembles with its steep potentials a hard-sphere mixture before also these pair potentials become softer at higher temperatures. The associated pair potentials (Fig. 5(b)), however, still reveal a high change rate among neighboring temperature sets. It is due to the fact that similar RDFs imply a similar Boltzmann weight factor leading to a nearly direct proportionality between $u$ and $T$ in the canonical ensemble average at overcritical temperatures (see inset). This behavior is in contrast to the temperature dependence of effective pair potentials of particles whose internal degrees of freedom were coarse-grained.Heinemann et al. (2014, 2015) As a result, the pair potential in units of $k_{\rm B}T$ (Subfig. (c)) becomes only slightly weaker with temperature and only values at $R\approx 0.9\sigma$ exhibit stronger weakening effects alongside the temperature progression (see inset). Given that the effective pair potential is a statistical potential, we interpret this weakening as a consequence of the kinetically-driven higher occurrence of the A-A-particle configuration with its Lennard-Jones well depth at $R\approx 1.1\sigma$ compared to the A-B configuration with its well depth at $R\approx 0.9\sigma$ (see Fig. 2). Our results infer that, since the Boltzmann weight factor is not exactly conserved along temperature, this effective monodisperse pair potential behaves similarly to other works in the standard ensembles where CG potentials vary linearly with temperature over a significant temperature range.Farah et al. (2011); Johnson, Head-Gordon, and Louis (2007); Lebold and Noid (2019); Lu and Voth (2011)

We now turn the focus also on the dynamics to assess diffusion-related phenomena when transforming from the bidisperse to the effective monodisperse system. In this regard, we consider the mean squared displacement defined by

[TABLE]

for the bidisperse and the effective monodisperse system in Fig. 6(a)-(b) for the temperature range $T^{*}=0.7,\dots,0.1$. As expected, we detect in the bidisperse system (Fig. 6(a)) a subdiffusive regime (i.e. $MSD\propto t^{a}$ with $a<1$ holds), where the MSD is almost approaching a constant value when cooling the system. This is reflected by the plateau separating ballistic motion characterized by $MSD(t)\propto t^{2}$, and diffusive motion, at which $MSD(t)=6Dt$ holds, with $D$ being the diffusion coefficient.

In the effective monodisperse system (Fig. 6(b)), we detect higher MSD-values which results from a lack of a significant subdiffusive regime, in which particles move by collectively rearranging its vicinity, but also from a lower packing fraction. Particles thus undergo collisions later and accordingly leave the ballistic regime on average later what can be seen when comparing the rectangular marked regions in Fig. 6(a,b). These inconsistencies in dynamics reflect the common failure of structure based coarse-graining.Ghosh and Faller (2007); Szukalo and Noid (2020)

III Temperature-independent RDF-matching approach

In an effort to retain subdiffusive dynamics, we attempt to do the RDF-matching at a reference temperature $T_{\rm ref}$ and intend to use that model –which we denote as $MT_{\rm ref}^{*}$-model with $T_{\rm ref}^{*}=k_{\rm B}T_{\rm ref}/\epsilon$ – for a broad temperature range. Such a non-intrinsic temperature dependent model is even more desired since also the original bidisperse model has per se no temperature dependence. Moreover, maintaining microscopic resolution in a temperature-fixed CG model might still result in a good temperature transferabilityIzvekov, Swanson, and Voth (2008) since the ideal gas part of the entropy is not affected. We chose $T^{*}_{\rm ref}=0.9$ and started all corresponding simulations from an isotropic configuration for which we took the equilibrated configuration of the effective monodisperse system at $T^{*}_{\rm ref}$. In order to asses states with fast crystallization, we also took ensemble averages in a so-called short sampling scheme covering the early time interval $[100~{}t.u.,1100~{}t.u.]$ after starting the simulation.

We observe in a so-formed M0.9-model towards low temperatures for the most part a more dominant pair potential in thermal units compared to those of the effective monodisperse model as depicted in Fig. 7. This is manifested in more attraction at $R\approx 1.05\sigma$ and more repulsion at $R\approx 1.5\sigma$. By more attracting direct neighbors but repelling particles slightly further away, we have induced a subdiffusive motion since particles tend to be trapped. Thus it seems possible to achieve a glass transition temperature and even tune its value towards the one of the original KA model with the proper choice for $T_{\rm ref}$. Other works have proposed temperature dependent renormalization approaches of the pair potential for achieving the desired dynamics which also involved corrections of the pair potential towards more cohesive strength at low temperatures.Xia et al. (2017, 2018, 2019) These even considered an enlargement of the contact distances (i.e. the packing fraction also increases) which additionally dampens diffusion of IBI models by reducing the excess entropy.Rondina, Böhm, and Müller-Plathe (2020); Dyre (2018)

With respect to structure, the M0.9-modeled system possesses at high temperatures the same/similar RDF like the bidisperse system (or equivalently effective monodisperse system) as displayed in Fig. 8(a). At low temperature, e.g., $T^{*}\approx 0.45$, the structural difference becomes more pronounced, what can be seen as an artifact of the increased potential well depth at $R\approx 1.05$ leading to a shift of particle distances towards that value.

By extracting from the RDF the running coordination number, defined as

[TABLE]

we can recognize that the RDF deviations of M0.9 with respect to the effective monodisperse system along $R$ lead each to an overcompensation of the coordination number deficit or surplus so far. That means the coordination number differences along $R$ act as diffusive barriers a particle has to overcome and reflects more hampered motion. As a consequence, MSD results for M0.9 (see Fig. 6(c)) in the short sampling scheme agree with the expected temperature progression of a glass forming fluid (cp. Fig. 6(c) with Fig. 6(a)). However, the M0.9’s MSD temperature progression in the long sampling scheme reveals some abrupt change between $0.55$ and $0.5$ due to crystallization, which has a faster pace at $T^{*}=0.5$ (MSD has no diffusive regime) than at $T^{*}=0.45$. Fortunately, at lower temperatures, no signs of crystallization have been detected which indicates the glassy state. Nonetheless, this fast crystallization is contrary to the KA system, which crystallizes at a very long time scale since the one particle species prevents the other species from crystallizing and these species might even segregate before crystallizing.Toxvaerd et al. (2009) So for industrial applications such as hot glue or inks for inkjet printing, which might make use of strong viscosity change around the glass transition temperature, one could avoid the vicinity around glass transition region, e.g., by rapid cooling or heating.

In accordance with the other models, we have also investigated the pressure for the M0.9 (and the M0.45) model in the prior introduced Fig. 4. Similar to the (intrinsically temperature-dependent) effective monodisperse model, we also recognize low-pressure profiles in both new models. Interestingly, the M0.9-modeled system exhibits no negative pressure for temperatures $T^{*}\gtrsim 0.3$ onwards, irrespective of the used sampling scheme. We, therefore, observe for this system in the temperature region above its glass transition–which can be estimated from the MSD analysis (Fig. 6(c))–no cavitation. Contrary to the short sampling scheme we observe in the long sampling scheme, close but above the glass transition temperature, that the pressure curve reveals for both new models a kink which peaks at about the highest temperature where crystallization has been observed. At that temperature(s) the dynamics is fast enough to establish a phase transition within the long sampling scheme towards crystallization as pointed out in the forthcoming investigation.

In order to evaluate the influence of using this non-intrinsic temperature-dependent approach on structure, we next investigate for all considered systems two-dimensional RDFsAndrienko, Marcon, and Kremer (2006) for the particle projections onto the x-y, y-z as well as z-x plane. These are defined as

[TABLE]

with $\mathbf{\hat{n}}$ being the normal vector for a considered plane. In Fig. 9(a-b), we present resulting curves of Eqs. (9a)–(9c) for the bidisperse and the effective monodisperse system at the temperatures $T^{*}=0.3,0.45,0.5,0.55,0.7$. These 2D-RDF curves show neither a dependence on the choice of orthogonal planes nor on the specific model, and only a weak dependence on temperature in accordance with the small temperature dependence of $g_{\rm KA}$ (see Fig. 5(a)).

As a consequence, the bidisperse and the effective monodisperse systems’ structure is isotropic and there is no sign of an abrupt structural change with temperature which would indicate a system-wide crystallization. For the $MT_{\rm ref}^{*}$ models, on the contrary, we expect some abnormality in shape since the underlying pair potentials are no longer a function of temperature. Corresponding distribution functions for the M0.9 model are shown in Fig. 9 for the long (c) and the short sampling scheme (d) and reveal that there is in fact a different structure with respect to all the curves in Subfigs. (a-b). Moreover, there is an abrupt structural change with temperature appearing in the long sampling scheme, especially at small projected inter-particle distances $w$. The shapes at $T^{*}=0.45,0.5$ clearly display the anisotropy resulting from the crystallization by having a pronounced periodic shape, whereas the system remains purely isotropic at shorter time scales (short sampling scheme used) as depicted in Subfig. (d). For the latter case, no fundamental structural change can be detected, neither as a function of temperature nor with respect to the plane in which we evaluate $g_{\bot}$, i.e. isotropy remains a meta-stable solution in M0.9 and no bond ordering occurs. Nonetheless, we detect in Subfig. (d) a deviation from the curves in Subfigs. (a,b) at $w\approx\sigma$. This results from the fact that both pair potentials only coincide at $T^{*}=0.9$ and reflects the expected structural change when choosing a non-intrinsically temperature-dependent model. With respect to the crystallization observed in Subfig. (c) at $T^{*}=0.45,0.5$, the structural change is even higher at $T^{*}=0.5$ than at $T^{*}=0.45$, which is an artifact stemming from not enough equilibration time given for completing crystallization, and as we see in the forthcoming investigation, reflects that this model is close to its glass transition temperature.

IV Assessing glassy aspects of the systems’ dynamics

IV.1

Diffusion and relaxation dynamics

This section deals with investigations about the dynamics of former-introduced particle systems in anticipation of detecting glassy aspects. For this purpose, we first present in Fig. 10 results for the diffusion coefficient $D$ as a function of temperature $T$ and can recognize that $D$ in the bidisperse system follows a path with a monotonous positive slope while having vanishing values below the critical temperature. The curve of the effective monodisperse system follows also a positive trend but is significantly higher than the bidisperse one and shows no clear inflection point throughout all temperatures. Consequently, a glass phase cannot be realized. If the bidisperse as well as the effective monodisperse system’s particle distribution were fully described by their equal RDF, the systems’ excess entropy, $S_{\rm exc}$ would depend exclusively on it and result in equal diffusion which depends exponentially on $S_{\rm exc}$ (see Refs. 61; 62). The higher packing fraction in the bidisperse system, however, leads to a decrease in the available positions for particles once a number of particle positions is fixed, i.e. the contribution for $S_{\rm exc}$ caused by higher order correlation functions is smaller (with respect to the effective monodisperse system) and leads to a comparatively low $S_{\rm exc}$ and thus a lower $D$.

Contrary to the effective monodisperse system, the $MT_{\rm ref}^{*}$-modeled systems reveal an inflection point that can be shifted depending on the choice of $T_{\rm ref}^{*}$. By closer inspecting their glass transition region (inset of Fig. 10), we denote for ${D\lesssim 2\cdot 10^{-3}\times\sigma\cdot(\epsilon/m_{p})^{1/2}}$ a sudden decay of the diffusion, whereas this is not observed in the short sampling scheme (dashed lines). This phenomenon is consistent with our prior observation of the unusual progression of the MSD in the same temperature region (see Fig. 6(c)) due to the crystallization leading to positional restraints and accordingly vanishing MSD. What is also obvious from Fig. 10 is the fact that $D$ of the $MT_{\rm ref}^{*}$-modeled systems grows at a steeper slope above their glass transition temperatures compared to the effective monodisperse system, which is caused by not increasing potential strength (acting as a diffusive barrier) with temperature in the $MT_{\rm ref}^{*}$-modeled systems.

A faster crystallization can thus be seen as a side effect of the increase in diffusion and might constitute an advantage or disadvantage depending on the application. A solution might consist in avoiding these temperatures (by rapid cooling/heating) or by establishing only a moderate increase in diffusion instead. In particular, one could continuously vary the $T_{\rm ref}$ values at each temperature above the glass transition temperature but this would further alter the structure since very high values for $T_{\rm ref}$ had to be chosen and would influence the critical exponents of the vitrification. Moreover, the manufacturing of this temperature dependence might be challenging. Another approach in preventing crystallization consists in overlaying the pair potential by time and temperature-dependent noise, which does not influence the structure on average. This, however, would alter the equations of motion, which is not our goal. Adding a time-independent noise to the pair potential, however, might represent another solution. The nature of the Newtonian equations of motion would remain untouched but one would further disturb the structure and it would be more difficult to properly simulate the system since a very short simulation time step would be required. Also, manufacturing of particles possessing a fuzzy overlay in their pair-potential might be difficult, if not impossible.

We next scrutinize the relaxation of the systems in anticipation of an unusual behavior in their subdiffusive regimes. For this purpose, we present in Fig. 11 results for the self-intermediate scattering function, that is

[TABLE]

This function characterizes the relaxation of a density mode with wave vector $\mathbf{k}$. We hereby focus on the spherical component and further set the wavenumber $k$ equal to $k_{0}$, being the wavenumber corresponding to the maximum of the structure factor in the bidisperse system at $T^{*}=0.45$, yielding $k_{0}\sigma\approx 7.3$.

As expected, $F_{s}$ in the bidisperse system (Fig. 11(a)) reveals the characteristic two-step relaxation behavior of supercooled fluids, which was already observed for each particle species separately.Kim and Saito (2013) For the effective monodisperse system (Fig. 11(b)), we only observe a simple decline for all temperatures, which coincides with the absence of a subdiffusive behavior in the MSD (see Fig. 6(b)). Fortunately, a two-step relaxation dynamics is recovered within our $MT_{\rm ref}^{*}$-modeled systems as exemplarily shown for the M0.9-modeled system in Fig. 11(c). In coincidence with the MSD, we also detect slightly below $T^{*}\lesssim 0.5$ a sudden strong influence of the equilibrium scheme on $F_{s}$. In fact, we observe in the latter mentioned temperature regime no more alignment between curves stemming from different equilibrium sampling schemes. We interpret this behavior as an artifact of crystallization (leading to non-vanishing relaxation). Moreover, $F_{s}$ at $T^{*}=0.5$ possesses higher values compared to $F_{s}$ at $T^{*}=0.45$, due to a further progressed crystallization process within the long sampling scheme at $T^{*}=0.5$, which is not finished yet at $T^{*}=0.45$. The short sampling scheme, however, leads to a rather “normal” temperature progression of $F_{s}$ (see next paragraph for more details), since the system is not given enough time to crystallize such that it remains supercooled as the bidisperse system.

In order to analyze the temperature progression of $F_{s}$ qualitatively, we next focus on its decay behavior and thereby quantify the $\alpha$-relaxation regime. Additionally, this also leads us to the characterization of the type of glass-former by focusing on the alpha relaxation time $\tau_{\alpha}$ being implicitly defined through

[TABLE]

By determining the $\tau_{\alpha}$ values at different temperatures for each model, we can see in Fig. 12(a) that in the bidisperse system, as well as the M$T^{*}_{\rm ref}$-modeled systems, the following relation:

[TABLE]

approximately holds if choosing $p=2$. This leads to the classification of a fragile glass-former (since $p>1$) and the relaxation behavior is called super-Arrhenius-like. The same relaxation behavior with $p=2$ is observed for the most part in the M$T^{*}_{\rm ref}$-modeled systems. Deviations from super-Arrhenius behavior were observed within the framework of the short sampling scheme and stem from a beginning crystallization in a small temperature interval short above the glass transition temperature (orange and purple dashed lines in Fig. 12(a)). Within this interval, the highest deviations occurred inside (rather than at its the edges) where the temperature is low enough and the mobility still high enough to enable initial crystallization during the short observation time. Results for the long sampling scheme are omitted in that interval due to very large $\tau_{\alpha}$ values which are out of this plotted range. Concerning the effective monodisperse system, we do not observe a super-Arrhenius law; instead, we could fit a curve (red line) with $p=0.2$. Accordingly, we yield a sub-Arrhenius behavior for this non-glass-forming model system. Hereby, the growth in structural relaxation time towards low temperatures is too weak which can be seen as a consequence of the decreasing pair potential strength as shown in Fig. 5(b).

We next provide an explanation for why $\tau_{\alpha}$ is higher in the M0.9(0.45) compared to the effective monodisperse system by focusing on the Adam-Gibbs theoryAdam and Gibbs (1965) stating that

[TABLE]

holds with $\tau_{0}$, $\Delta\mu$ being material constants and $z$ is the number of particles in a hypothetical cooperative spatial region.Dudowicz, Freed, and Douglas (2005a, b, c, 2006) If assuming the same material constants among the M0.9(0.45)-modeled and the effective monodisperse system, the difference has to be primarily in $z$, which is influenced by the next neighbor structure. Since the running coordination number in M0.9(0.45) is more wavy at low temperatures compared to the effective monodisperse one (see Fig. 8(b)), the particles have to overcome barriers which require the local environment to restructure more and explains the higher $z$ and thus the higher $\tau_{\alpha}$ or lower $D$ accordingly.

In order to quantify the glass transition temperature aside from estimates based on values of the diffusion constant, Kob and AndersonKob and Andersen (1995) used the power-law fit (as predicted in the mode-coupling theory (MCT))

[TABLE]

with $T_{c}^{*}$ being the critical temperature and $\gamma^{\tau}$ the corresponding critical exponent. In Fig. 12(b) we present the $\tau_{\alpha}$–$(T^{*}\!-\!T_{c}^{*})$ relationship via a log-log plot for the bidisperse and the M0.9(0.45) model. With respect to the bidisperse system, we observe with $\gamma\approx 2.4$ when approaching $T_{c}^{*}=0.435$ the same power-law behavior as presented in the work of Kim and SaitoKim and Saito (2013) resulting from an $F_{s}$ defined only among the A-type particles (original work:Kob and Andersen (1994, 1995) $\gamma\approx 2.5(2.6)$ for A(B) particles). From that plot, we can see that the power-law of Eq. (14) also holds for our M$T^{*}_{\rm ref}$ models close to their respective glass transition temperature as well and their exponents exhibit similar values, i.e. $\gamma^{\tau}_{\mathrm{KA}}\approx\gamma^{\tau}_{\mathrm{CG}}$. In Appendix A, we provide reasoning for such similarities. With respect to the effective monodisperse system, no critical temperature can be extrapolated, which is why it is not displayed.

IV.2

Stokes-Einstein relation

We next investigate how the diffusion coefficient $D$ depends on the alpha relaxation time $\tau_{\alpha}$. If, for example, we consider one large particle that moves slowly with respect to a viscous medium, the diffusion coefficient $D$ can be expressed via the Stokes-Einstein (SE) relation, that is

[TABLE]

with $a$ and $\eta$ being the Stokes radius (the radius of the considered particle) and the shear viscosity, respectively. Literature devoted to the violation of the SE relation when supercooling a liquid can be found in Refs. 70; 71; 72; 73; 74; 75; 76. With respect to $\eta$, Debye provided an explicit expression, that isGrant (1957)

[TABLE]

This expression simplifies Eq. (15) towards $D\propto\tau_{\alpha}^{-1}$ under the assumption of a fixed (i.e. temperature invariant) Stokes radius $a$. However, this relationship breaks down when supercooling the fluid.Yamamoto and Onuki (1998); Tarjus and Kivelson (1995) Other works investigating glassy systemsKim, Park, and Jung (2016); Park, Kim, and Jung (2015); Jeong et al. (2010); Jung, Garrahan, and Chandler (2004) used a similar but fractional version of this expression, that is $D\propto(\tau_{\alpha})^{-\xi}$ to cover the dynamics in the highly supercooled regime.

In the following, we look for a power-law between $D$ and $\tau_{\alpha}$ for our modeled systems in Fig. 13(a) by using a log-log plot.

We find that $\xi$ seems to differ from the ideal value of $1$ for almost all of our models, especially in the limit of a low $\tau_{\alpha}$ (see inset in Fig 13(a)) as well as high $\tau_{\alpha}$. These limits make the curves appear in a slightly convex shape. With respect to the effective monodisperse system, we detect with $\xi\approx 1.05$ an almost ideal SE behavior, although the glass phase cannot be reached. In the bidisperse system, however, we yield at large $\tau_{\alpha}$ (i.e. $\tau_{\alpha}^{*}\gtrsim 10$) or equivalently low $T$ (i.e. $T^{*}\lesssim 0.6$) with $\xi\approx 0.78$ a value significantly smaller than $1$. This is expected since we already know from the A- as well as the B-particles’ diffusive motion that a break-down of the SE relation towards a fractional one with $\xi\approx 0.8$ for A-particles, and $\xi\approx 0.65$ for B-particles occurs in the deeply supercooled regime (Ref. 24).

A further significant point in this power law investigation is the similarity of the $\xi$-values among the systems formed by the CG models M0.9/M0.45 and the bidisperse model. In particular, the M0.45/M0.9 systems’ $\xi$ coefficients possess in their high temperature regimes with $\xi_{M0.45}(T^{*}>0.45)=1.35$ and $\xi_{M0.9}(T^{*}>0.7)=1.22$, respectively, an only slightly higher as well as similar value compared to the one of bidisperse reference system. For temperatures $T\rightarrow T_{c}$, the relation $\xi<1$ also holds for our M$T_{\rm ref}^{*}$-modeled systems with $\xi_{M0.45}(T^{*}<0.35)=0.79$, $\xi_{M0.9}(T^{*}<0.6)=0.8$ being close to $0.78$ as detected among all particles in the bidisperse system. Upon this observation, we provide in Appendix B a theoretical explanation by showing similarity in the $\gamma^{D}=\epsilon\cdot\gamma^{\tau}$ values, i.e. $\gamma^{D}_{\mathrm{KA}}=\gamma^{D}_{\mathrm{CG}}$ (similarity for $\gamma^{\tau}$’s already shown in Appendix A). Therefore, the following relation approximately holds for the bidisperse KA and the M$T_{\rm ref}^{*}$-modeled (CG) systems ($\alpha^{D}$ values mark scaling parameters).

[TABLE]

In order to investigate more of the entire $D$-$\tau_{\alpha}$-T relationship, we present in Fig. 13(b) a “$D\tau_{\alpha}$” versus “$T$” plot, in which a strong $T$ or $\tau_{\alpha}$ dependence of $D\times\tau_{\alpha}$ at $T^{*}\lesssim 0.6$ for the bidisperse system (black-solid line) indicates the observed SE breakdown. Nonetheless, consistency for that model with the SE relation is found at an intermediate temperature regime ($0.7\lesssim T^{*}\lesssim 1$) directly beneath the melting point beyond which $\xi=1.19$ holds and a linear growth of $D\times\tau_{\alpha}$ with $T$ is found (inset of Fig. 13(b)). Such change in dynamics is also observed in the M$T^{*}_{\rm ref}$ models at temperatures at which the relaxation time is equal to the value observed in the bidisperse system at its melting temperature $T^{*}\approx 1$ (see Fig. 12). The observed linear relationship can be obtained from Eq. (15) if the shear viscosity $\eta$ is related to $\tau_{\alpha}$ as

[TABLE]

and was already observedChen et al. (2006); Kumar et al. (2007) even in fractional form.Xu et al. (2009); Becker, Poole, and Starr (2006) The same dynamic characteristic is also found in our M0.45 and M0.9-modeled systems. Investigations involving these two versions of the SE relation (i.e. Eq. (15) with Eq. (16) or Eq. (18)) can be found in Refs. 88; 89 whereas generalized versions of this relation are investigated in Refs. 87; 90.

IV.3

Statistical analysis of dynamic heterogeneity

More insights into the dynamics can be obtained by analyzing the particle’s activity via excitation events which are characterized by persistence and exchange times defined as follows:Jeong et al. (2010) the persistence time $\tau_{p}$ for a particle $i$ is the minimal waiting time $t_{1}$ to undergo its first excitation such that $\left\|\mathbf{r}_{i}(t_{1})-\mathbf{r}_{i}(0)\right\|\geq d$ holds whereas the exchange time $\tau_{x}$ represents the time between subsequent excitation events, i.e. $\left\|\mathbf{r}_{i}(t_{2}+t_{1})-\mathbf{r}_{i}(t_{1})\right\|\geq d$, $\left\|\mathbf{r}_{i}(t_{3}+t_{2}+t_{1})-\mathbf{r}_{i}(t_{1}+t_{2})\right\|\geq d$, etc. The parameter $d$ represents the critical displacement for which we chose $d=\sigma$ for all considered models. A full alignment between the $\tau_{p}$ and $\tau_{x}$ distributions means that the probability distribution of the time length, an arbitrarily chosen particle at an arbitrarily chosen simulation time needs to undergo a full excitation event, is equal to those of a particle that just went through such event. This seems plausible at high temperatures at which there is a very short structural relaxation time $\tau_{\alpha}$ leading quickly to an equal mobility characteristic throughout the system. If, however, $\tau_{\alpha}$ becomes more dominant in relation to $\left<\tau_{p}\right>:=\mathbb{E}(P(\tau_{p}))$, local structures in configuration and momentum space, that have favored a prior excitation event, might still exist and more easily tend to trigger another event. As a consequence of this behavior, the system can exhibit dynamic heterogeneity (spatial partition of the system in dynamically active and inactive regions).

In Fig. 14, the distributions $P(\tau_{p})$ and $P(\tau_{x})$ and their first moment ratios $\left<\tau_{p}\right>/\left<\tau_{x}\right>$ are depicted.

With respect to the bidisperse system (Fig. 14(a,b)), alignment between $P(\tau_{p})$ and $P(\tau_{x})$ is found at high temperatures, i.e. $T^{*}\gtrsim 0.7$. However, alignment quickly disappears upon cooling leading to persistence times far above excitation times. As a result, the first moment ratio significantly increases with respect to the ideal value of $1$ upon cooling as depicted in Subfig. 14(b). In the effective monodisperse system (Fig. 14(c,d)), the distribution alignment holds even at lower temperatures and thus preventing the ratio $\left<\tau_{p}\right>/\left<\tau_{x}\right>$ to diverge. The M0.9-modeled system (Fig. 14(e,f)), on the contrary, possesses a similar characteristic as those of the bidisperse system as long as we constrain the analysis to the short sampling scheme, at which we have not a progressed crystallization process. In Fig. 14(f) we provide also the curve resulting from the long sampling scheme. From that graph, we can see that from a temperature of $T^{*}\approx 0.6$ onwards, we have full alignment between both equilibrium sampling schemes, which is expected from our upper investigations.

In order to quantify dynamic heterogeneity for the goal of investigating a corresponding critical exponent, we focus in Fig. 15 on the rescaled variance of the $\hat{F}_{s}$ function defined in Eq. (10), which is also called four-point dynamic susceptibility function,Toninelli et al. (2005) and defined through

[TABLE]

Its maximum value is by definition located where the dynamic fluctuation of $\hat{F}_{s}$ is maximized and the positions of the maxima follow the alpha relaxation time,Szamel and Flenner (2006) which we have displayed in Fig. 12. With respect to the bidisperse system (Fig. 15(a)), we observe that the position of the maximum of $\chi_{4}$ is growing by lowering the temperature. This is expected since we already know that a similar result was found in a former investigationKim and Saito (2013) of the dynamics of the A-particles in the bidisperse KA system. In the effective monodisperse system (Fig. 15(b)), we also observe a growth of the maximum of $\chi_{4}$ but the maxima of neighboring curves only slightly differ, meaning that dynamic heterogeneity is not a dominant artifact being observed in this low-fragility ($p=0.2$) fluid while cooling. This is consistent with a prior investigation,Kim and Saito (2013) in which it was observed that more fragility leads to more dynamic heterogeneity. In Subfig. (c), we present $\chi_{4}$-results for the M0.9 model.

Clearly and in a similar fashion to the bidisperse system, we can detect a significant growth of $\chi_{4}$’s maximum corresponding to strongly increasing dynamic heterogeneity upon cooling when considering the short sampling scheme. As expected, we also detect a strong alignment between curves of both sampling schemes at temperatures above the glass transition temperature at which no crystallization occurred within the framework of our investigation (i.e. $T^{*}\gtrsim 0.5$). The growth of $\chi_{4}$ at the alpha relaxation time, when approaching the critical temperature of the respective model, is further depicted in Fig. 15(d), where we can roughly recognize in all models exhibiting a glass transition a $\chi_{4}(\tau_{\alpha})\propto(T^{*}-T^{*}_{c})^{-\nu}$ behavior with a similar critical exponent of about $\nu\approx 1.2$ as observed for the pure A particle subsystem in Ref. 28. With respect to the considered temperature regimes in Fig. 15(d), we restrict to values not too close to the corresponding critical temperature to avoid artifacts, like roughness or crystallization, due to the requirement of a large sampling time for this variance $\chi_{4}$. In order to compare with theoretic predictions and understand why we observe in Fig. 15(d) similar $\nu$ values, we refer to the inhomogeneous MCT, which predicted $\nu=1$ (Refs. 28; 92; 93). The critical exponents of the diffusion coefficient $D$, alpha-relaxation time $\tau_{\alpha}$ and 4-point susceptibility $\chi_{4}$ thus seem to coincide among the coarse-grained M$T_{\rm ref}$-modeled systems and the bidisperse system.

V Conclusions

In this article, we presented a coarse-graining (CG) method for mixtures for achieving on the one hand similar structure and higher particle diffusion, which is still based on the Newtonian equations of motion, while satisfying similar dynamic features on the other. In particular, we aimed at conserving the glass transition temperature, but not all dynamic aspects, which would result in a shift of complexity from the static modeling towards more complexity in the equations of motion. We, however, realized our goals by solely focusing on effective pair potentials and think that such strategy may contribute to real-world and in silico applications as pointed out further below. The fundamental concept is based on a reduction of the number of particle species instead of particles itself. Hereby the well-known bidisperse Kob-Andersen mixtureKob and Andersen (1994) was transformed into an effective monodisperse system at the same particle density and coarse-grained by conserving in a first attempt the all-particle radial distribution function (RDF).

As a result, we yield in comparison to the original Kob-Andersen mixture quite high values for the diffusion coefficient, which is desired as it acts to increase dynamics, however, this coefficient is unfortunately too high at low temperatures (see Fig. 10) leading to a complete absence of the glass transition, which is a dynamic feature aimed to conserve within this method. We did not anticipate an absence, although, we expected a shift in the glass transition temperature (towards lower temperatures) since coarse-graining smoothens the energy landscape and might thus lead to higher diffusion in formerly glassy temperature regions. This high diffusion is also reflected in the weak growth of potential strength in thermal units among effective particle pairs towards lower temperatures (see dashed lines in Fig. 7).

In order to reintroduce the ability to appear in a glassy state, we softened our structural constraint of an identical RDF at all temperatures by introducing two effective models having no intrinsic temperature dependence. These models each correspond to the temperature-dependent effective monodisperse model at a specific dimensionless temperature $T_{\rm ref}^{*}$ and were denoted with M0.45 and M0.9 (corresponding to $T_{\rm ref}^{*}=0.45,0.9$). By using this type of approach it was possible to obtain a stronger growth of potential strength in thermal units towards lower temperature (see solid lines in Fig. 7) and thus achieving the glass phase. The stability of the supercooled region, however, has decreased, such that right above the glass transition temperature a shorter lifetime of the system in the isotropic phase results before crystallization sets in. Nonetheless, the lifetime is still long enough to reach a two-step relaxation behavior for a short amount of time (see Fig. 11(c) displaying the self-intermediate scattering function for a short and a long sampling interval). This faster crystallization results from the increased diffusion which compensates low kinetic energies and thus facilitates the process to overcome potential hurdles for crystallization. However, this is only an issue at temperatures slightly above the glass transition temperature in the supercooled regime, which is of less extension with respect to the temperature range of the supercooled phase of the mixture. The general nature of the glass phase, however, is not affected by this phenomenon. We even observed similar dynamical critical exponents. Furthermore, we identified that the glass transition temperature can be changed by taking the right choice for $T_{\rm ref}^{*}$, which for $T_{\rm ref}^{*}=0.9$ satisfied a critical temperature of about $0.41$ being close to the value $0.435$ of the original work.Kob and Andersen (1994)

At this point, we would also like to mention that the underlying methodology can be used for reducing an A,B,C,…-species mixture towards a monodisperse one as well. More complex cases, for which the outcoming fluid is still a mixture, can be treated in a similar fashion, however, more assumptions might be required and would accordingly constitute another topic of research.

Besides studying physical aspects of this method, our findings suggest a possible use for designing or manufacturing a less or even one-component type of colloidal model fluid out of complex fluids such as paints, lubricants, inks for jet printing, or lacquers which then become disordered solids beneath the same critical (glass transition) temperature as the original mixture. One could also imagine the development of hot glue possessing less starch molecule variety which–due to the higher diffusion at warmer temperatures–better lubricates small pores distributed on surfaces. As a consequence, such glue could provide increased adhesion strength between the materials intended to be glued together. With respect to machinery, we could imagine the further development of improved motor oil of less added adhesions (or less variety among the often not so uniform nanoparticles) which, e.g., better reduces friction exhibited in the gap between piston and cylinder or even circumvents patents on certain adhesions. The one-component fluids resulting from prior applications might be even capable to remain in the liquid state for a broad parameter range, i.e. also in regions where an unwanted eutectic phase separation would occur in the corresponding polydisperse system. However, designing such effective particles might be difficult or even impossible for certain systems since the effective interactions might be too complex for manufacturing the corresponding colloids or nanoparticles.

There are also in silico applications since a faster diffusion leads not only to shorter relaxation times, which are useful to quickly equilibrate a fluid, there is an additional increase in achievable simulation time since the interaction range might get shorter (but not necessarily at high densities). An exemplary system for an in silico application is given through ionic liquids which in their simplest form are described through bidisperse fluids possessing glassy dynamics.Jeong et al. (2011) Hereby, each ionic particle consists of one or more charged atoms or coarse-grained sites.Jeong et al. (2011); Wang et al. (2013) Computer simulations in such a system normally have the drawback that evolving the equations of motion is quite involved due to long-range electrostatic interactions among numerous charges. Established bottom-up CG methods dealing with this issue propose a reduction of the atomic detail towards a CG site descriptionWang et al. (2006); Marrink, de Vries, and Mark (2004) leading to a reduced electrostatic range due to charge compensation effects within and between CG sites. For the case of only one CG site per ionic particle, however, one cannot group more atoms into larger CG sites anymore. For this purpose, there exists already a CG method in the interaction space which foots on a usual force matching method under the constraint that the effective force is created only by particles within a certain range.Izvekov, Swanson, and Voth (2008); Shi, Liu, and Voth (2008) That method is thus more advanced than a simple cut-off scheme. Our method would also coarse-grain the interaction space, but contrary to the latter one, by introducing symmetries in it, which might effetively account for the electrostatic screening effect. In a fully mixed ionic liquid, one can thus also obtain a short range pair potential for the electrostatic part, and accordingly reduce computational costs. However, more important is the fact that one might design a (perhaps manufacturable) monodisperse model imitating glassy aspects of the two-particle ionic counterpart.

Based on the results of this investigation, we could imagine the following directions of study. One could construct an efficient algorithm to approach $T_{\rm ref}$, which is the most important parameter in our work for developing an effective monodisperse model with the same glass transition temperature. We could also imagine studying the possibility of altering or tuning the chemical structure of a polymer to account for structural influences caused by other similar polymers in a mixture. From the methodical point of view, it would be interesting to study the combination of this type of coarse-graining with the classical coarse-graining involving particle reduction. Another idea concerns the reversion of this method, i.e. one could increase the number of species while having a similar structure to achieve lower diffusion of particles. Finally, we would like to point out that the investigation, so far, considers only the inter-particle distance as a variable to describe a pair configuration. Perhaps, one might extend such an investigation towards more complex pair configuration descriptions.Heinemann and Klapp (2017)

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgements.

This work was supported by the BK21 program and the Samsung Science and Technology Foundation under Project Number SSTF-BA1601-11. Furthermore, we would like to thank Sang-Won Park, Gyehyun Park and Inrok Oh for fruitful discussions.

Appendix A Similar critical exponents $\gamma^{\tau}_{\mathrm{KA}}\approx\gamma^{\tau}_{\mathrm{CG}}$ in Eq. (14)

We present in this section an explanation for why the critical exponent $\gamma^{\tau}$ in the scaling law for the alpha relaxation time (see Eq. (14)) reveals similar values for the bidisperse KA and the MTref-modeled CG systems as depicted in Fig. 12(b). As a first step, $\gamma^{\tau}$ can be decomposed into two other coefficients within the framework of the mode coupling theory (MCT), denoted as $a$ and $b$, that isKob and Andersen (1995)

[TABLE]

These coefficients ($a\in(0,0.4]$ - exponent of critical decay; $b\in(0,1]$ - von Schweidler exponent) can be calculated from the exponent parameter $\lambda$ of the bidisperse system through

[TABLE]

A graphical representation of these parameters is depicted in Fig. 16(a) and reveals the fact that for a given $\lambda$ (with $0.5\leq\lambda<1$), the coefficients $a$ and $b$ can be uniquely determined within the considered range. The exponent parameter can be approximately related to the structure factor $S$ viaBengtzelius, Gotze, and Sjolander (1984)

[TABLE]

with $k_{0}$ and $A$ being the wavenumber corresponding to the major peak in $S$ and its peak area. In Fig. 16(b) we provide the structure factor $S$ for the bidisperse KA and the CG model M0.9. At $T^{*}=0.9$, $S_{\rm KA}$ and $S_{\rm M0.9}$ coincide per definition (due to the equal RDF), whereas at $T^{*}=0.55$ (close but not too close to $T^{*}_{c}$), the peaks of $S_{\rm M0.9}$ are significantly higher with respect to $S_{\rm KA}$. Contrary to the relative growth of $S_{\rm M0.9}$ at $k_{0}$, we detect a decline around its minima leading to compensation effects for the peak area $A$ and accordingly to a compensation for the growth of $\lambda_{\rm M0.9}$ with respect to $\lambda_{\rm KA}$. In this respect, we can assume a similar $\lambda$-value resulting in $a_{\rm CG}\approx a_{\rm KA}$ and $b_{\rm CG}\approx b_{\rm KA}$. Following these thoughts, we finally yield

[TABLE]

We also want to point out that the upper investigation can be applied to the M0.45 CG model analogously.

Appendix B Similar critical exponents $\gamma^{D}_{\mathrm{CG}}\approx\gamma^{D}_{\mathrm{KA}}$ in Eq. (17)

We next provide a theoretical explanation for the similar critical exponent $\gamma^{D}$ of the diffusion coefficient (see Eq. (17)) between the bidisperse KA and the CG modeled system M0.9 (M0.45 in a similar way). Following that investigation and using the connection $\xi=\gamma^{D}/\gamma^{\tau}$ can explain why our observation from Sec. IV.2, stating that $\xi_{\mathrm{C}G}\approx\xi_{\mathrm{K}A}$, holds (similarity for $\gamma^{\tau}$ shown in Appendix A). For this purpose, we first consider an arbitrary complex fluid system whose dynamics on a large (“diffusive”) time scale is given through the diffusion equation, that is

[TABLE]

with $\rho(\vec{r},t)$ being the local density (in space and time) in the continuum system or a probability density. Using this time scale, we can investigate applying the linear response approach, which connects the pressure growth with those of the density in a linear fashion, yielding

[TABLE]

Hereby, $\chi_{T}$ denotes the corresponding system-wide susceptibility, which is the isothermal compressibility. As a result of this linear response approach, the diffusion equation (24) turns into

[TABLE]

With this equation we can connect static and dynamic properties since $\chi_{T}$ can be written in terms of the RDF g through the compressibility equation

[TABLE]

The compressibilities of all considered models are displayed in Fig. 16(c) at various temperatures and reveal only slightly positive values above $T^{*}\approx 0.51$. The fluids are thus less compressible. Technically, we would expect a perfect match between curves of the bidisperse and the effective monodisperse model since the compressibility explicitly incorporates the RDF as its only structural information being equal in both models according to the IBI scheme. However, there are a few incidences at which the RDF has become rather unstable when equilibrating after our IBI procedure. We have observed such incidences at $T^{*}=0.3$ and $T^{*}=0.4$ at which the compressibility exhibits large positive values due to minor changes in the RDF as exemplarily displayed in the inset in Fig. 16(c) at $T^{*}=0.3$.

Despite differences close to the glass transition temperature, we can assume a similar compressibility between the M0.9 model (currently referred to as the CG model) and the bidisperse KA model, i.e.

[TABLE]

Consequently, Eq. (26) can be applied to the latter identity leading to

[TABLE]

We hereby want to point out that we are rather interested in a similarity in their logarithmized values as proposed later. In particular, we mean a similarity even though the computer experiment is not providing a good match. With respect to the diffusion coefficients $D_{\mathrm{CG}}$ and $D_{\mathrm{KA}}$, it follows from MCT that at temperatures slightly above the corresponding critical temperature $T_{c}$, the following scaling laws hold:

[TABLE]

Since we chose the temperature $T_{\mathrm{ref}}^{*}$ in the MT${}_{\rm ref}^{*}$ CG model M0.9 to satisfy $T_{c,\mathrm{CG}}^{*}\approx T_{c,\mathrm{KA}}^{*}=T^{*}_{c}$, we can focus on the same temperature range $T^{*}_{\mathrm{CG}}=T^{*}_{\mathrm{KA}}=T^{*}>0.51>T^{*}_{c}$. By taking the logarithm of Eq. (29), which incorporates Eqs. (30a) and (30b), we yield

[TABLE]

By, then, applying the limit $T^{*}\rightarrow[T^{*}_{c}]_{+}$ (and neglecting that the compressibility then slightly turns negative in the M0.9-modeled system close to $T^{*}_{c}$), we yield similar critical exponents in the diffusive limit among the bidisperse (KA) and the coarse-grained system, i.e.

[TABLE]

The latter similarity can be motivated for the our other CG model M0.45 as well if assuming for Eq. (28) the expression $\chi_{T}^{\mathrm{CG}}(T^{*})\stackrel{{\scriptstyle!}}{{\approx}}\chi_{T}^{\mathrm{KA}}(T^{*}-T^{*}_{c,\mathrm{CG}}+T^{*}_{c,\mathrm{KA}})$.