Obtaining equilibrium states in ultrasoft cluster forming systems using a combined thermo- and barostat

In what follows we concentrate at a fixed target temperature ${{T}_{\text{t}}}=0.4$ and present data obtained for three different kinds of computer experiments: (i) a compression run, using a pressure increment of $\Delta P=1$, starting from a random (fluid) configuration at ${{P}_{\text{t}}}=5$; (ii) three expansion runs, applying a pressure increment of $\Delta P=-1$, starting from a particle arrangement during the preceding compression run and assuming initial pressure values of ${{P}_{\text{m}}}=50$, 80 and 95, respectively; and (iii) two 'new' compression runs, applying a compression rate of $\Delta P=1$, starting from an initial configuration obtained during the expansion runs with an initial pressure value of ${{P}_{\text{m}}}\sim 45$. In these experiments we used a rather small value for the pressure increment $\Delta P$, which ensured a better sampling along the ${{P}_{\text{m}}}$-axis and consequently to better statistics of the data displayed in the following in the $\left(\rho ,{{P}_{\text{m}}}\right)$-plane.

In figure 1 we display the measured data for the density of the system as a function of ${{P}_{\text{m}}}$ along the aforementioned compression, expansion and 'new' compression runs. The main panel shows the data of the compression (red) and of the three expansions runs (blue); the initial pressure values for the latter ones are marked by big black circles. From a visual inspection we identify the transition between the liquid and the ordered phase to be located at ${{P}_{\text{m}}}\sim 15$; we note that in the liquid phase (i.e. for ${{P}_{\text{m}}}\lesssim 15$) the density data obtained along compression and expansion runs coincide, i.e. ${{\rho}_{\text{comp}}}(P)={{\rho}_{\text{exp}}}(P)$. This is definitely not the case for the ordered phase, as can be seen by comparing the different sets of data shown in figure 1: starting from different ${{P}_{\text{m}}}$-values, the density data along the expansion curves first follow for ${{P}_{\text{m}}}\gtrsim 30$ parallel trajectories in the ($\rho ,{{P}_{\text{m}}}$)-plane, which differ distinctively from the density data accumulated along the compression curve. As the pressure is further decreased, a second regime can be identified for $15\lesssim {{P}_{\text{m}}}\lesssim 30$, where the three density curves become as functions of ${{P}_{\text{m}}}$ steeper and eventually merge: they finally join for ${{P}_{\text{m}}}\lesssim 15$ the aforementioned branch of data for the liquid phase. The diverging data along the different expansion runs provide unambiguous evidence that the system is not in equilibrium, neither during during the compression nor during the expansion runs. Thus, the measured data for the density, i.e. $\rho =\rho \left({{P}_{\text{m}}},{{T}_{\text{t}}}=0.4\right)$ do not correspond to data that one would obtain from the as yet unknown equation-of-state of the system. This issue will be addressed in the subsequent subsection.

Before we focus on this aspect, we try to shed some light onto those processes that are characteristic when compressing and expanding our system. From the results discussed in the previous paragraph one might suspect that the system is suffering from dynamic lag and that—because the measured values for the density differ along the compression and the expansion runs—our observations are related to a rate-dependent hysteresis effect. To put these hypotheses to a thorough test we select two particle configurations that were realised along two of the expansion runs and start to compress them again from a particular pressure value onwards (using $\Delta P=+1$). The results for $\rho \left({{P}_{\text{m}}}\right)$ as measured during these new compression runs are shown in the inset of figure 1, where the initial pressure values are marked by big black triangles. If, according to our assumption, the system indeed suffers from rate-dependent hysteresis effects, then we would expect that the expansion runs contribute to healing those lattice defects that are created during the first compression run; as a consequence, we would expect that the density data obtained during the new compression runs would be located along a trajectory that is parallel to the data obtained during the first compression run, but shifted to higher density values. On the contrary, the curve of the new compression run follows the trajectory of the expansion run, and when the point where the data of the expansion and the initial compression run meet is reached, i.e. the $\rho \left({{P}_{\text{m}}}\right)$-value at which the expansion run was launched, the slope of the new compression curve changes and it follows the initial compression one.

In an effort to definitely exclude the occurrence of dynamic lags during the compression runs we compare in figure 2 the measured data for $\rho \left({{P}_{\text{m}}},{{T}_{\text{m}}}\right)$ with ${{T}_{\text{t}}}=0.4$, obtained during compression runs applying different pressure increments: $\Delta P=1,5,10$ and 25. Different symbols are used for the different $\Delta P$-values (as labeled in the figure caption), in addition each symbol is colour coded according to the respective ${{T}_{\text{m}}}$-value (see equation (3)). In panel (a) of this figure we display results for $\rho \left({{P}_{\text{m}}},{{T}_{\text{m}}}\right)$, retaining only those data with $\delta T<0.02$. We observe that, even though the measured density values are roughly located on the same curve, there is a considerable scattering of the data. This spreading of the results can be decreased by reducing the tolerance parameter $\delta T$: the corresponding data are shown in the other two panels, where we have imposed via a smaller $\delta T$-value a harder criterion on the data: $\delta T<0.002$ (panel (b)) and $\delta T<0.0002$ (panel (c)). The results presented in these two panels provide unambiguous evidence that—provided we have chosen a sufficiently small $\delta T$-parameter—the measured data for the density of the system are independent of the actual value of $\Delta P$; thus we can definitely exclude the occurrence of dynamic lag effects in our computer experiments.

Zoom In Zoom Out Reset image size

We come back to figure 1 and focus now on the observation that curves in the ($\rho ,{{P}_{\text{m}}}$)-plane that collect data along different runs have different slopes. In an effort to provide a deeper insight into the mechanisms behind those processes that occur along these branches, we recapitulate the two principle mechanisms of how a cluster crystal reacts to a change in volume: (i) shrinkage/growth of the lattice constant, a, which is accompanied by a shift of the position of the main peak of the radial distribution function, ${{d}_{\text{nn}}}$, of the system towards a lower/higher r-value; (ii) deletion/creation of new lattice positions, reflected by an increase/decrease of the average cluster occupation, $\langle {{n}_{\text{occ}}}\rangle$; this quantity is obtained via a cluster analysis of the particle configurations, as specified in the appendix of [18].

We have displayed the data for $\langle {{n}_{\text{occ}}}\rangle$ as a function of $\rho$ along the compression, the expansion and the 'new' compression runs in figure 3 in two different representations: (i) in the main panel the $\langle {{n}_{\text{occ}}}\rangle$-data are colour coded according to their corresponding ${{d}_{\text{nn}}}$-value; (ii) $\langle {{n}_{\text{occ}}}\rangle$-results in the inset are coloured according to the type of runs, i.e. data from compression and 'new' compression runs are displayed as red symbols, while the results from expansion runs are coloured in blue, thereby allowing the reader to establish a connection from these data to those in figure 1.

Zoom In Zoom Out Reset image size

Surprisingly, the reactions of $\langle {{n}_{\text{occ}}}\rangle$ and of ${{d}_{\text{nn}}}$ on a compression or an expansion of the system help us to classify the different branches of data as functions of $\rho$ into three categories:

• 1.  
  Compression branches (dashed line): along these branches $\langle {{n}_{\text{occ}}}\rangle$ increases linearly with the density. After an initial decrease for small densities, ${{d}_{\text{nn}}}$ assumes an essentially constant value of ${{d}_{\text{nn}}}\sim 1.28$ up to high densities. This value corresponds to the minimum value that was observed in [18] in the compression experiments presented there;
• 2.  
  Expansion branches (dash–dotted line): these data were extracted during the final parts of the expansion runs, i.e. (${{P}_{\text{m}}}\gtrsim 30$). Along these branches $\langle {{n}_{\text{occ}}}\rangle$ shows again a linear dependence on the density, although with a larger slope than along the compression branches. Similar to along the latter branches, we observe a small variation of ${{d}_{\text{nn}}}$ at small densities; a similar observation is made at rather high densities; however, over a large density range ${{d}_{\text{nn}}}$ assumes a constant value of ${{d}_{\text{nn}}}\sim 1.48=d_{\text{nn}}^{\text{max}}$;
• 3.  
  Intermediate branches (solid lines) data along these branches stem either from the 'new' compression runs or from the initial part of the expansion runs (i.e. ${{P}_{\text{m}}}\gtrsim 30$). As $\langle {{n}_{\text{occ}}}\rangle$ remains constant along these branches no cluster merging or splitting events occur. These sets of data connect the compression and the expansion branches, thus ${{d}_{\text{nn}}}$ varies between $d_{\text{nn}}^{\text{max}}=1.28$ and $d_{\text{nn}}^{\text{min}}=1.48$.

We summarise our observations as follows: if an ordered system of GEM-4 particles is compressed or expanded, the two fundamental reaction mechanisms of the system are completely decoupled and they obey the following scheme: (i) as an immediate reaction on compression/expansion the lattice spacing a of a cluster crystal is reduced/enhanced; it varies within a range $\left[d_{\text{nn}}^{\text{max}},d_{\text{nn}}^{\text{min}}\right]$ while $\langle {{n}_{\text{occ}}}\rangle$ remains constant; according to our above nomenclature, the system 'moves' along an intermediate branch; (ii) once the system has reached a state where $a=d_{\text{nn}}^{\text{min}}$, further compression will induce cluster merging processes [18]: the required volume reduction is effectuated by deleting lattice positions via cluster merging events (i.e. $\langle {{n}_{\text{occ}}}\rangle$ increases), while keeping the lattice spacing constant fixed at $a=d_{\text{nn}}^{\text{min}}$; for this part of the compression process the system 'moves' along a compression branch; (iii) alternatively, if the system—once it has reached a state where $a=d_{\text{nn}}^{\text{max}}$—is further expanded, the considerably increased inter-cluster distance will reduce the repulsion among the clusters, and consequently the aggregates will start to split: the larger available volume will be filled by an increased number of newly created, but smaller clusters (i.e. $\langle {{n}_{\text{occ}}}\rangle$ decreases), while keeping the lattice spacing constant; for this part of the expansion process the system 'moves' along an expansion branch.

We now return to the question of how we can obtain in our experiments the density value that corresponds at a given pressure to the density that one would obtain from the equation-of-state of the system (i.e. the equilibrium state). So far, we can summarise our observations as follows: we obtain—depending on the protocol of our compression and expansion processes—different results for $\rho$, a and $\langle {{n}_{\text{occ}}}\rangle$ for a given ${{P}_{\text{m}}}$-value; thus, we can conclude that the system is driven by the different processes into metastable states that are separated from the equilibrium state by large energetic barriers. To surmount these, we launch several annealing processes at constant ${{P}_{\text{m}}}$, hoping we will recover in this manner an equilibrium state in the system. The starting configurations for the different annealing runs are selected from state points along the compression, expansion and intermediate branches.

During the annealing processes we keep the target pressure constant and increase the temperature until the system melts; then we cool the system down to the desired target temperature ${{T}_{\text{t}}}=0.4$. We hope that the energy that is introduced into the particles via the heating process allows the system to overcome the energetic barrier that separates the metastable states from the equilibrium state. If these annealing runs indeed lead to the equilibrium state, we expect all these processes (that start for a given ${{P}_{\text{m}}}$-value from different $\rho$-values) to end up at the same density, which we identify as the equilibrium density at this pressure value, i.e. ${{\rho}_{\text{equ}}}\left({{P}_{\text{m}}}\right)$. In total we have performed annealing runs at ${{P}_{\text{m}}}=30,40,53,60$ and 80.

In figure 4 we show the data for the measured density as functions of ${{T}_{\text{m}}}$ for three annealing processes: they have been launched from state points that are characterised by ${{P}_{\text{m}}}=70$ and are located (i) at the compression branch (red symbols, initial density value of ${{\rho}_{\text{init}}}\sim 6.4$), (ii) the middle intermediate branch (green symbols, ${{\rho}_{\text{init}}}\sim 6.6$) and (iii) the upper intermediate branch (blue symbols, ${{\rho}_{\text{init}}}\sim 7.25$). The ${{\rho}_{\text{init}}}$-values are marked by grey rectangles and the final value of the density after the annealing process ${{\rho}_{\text{final}}}$ is marked as a purple circle.

Zoom In Zoom Out Reset image size

During the initial heating, the sets of data of the compression branch and of the middle intermediate branch merge at ${{T}_{\text{m}}}=0.6$ and follow for the rest of the annealing process the same $\rho \left({{T}_{\text{m}}}\right)$-curve. The results from the upper intermediate branch merge with the other sets of data at ${{T}_{\text{m}}}=1.1$; by a visual inspection of the particle configurations we can conclude that at this temperature the nanocrystal begins to melt. These annealing processes are stopped as the system reaches a temperature of ${{T}_{\text{m}}}=1.2$, a typical snapshot of a particle configuration taken at that instant is shown in the right panel of figure 4, providing evidence that the system has completely lost its ordered structure. Now the increment in temperature, $\Delta T$, is reversed from $\Delta T=0.1$ to $\Delta T=-0.05$ and the three sets of data follow during the subsequent cooling run essentially the same curve; eventually they reach a final, common state whose density is marked in figure 4 as a violet circle; obviously this value corresponds to the density of the equilibrium state, ${{\rho}_{\text{final}}}={{\rho}_{\text{equ}}}\left({{P}_{\text{m}}}=70\right)$.

This protocol of annealing processes was repeated for different pressure values, namely for ${{P}_{\text{m}}}=30,40,53,60$ and 80; the corresponding sets of data, displayed in the ($\rho ,{{T}_{\text{m}}}$)-plane are similar to those obtained for ${{P}_{\text{m}}}=70$ (and shown in figure 4). Each of these annealing runs carried out for a given ${{P}_{\text{m}}}$ leads to the desired ${{\rho}_{\text{equ}}}\left({{P}_{\text{m}}}\right)$-value. These density data (corresponding to ${{\rho}_{\text{equ}}}\left({{P}_{\text{m}}}\right)$, shown in purple) are summarised in figure 5, where, in addition, the original data for the density, $\rho \left({{P}_{\text{m}}}\right)$, as obtained from the compression (light-red symbols) and from the expansion (light-blue symbols) runs are shown (see figure 1). The respective ${{\rho}_{\text{init}}}$-values are plotted with dark-grey asterisks. The curve that connects as a function of ${{P}_{\text{m}}}$ the ${{\rho}_{\text{equ}}}\left({{P}_{\text{m}}}\right)$-data represents at ${{T}_{\text{t}}}=0.4$ the equation-of-state, i.e. ${{\rho}_{\text{equ}}}\left({{P}_{\text{m}}},{{T}_{\text{t}}}=0.4\right)$.

Zoom In Zoom Out Reset image size

Of course, the final answer if the obtained state indeed corresponds to an equilibrium state can only be given via additional free energy calculations. However, due to the high numerical costs for our different types of compression, heating and annealing runs we postponed additional free energy calculations (which themselves would have been rather time consuming) to future investigations on this system.

Our conclusions are supported by a structural analysis of the systems, putting again focus on ${{T}_{\text{t}}}=0.4$: we have investigated the correlation between the average distance of two neighouring clusters, ${{d}_{\text{nn}}}$, and the hexagonal order parameters, ${{\Psi}_{6}}$. The latter one is defined as (see [21])

$$\begin{eqnarray}&&{{\Psi}_{6}}=\frac{1}{N}\underset{j}{\overset{N}{\sum}}\,\frac{1}{{{n}_{j}}}\underset{k}{\overset{{{n}_{j}}}{\sum}}\,\text{exp}\left(\text{i}6{{\Theta}_{jk}}\right),\end{eqnarray} \tag{ 4 }$$

where n_j denotes the number of nearest neighbours of cluster j, and ${{\Theta}_{jk}}$ the angle between the line connecting cluster j and its nearest neighbour k with a fixed axis. By definition, ${{\Psi}_{6}}$ attains a value of 1 for a perfect hexagonal arrangement, while we typically found ${{\Psi}_{6}}\approx 0.4$ for the disordered liquid state.

In figure 6 we show this correlation between ${{d}_{\text{nn}}}$ and ${{\Psi}_{6}}$, considering all the expansion (red symbols), compression (blue symbols) and annealing (violet symbols) runs, as well as all pressure values P_rmm. We observe that during the compression and the expansion runs the maximum value of the order parameter obtains values of ${{\Psi}_{6}}\sim 0.975$ while after the annealing processes we obtain ${{\Psi}_{6}}\sim 0.985$: thus, we conclude that annealing is able to heal out lattice defects in the cluster crystal. The maximum value obtained for ${{\Psi}_{6}}$ corresponds to a nearest neighbour distance of ${{d}_{\text{nn}}}\sim 1.37$ which we thus consider the lattice spacing of the equilibrated state. This result differs from the corresponding value predicted density functional theory calculations (i.e. ${{d}_{\text{nn}}}=1.42$); we attribute this difference to the fact that we are investigating a nanodroplet of finite size and not a bulk system.

Zoom In Zoom Out Reset image size