Postponing the dynamical transition density using competing interactions

TL;DR

This study investigates whether adding a longer-ranged repulsive interaction to short-ranged attractive spheres delays the dynamical arrest, testing a mean-field prediction through numerical simulations.

Contribution

The paper provides a comprehensive numerical test of mean-field predictions on delaying dynamical arrest with competing interactions in three-dimensional systems.

Findings

No significant delay observed beyond short-range attraction effects

Simulations do not confirm the mean-field prediction of delayed arrest with added long-range repulsion

Discussion of possible reasons for the discrepancy between theory and simulation

Abstract

Systems of dense spheres interacting through very short-ranged attraction are known from theory, simulations and colloidal experiments to exhibit dynamical reentrance. The liquid state can thus be fluidized to higher densities than otherwise possible with interactions that are purely repulsive or long-ranged attractive. A recent mean-field, infinite-dimensional calculation predicts that the dynamical arrest of the fluid can be further delayed by adding a longer-ranged repulsive contribution to the short-ranged attraction. We examine this proposal by performing extensive numerical simulations in a three-dimensional system. We first find the short-ranged attraction parameters necessary to achieve the densest liquid state, and then explore the parameters space for an additional longer-ranged repulsion that could enhance the effect. In the family of systems studied, no significant (within numerical accuracy) delay of the dynamical arrest is observed beyond what is already achieved by the short-ranged attraction. Possible explanations are discussed.

Equations

$$V_{ij}=\{\begin{array}{ll}\infty & r_{ij}\le {\sigma }_{ij}\\ -U_0& {\sigma }_{ij}<r_{ij}<{\sigma }_{ij}+{\Delta }_{ij}^0\\ U_1& {\sigma }_{ij}+{\Delta }_{ij}^0<r_{ij}<{\sigma }_{ij}+{\Delta }_{ij}^0+x{\Delta }_{ij}^1\\ 0& {\sigma }_{ij}+{\Delta }_{ij}^0+x{\Delta }_{ij}^1<r_{ij}\end{array}$$

$$Q(t)=\frac{1}{N}\sum _{i=0}^N\Theta (a-|r_i(t)-r_i(0)|),$$

Full text

1 Introduction

Particles with short-ranged attractive and long-ranged repulsive (SALR) interactions can form fairly elaborate structures Stradner et al. (2004); Mani et al. (2014); Zhuang and Charbonneau (2016); Ciach et al. (2013); Lindquist et al. (2016); Bollinger and Truskett (2016a, b); Jadrich et al. (2017); Cao et al. (2017). Despite the spherical symmetry of their pair interaction potential, at low temperatures these models assemble into exotic ordered and disordered mesophases, and their structural richness has clear dynamical consequences, even in the disordered regime Zhuang and Charbonneau (2017); Coniglio et al. (2004); Abete et al. (2007). A recent theoretical proposal suggests that certain SALR models exhibit unusual dynamical features in the very dense fluid regime Maimbourg et al. (2018) as well. Maimbourg et al. Maimbourg et al. (2018)’s extension of a high-dimensional treatment of the glass transition Charbonneau et al. (2017); Sellitto and Zamponi (2013a, b) suggests that certain models should display a very pronounced dynamical reentrance upon changing temperature. More precisely, the theoretical analysis suggests that a carefully chosen high-density SALR system that is glassy at low temperature should, upon heating, first melt and then get dynamically arrested once again, all while remaining completely disordered, i.e., without crystallizing.

On its own, such reentrance is not exceptional. The phase behavior of systems with core-softened interactions can exhibit multiple dynamically arrested phases leading to high-order singularities, as first proposed by mode-coupling theory  Fabbian et al. (1999); Bergenholtz and Fuchs (1999); Dawson et al. (2000); Götze and Sperl (2002); qiang Chu et al. (2010), and then verified by both experiments  Pham et al. (2002); Chen et al. (2003); Eckert and Bartsch (2002); Mallamace et al. (2000); Lu et al. (2008); qiang Chu et al. (2010) and numerical simulations Gnan et al. (2014); Sciortino et al. (2003); Zaccarelli et al. (2002); Moreno and Colmenero (2006, 2007); Foffi et al. (2002); Puertas et al. (2002); Charbonneau and Reichman (2007). Dynamical quantities, such as the density-density correlator, then exhibit a logarithmic decay instead of a typical two-step relaxation, and the mean-squared displacement grows sub-diffusively instead of plateauing at intermediate times. A common physical interpretation of this effect is that introducing short-ranged attraction leads to an interplay between two localization mechanisms: caging from the hard-core repulsion and interparticle bonding. As a result liquids with a higher packing fraction than is possible from either mechanism can then be stabilized Sciortino (2002). Adding a supplementary, longer-ranged repulsion is understood as effectively deepening the well created by the short-range attraction, and thus leads to a slightly more efficient packing of neighboring spheres in the liquid state Maimbourg et al. (2018). In the mean-field description, the nonergodicity transition to a glass phase is then pushed to even higher densities although only over a very narrow temperature window Maimbourg et al. (2018). Even though this improvement over a system with purely short-ranged attraction is predicted to be about $3\%$ in the $d\rightarrow\infty$ limit, the effect should be large enough to be numerically distinguishable if it indeed persists in a finite-dimensional system. An additional methodological challenge, however, is that this transition is only a crossover away from the $d\rightarrow\infty$ limit  Charbonneau et al. (2017).

In this article, we attempt to test this prediction in three dimensions via extensive numerical simulations. First, we tune the attraction range of a system of particles interacting via a hard core followed by a short-ranged square-well attraction (SW) to maximize the high-density extension of the liquid phase. We then optimize the interaction parameters of a system with an additional larger-ranged square-shoulder repulsion (SW+SS) in an attempt to push the dynamical arrest to even higher densities. The plan for the rest of this article is as follows. In Sec. 2 we describe the model, the simulation method and the observables of interest. In Sec. 3, we present the simulation results, and we briefly conclude in Sec. 4.

2 Models and Simulation Method

We study $50\%-50\%$ binary (A-B) mixtures of $N=1000$ spherical particles interacting via two potentials: (i) a simple square-well (SW) interaction, and (ii) a SALR square-well plus square-shoulder (SW+SS) interaction. The hard core diameter ratio of the two particle types, $\sigma_{\rm A}/\sigma_{\rm B}=1.2$, with an additive hard-core interaction, i.e., $\sigma_{ij}=(\sigma_{i}+\sigma_{j})/2$ $\forall ij$, is chosen so as to strongly suppress crystallization. The interaction potential can then be generically expressed as

[TABLE]

where $\Delta^{0}_{ij}=\lambda_{0}\sigma_{ij}$ and $U_{0}$ are the width and depth, respectively, of the square well, and $\Delta^{1}_{ij}=\lambda_{1}\sigma_{ij}$ and $U_{1}$ are the corresponding parameters for the square shoulder. Model (i) has $x=0$, while model (ii) has $x=1$, and in both cases temperature, $T$, is expressed in reduced units of $U_{0}$ with Boltzmann constant, $k_{\mathrm{B}}$, set to unity. Hence, model (i) has a single tuning parameter, $\lambda_{0}$, while model (ii) has three: $\lambda_{0}$, $\lambda_{1}$, and $U_{1}$. We consider the dynamics of these systems at constant $N$, volume $V$ and $T$ using a Monte Carlo dynamics that only consists of $N$ single-particle translations per unit time, $t$. These translation are taken along a vector randomly drawn within a three-dimensional cube of side $\delta\ell$, such that the relaxation time is minimum at a packing fraction close to the dynamical transition. The results of such Monte Carlo dynamics are known to be similar to those of other local dynamics in the dense fluid regime which is the regime of interest for this work Berthier (2007); Berthier and Kob (2007); Rutkai and Kristóf (2010).

Equilibration of the initial system is ensured by running Monte Carlo dynamics for at least ten structural relaxation times, $\tau_{\alpha}$, defined from the characteristic decay, $Q(\tau_{\alpha})\equiv e^{-1}$, of the self-part of the particle-scale overlap function

[TABLE]

where $\Theta$ is a step function and $a=0.3{\sigma}_{\mathrm{B}}$ is a microscopic length chosen to be close to the typical particle cage size. This function therefore represents the fraction of particles having moved a distance smaller than $a$ by time $t$.

The equilibrium $Q(t)$ for the liquid is averaged over the trajectory that begins after equilibration. Typical plots for the relaxation time as a function of the packing fraction are shown in Fig. 1 for different temperatures. At fixed $T$, $\tau_{\alpha}(\varphi;T)$ is used to estimate the (avoided) dynamical transition density, $\varphi_{\mathrm{d}}(T)$, by fitting to the critical scaling form, $\tau_{\alpha}(\varphi;T)=A(\varphi_{\mathrm{d}}(T)-\varphi)^{-\gamma}$– see inset of Fig. 1. Because of the presence of activated processes in finite dimensions, this power-law scaling persists for at most a couple of decades Charbonneau et al. (2017), but this range suffices to provide a fairly robust estimate of $\varphi_{\mathrm{d}}$. Estimation of $\varphi_{\rm d}(T)$ provides the dynamical diagram in the $\varphi$-$T$ plane.

3 Results and Discussion

We first tune the interaction range of the simple SW system in order to maximize the depth of the fluid pocket. To the best of our knowledge this optimization had not been previously attempted in simulations. Most previous studies considered models with $\lambda_{0}=0.03$, following the MCT prediction for the existence of an anomalous glassy regime for that interaction range. The dynamical diagrams for different $\lambda_{0}$ around 3% are shown in Fig. 2; the dynamical reentrance of the liquid is clearly visible. The maximum accessible liquid density, $\varphi^{\ast}_{\rm d}$, is however, not attained with $\lambda_{0}=0.03$, but rather with one of $\lambda^{\ast}_{0}\approx 0.019$. Although our result for $\lambda^{\ast}_{0}\approx 0.019$ is in the vicinity of the infinite-dimensional theoretical prediction for this optimization ($\lambda^{\ast}_{0}\approx 0.029$) Maimbourg et al. (2018), it is nonetheless significantly different from it. Because the intricate liquid structure of finite-dimensional systems is neglected in the analytical study, this discrepancy is not particularly surprising. In three dimensions, the nearest-neighbor shell structure is indeed much tighter than what is theoretically assumed. A possible explanation for the discrepancy is therefore that a smaller attraction range suffices in simulations to obtain an energetic stabilization comparable to what is expected in the $d\rightarrow\infty$ limit.

We next explore whether adding a suitably tuned repulsive part to the potential can further delay the dynamical arrest. In this case, three parameters are to be optimized: $\lambda_{0}$, $\lambda_{1}$, and $U_{1}$. We expect the three-dimensional parameter space ($\lambda_{0}$, $\lambda_{1}$, $U_{1}$) for the SW+SS system to be simple with a single minimum (corresponding to the densest liquid configuration) connected to the minimum of the SW system ($U_{1}=0,\lambda_{1}=0$) through a path without large barriers. To nonetheless ensure that our optimization does not miss its target, we explore a wide range of parameter values. We search for an optimum over $\lambda_{0}\in(0.010,0.060)$, $\lambda_{1}\in(0.5,5.0)$, and $U_{1}\in(0.0,0.40)$ by gridding the parameter space, and compute $\varphi_{\mathrm{d}}$ for a few temperatures around the reentrance regime in the dynamical diagram for each grid point to estimate $\varphi_{\mathrm{d}}^{\ast}$. From this scheme we identified the set of parameters that pushes the dynamical transition furthest as $\lambda^{\ast}_{0}=0.019\pm 0.004$, $U^{\ast}_{1}\leq 0.10$, and $0.8\leq\lambda^{\ast}_{1}\leq 3.0$. All directions away from this optimum lead to lower or comparable values of $\varphi_{\rm d}^{\ast}$. The optimal parameters identified are in qualitative agreement with the theoretical prediction that the repulsion should be much weaker and longer ranged than the attraction and that the attraction range does not markedly broaden in going from a SW to a SW+SS model. In our case, however, the attraction range barely changes, while the theoretical prediction has $\lambda^{\ast}_{0}$ increase from $0.029\to 0.054$. Here again, the tightness of the finite-dimensional neighbor shell is likely at play.

A more significant difference is that while the dynamics (and thus $\varphi_{\mathrm{d}}$) is somewhat sensitive to the repulsion strength $U_{1}$, its dependence on the repulsion range $\lambda_{1}$ is much weaker over the parameter window considered. The dynamics of all models with $\lambda_{1}\in[0.5,2.5]$ indeed roughly coincides (see Fig. 3 inset). Once more, the pronounced shell structure of three-dimensional dense fluid is likely at play. While in the $d\rightarrow\infty$ limit, $g(r)=\exp(-\beta V_{ij}(r))$, and hence self-solvation can be strongly impacted by the interaction potential; in three dimensions, the influence of the hard core-repulsion is felt much more strongly (see Fig. 3). As a result, adding a weak repulsive contribution to the interaction potential results in a much weaker dynamical effect. For a repulsion range that falls within the intershell depletion regime, no notable effect on the dynamics are thus observed.

Given the relative insensitivity of the optimization to $\lambda_{1}$, we can concentrate on the two-dimensional parameter space, $\lambda_{0}-U_{1}$, for $\varphi_{\mathrm{d}}^{\ast}$. Fig. 4 shows the maximum fluid packing fraction $\varphi_{\rm d}^{\ast}$ in the space of $U_{0}$ and $\lambda_{1}$ where $\lambda_{1}=2.5$. Interestingly, the optimization landscape is relatively flat along $U_{1}$. The SW optimum is therefore connected by a fairly soft mode to the SW+SS optimum. Along $\lambda_{1}$, by contrast, $\varphi_{\rm d}^{\ast}$ changes much more rapidly. This landscape projection is therefore consistent with the above discussion.

The resulting dynamical diagrams for the SW and the SW+SS optima are compared in Fig. 5. The results show that the corresponding $\varphi_{\rm d}(T)$ values are not significantly different (within numerical uncertainty) from one another. If any enhancement of the reentrance pocket is present, it is therefore much smaller than the $3\%$, predicted by the $d=\infty$ calculation.

4 Conclusion

Motivated by a recent mean-field prediction that the dynamically sluggish fluid regime for models with SALR interactions can be pushed to higher densities than for models with purely short-ranged attraction, we have performed extensive Monte Carlo simulations of a family of SW and SW+SS models. Our exploration of model parameters did not identify (within numerical uncertainty) any SALR model that pushes the dynamical transition significantly beyond the densest packing achievable by only short-ranged attraction. We did, however, identify a branch of parameters over which the optimum extends. This nontrivial feature could be a finite dimensional echo of the $d\to\infty$ prediction. The theoretical prediction that further tuning the interaction potential could engender additional (smaller) gains in $\varphi_{\mathrm{d}}^{\ast}$ Maimbourg et al. (2018) is nevertheless unlikely to be verifiable in three-dimensional systems.

Data associated with this work are available from the Duke Digital Repository at “will be added".

Acknowledgements.

This paper is dedicated to the late Bob Behringer, who has always been warm, wise and supportive to this junior colleague (PC). He will be sorely missed. We acknowledge funding from the Simons Foundation (Grant # 454937 to PC) and computer time of Duke Compute Cluster (DCC) and Extreme Science and Engineering Discovery Environment (XSEDE), supported by National Science Foundation grant number ACI-1548562.