Cluster glasses of ultrasoft particles

We present molecular dynamics (MD) simulations results for dense fluids of ultrasoft, fully-penetrable particles. These are a binary mixture and a polydisperse system of particles interacting via the generalized exponential model, which is known to yield cluster crystal phases for the corresponding monodisperse systems. Because of the dispersity in the particle size, the systems investigated in this work do not crystallize and form disordered cluster phases. The clustering transition appears as a smooth crossover to a regime in which particles are mostly located in clusters, isolated particles being infrequent. The analysis of the internal cluster structure reveals microsegregation of the big and small particles, with a strong homo-coordination in the binary mixture. Upon further lowering the temperature below the clustering transition, the motion of the clusters' centers-of-mass slows down dramatically, giving way to a cluster glass transition. In the cluster glass, the diffusivities remain finite and display an activated temperature dependence, indicating that relaxation in the cluster glass occurs via particle hopping in a nearly arrested matrix of clusters. Finally we discuss the influence of the microscopic dynamics on the transport properties by comparing the MD results with Monte Carlo simulations.


I. INTRODUCTION
A common strategy to facilitate the study of the physical properties of complex macromolecular aggregates is to coarse-grain the intramolecular degrees of freedom [1,2]. By using standard statistical mechanical tools, it is possible to represent each macromolecule as a single point particle and to obtain an effective pair potential accounting for the free energy of interaction between two such macromolecules. For several macromolecular architectures, including linear chains [3,4], rings [5,6], stars [7][8][9], dendrimers [10][11][12] or microgels [13,14], the so-derived effective potentials are "ultrasoft", i.e., the centers of mass of the macromolecules can coincide at a modest energetic cost (of order k B T ) without violating excluded volume interactions between monomers.
Ultrasoft particles exhibit a more complex phase behavior than that of hard ones. This originates from the interplay between entropy, which governs the structural properties of hard-sphere solutions, and energetic contributions arising from the fully-penetrable character of the ultrasoft particles. The topology of the phase diagram in the temperature-density plane can be classified in two classes: reentrant or monotonic behavior. The behavior of the Fourier components of the ultrasoft bounded potential provides a necessary and sufficient condition for observing one or the other class [15,16]. Thus, if all the Fourier components are positive the crystallization lines are reentrant. A complex cascade of crystalline phases is found on increasing the density and these depend on the specific ultrasoft potential [17][18][19][20]. On the * Email: daniele.coslovich@univ-montp2.fr contrary, if the Fourier transform of the potential shows negative values the crystallization lines are monotonic in the temperature-density plane. The corresponding crystalline phases of this class are non-conventional: the ultrasoft particles form a cluster crystal [21]. This crystal consists of clusters of particles located in the nodes of the lattice. Another particular feature of this phase is that the lattice constant is density independent. A direct consequence of this property is that the cluster population is directly proportional to the density of the fluid [16]. The generalized exponential model (GEM) [21] is a well-known example of ultrasoft bounded potential leading to the two former scenarios, depending on the specific parameters of the model (see Section II). The cluster crystal scenario of the GEM has been confirmed in a series of computational investigations [16,[21][22][23][24][25]. Detailed investigations of the phase behavior have revealed an extremely complex map of cluster crystal structures [26]. Some of these works [22][23][24][25] have focused on the dynamic aspects of cluster crystals, revealing interesting properties. The stability of the lattice, which has a noninteger average cluster population, is maintained by incessant hopping of all the particles between the clusters. In contrast to the usual observation in glass-forming liquids [27,28], a comparison between Newtonian, Brownian and Monte Carlo (MC) simulations reveals a significant role of the microscopic dynamics on the long-time dynamics [25]. In particular the hopping dynamics in Brownian and Monte Carlo simulations is characterized by short-range jumps, and the long-range, highly directional jumps found in Newtonian dynamics are strongly suppressed. Recent non-equilibrium simulations of cluster crystals reveal novel features for their rheological response [29].
It is worth mentioning that potentials of the clustercrystal class have been derived for some macromolecules [5,30]. However, though a certain degree of clustering was found in concentrate solutions of such macromolecules, cluster crystal formation has not been observed yet [5,31]. The reason is that, because of increasing many-body effects at high densities, the obtained effective potentials were no longer valid at the densities for which crystallization was predicted. Whether one can design specific macromolecules that can form cluster crystal phases at high concentrations is still an open issue.
As usual in colloidal systems crystallization may be avoided in some situations, leading to amorphous states of ultrasoft particles. This allows to investigate the formation of glassy states for sufficiently high densities or low temperatures. To the best of our knowledge, all investigations on this issue have been performed in systems of ultrasoft particles showing the reentrant scenario for crystallization. The counterpart of this phenomenon in the amorphous case is a reentrant glass transition. This feature was predicted by the Mode Coupling Theory in a coarse-grained model of star polymers [32], though arrested states could not be investigated because of crystallization. A natural way of preventing crystallization is to introduce dispersity in the particle size. This allowed to investigate glassy states of Hertzian spheres [33], confirming the presence of a reentrant glass transition. Interestingly, it has been recently shown that Gaussian spheres can form glasses at high density even in the absence of size dispersity [34] and that these have a strong mean-field character [35].
In this work we aim to get further insight in the dynamics of ultrasoft particles by investigating the glassy behavior for the class of cluster crystal-forming systems. We present extensive computer simulations of a binary mixture and a polydisperse system of ultrasoft particles interacting through the generalized exponential potential. We investigate structural and dynamic properties around and below the clustering transition. The introduced size dispersity is sufficient to prevent crystallization and to produce a disordered arrangement of the clusters' centers-of-mass. We observe the signatures of an incoming glass transition, leading to a state that we denote as "cluster glass", akin to the dynamically arrested states observed in colloidal systems with competing interactions [36][37][38]. Note, however, that the GEM potential is purely repulsive and have no minima, i.e., for the system investigated in this work clustering is found in the absence of attractive interactions. A detailed analysis of the dynamics reveals a progressive arrest of the clusters' centers-of-mass on decreasing temperature, with the relaxation of the particles taking place by hopping between the nearly arrested clusters. Finally we provide indications that the role of the microscopic dynamics (Newtonian or stochastic) on the long-time dynamics may be less important in cluster glasses than in cluster crystals.
The article is organized as follows. In Section II we de-scribe the investigated model and give simulation details. In Section III we present the simulation results and discuss structural and thermodynamic properties (IIIA), as well as dynamic properties (IIIB). We discuss the dependence of our results on thermal history and microscopic dynamics in subsection IIIC. Conclusions are give in Section IV.

II. METHODS
We investigate the dynamics of ultrasoft fullypenetrable particles by means of extensive computer simulations of a generalized exponential model of index n (GEM-n) [21]. In this model the interactions between the two particles are given by the bounded potential For exponents n > 2 the Fourier transform of the potential has negative components, and hence the monodisperse system can form cluster crystal phases. In this work we focus on two different values of the exponent n in Eq. (1): a binary mixture with n = 4 and a polydisperse model with n = 8. Binary mixture.-The system is composed of a mixture of two species {1, 2} of particles interacting via a GEM-4 potential. The potential is cut and quadratically shifted at a distance r c = 2σ αβ where α, β ∈ {1, 2}. The ratio between the particles' diameters is σ 22 /σ 11 = 1.3, and the cross-diameter σ 12 = (σ 11 + σ 22 )/2 = σ = 1 is set as the unit of length. In this work we focus on the case of an equimolar mixture, i.e., with the same number of particles for both species 1 and 2.
Polydisperse model.-The system is composed of N polydisperse particles of diameter σ i . Polydispersity is introduced by means of a flat distribution of the variable σ i . The distribution is centered at σ = 1 and the minimum and maximum values are σ min = 0.826 and σ max = 1.164, respectively. Particles interact through the pair potential in Eq. (1) with n = 8. The interaction is cut off at a distance r c = 1.5σ ij , where σ ij = (σ i + σ j )/2, σ min ≤ σ ij ≤ σ max , and i, j ∈ {1, .., N }. In order to discriminate between particles of different sizes, we introduce three subpopulations of particles labelled by α = 1, 2, and 3. If we sort the particles by increasing value of the diameter σ i , we say that the particle i belongs to the species α if i ∈ [1 + (α − 1)∆, α∆], with ∆ = N/3. Thus α increases with increasing average size of the particles.
In both the mixture and the polydisperse system we use a common energy scale αβ = 1 and particle mass m = 1. The particles are placed in a cubic box with periodic boundary conditions. The static and dynamic properties of these models were investigated by means of Molecular Dynamics (MD) and MC simulations, performed over a wide range of temperatures T and of densities ρ = N/V , with N the number of particles and V the volume of the simulation box. Namely we used Production MD runs were performed in the microcanonical ensemble (Newtonian MD). Newton's equations of motion were integrated by means of the velocity Verlet algorithm [39]. For the polydisperse model the time step δt ranged from 10 −3 at high temperature to 10 −2 at low temperature. For the binary mixture the time step δt was 0.02 independent of temperature and density. These values of δt allowed to keep the degree of energy conservation, determined from the ratio between the root mean square deviations of the total and potential energy, to less than 2% at all the investigated state points. Thermalization in the equilibration runs was achieved by periodic velocity rescaling in the polydisperse model and by means of the Berendsen thermostat [40] using a time constant t T = δt/0.1 in the binary mixture.
For comparison with the MD results we also performed Monte Carlo (MC) simulations in the canonical ensemble. Propagation of the particles was implemented according to the standard Metropolis algorithm [39]. The trial moves performed during the MC simulations involved random particle displacements, generated over a cube of side 0.1. We observed that the acceptance ratio varied between a 80% (high T ) and a 50% (low T ).
To test the reliability of our results different thermalization criteria were adopted and compared. The key quantity to assess equilibration was the root mean squared displacement (RMSD) evaluated at the end of the simulation. In the case of the binary mixture, the duration t eq of the equilibration run at any temperature was such that R(t eq ) > 8. The duration of corresponding production runs, t prod , was typically four times longer than t eq , so that R(t prod ) ≈ 2R(t eq ). At the end of the production run, the system was cooled to a lower temperature, and the procedure was reiterated. To ensure that all the single-particle degrees of freedom, i.e., those corresponding to both small and large particles, were equilibrated, we also implemented an analogous equilibration criterion based on the partial RMSD where the sum is restricted to particles of species α. A small, systematic dependence on the target value R α (t eq ) was observed at the lowest temperatures and will be discussed in Sec. III C.
In the case of the polydisperse model, equilibration and production runs were such that R(t eq ) and R(t prod ) typically exceeded three interparticle diameters. In the runs at the lowest investigated temperatures, which covered about 10 8 time steps, the target RMSD reached values of about one interparticle diameter. Even in these runs at very low temperature, however, no drift in the time dependence of the potential energy and pressure during the production runs could be observed. In addition to the same gradual cooling protocol used for the binary mixture, we also implemented an infinite-quench rate protocol, whereby the initial configuration was always prepared by placing the particles randomly in the simulation box and then performing an infinite-rate quench to the target temperature T , which was subsequently equilibrated by monitoring the value of the potential energy and pressure as a function of time. We made sure that no drift was observed during the production runs. We found that the static and dynamic properties of the polydisperse model did not depend appreciably on the quenching protocol employed.
Temperature T , time t, distance r, wave vector q and density ρ are given respectively in units of /k B (with k B the Boltzmann constant), σ(m/ ) 1/2 , σ, σ −1 and σ −3 . Unless otherwise specified, in the following the presented results will correspond to the Newtonian MD simulations. The comparison between MD and MC results will be discussed at the end of Section III.

III. RESULTS AND DISCUSSION
Unless specified otherwise, in the following we will present simulation results for two selected densities: ρ = 4.0 for the binary mixture and ρ = 5.0 for the polydisperse model. As we will show below, the dynamic properties of the models investigated here exhibit, at sufficiently low temperatures, the ρ/T -scaling found in the cluster crystal phase of the monodisperse system [22,23]. The scaling becomes effective for ρ ≥ 3.0 in the binary mixture and for ρ ≥ 5.0 in the polydisperse model (see below). The selected isochores are therefore representative of the behavior in the high-density scaling regime.

A. Structure and thermodynamics
We start our discussion by analyzing the static pair correlations. Figures 1 and 2 show the temperature variation of the partial radial distribution functions g αβ (r) for the binary mixture and the polydisperse model, respectively. Both models follow a very similar structural evolution along the selected isochores: a prominent peak located around r ≈ 0 builds up as the temperature decreases, indicating an increased interpenetration of the particles. At the same time, the first minimum in g αβ (r) becomes deeper, suggesting the formation of well-defined clusters of typical maximum size ≈ 0.7. This value can be read off from the positions of the first minima of the radial distribution functions. The ability of the particles to interpenetrate stems of course from the ultrasoft and bounded character of the interaction potential, Eq. (1), but formation of clusters in purely repulsive models is  a non-trivial collective phenomenon, that arises only at sufficiently high density [16]. A close inspection of the g αβ (r) in the binary mixture shows that the peaks around r ≈ 0 are significantly higher for like correlations. This suggests that in this latter model clusters are mostly populated by particles of the same species-a phenomenon that can be described as "chemical segregation". The data for the polydisperse model indicate a similar tendency towards homo-coordination, although the effect is significantly weaker than in the binary mixture.
To characterize these effects more precisely, we performed a simple cluster analysis. A particle belongs to a given cluster if its distance to at least one of the other particles of that cluster is smaller than a preselected cut off r cut . When using a fixed value of r cut , it is sometimes difficult to unambiguously identify the cluster to which a particle belongs, due to the continuous flow of particles from cluster to cluster [25], in particular at high temperature. In practice, however, we do not observe any major artifacts for the systems at hand. In particular, our data indicate that merging of neighboring clusters [25] is not a serious issue in our models at sufficiently low temperature (see Fig. 4). In view of the typical widths of the first peaks of g αβ (r) at low temperature, a reasonable choice of the cut-off distance is r cut ≈ 0.4 for the binary mixture and r cut ≈ 0.35 for the polydisperse model. Small changes of this parameter (±30%) did not affect qualitatively our analysis. Let us first analyze the distribution P (n cl ) of cluster population numbers n cl . The temperature variation of P (n cl ) is shown in Fig. 3 for the two models. At high temperature the distribution consists of a sharp peak around n cl = 1 (isolated particles), and a featureless background of more populated clusters. Thus in the high-temperature fluid phase the system is mostly "dissociated", though particles may temporarily overlap due to collisions. An analysis of the cluster lifetimes (not shown) confirms this picture. Moreover, at such high temperatures large values of n cl , leading to the observed background in P (n cl ), may result from the fluid structure being spatially more homogeneous, making the definition of clusters meaningless. As the temperature is lowered, the height of the peak at n cl = 1 decreases, larger clusters become more frequent and P (n cl ) approaches a welldefined ultimate profile. In the polydisperse model at ρ = 5.0, a clear peak is visible around n cl ≈ 8. In the binary mixture at ρ = 4.0 the distribution is clearly bimodal, indicating the formation of two distinct types of clusters with average populations of n cl ≈ 6 and 10, respectively. We have performed a similar analysis (not shown) at the other simulated densities. As expected, we find that the "preferred" values of n cl are density dependent, namely they tend to increase with increasing ρ.
A closer inspection of the chemical composition of the clusters reveals that the double peak structure of P (n cl ) in the binary mixture reflects the same chemical segregation indicated by the radial distribution functions: "small" clusters (n cl ≈ 6) are mostly formed by small particles (species 1), whereas "big" clusters (n cl = 10) are mostly formed by big particles (species 2). This effect is nicely illustrated in Fig. 4. The population number of each cluster is first decomposed as n cl = n (1) cl +n (2) cl , where n (1) cl and n (2) cl indicate the number of particles of species 1 and 2 composing the cluster, respectively. Then, the two-dimensional map of the distribution of clusters with "chemical composition" (n (1) cl , n (2) cl ) is constructed. The distributions shown in the figure correspond to thermodynamic states located in the temperature regime where strong clustering is observed. For the polydisperse model, the analysis is performed using (n (1) cl , n cl ) and averaging over all possible values of n (2) cl . In the binary mixture, the tendency towards chemical segregation is rather evident. It is worth mentioning that at these low temperatures "merged clusters" are extremely rare. Cluster merging in the binary mixture would lead to the appearance of chemical compositions such as (n (1) cl , n (2) cl ) = (12, 0) or (0, 20) (i.e., union of preferred clusters composed of small and big particles, respectively), or (6, 10) (i.e., union of two different preferred clusters). However, the distribution shown in Fig. 4a has no contributions around (n (1) cl , n (2) cl ) = (6, 10) and (0, 20). Only a few instances of (12,0) and (12,1) clusters are visible, the latter arising from two (6,0)-clusters merged by one big particle, but their fraction is negligible (< 10 −5 ). Thus, we conclude that the binary mixture of GEM-4 particles self- assemble at low temperatures into a binary mixture of clusters. The results for the polydisperse model indicate an anti-correlation between subpopulations of small and large particles, but the effect of chemical segregation is much weaker than in the binary mixture (see Fig. 4b).
Thus, the clusters in the polydisperse model remain intrinsically polydisperse in character.
In the following, the crossover from the dissociated fluid phase at high temperature to the cluster-dominated regime at low temperature will be identified, without too much rigor, as a "clustering transition". Though we have no evidence of the possible thermodynamic nature of this phenomenon, inspection of the excess spe- dT indicates a clear operational definition of the transition. The total potential energy per particle, U (T ), and C V (T ) are shown in Figs. 5 and 6 for the binary mixture and polydisperse system, respectively. Both systems display a smooth decay in U (T ) on decreasing T , and a broad peak in C V (T ). The maximum of this peak, at temperature T * , marks the clustering transition. T * takes value 0.52 and 1.3 for the binary mixture and polydisperse model, respectively. An analogous definition of the clustering transition has been used in the study of the microphase separation of a model with short-range attraction and long-range repulsion [41,42]. For comparison, we have included in Figs. 5-a and 6-a the respective results for the monodisperse systems (GEM-4 and GEM-8, respectively). These systems were prepared in their equilibrium cluster crystal phases and subsequently heated up to high temperatures. The density was ρ = 5.0 and ρ = 4.097 for the GEM-8 and GEM-4 models, respectively. The latter value of the density was chosen so as to match precisely the effective packing of the binary mixture within an "effective one-component" description [43]. This adjustment was necessary to achieve the expected full collapse of the potential energies of the two GEM-4 models at high temperatures. The abrupt change in the potential energy of the monodisperse systems indicates the transition from the fluid to the cluster crystal phase, which is a true thermodynamic transition. In the following we will denote the corresponding melting temperature as T m . Interestingly, T * of the polydisperse system is very close to the corresponding T m of the monodisperse model, whereas the difference is more pronounced for the binary mixture.
To correlate the variation of the former thermodynamic quantities to the formation of clusters in the binary Concomitantly, the fractions of the preferred clusters increase, but in a rather smooth fashion. Therefore, although a clear signature of the clustering transition can be found in the thermodynamic properties, one should bear in mind that it may well represent a crossover and not a sharp, thermodynamic transition. Around T * , clusters and isolated particles are indeed in continuous and dynamic exchange-a sort of "chemical equilibrium" picture. Figs. 7 and 8 show typical snapshots of the binary mixture and polydisperse system above and below the clustering transition T * . The densities (ρ = 4.0 and 5.0 for the mixture and polydisperse system, respectively) are the same for which structural and thermodynamic observables have been presented in previous figures. The snapshots confirm visually the general structural features discussed above. As temperature decreases well below the clustering temperature T * , the fraction of isolated particles becomes negligible ( 10%) and the system enters in a "cluster phase", i.e., a regime where almost all particles form tightly bounded clusters and only rare jumps allow particles to be transferred from one cluster to another. What is the structure of the clusters of such a cluster phase? To address this point, we analyze the radial distribution functions g CM αβ (r) of the clusters' centers of mass (CM). The index α indicates here the "chemical composition" of the clusters. Namely, α-clusters are defined as those composed by a majority of particles of the species α. Isolated particles are excluded from this analysis. In Figs. 9 and 10 we show, for temperatures below the clustering transition, the functions g CM αβ (r) of the binary mixture and polydisperse model, respectively. We find that the cluster structure changes only slightly upon cooling, without any major transformation. In particular, the cluster structure of the polydisperse model is rather insensitive to temperature variation and clearly amorphous. The temperature dependence of the radial distribution functions of the binary mixture is somewhat stronger than that of the polydisperse model. Moreover these functions exhibit a more marked splitting of the second peak, especially for 2-2 correlations. The typical distances between neighboring clusters can be read off from the location of the first peaks and reflect the different "chemical" composition: clusters composed of big particles tend to have larger distances to their first neighbors. In Figs. 11 and 12 we include the corresponding data for the static structure factors of the clusters' CMs, S CM αβ (k). The maxima of the different structure factors take moderate values, S CM αα (k) ∼ 2 for correlations between same species and S CM αβ (k) ∼ 1 for distinct species. No signature of Bragg peaks, which would be present in cluster crystals, is found, confirming the amorphous character of the arrangement of the clusters' CMs. The significant anticorrelation, S CM αβ (k) ∼ −0.5, of distinct species in the low-k regime suggests a certain segregation between clusters of big and small particles.
In summary, from the former analysis we conclude that the structure of the clusters' CMs below the clustering transition is essentially amorphous, at least down to the lowest simulated temperature, but differences are visible between the two investigated models (binary mixture and polydisperse system). We will study the dynamic character of these amorphous cluster phases (whether they consist of fluids or glasses of clusters) in the next subsec- tion.

B. Dynamics
We now turn our attention to the dynamics of the two investigated models. In doing so, we will better characterize the nature of the cluster phases identified in the previous section, and we will highlight the differences and similarities with respect to the dynamics of ultrasoft particles in cluster crystals [22,25].
Let us start with the analysis of the temperature dependences of the diffusion coefficients. These have been extracted from the long time limit of R 2 α (t) /6t, where the partial mean square displacements are defined as: The index i runs over all particles of species α. In Fig. 13 we show the species-dependent diffusion coefficients D α of the binary mixture (ρ = 4.0, panel (a)) and the total diffusion coefficient D of the polydisperse model (ρ = 5.0, panel (b)). The very high temperature regime is omitted for clarity and will be shown in Fig. 23. An Arrhenius representation is used to highlight the development of slow dynamics upon decreasing the temperature. A first portion of the data, covering temperatures for which the diffusion coefficients are between ≈ 1 and ≈ 10 −2 , can be reasonably well described by a Vogel-Fulcher-Tamman (VFT) law The fitted values of the strength parameters A (α) depend on both species and models under consideration. We observe that the clustering temperature T * lies within this first portion of data. Thus no dynamic signature of the thermodynamically defined clustering transition can be evidenced from this representation of the data (see also Fig. 23). At lower temperatures the diffusion coefficients undergo a crossover to a milder temperature dependence, which can be well described by an Arrhenius law The activation energies are species-dependent and tend to be higher the bigger the particles. We interpret this "fragile-to-strong" crossover (i.e., from VFT to Arrhenius temperature dependence) as a signature of a change in the microscopic transport mechanisms, from the high-T collective dynamics to activated, single-particle hopping taking place in a nearly frozen cluster structure. We re-  [22]. Empty symbols are data for the polydisperse system of this work. The thick dashed lines indicate Arrhenius-like behavior in the cluster crystal and cluster glass of the monodisperse and polydisperse systems respectively. mark that this crossover takes place below T * , in a regime where clustering is nearly complete and the fraction of isolated particles P 1 is typically lower than 15-20%. Figure 14 shows the diffusivities of the polydisperse system versus the variable ρ/T at the different investigated densities. At low temperatures, data for ρ ≥ 5.0 collapse to a common Arrhenius law, i.e., with an activation energy E act ∝ ρ. Low-T data for ρ = 3.0 do not collapse though still are close and parallel to the data for ρ ≥ 5.0. This scaling behavior of the diffusivities in the cluster phase of polydisperse GEM-8 particles is analogous to that observed in the fcc cluster crystal phase of the monodisperse counterpart. The data of the fcc cluster crystal (taken from Ref. [22]) are included in Fig. 14 for comparison. In analogy with the observation for the cluster crystals (see discussion in Ref. [22]), the Arrhenius temperature dependence suggests that the particles perform hopping dynamics on the sites of an almost frozen matrix of clusters. Though the same qualitative ρ/Tscaling is found, quantitative differences are observed. In particular the activation free energy in amorphous cluster phases is lower than in the cluster crystal. This can be tentatively assigned to structural disorder leading to a larger number of available pathways ("entropic" contribution) or to a smoother energy landscape probed by the single particles in the amorphous cluster phase, in contrast to the ordered structure of local minima separated by high barriers in the cluster crystal [23].
The ρ/T -scaling of the polydisperse system breaks down at sufficiently high temperatures, i.e., above the fragile-to-strong crossover. This suggests a more complex, cooperative transport mechanism above this crossover, in analogy with the mildly supercooled regime of glass-forming liquids. We will show below that indeed the systems investigated here exhibit, in this regime, characteristic dynamic features of that scenario. Finally, it is worth mentioning that the data for the polydisperse system and its monodisperse counterpart show a perfect overlap in the high-temperature fluid phase, i.e, above the crystallization and clustering transition for the monodisperse and polydisperse system respectively. Thus, the effects of the size dispersity have no relevance for the dynamics at such high temperatures. Consistently with the observations for the potential energy and specific heat (Figs. 5 and 6) the diffusivities of the polydisperse system display a smooth variation around T * , in contrast to those of the monodisperse system around T m . All the previous qualitative observations are also found for the ρ/T -scaling (obeyed for ρ ≥ 3.0) of the diffusivities in the binary mixture (not shown).
To characterize the dynamics in more detail, we study the temperature evolution of both incoherent and coherent intermediate scattering functions (Figs. 15, and 16). In the case of the binary mixture, we report separately the scattering functions for small and large particles, whereas an average over all particles is performed for the polydisperse model. The selected wave vectors reflect the positions of the respective first peaks in the static structure factors. For the binary mixture the values are k * = 6.0 (for α = 1) and k * = 5.0 (for α = 2), whereas for the polydisperse model k * = 5.8.
We first discuss the incoherent scattering functions F s (k, t). These are calculated as where the sum is performed over the positions, r j , of the N α particles of the species α. At high temperature the correlation functions decay rapidly to zero in a simple exponential fashion. Around and below T * , marked oscillations appear at short times t ∼ 1, which can be attributed to the vibrations of individual particles within the clus- ters. Similar features have been observed in the cluster crystal phases of similar ultrasoft particles [25], and are known to be associated to single-particle vibrational modes [44]. As the temperature is further decreased, the damping of these oscillations becomes weaker and the amplitude gets larger, consistently with an increased stability of the clusters. The oscillations are also visible in the polydisperse model (see Fig. 16a), though they are less pronounced due to the average over all particles' sizes.
The appearance of oscillations in F s (k, t) is accompanied by evident signatures of glassy dynamics. A plateau in F s (k, t) develops at intermediate times and apparently grows in a continuous fashion from zero to finite values. Such an increase of the plateau height indicates a dynamical slowing down of the continuous type. Interestingly, this feature resembles the type-A transition scenario predicted by the Mode Coupling Theory for certain classes of glassy systems [45,46]. It should be noted that the emergence of the plateau occurs clearly above T * , i.e., the onset of slow dynamics already takes place prior to the thermodynamically defined clustering transition, both for the binary mixture and for the polydisperse system. This is consistent with the smoothness of the transition (see Figs. 5 and 6). Concomitant with the observed decrease of the diffusivities observed in Fig. 13, the plateau extends over longer time scales and is followed by a slow decay with increasing relaxation time as temperature decreases.
The coherent scattering functions are calculated as where . Inspection of the coherent functions F (k * , t) of the binary mixture ( Fig. 15c-d) and the polydisperse model (Fig. 16b) reveals similar features of glassy dynamics. Since coherent functions probe collective correlations, the data of Figs. 15c-d and 16b reflect a progressive dynamic arrest of the cluster structure. At very low temperature, F (k * , t) does not relax to zero within our observation times and the cluster structure is effectively frozen. The highest temperature at which F (k * , t) does not completely relax to zero over the available time scale defines the cluster glass transition and will be denoted as T g . We find T g ≈ 0.48 and T g ≈ 1.02 for the binary mixture and the polydisperse model, respectively [47].
The slowing down of F (k * , t) above T g follows the two-step scenario typically observed in supercooled liquids. This allows one to define, as usual, the corresponding structural relaxation times τ from the condition F (k * , τ ) = 1/e. Figure 17 shows the temperature variation of the species-dependent relaxation times τ 1 and τ 2 for the binary mixture and the total relaxation time τ for the polydisperse model. The data have been tentatively fitted using the VFT function, yielding values of T 0 in the range 0.42-0.45 for the binary mixture and 0.65-0.7 for the polydisperse model. These temperatures can be taken as lower bounds for the cluster glass transition, corresponding to infinitely slow cooling rates. Furthermore, the separation between T * and T 0 is very pronounced in the polydisperse model, which indicates that the cluster glass formed by the latter is "stronger", in the Angell's terminology [48], than the one formed by the binary mixture. All this clearly shows that the clustering and cluster glass transitions are distinct and different in nature: the former is an equilibrium crossover, the latter an ergodicity breaking transition that would occur at different temperatures depending on observation times.
Let us summarize the scenario emerging from our simulations. Consistent with the structural features discussed in subsection IIIA, the clustering transition is associated to a crossover from a fluid of mostly dissociated particles to a fluid of clusters. The dispersity of the clusters evidenced in IIIA frustrates crystallization of the clusters CMs and allows quenching the fluid to lower temperatures. In the temperature range T g T T * the underlying transport processes are rather complex: by visual inspection of animated particles' trajectories, we observed diffusive motion of clusters as a whole, branching of clusters as well as random walks of isolated particles. This results in non-Arrhenius temperature dependence of the transport coefficients and a non-trivial behavior of the intermediate scattering functions. Below T g the structure of the clusters' centers-of-mass becomes practically arrested on the time scale of our simulation (see below) and the system enters in the cluster glass regime. Having noted this, the incoherent scattering functions still decay to zero in the time scale of the simulation at some temperatures T < T g (see Figs. 15a-b and 16a). In other words, single-particle dynamics is still effectively ergodic in the cluster glass and is driven by hopping between the effectively arrested clusters. This hopping dynamics is manifested by the Arrhenius-like behavior observed in the diffusivities (see Fig. 13). Therefore, as in cluster crystals [22], single-particle and collective dynamics in the cluster glass are strongly decoupled. Representative particles' displacements illustrating the difference in transport below and above T g are depicted in Fig. 18a-b and Fig. 18c-d, respectively. Panels (a) and (b) show typical single-particle hopping processes at low T , involving jumps over a few neighbor distances. Note that, during the jumps, the motion of the particle is nearly ballistic, as in cluster crystals [25], and may involve passing through different clusters. Panels (c) and (d) illustrate the more complex motions occurring above T g , which involve both slow diffusion of particles residing within clusters (panel (c)) and diffusion of isolated particles (panel (d)).
To conclude this section, we provide now explicit evidence that the cluster glass transition corresponds to the loss of ergodicity in the degrees of freedom associated to the clusters' centers of mass. At each time t we identify the positions R j (t) of the clusters' centers of mass, where 1 ≤ j ≤ N cl (t) and N cl (t) is the total number of clusters in the system at time t. We then calculate the coherent scattering functions of the clusters' CMs as where and the sum is done over the N α cl (t) clusters of the species α. The chemical species α of a given cluster is again defined as the majority species in the cluster (see above). Single clusters (i.e., isolated particles) are excluded from this analysis. It is worth mentioning that the averages in Eq. (4) are performed in a grand canonical ensemble, since the number of clusters in the sample fluctuates in time. A comparison between F αβ (k, t) and F CM αβ (k, t) of the binary mixture is shown in Fig. 19(a-b), for wave vectors at the first peak in the corresponding static structure factors, and at different temperatures. Both data sets follow closely each other. When looking at the former scattering functions at some fixed, low temperature, but changing the wave-vector (Fig. 19(c-d)) the agreement remains overall good, especially between big particles and big clusters (2-2 correlations). With this, we conclude that the col-lective slowing down of the fluid is indeed driven by the progressive arrest of the clusters' structure.

C. Dependence on thermal history and microscopic dynamics
To corroborate our interpretation of T g as a glass transition occurring at the clusters' level, we now study how the clusters' structure and dynamics depend on the quench history obtained from three different thermal histories for the binary mixture. The difference between the data sets considered in this section lies in the time t eq allowed for equilibration: the equilibration runs are such that the partial RMSD R α (t eq ), evaluated at the end of the run, is always larger than one, two, and three, respectively. Note that this equilibration criterion is slightly different from the one employed in the rest of the work (indicated in Fig. 20 for comparison), which is based on the total RMSD of particles. Production runs are typically 4 times longer than the corresponding equilibration runs [49]. In Fig. 20 we analyze the influence of the different quench rates on the T -dependence of the diffusion coefficients. We remark that, since the equilibration times are adjusted at each thermodynamic state to obtain the target RMSD, a measure of the average quench rate is not appropriate. In fact, as shown in the insets of Fig. 20, the equilibration times increase exponentially with decreasing T . Within statistics the three data sets of Fig. 20 show a perfect overlap for T > T g . However, small but discernible differences are visible below T g (error bars are of the order of the symbol size). Thus, longer equilibration runs lead to slightly lower diffusion coefficients. A similar and more pronounced effect is visible in the relaxation of the coherent intermediate scattering functions, which are shown in Fig. 21 at T = 0.48 = T g for two different quench rates. Not only τ increases with increasing R α (t eq ), but also the extent to which F α (k, t) has relaxed during the simulation increases with increasing R α (t eq ), hence with the total length of the production run. This, in turn, shifts the cluster glass transition to lower temperatures.
Whereas at the superficial level these observations resemble the typical aging process in structural glasses, the microscopic origin of the quench rate dependence of the diffusion coefficients is not trivial. In fact, at least two factors may affect the single particle dynamics below the cluster glass transition. First, the structure of the clusters' CM relaxes faster upon decreasing equili- bration times. The single particle dynamics, in turn, is affected by the relaxation of the cluster structure, since some relaxation channels (e.g., diffusion of a cluster as a whole), which would be suppressed upon longer annealing, remain active and increase the diffusivity.
One may also speculate that a further contribution to the single particles' mobility arises from differences in the topology of the cluster glass, which determines in turn the available pathways for particles' hopping. To assess whether this is a possible scenario, we analyze the radial distribution functions g CM 11 (R) of the clusters' CM in the binary mixture as a function of the equilibration criterion (see Fig. 22). We see that the cluster structure of the most slowly annealed sample differs appreciably from the others. Though the positions of the peaks are unchanged, their amplitude is increased, indicating a more pronounced ordering. These somewhat large differences can be explained by the fact that, as evidenced in the inset of Fig. 20, the average quench rate does not change linearly with the value of R α (t eq ). [50] From the comparison of cluster crystals and cluster glasses (see Fig. 14), we see that hopping in a more disordered matrix decreases the activation barrier for diffusion, and hence increases the particles' mobility. The actual dependence of D on the quench rate shown in Fig. 20 suggests therefore that a similar effect might be at play in this model. More systematic investigations should be performed to assess how the specific structure of the disordered matrix affects the transport properties at the single particle level.
Finally, we briefly analyze the role of the microscopic dynamics on the transport properties. Following [25] we compare the diffusion coefficients obtained from Newtonian and Monte Carlo dynamics for the polydisperse model. The results of this analysis in the polydisperse system are shown in Fig. 23. For temperatures above the clustering transition the diffusion coefficients obtained from the two methods differ qualitatively above T * , in that the MC data saturate at high temperature, whereas the MD data show a monotonic increase. Such differences are expected, since at high T the GEM-n interactions will be negligible (random walk on a nearly flat energy landscape). Therefore in that limit the MC dynamics will become purely random and T -independent. On the contrary, in Newtonian dynamics the onset time of the crossover from ballistic to diffusive motion in the mean square displacement increases systematically with increasing temperature, and so does the diffusivity. Interestingly, in this regime, the temperature dependence of D is approximately given by a power law T ν , with ν ≈ 2.2. This value of ν is higher than the ones expected for both ideal gas (ν = 0.5) and Brownian motion (ν = 1), and remains to be understood.
Only below T g (i.e., in the cluster glass) the data sets obtained from the two dynamics can be collapsed reasonably well onto a unique master curve, by appropriately rescaling the time units. With this, the microscopic dynamics does not seem to play a relevant role in the transport properties of the cluster glass. This result is rather different from the observation in cluster crystal phases [25]. In such systems the Newtonian dynamics leads to long persistent jumps over the lattice sites. These highly directional motions are suppressed in the MC dynamics, and the corresponding MD and MC diffusivities in the cluster crystal exhibit a qualitatively different T -dependence, i.e, they cannot be mutually rescaled by a T -independent time factor [25]. The fact that, instead, scaling of Newtonian and MC dynamics works reasonably well in the cluster glass suggests that the disordered structure of the clusters' centers-of-mass partly suppresses the mentioned long jumps in Newtonian dynamics. A more systematic analysis is required to completely settle this issue. Work in this direction is in progress.

IV. CONCLUSIONS
We have presented a computational investigation of the structural, thermodynamic and dynamic features of dense fluids of ultrasoft fully-penetrable particles. The investigated systems are a binary mixture and a polydisperse system of particles interacting via the generalized exponential model, for which equilibrium cluster crystal phases exist in the monodisperse case. Because of the introduced dispersity, the systems investigated in this work form disordered cluster phases.
We have characterized the structure of the clusters. The analysis reveals a microsegregation of the big and small particles, and a strong homo-coordination in the case of the binary mixture. The analysis of the thermodynamic observables does not provide yet evidence for a thermodynamic transition associated to the clustering transition. Instead, the clustering transition appears as a smooth crossover to a regime in which most of the particles are located in clusters, isolated particles being infrequent. The structure of the clusters' centers-of-mass mirrors that of the dissociated fluid: the binary mixture effectively self-assemble into a binary mixture of clusters, whereas in the polydisperse model the clusters remain polydisperse in population and size.
The dispersity of the clusters drives a progressive dynamical arrest at the level of the clusters' centers-ofmass on approaching the "cluster glass transition", T g . The latter is operationally defined by inspection of the coherent scattering functions, which exhibit characteristic features of glass-forming liquids above T g and do not completely relax to zero below T g . The diffusivities exhibit a fragile-to-strong crossover upon decreasing temperature. The onset of strong (Arrhenius-like) behavior occurs below the clustering transition and signals a change in the transport mechanisms. Relaxation below the cluster glass transition is driven by particle hopping between the nearly arrested clusters. The analysis of the dependence of dynamic and structural properties on the equilibration times confirms the glassy nature of the systems at low temperatures. Finally, we have investigated the role of the microscopic dynamics in the transport properties, by comparing results from Newtonian and Monte Carlo simulations. At low temperature, the diffusivities obtained by both methods can be related reasonably well by a single scaling factor. This suggests that the microscopic dynamics might play a less significant role for cluster glasses than for cluster crystals. In view of the recent observation of clustering in fully atomistic models of interpenetrable colloidal particles [31], we believe disordered cluster phases of ultrasoft particles deserve deeper and more systematic investigations.