Multiple character of non-monotonic size-dependence for relaxation dynamics in polymer-particle and binary mixtures

Adding plasticizers is a well-known procedure to reduce the glass transition temperature in polymers. It has been recently shown that this effect shows a non-monotonic dependence on the size of additive molecules [The Journal of Chemical Physics 150 (2019) 024903]. In this work, we demonstrate that, as the size of the additive molecules is changed at fixed concentration, multiple extrema emerge in the dependence of the system's relaxation time on the size ratio. The effect occurs on all relevant length scales including single monomer dynamics, decay of Rouse modes and relaxation of the chain's end-to-end vector. A qualitatively similar trend is found within mode-coupling theoretical results for a binary hard-sphere (HS) mixture. An interpretation of the effect in terms of local packing efficiency and coupling between the dynamics of minority and majority species is provided.


Introduction
To identify the physics underlying the glass formation in amorphous materials has been one of the challenging problems in the field of soft condensed matter physics over the last several decades, and still demands further studies [1]. During the glass transition in a material, a sharp increase of relaxation times occurs in a narrow temperature interval, while the static glassy structure undergoes no distinctive change. This feature has been investigated by a myriad of works in experiments, simulations and theoretical studies [2][3][4][5][6][7][8].
An interesting feature of some glass forming systems is the emergence of a nonmonotonic dependence of the relaxation dynamics (see [9] and references therein) by changing only one control parameter. For instance, the mode-coupling theory revealed a reentrant glass-transition line in colloidal systems with short-ranged attractive interactions [10]. This prediction was later confirmed by experiments on colloid-polymer mixtures [11][12][13][14][15][16] as well as computer simulations [8]. There have been also other studies in the literature with regard to the reentrant glass transition with different mechanisms such as quantum effects [17] and confinement [18][19][20].
In recent works [16,21], distinct scenarios are reported for the glass transition in asymmetric binary (colloidal) mixtures. In addition to the possibility of gel formation [16], the larger component can undergo a glass transition while the smaller one remains ergodic. Another path is marked by simultaneous loss of ergodicity in both components and is therefore termed as a double glass transition. While the latter is marked by a strong coupling between the dynamics of the two species, the single glass state is ascribed to a decoupling of self-dynamics from collective dynamics at short length scales [16,21].
The above studies explore the available states in binary mixtures by varying the packing fraction at a fixed size ratio. There have been few studies focusing on the effects of size ratio in bi-disperse systems of fixed composition [22,23]. A central result of these studies is the observation of a non-monotonic dependence of the glass transition on the size-ratio of the constituents. The systems investigated include binary mixture of soft spheres [22] and polymer-additive blends [23].
A non-monotonic size-dependence of glass transition is also predicted by a modecoupling theory (MCT) for binary hard-sphere mixtures [24,25] and by generalizations of nonlinear Langevin equations [26,27]. A detailed survey of MCT calculations reveals the existence of multiple extrema in relaxation time versus concentration [25]. This means that, at a fixed size ratio, varying concentration of the minority component leads to multiple passages between the liquid-like and glassy states (liquid-glass-liquid-glass-...). A question of interest here is whether this multiple re-entrance can occur at a constant concentration but by varying the size-ratio of the two components, and whether it can rationalize the non-monotonic effect that was recently found for the plasticization of polymers by addition of small molecules [23].
In the current work, we explore this issue by combining MD simulations of a polymer/additive system and MCT analysis of a colloidal hard sphere mixture. In both cases, the concentration of additive/small component is kept constant but the size-ratio small/large particles is varied. It will be shown that, indeed, additional extrema exist in the dependence of structural relaxation time on size-ratio. The fact that qualitatively this phenomenon occurs both in entangled polymers containing small spherical additives and in binary mixture of colloidal articles calls for an interpretation in terms of local competing packing effects. A discussion of this issue is provided.

Polymer-additive simulation model
The polymer is modeled as a linear chain of spherical particles connected via a finite extensible nonlinear elastic (FENE) force. The corresponding interaction potential reads [28,29], where k = 30 pp /σ 2 pp = 30 is the strength factor and R 0 = 1.5 the breaking limit of covalent bonds. Throughout this paper, the index p refers to polymeric beads (monomers) and s to (small) spherical particles added to the system. In addition to the FENE-potential which acts only between adjacent monomers along the chains' backbone, all particle pairs interact via a Lennard-Jones (LJ) potential, In Eq. (2), α, β ∈ {p, s} and r αβ is a short hand notation for the distance between a particle, i, of type α and another one, j, of type β: r αβ = |r i,α − r j,β |. The LJ potential is truncated at a cutoff radius of r c,αβ = 2 × 2 1/6 σ αβ . The monomer diameter, σ pp , is kept constant throughout the simulation and defines the unit of length (a convenient way to achieve this is to set σ pp ≡ 1). The size disparity is defined via δ = σ ss /σ pp . The parameter σ sp is chosen as the arithmetic mean, σ sp = 0.5(σ pp + σ ss ). For simplicity, the energy scale of the LJ potential is set to unity regardless of the particle type, i.e. pp = ss = sp = 1. The mass of an additive particle is set to be equal to that of a monomer, m s = m p = 1.
Temperature is measured in units of pp /k B with the Boltzmann constant k B . All other quantities are given as a combination of the above described units. The unit of time, for example, is given by τ LJ = (mσ 2 pp / pp ) 1/2 and that of pressure is pp /σ 3 pp . All quantities addressed here are expressed in this set of reduced LJ units. The equations of motion are integrated using the Velocity-Verlet algorithm with a time step of δt = 0.003. All the simulations are performed by the open source molecular dynamics simulator LAMMPS [30].
The number of monomers per chain is N p = 10. This choice is motivated by the need to equilibrate the system at all investigated temperatures, while keeping the computational cost acceptable. The total number of particles is N = 4000. Simulations are first prepared in the N pT -ensemble at a pressure of p = 1. We then switch to the N V T -ensemble for dynamics measurements. This way, we avoid undesirable effects on dynamics which could originate from the fluctuations of the simulation cell size [31]. The number-concentration of additive molecules is kept constant at 20%. In contrast to this, the size disparity is varied from 0.3 to 1 by adjusting the size of the additive particles, σ ss .

MCT for binary colloid mixtures
The central quantity for MCT is the matrix of time-dependent density autocorrelation are the microscopic density fluctuations to wave vector q associated to the N α particles of type α. Its equal-time limit is just the static structure factor, Φ(q, 0) = S(q). Following a projection-operator formalism [32], one derives a time-evolution equation of the form (in matrix notation) The J αβ (q) = q 2 (k B T /m α )x α δ αβ set the time scales for the short-time motion, chosen here to fix the unit of time through the thermal velocity of the bigger particles. For simplicity, we set both masses equal (as in the simulation); the long-time dynamics predicted by MCT does not crucially depend on this choice as long as the mass ratio is not too extreme [33]. The memory kernel M (q, t) describes the slow structuralrelaxation dynamics, and is approximated by MCT as a bilinear functional of the density-correlation functions [24,32], where p = |q − k|, and with vertices that are determined entirely by the static structure of the system. In a common approximation that simplifies triplet correlations, The static structure factors are assumed to be known; in the present work, we calculate them from the Percus-Yevick approximation for hard-sphere mixtures [34], to obtain a closed fully analytic theory.
The long-time limits of the density correlation functions, F (q) = lim t→∞ Φ(q, t) are called the glass form factors. They are identically zero in the liquid, and non-vanishing in the ideal glass. From Eqs. (3) and (4) one obtains an algebraic equation is the MCT memory kernel evaluated with the form factors. The bifurcation points of this equation identify the ideal glass-transition points of MCT and can be found via a simple iteration scheme.
The long-time diffusion coefficients D α of the individual particles of type α are obtained from the respective tagged-particle density correlation functions, ] is the single-particle microscopic density fluctuation. In the limit q → 0, one obtains an equation for the mean-squared displacement (MSD) of the tagged particle, with the thermal velocity v th,α = k B T /m α . The analysis of the long-time behavior of this equation yields which is finite as long as the tagged-particle memory kernelm s α (t) decays to zero for long times. This is always the case in the liquid, but if the host system is in an ideal-glass state, weakly coupled tagged particles can still retain a finite diffusivity in this glass.

Particle size effects in polymer-additive systems
In this section, we analyze the relaxation dynamics in the polymer/additive system for a fixed concentration of the additive component with a focus on the occurrence of multiple extrema as the size of added particles (minority species) is varied.

Overall Dynamics
The main observation for the case of the polymer-additive system is illustrated in Fig. 1, where both a polymer-specific quantity -the autocorrelation function of the chain's end-to-end vector -and a quantity which stands for single particle dynamics -meansquared displacement (MSD) of a monomer, averaged over all monomers in the systemare shown. The autocorrelation function of the end-to-end vector is obtained from EtE- denotes the vector connecting the first and the last monomers of the i-th chain at time t and the sum runs over all N c chains. The all-monomer mean-squared displacement is determined via with decreasing the size of the additive from δ = 1 to δ ≈ 0.85; then it speeds up upon lowering the additive size down to δ ≈ 0.5; and finally, it strongly slows down again at even smaller δ. The dependence of typical relaxation times on δ thus shows multiple extrema (insets of Fig. 1). The fact that both the end-to-end vector ACF and the monomeric MSD show the same qualitative features, suggests that the observed effect is not a consequence of a coupling between chain connectivity and dynamics of additive molecules. Rather, it seems to reflect the influence of added molecules on the local packing of monomers, regardless of their connectivity along the chains' backbone.
As will be shown later below, this interpretation is further corroborated by the fact that a binary mixture of hard sphere particles is predicted within MCT to also exhibit qualitatively the same features in terms of size disparity, δ.
To provide further quantitative evidence for the existence of multiple extrema, we have performed an extensive set of simulations and have sampled the parameter range for size disparity more densely. Figure 2 shows the thus obtained results on the structural relaxation times versus δ. The data are extracted both from the relaxation dynamics of chains' end-to-end vector and from the all-monomer mean-squared displacements.
It is noteworthy that, in contrast to the main (global) minimum, which is associated with a relatively large change in the relaxation time and is thus rather easy to detect, resolving the existence of additional extrema requires significant computational effort to reach a sufficiently high signal to noise ratio. As can be inferred from the error bars in the relaxation data, this goal has been achieved convincingly in our simulations. The effect is demonstrated for two temperatures close to T c and does not qualitatively depend on temperature in this regime. This suggests that it is an underlying change in the glass-transition point T c itself as a function of δ that explains the observed nonmonotonic dependence of the relaxation time on the size ratio. In principle, one has to be careful in distinguishing possible pre-asymptotic effects on the dynamics (that may have a different size-ratio dependence) from changes of the glass-transition point [33], but our MCT calculations discussed below support this interpretation.
Figure 3a underlines these findings further by showing polymer diffusion coefficient versus size disparity. Here, diffusion coefficient is extracted from the long time limit of the chains' center of mass displacements, D p = lim t→∞ (r cm (t) − r cm (0) 2 /6t. Noteworthy, as shown in Fig.3b, the mobility of additive molecules increases rapidly and in a monotonic way as size disparity decreases (δ = σ ss /σ pp → 0). In contrast to this, the enhancement of polymer dynamics is non-monotonic ( Fig. 2 and Fig. 3a). A possible way to rationalize this observation is as follows. As the size of additive molecules decreases, they can explore the small available free volume while at the same time experiencing fewer collisions with monomer beads. The coupling between the motion of polymer chains and that of additive molecules thus decreases with smaller size ratios. This competition of increasing mobility and decreasing coupling strength is one of the main reasons for the occurrence of the observed non-monotonic effect [23].

Dynamics of Rouse modes
The above data show that multiple extrema occur both in a polymer specific quantity (decay of the end-to-end vector) and in single particle dynamics (MSD). Here we analyze this issue in more details with a focus on Rouse modes [35], defined via X p (t) = 1 Np Np n=1 r n (t) cos( (n−0.5)pπ Np ), p = 0, 1, . . . , N p − 1 [36]. In this definition, r n is the position of the n-th monomer of a chain. The parameter p plays qualitatively the role of a wave number in Fourier series. Rouse modes thus provide information on chain relaxation dynamics at different length scales, from the entire chain (p = 0) down to a single monomer (p = N p − 1).
Taking advantage of this flexibility of the Rouse modes in focusing on different length scales, we investigate their relaxation behavior as a function of size ratio. Before doing this, and as a benchmark of our simulations, we first check how the   relaxation time of Rouse modes depends on the mode index, p. Within the Rouse model, where monomers are subject to random forces acting on them and are at the same time connected via harmonic springs to their neighboring beads, one expects that τ Rouse ∝ 1/p 2 [35,36]. As shown in Fig. 4, for small p (which correspond to large lengths close to the size of a chain), this prediction provides a good approximation to the behavior of Rouse modes within our model system. Deviations from Rouse prediction at large p (small length scales) reflect excluded-volume interactions, which correlate particle motion at short lengths and thus slow down the structural relaxation at these length scales as compared to the predictions of ideal Rouse model, which treats monomers as point particles.
Results on the dependence of relaxation times of Rouse modes on size ratio δ are shown in Fig. 5. The non-monotonic effect and the existence of multiple extrema are clearly visible in this figure for all Rouse modes investigated. This observation supports the view that connectivity along the chain's backbone has no major effect on the coupling between dynamics of additive molecules and polymer chains. Interestingly, also the size ratios for which different extrema occur, seem to be largely independent of the mode index p. This suggests that the phenomenon can be rationalized letting aside polymer specific features. Following this idea, we provide below a study of the same effect in a binary mixture of hard spheres. It is shown that local packing effects are responsible for the non-monotonic effect and the occurrence of multiple extrema.

Binary hard-sphere (HS) mixtures
The simplest model to capture excluded-volume interactions and specifically the ensuing mixing effects of particles with different sizes, is the binary hard-sphere mixture. We now turn to a discussion of MCT predictions for this system to support the interpretation given above for the dynamics of the polymer/additive system. The overall scenario predicted by MCT for the binary HS model has been discussed at length [24,25], so we focus on the specific aspects that are relevant for the present discussion.
In the binary HS system, the role of temperature is replaced by that of the overall packing fraction ϕ; in the monodisperse HS system with the Percus-Yevick structure factor, MCT predicts a glass-transition point ϕ c ≈ 0.5159, thus we choose packing fractions ϕ = 0.514 and 0.515 as proxies for the temperatures close to T c that have been studied in the simulation. Essentially, as one enters the asymptotic regime of small |ϕ − ϕ c |, the effects to be discussed below do not change qualitatively, but become more pronounced as one approaches ϕ c . For simplicity, we keep the overall packing fraction of the binary HS system constant as well as the number concentration of small particles, and vary the size ratio δ as the single control parameter.  To quantify the change in relaxation of the collective dynamics, we show in Fig. 6a the structural relaxation times τ extracted from the density correlation function of the large particles at a wave number q = 7, corresponding to the main peak of the static structure factor in the monodisperse system. As is evident already for ϕ = 0.514 and more pronounced for ϕ = 0.515, this relaxation time displays a maximum around δ = 0.8, and a minimum around δ = 0.4, for the chosen composition. There is an additional slight maximum around δ = 1.05. This is in qualitative agreement with the non-monotonic variation discussed in Fig. 2, where a maximum around δ = 0.85 and a minimum around δ = 0.6 are seen; also there, a slight maximum is visible around δ = 1. We attribute the small differences in δ values to the differences in the models: In particular, for the HS mixture that case δ = 1 implies that there is no difference between large and small particles, while for the polymer/additive system even at δ = 1 there are structural and dynamical differences between the monomers of a polymer and the same-sized additive particle.
Also shown in Fig. 6b are the critical packing fractions ϕ c found for the various size ratios. The ϕ c -versus-δ curve shows a clear maximum around δ = 0.4, and a minimum around δ = 0.8. These two extrema thus explain the pronounced corresponding minima and maxima of the relaxation times -as the overall packing fraction of particles is kept constant, one moves further away from, respectively, closer to the glass transition point. These extrema are also robust against changing the small-particle concentration within a certain range, as can be seen by noticing that the ϕ c -versus-δ curves shown in Fig. 6b for two different concentrations display qualitatively the same variation. The additional extrema seen in the relaxation-time curve do not have corresponding extrema in the variation of the glass-transition point. They are thus more subtle features of the asymptotic approach to the glass transition.
The diffusion coefficients predicted by MCT for the binary HS mixture, Fig. 7, likewise capture qualitative trends that are visible in the results for the polymer/additive system, Fig. 3. At small δ, a decoupling of the diffusivity of the smaller particles from the diffusivity of the rest of the system becomes apparent; this is a precursor to the transition to a partially arrested glass that MCT predicts in this range of size ratios. As a consequence, the small-particle diffusivity for small δ has a weaker dependence on ϕ than the large-particle diffusivity, because it is not influenced by the vicinity of ϕ c . This effect is qualitatively similar to the observed decoupling in the polymeric system [23]. It, in particular, explains the weak temperature dependence of the additive diffusivity seen in the simulations for small δ. At the same time, the large-particle diffusivity shows a maximum around δ = 0.4 as expected from the collective relaxation dynamics, and this maximum qualitatively agrees with the one occurring around δ = 0.6 in the polymer/additive system.
At large δ, differences appear that are again due to the simplification made in the binary HS model: here, necessarily the ratio of diffusion coefficients approaches unity as δ → 1; this is not the case in the polymer/additive system, where equalsized free particles remain faster diffusing than their monomeric counterparts in the polymer chains. This originates from chain connectivity, which slows down the dynamics of a monomer as compared to a single non-bonded particle. It is, therefore, quite interesting that, despite this asymmetry, all the qualitative features remain similar between polymer-additive system on the one hand and binary HS mixture on the other hand. This again underlines the dominance of local packing effects in the phenomenon addressed in this work.

Summary and outlook
In this study, the effect of size disparity on the glass transition is explored for two different systems, namely polymer-additive blends and hard-sphere mixtures. The first system is studied by molecular dynamics simulations and the latter via modecoupling theory. In the polymeric systems, monomers act as larger components which are connected to each other forming linear chains. The smaller component (additive molecule) is modeled as spherical LJ particles. Relaxation time is found to exhibit multiple extrema upon a variation of the additive/monomer size-ratio at constant additive's number concentration. The trend is robust and shows itself in structural relaxation at various length scales ranging from a monomer diameter to the chain's end-to-end vector. Interestingly, a similar trend is found also in mode coupling theoretical calculations of a binary hard sphere mixture. These observations suggest an interpretation of the additive's effect in terms of local packing effects and the related coupling between additive molecules/minority species and polymer/majority species. Essentially, the small additives are most effective in speeding up the relaxation dynamics when they are sufficiently different in size from the majority species, but also not too small to merely diffuse inside the structure set by the larger particles. This leads to a strong acceleration in relaxation around δ ≈ 0.4 (δ ≈ 0.6 in the polymeric system). If the additives are just slightly smaller than the majority species (δ ≈ 0.8), their smaller size increases free volume, but does not yet disrupt the overall structure, so that in the HS system, a slightly smaller overall packing fraction is sufficient to achieve the same slow relaxation as in the monodisperse system.
The interpretation of the additive effect given here in terms of local packing effects between monomers and additive molecules is also borne out in the application of MCT to polymer melts [37,38]. It would be interesting to see if the inclusion of specific intrachain interactions in that theory will be able to improve the qualitative agreement that we found here also for large δ. A further interesting question concerns ways to control this non-monotonic effect via additional mechanisms. One possibility would be to tune interaction energies in order to make the interactions between the two species more repulsive.