We investigate the phase ordering (pattern formation) of systems of two-dimensional core–shell particles using Monte Carlo (MC) computer simulations and classical density functional theory (DFT). The particles interact via a pair potential having a hard core and a repulsive square shoulder. Our simulations show that on cooling, the liquid state structure becomes increasingly characterized by long wavelength density modulations and on further cooling forms a variety of other phases, including clustered, striped, and other patterned phases. In DFT, the hard core part of the potential is treated using either fundamental measure theory or a simple local density approximation, whereas the soft shoulder is treated using the random phase approximation. The different DFTs are benchmarked using large-scale grand-canonical-MC and Gibbs-ensemble-MC simulations, demonstrating their predictive capabilities and shortcomings. We find that having the liquid state static structure factor S(k) for wavenumber k is sufficient to identify the Fourier modes governing both the liquid and solid phases. This allows us to identify from easier-to-obtain liquid state data the wavenumbers relevant to the periodic phases and to predict roughly where in the phase diagram these patterned phases arise.
I. INTRODUCTION
Systems of particles interacting pairwise via purely repulsive interaction potentials can exhibit a surprisingly rich phase behavior. The particles that we consider here have hard impenetrable cores, surrounded by a penetrable corona. The resulting pair potential between the particles exhibits a hard core of diameter σ, beyond which is a repulsive shoulder of range λσ, with energy penalty ϵ > 0 for overlapping coronas. The range parameter λ > 1 depends on the extent of the coronas. Here, we treat the system as being two-dimensional (2D); the three-dimensional analog is also very interesting, but beyond the scope of the present study. An experimental realization of the 2D model is the system studied in Ref. 1, consisting of metallic nanoparticles with a corona of attached polymers. These particles are then adsorbed at a water–air interface, rendering the system effectively 2D. Other such systems include microgels2–4 or various types of colloidal spheres decorated with polymer-chains.5–7 These core–shoulder systems are of fundamental importance because understanding and predicting the phase behavior provide a crucial benchmark challenge for computer simulations and theories for interacting classical many-body systems. They are also of practical importance because the colloidal structures that self-assemble promise useful optical and/or material properties.5
Core–shoulder systems do not exhibit gas–liquid phase separation due to the lack of attractive interactions. Nonetheless, complex phase ordering can occur depending on the range of the shoulder potential. Typically, they form cluster, stripe, and hole phases at intermediate densities, as long as the shoulder repulsion is strong enough (or, equivalently, when the temperature is low enough),8 which at first sight can appear like the generic pattern formation exhibited by systems with competing interactions.9 However, at higher densities and lower temperatures, a plethora of different ordered phases can arise, depending on the density, temperature, and range of the shoulder repulsion λ, often with a sensitive dependence on the precise model parameter values. Belying the simplicity and isotropic nature of the interactions, the structures that can form range from the clusters and stripes already mentioned to various anisotropic structures6,10–14 and even quasicrystals of various different symmetries.15–17
The notion that softening the repulsive part of the interparticle pair potential can lead to the appearance of additional phase transitions goes back to Hemmer and Stell,18 who investigated systems of hard-core particles “softened” with a weakly attractive potential. Later, Jagla introduced a more general class of two length-scale models19,20 with a rich phase behavior, popularizing the topic. Such potentials have even been used to model liquid–liquid transitions in water.21 Here, we report results for systems having a shoulder potential that is step-like. However, various soft-shoulder analogs have been considered, such as exponential and tanh-shaped shoulders or ramp-like potentials.19,20,22 Our general approach is relevant to this whole class of models.
Our focus here is to identify features of the liquid state that are indicators or precursors of the complex phase behavior that arises at lower temperatures. In other words, to allow for easy navigation of the phase diagram, we identify the features to look for in the liquid state that are associated with the anisotropic phases. To describe the structural properties of our system, we have on the one hand performed extensive Monte Carlo (MC) based simulations, including both grand-canonical-ensemble MC (GCMC) and Gibbs-ensemble MC (GEMC),23,24 and on the other hand, we have used various complementary theoretical approaches, all based on classical density functional theory (DFT).25–27
A key part of our approach is to use DFT to determine the major density modes (i.e., the characteristic wavelengths) of the liquid state structure. We use two different DFT approximations. The first is a simple DFT that has the advantage of providing simple analytic expressions for many of the quantities of interest, while the second is a more sophisticated (fundamental measure theory26 based) DFT, which still provides semi-analytic results for some quantities. We compare results from the two DFTs and benchmark all predictions against the results from the MC simulations. Specifically, we compare results for the static structure factor S(k) and the radial distribution function g(r)26 obtained from two different DFT approximations and compare with simulation results in order to identify strengths and shortcomings of the DFT theories. We identify two types of errors arising in the DFT results for S(k), specifically a small systematic error in the determination of peak position and an overestimation of peak height. However, by comparing with the MC simulations, we can understand the scale of the shift and so compensate for the errors in order to use the DFT results in order to accurately predict the location of the peaks in S(k). This, therefore, allows us to identify the characteristic wavenumbers k of the density modes favored by the system. This analysis also allows us to understand how these errors manifest in DFT calculations for the density profiles in various structured phases. In looking for hallmarks of structured phases in our MC computer simulations of the liquid phase, in addition to calculating S(k) and g(r), we also calculate bond orientational order parameters and identify how they vary on approaching the structured phase from liquid state side of the phase diagram.
This paper is structured as follows: In Sec. II, we briefly explain the pair potential defining our model system and the parameters on which the state of the system depends. In Sec. III, we first explain various MC simulation approaches developed, including both the grand-canonical and Gibbs ensemble MC schemes we use. Then, in Sec. III B, we describe the two DFTs that we have applied and how to obtain the correlation functions. In Sec. IV, we present our results, structured as follows: In Sec. IV A, we consider the liquid–cluster transition, as an example to illustrate the changes in structural properties across a phase transition using GCMC simulations. By direct simulation of the phase coexistence using GEMC, we demonstrate that the dominant density modes in the liquid can be associated with the wavevectors shaping the structured phases. In Sec. IV B, we present results from a stability analysis of a HCSS bulk liquid. We determine the relevant specific density modes and corresponding wavevectors introduced by the shoulder part of the pair-interaction, elucidating how the ordering arises from the interplay between these and those from the entropic volume-exclusion effects of the hard cores. The influence of the shoulder-interaction on the freezing transition of a pure hard-core fluid is also discussed. In Sec. IV C, we compare the liquid state static structure factors S(k) obtained using GCMC with the analytic results using DFT. This allows us to highlight shortcomings and capabilities of two different DFTs and quantify errors made in the determination of the main structural properties of the liquid. In Sec. IV D, the effect of the latter errors on the predicted equilibrium density profiles ρ(r) is illustrated by comparison to GCMC results. Furthermore, we demonstrate that the analytic predictions from DFT for the bulk liquid suffice to accurately determine the significant wavevectors shaping striped, clustered, and holed phases of the HCSS system. Finally, in Sec. V, we finish with a few concluding remarks.
II. SYSTEM
III. SIMULATION AND THEORY
A. Simulation methods
1. MC simulations in the grand-canonical ensemble
We have performed large-scale MC simulations in the grand-canonical ensemble (GCMC), imposing particular values for the temperature T, the volume V, and the chemical potential μ of the system.23,24 The number of particles N can vary within the system but has been limited—for practical reasons—to 20 000 particles. The ensemble is confined within a square box for the results presented in Secs. IV A–IV C, whereas we use rectangular boxes with aspect ratio for the solid state simulations presented in Sec. IV D, in order to minimize frustration from straining the crystal via the applied periodic boundary conditions. Three types of MC “moves” are implemented in the course of GCMC simulations: translational moves, particle insertions, and particle deletions, with respective probabilities αt, αi, and αd and suitably adapted acceptance-/rejection-criteria.23,24 In our simulations, we have assumed that αt = αi = αd = 1/3 and have applied the three types of moves in a random consecutive manner. Typical simulations extend in total over 15–20 × 109 of such MC-moves. Equilibration of the energy was typically observed after 4–5 × 109 MC-moves (see subsection 2 of the Appendix), and ensemble averages were performed using 150–250 states drawn from the end of the simulation separated by 4 × 108 MC-moves.
2. Gibbs ensemble MC simulations
The Gibbs ensemble represents a very particular ensemble within the framework of statistical mechanics.23,24,28 The entire system is considered in the canonical (i.e., NVT) ensemble, which is subdivided into two subsystems (confined in sub-boxes), which are assumed to be in phase coexistence, meaning that the two subsystems are at equal temperature T, equal pressure P, and equal chemical potential μ. This is achieved by introducing three types of MC “moves”: translational moves of particles inside each of the boxes, particle exchange between the two boxes, and volume exchange between the two subsystems while maintaining the total volume of the system. Related acceptance/rejection criteria for each type of moves can be found in the literature.23,24 Periodic boundary conditions are applied to each sub-box.
In our simulation, we start from two square sub-boxes (indices “1” and “2”) with N1 = N2 = 1500 particles, the total number of particles N1 + N2 = 3000 being fixed. The initial box size is chosen such that V1 = V2 = V/2, where the total volume V = N/ρb, where ρb is the starting total average density and N = N1 + N2. The positions of the particles are initialized in hexagonal lattices.
The simulations are performed via consecutive blocks where each block represents a sequence of the following consecutive steps: (i) ntrans = 200 attempts for translational MC-moves of randomly chosen particles, i.e., 100 moves per sub-box; (ii) a single nvolume = 1 volume change, changing the volume of one sub-box by ±ΔV at the cost of the other, where ΔV is drawn uniformly from the range ΔV ∈ [0, 4]; and (iii) nswap = 500 attempts to swap particles from one sub-box to the other, where in each step the box in which a particle is deleted is selected randomly. The position of a newly created particle is chosen uniformly from within the hosting sub-box; the move is rejected if an overlap of hard cores with the existing particles occurs. The simulation is essentially a sequence of repeated blocks, being continued until a total number of MC-move attempts nmoves = nblocks × (ntrans + nvolume + nswap) ≥ 8 × 109 is attained. For the ensemble averages (e.g., to calculate the structure factor—see below), each box is treated as a μVT-ensemble of its own, analogous to the GCMC calculations.
B. DFT for the HCSS system
1. DFT formalism and thermodynamics
The key features of the functional Ω[ρ] are that (i) Ω[ρ] is minimized by the equilibrium one particle density profile of a system, ρeq(r), and (ii) for ρeq(r), the functional Ω[ρ] takes the value of the thermodynamic grand potential of the system, i.e., Ω = Ω[ρeq]. In a homogeneous liquid, the equilibrium density profile ρ(r) (henceforward, we drop the subscript “eq”) is equal to the bulk density ρb. However, in the cluster, stripe, crystal, and other such phases, the equilibrium density profile ρ(r) is nonuniform.
For the contribution to the free energy , which originates from the hard-core interactions, we consider here two different approximations: (i) a very simple approximation, namely, a local density approximation (LDA),26 and (ii) a fundamental measure theory (FMT) based functional.31,32 Both functionals are briefly be introduced in the following.
2. Structure
With the LDA- and the FMT-functionals at hand—see Eqs. (9) and (12), respectively—one can calculate in a straightforward manner (albeit somewhat lengthy in the case of the FMT) the function .
Since the static structure factor S(k) can easily be calculated from via Eq. (17) and since in both approximations is given by analytic expressions, one can—in principle—write down the expression for S(k).
3. Stability analysis
Identification of the regions of the phase diagram where the uniform bulk liquid becomes unstable simultaneously enables identification of the regions where the formation of ordered phases can be expected. In the regions where the bulk liquid is unstable, it is because it becomes unstable against the growth of periodic modulations of the density, characterized by a specific wavevector. Here, we demonstrate that the specific wavevectors that play a significant role in describing the onset of instability of the liquid are also sufficient to predict the length-scales that one can expect in (possibly highly complex) ordered structures.
When considering the stability of a liquid, one can take either a thermodynamic/structural or a dynamical perspective. Within the DFT framework, these are closely related since both originate from the free energy functional. A dynamical analysis is based on dynamical density functional theory.26,40–42 To identify the limit of linear stability of the bulk liquid, one approach is to identify the locus in the (T, ρ)-plane where the first peak of the static structure factor (located at kc > 0)—see (17)—is diverging, i.e., S(kc) → ∞. This marginal stability limit of the liquid against the critical wavenumber kc > 0 is referred to here as the λ-line, following the terminology of Ref. 43. This should not to be confused with our use of λ as the parameter defining the range of the pair potential (1). From the dynamical perspective, inside the λ-line, the uniform liquid is dynamically unstable against periodic density modulations with this unstable wavenumber, where (at least within the linearized regime) these modes grow over time without any (free) energetic barrier to surmount.44 Either way, one can formulate the necessary conditions for locating the λ-line based on properties of S(k)−1 [see Eq. (17)] via the following criteria: (i) the extremum condition, i.e., , and (ii) the divergence condition, i.e., .
IV. RESULTS
A. Structural properties across the liquid–cluster phase transition
To identify signatures of incipient self-assembled structures in the liquid phase and to elucidate the changes in the structural properties as the liquid becomes unstable to periodic phases, we first focus on the liquid–cluster phase transition. This has been well studied for the different but somewhat related “short-range-attraction, longer-range-repulsion” (SALR) system.45 We find that at low densities, the phase ordering of the SALR and HCSS systems is similar, but at higher densities, they can be different, depending on the value of λ. The low-temperature clustered phase is characterized by two distinct length-scales, which shape the self-assembled pattern: (i) The larger one is the typical cluster spacing, ac = 2π/kc, which governs both the area occupied by a single cluster as well as the self-assembled hexagonal ordering pattern of the clusters. Below, we show that ac is related to the wavevector k1, corresponding to the first peak in S(k). (ii) The shorter length-scale is the length-scale of freezing of particles within clusters, af = 2π/kf (where kf is the freezing wavevector), governing the spacing and the arrangement of individual particles within a given cluster. Below, we show that af is related to the wavevector of one of the higher-k maxima in S(k). For example, for the case where λ = 3.7 considered below, this corresponds to the fourth maximum, at k4. While ac is predominantly determined by the shoulder width λ, af is in general associated with the volume-exclusion effects of the hard core. In fact, as the density increases, kf converges toward in the limit of a closest hexagonal packing of hard cores, occurring at η ≈ 0.906.
In order to obtain deeper insights into the implications of the DFT-based results in Sec. III B 3 and to determine the relevance of the wavevectors corresponding to peaks in S(k) on the structure, we first need to verify that these are indeed the wavevectors governing the density fluctuations in a marginally stable liquid and also are the ones that ultimately shape the structures that emerge as the liquid becomes unstable. We emphasize that the investigations discussed here are entirely based on GCMC simulations. For this purpose, we focus on one system, characterized by λ = 3.7. This choice is in some sense arbitrary (as we expect similar results for other such intermediate values of λ) and is justified by the fact that for shoulder widths around this value, the relevant wavevectors can be nicely displayed in graphical representations. Furthermore, we fix the values of the chemical potential μ (as specified below) and the temperature to be located in the vicinity of the liquid–cluster phase transition—see the Appendix.
In Fig. 1, we present results obtained from GCMC simulations. These are all state points where the DFT predicts the uniform liquid to be unstable, i.e., are inside the λ-line. In Fig. 1(a), we display the static structure factor S(k) at different temperatures, for states located in the vicinity of the liquid–cluster transition and for the fixed chemical potential value μ(η = 0.2), as obtained from Eq. (15). For this value of μ, the liquid–cluster transition occurs at Tc ≈ 0.36. This is also identified by GEMC—see Fig. 2—but not for exactly the same value of the chemical potential considered here. Additional evidence for this value of the transition temperature at μ(η = 0.2) also comes from inspecting the average packing fraction and bond-order-parameter as a function of temperature; further details are given in the Appendix (see especially Fig. 8). S(k) is obtained via the inverse Fourier transform of Eq. (18) together with the radial distribution function g(r), which is displayed in Fig. 1(b) (and discussed further below). g(r) is calculated in GCMC simulations, using a simulation box of size Lx × Ly = 100 × 100.
(a) The radially symmetric structure factor S(k) as a function of the wavevector k for a HCSS system with λ = 3.7, evaluated for different temperatures T (as labeled) and for constant chemical potential μ(η = 0.2); see Eq. (15). S(k) is obtained via Eq. (18) from the radial distribution function g(r) [see (b)], obtained in GCMC simulations. (b) The radial distribution function g(r) as a function of r, obtained in GCMC simulations for the different temperatures indicated in the key in (a). At the top of (b) are displayed three characteristic snapshots of the liquid and of the cluster phase. The dashed vertical line indicates the distance r = 1.5, relevant to identify nearest neighbors of a tagged particle, while the solid vertical line indicates the range of the shoulder λ = 3.7. (c) The bond orientational order parameter BOO(m) for a selection of m-values (those indicated on the axis) and for different temperatures [as labeled in (a)]. The data for an additional set of BOO(m) have been added, which were evaluated for T = 0.25 (shown in magenta). The dashed lines are BOO(m) evaluated for just nearest neighbors (with distances r ≤ 1.5), while the solid lines are for all neighbors separated by distances up to the shoulder width, i.e., for r ≤ λ.
(a) The radially symmetric structure factor S(k) as a function of the wavevector k for a HCSS system with λ = 3.7, evaluated for different temperatures T (as labeled) and for constant chemical potential μ(η = 0.2); see Eq. (15). S(k) is obtained via Eq. (18) from the radial distribution function g(r) [see (b)], obtained in GCMC simulations. (b) The radial distribution function g(r) as a function of r, obtained in GCMC simulations for the different temperatures indicated in the key in (a). At the top of (b) are displayed three characteristic snapshots of the liquid and of the cluster phase. The dashed vertical line indicates the distance r = 1.5, relevant to identify nearest neighbors of a tagged particle, while the solid vertical line indicates the range of the shoulder λ = 3.7. (c) The bond orientational order parameter BOO(m) for a selection of m-values (those indicated on the axis) and for different temperatures [as labeled in (a)]. The data for an additional set of BOO(m) have been added, which were evaluated for T = 0.25 (shown in magenta). The dashed lines are BOO(m) evaluated for just nearest neighbors (with distances r ≤ 1.5), while the solid lines are for all neighbors separated by distances up to the shoulder width, i.e., for r ≤ λ.
(a) and (b) GEMC snapshots of a HCSS system with λ = 3.7 and T = 0.3663; these states are at the transition from the liquid (a) to the cluster phase (b), having average densities ⟨ρliq⟩ = 0.0902 and ⟨ρcl⟩ = 0.0898, respectively. (c) and (d) Heatmap plots (see color bars) of the static structure factor S(k) as a function of k = (kx, ky). The radii of the white dotted concentric circles in both (c) and (d) are the wavevectors k1, …, k4 of the radially averaged structure factor S(k) in the liquid phase, thus corresponding to the four main peaks in S(k) in the liquid phase. The radius of the white dashed circle is (see text). The small green circles in (d) highlight the Bragg-peaks associated with the hexagonal patterns observed between and within clusters (see text).
(a) and (b) GEMC snapshots of a HCSS system with λ = 3.7 and T = 0.3663; these states are at the transition from the liquid (a) to the cluster phase (b), having average densities ⟨ρliq⟩ = 0.0902 and ⟨ρcl⟩ = 0.0898, respectively. (c) and (d) Heatmap plots (see color bars) of the static structure factor S(k) as a function of k = (kx, ky). The radii of the white dotted concentric circles in both (c) and (d) are the wavevectors k1, …, k4 of the radially averaged structure factor S(k) in the liquid phase, thus corresponding to the four main peaks in S(k) in the liquid phase. The radius of the white dashed circle is (see text). The small green circles in (d) highlight the Bragg-peaks associated with the hexagonal patterns observed between and within clusters (see text).
The plot of S(k) in Fig. 1(a) shows four distinct peaks with values S(ki) = Si, located at wavevectors ki, i = 1, …, 4. These four are those that are within the displayed physically relevant range of k-values, i.e., 0 < ki ≲ khex (i.e., there are other maxima at larger-k, and there can also be a local maximum at k = 0, which are not of interest here). In this plot, we observe that the peak heights S1 and S4 increase significantly as the temperature is decreased. It is interesting to enquire if the so-called Hansen–Verlet (HV) criterion (generalized to 2D systems) applies in the present case.46–49 HV observed empirically for the 3D Lennard-Jones (LJ) system that as the main peak of S(k) passes the threshold S1 ≈ 2.85, the LJ system is prone to freeze. This criterion has been found to roughly apply to a wide variety of different systems, albeit for 2D systems the HV threshold is Si ≈ 5.48,49 For the two higher temperatures displayed in Fig. 1 (i.e., for T = 0.5119 and T = 0.4048), the HCSS system is still in the liquid phase, which agrees with what the 2D HV criterion suggests.
At the lower temperature T = 0.3571, which is close to the transition temperature obtained from the MC simulations (see also the Appendix), Tc ≈ 0.36, where the phase transition to the cluster phase (consisting of a hexagonal lattice of clusters) occurs for the imposed value of μ, S(k) exhibits several pronounced peaks. This is displayed in Fig. 1(a). The first peak (at k1) has height S1 ≈ 4, which is a little below the 2D HV threshold value. In addition, the fourth peak S4 is also prominent, but some way short of the HV threshold. It is possible that this state is located in the two phase coexistence region between the liquid and the cluster phase.
At the lowest temperature T = 0.3334, the system is definitely in the cluster phase, having pronounced peaks with large values of both S1 and S4 (with S1 well above the HV threshold). In this case, S3 assumes a relatively modest value. As a consequence, the length-scale associated with k3 turns out to play only a minor role in the cluster phase. It is also worth noting that the values of k1 and k4 remain essentially unchanged as we vary the temperature, while one can observe a notable temperature-dependent variation of the k2- and k3-values. As discussed below, the large peaks S1 and S4 and the corresponding wavevectors impact the shape of the radial distribution function.
The corresponding radial distribution functions, g(r), are displayed in Fig. 1(b). These display the characteristic features of g(r) for a HCSS-system, namely, the two distinct discontinuities at r = 1 (=σ) and r = λ. These two discontinuities are to be expected and arise essentially as contributions of the core and of the shoulder to the overall pressure in the system. The magnitudes of the discontinuities are related to the bulk pressure of the system.26,50 Furthermore, g(r) is characterized by oscillations on both long- and short-range scales, with wavelengths that can be associated with k1 and k4, respectively. For the two higher temperature liquid states, both types of oscillations decay relatively fast. For the temperature near the transition (T = 0.3571), the long wavelength oscillations in g(r) are of large amplitude, indicating strong inter-cluster correlations. At the lowest temperature displayed (in the cluster phase), the short wavelength oscillations also decay rather slowly and with large amplitude, indicating significant correlations of particle positions both within and between the clusters.
Thus, in the case of an intermediate or large shoulder width and at low densities, where particles have sufficient space to move around freely inside a cluster, we can identify the following two stages of cluster formation: (i) below Tc, but still at sufficiently high temperatures (where S1 is large and well above the HV threshold, but S4 is not), clusters populate a hexagonal lattice, while inside each cluster, the particles form a disordered, liquid-like structure; (ii) at even lower temperatures (where we observe S4 approaches or exceeds the HV threshold, along with S1), we find that the particles freeze inside the clusters into hexagonal particle arrangements (with close core-contact). Furthermore, the clusters themselves are oriented in the same direction. For visualization of this trend, typical snapshots are shown along the top of Fig. 1(b). For the liquid–cluster phase transition at higher densities—which we investigate using GEMC simulations below—the stages (i) and (ii) occur simultaneously, with a large value of S1, above the HV threshold. The particles are constrained to the area of the clusters, which at higher densities immediately leads to freezing due to packing of the cores. Only at lower densities and large λ-values are clusters of liquid-like droplets observed.
The above analysis of the (radially symmetric) structure factor S(k) shows that for the system at hand the transition from the liquid into the cluster phase is characterized by a significant increase of S1 and S4. In the following, we provide evidence that the wavevectors associated with these peaks in S(k) in the liquid state are identical to the wavevectors that are responsible for forming the cluster phase. To do this, in Fig. 2, we now analyze in detail the full angular dependent (i.e., vector k-dependent) structure factor S(k) for the system at hand (i.e., for the cluster phase that is at coexistence with the liquid, which have average densities ⟨ρcl⟩ = 0.0898 and ⟨ρliq⟩ = 0.0902, respectively). S(k) is obtained in GEMC simulations [see the Appendix, Eq. (A6)], directly from the positions of the particles as input. In Figs. 2(a) and 2(b), we display typical snapshots of the two phases involved, while Figs. 2(c) and 2(d) show heatmap plots of S(k) as functions of k = (kx, ky). Figure 2(c) shows the typical radially symmetric result for a fluid state, while Fig. 2(d) exhibits the typical Bragg-peaks of a crystalline state. On top of these, we superimpose four white dotted concentric circles, whose radii correspond—from the center to the periphery—to the aforementioned wavevectors k1, …, k4, which correspond to the four main peaks (S1, …, S4) of the radially averaged structure factor—see Fig. 1. An additional white dashed circle of radius is also displayed. In the case of Fig. 2(d), we observe the following: (i) The circles with radii k1 and k4 match very nicely the brightest Bragg-peaks in S(k) (circled in green) in the cluster phase. Thus, the peaks associated with hexagonal patterns observed in S(k) are located with high accuracy on the circles corresponding to wavevectors k1 and k4. (ii) The Bragg-peaks near k3 are found to have a somewhat bigger k-value in the cluster phase as compared to the liquid phase (i.e., these Bragg-peaks do not lie on the k3 circle), an observation which is in nice agreement with the findings from the GCMC simulations discussed above in relation to Fig. 1. With all this in mind, we can conclude that the cluster phase that emerges as the liquid becomes unstable is shaped by those wavevectors (namely, k1 and k4) that predominantly govern the density fluctuations in the liquid. This points to the possibility that once identified, one can tune these wavevectors, gaining insights into the complex phase behavior of the HCSS-system and be able to design self-assembled phases using specific combinations of wavevectors. We flesh this idea out further in Secs. IV B–IV D, but pursue the “design” aspect in a future publication.
In the following subsections, we test whether the two different DFTs introduced in Sec. III B are able to predict (in a quantitative manner) the values of the relevant wavevectors.
B. DFT-based stability analysis of the HCSS liquid
Having shown which wavevectors characterize the structure of the liquid and are also involved in the structure of the cluster phase and also shown that monitoring the values of S1 and S4 allows us to predict the emergence of ordered structures, we now move on to discuss how DFT fares in predicting these wavevectors. In Fig. 3, we present results for the negative reciprocal of the structure factor −1/S(k) and for the location in the phase diagram of the λ-line. As mentioned already, DFT predicts the liquid to be linearly unstable (i.e., within the λ-line) at temperatures higher than where the cluster formation actually occurs. This is similar to the case with the SALR class of models, where the connection between the location of the λ-line and the cluster formation has been investigated in some detail.43,45,53–56 Figure 3(a) illustrates how the presence of the square-shoulder interaction modifies the structure factor of the pure HC system. Note that we plot −1/SFMT(k), the structure factor from the FMT–DFT, as a function of k rather than SFMT(k) itself because the exercise of investigating where and how peaks in S(k) become large (and within the DFT treatment diverge S(ki) → ∞) is entirely equivalent to investigating where and how peaks in −1/SFMT(k) approach zero from below.
(a) Plots of the negative reciprocal of the structure factor −1/SFMT(k) as a function of k, for different values of the packing fraction η (as labeled). The solid lines are for the ϵ = 0 pure HC system as obtained via FMT, from Eq. (23). The dotted and dashed lines are related results for the HCSS-system, with ϵ = 1 and λ = 1.3 (dotted) or λ = 3.7 (dashed), for the same range of values of η. The black vertical lines indicate the wavenumber for HC contact kcon = 2π (dotted), the wavenumber associated with closest hexagonal packing (dashed), and the predicted critical wavenumber for freezing of the pure HC fluid kf = 6.616 (full). (b) Scaled temperature location of the λ-line as functions of η, calculated numerically within FMT (solid lines) and within LDA (dashed line) for different values of λ (as labeled). For clarity, the curves were multiplied by a factor λ−2 in order to scale out the shoulder dependency of the LDA results—see Eq. (26). Note that when displayed this way, the LDA λ-line corresponds to the same curve for all values of λ. The horizontal dotted lines (and crosses) indicate the temperatures of the GCMC simulations carried out to investigate specific liquid states (see discussion in Sec. IV C). (c) Critical wavenumbers kc as functions of η, calculated within FMT (solid lines) and within LDA (dashed line) for different values of λ, as labeled. For clarity, the curves were re-scaled by the density-independent, critical wavevector predicted by the LDA kc = j2,1/λ—see Eq. (25). The black solid vertical line [which extends over both (b) and (c)] indicates the packing fraction at which entropic freezing is predicted by the FMT for the pure HC fluid ηf = 0.744.
(a) Plots of the negative reciprocal of the structure factor −1/SFMT(k) as a function of k, for different values of the packing fraction η (as labeled). The solid lines are for the ϵ = 0 pure HC system as obtained via FMT, from Eq. (23). The dotted and dashed lines are related results for the HCSS-system, with ϵ = 1 and λ = 1.3 (dotted) or λ = 3.7 (dashed), for the same range of values of η. The black vertical lines indicate the wavenumber for HC contact kcon = 2π (dotted), the wavenumber associated with closest hexagonal packing (dashed), and the predicted critical wavenumber for freezing of the pure HC fluid kf = 6.616 (full). (b) Scaled temperature location of the λ-line as functions of η, calculated numerically within FMT (solid lines) and within LDA (dashed line) for different values of λ (as labeled). For clarity, the curves were multiplied by a factor λ−2 in order to scale out the shoulder dependency of the LDA results—see Eq. (26). Note that when displayed this way, the LDA λ-line corresponds to the same curve for all values of λ. The horizontal dotted lines (and crosses) indicate the temperatures of the GCMC simulations carried out to investigate specific liquid states (see discussion in Sec. IV C). (c) Critical wavenumbers kc as functions of η, calculated within FMT (solid lines) and within LDA (dashed line) for different values of λ, as labeled. For clarity, the curves were re-scaled by the density-independent, critical wavevector predicted by the LDA kc = j2,1/λ—see Eq. (25). The black solid vertical line [which extends over both (b) and (c)] indicates the packing fraction at which entropic freezing is predicted by the FMT for the pure HC fluid ηf = 0.744.
In Fig. 3(a), we compare −1/SFMT(k) for the pure hard-core (HC) liquid with that for the HCSS liquid with T = 1 (i.e., ϵ = 1) and either λ = 1.3 or λ = 3.7. The pure HC system becomes marginally unstable at the packing fraction ηf ≈ 0.744, i.e., where −1/SFMT(kf) = 0, with a critical freezing kf ≈ 6.616. This freezing is purely entropic in nature. At very high densities, this peak is located at , representing the wavenumber that is responsible for a closest hexagonal packing, occurring at . Note that the packing fraction at which the freezing transition occurs as predicted by minimization of the FMT functional in a 2D domain is at packing fractions somewhat lower than this: The FMT predicts the liquid–solid coexistence packing fractions to be ηliq = 0.711 and ηsol = 0.732,57 for the liquid and solid, respectively. For 0 < η < ηf, the liquid is linearly stable, but only for η ≤ ηliq, is it the global free energy minimum and therefore the thermodynamic equilibrium state.
The location of the peak in S(k) for the high density liquid predicts fairly well the wavenumber determining the length scale of the particle ordering in the pure HC crystal. Since this freezing is entropic in origin (in contrast to the cluster freezing that is driven by the energetics of the shoulder part of the potential), we refer to this wavenumber as the entropically favored wavenumber in our discussions here. For the HCSS-liquid (with results shown for λ = 1.3 and λ = 3.7 in Fig. 3), we observe a structure factor that shows a similar behavior to that of the pure HC system for k ≈ kf. However, at smaller k, the smooth curve of SHC(k) is superposed by modulations with maxima occurring at , with integer n = 1, 2, 3, … [a better approximation for the locations of these maxima comes from the locations of the maxima in Eq. (20), i.e., Eq. (25)]. Each of these peaks correspond to wavenumbers of density modes that may or may not be favored by the nonlinear couplings in the free energy functional. At intermediate to high densities, the FMT-based HC contribution of the excess free energy functional dominates, leading to a favoring of modes in the vicinity of the peak of the structure factor of the pure HC fluid. However, at lower densities, modes at the other peaks can be favored. For very high densities (i.e., for η > 0.7), the modulations in −1/SFMT(k) are no longer visible and only one single peak, which is associated with entropic freezing, emerges in the region kcon < k < khex, where kcon = 2π is the wavevector associated with core contact. Hence—and as already discussed in the context of Eq. (25)—longer ranged shoulder interactions lead to the introduction of additional peaks in the structure factor and, consequently, other preferred wavevectors that can potentially play a role in the formation of structured phases. We discuss further the role of these additional wavevectors in the liquid phase in Sec. IV C.
In Fig. 3(b), we display the location in the phase diagram of the λ-line as a function of η, for different values of λ. By virtue of Eq. (26), the λ-line scales within the LDA–DFT treatment as T ∝ λ2. Therefore, we display the λ-line in a scaled form, namely, by plotting in the scaled temperature vs η plane. This helps us to compare systems with different shoulder widths λ in a more transparent manner. For small values of η, the LDA–DFT and FMT–DFT results are in a reasonably good agreement, but rather pronounced differences between the two can be observed at intermediate and high densities. These discrepancies are largest for the smallest λ-value investigated, while for the larger λ’s, the FMT-based curves essentially converge toward the LDA–DFT results. Eventually, as the shoulder becomes very short ranged, i.e., as λ approaches unity, the first peak in the static structure factor is roughly located at kc ≈ j2,1 ≈ 5.135—see Eq. (25)—where it is close to the (η-dependent) position of the main peak of the HC structure factor at intermediate densities (see Fig. 3). This leads to a strong η-dependence of both the λ-line and kc for smaller λ-values. The predicted freezing packing fraction ηf of a HC fluid is indicated by the solid vertical line.
The value of the critical wavevector kc related to the λ-line is plotted as a function of η in Fig. 3(c), for various values of λ. To highlight the differences between the FMT- and the LDA-results for this quantity, in this plot, we have divided the FMT kc by the corresponding LDA result, given in Eq. (25). Again, we observe that for larger λ-values (i.e., for λ > 2), kc only weakly depends on η and the respective FMT-results are very close to the data originating from the LDA. For the smallest λ-value investigated, the results from the FMT deviate strongly from the LDA-results. Furthermore, we observe a pronounced η-dependence of kc. The origin of this η-dependency lies in the interplay between the purely entropic HC and the energetic soft shoulder contributions to the excess free energy.
To shed further light on how repulsive shoulder affects the HC freezing (with corresponding associated wavenumber kf), it is instructive to take a closer look on the DFT-predictions for the structure of the HCSS system in comparison to the pure HC liquid. In Fig. 4(a), we present a heatmap plot of −1/SFMT(k) for the pure HC system (calculated with FMT) as a function of both η and k. Cuts through this are also displayed as line plots in Fig. 3(a). We have chosen a cutoff of −1/S(k) ≥ −4 to emphasize the dominant wavevectors, which at high η determine the structure (lattice spacing) of the pure HC crystal, i.e., which are “entropically favored.” In Fig. 4(a), the red solid line denotes the position of the first maximum of S(k), while the dashed red lines mark the inflection points of this function, thereby defining the width of the peak in S(k) as the distance between these. The blue star marks the η-value for freezing (ηf = 0.744—blue vertical line) of the system into a hexagonal crystal, which is governed by kf = 6.616 (blue horizontal line). The HC liquid is marginally stable for kf ≈ 6.616 at ηf. As the density is increased, the distances in the hexagonal crystal shrink (with a corresponding increase of kf) until particles reach the closest hexagonal packing, occurring at η ≈ 0.906 with (see the red star in Fig. 4). We should emphasize that for the states bounded by the green curve in Fig. 4(a), the structure factor of the fluid attains nonphysical negative values. Of course, this is a sign that the system is not a fluid at these state points and is in fact in a solid state. Note in Fig. 4(a) that even at intermediate values of the packing fraction (i.e., for η ≈ 0.5), the range of favored wavevectors narrows down, approaching a fairly narrow range in the vicinity of the maximum (red solid line).
(a) A heatmap plot of the negative reciprocal of the structure factor −1/S(k) for the pure HC fluid as a function of η and k, as obtained from FMT. We employ a cutoff in the range displayed, −1/S(k) ≥ −4, in order to emphasize those wavevectors that predominantly govern the density modes in the pure HC fluid and determine the structure of the HC crystal. The positions of the first peak of −1/S(k) [corresponding to the first peak in S(k)] are indicated by the red solid line. The red dashed lines mark the inflection points of the first peak of S(k), i.e., where . The distance between these dashed lines represents a measure for the width of the peak of S(k). The blue star indicates the point in the (η, k)-plane where the HC liquid becomes marginally stable against kf, i.e., where −1/S(kf) → 0 (see text). The red star marks the closest hexagonal packing at η = 0.906 and . The green lines indicate the locus of the points where −1/S(k) = 0. (b) A plot of the freezing kf as a function of λ to show the effect of the shoulder interaction (in terms of the shoulder width λ) (black and left vertical axis) together with a plot of the freezing packing fraction ηf(λ) (blue and right vertical axis). These are for the scaled shoulder heights ϵ* = ϵλ [see Eq. (24)] of ϵ* = 0.2 (solid), ϵ* = 0.4 (dashed), ϵ* = 0.6 (dotted), and ϵ* = 0.8 (dashed-dotted). The black horizontal line indicates the pure HC case. If we assume that at the freezing transition the particles form a perfect hexagonal lattice (governed by kf), one can derive an ideal packing fraction of such a lattice, given by . The relative error between the predicted ηf and this “ideal” value ηid(kf) is plotted in (c) for the different ϵ* values specified above.
(a) A heatmap plot of the negative reciprocal of the structure factor −1/S(k) for the pure HC fluid as a function of η and k, as obtained from FMT. We employ a cutoff in the range displayed, −1/S(k) ≥ −4, in order to emphasize those wavevectors that predominantly govern the density modes in the pure HC fluid and determine the structure of the HC crystal. The positions of the first peak of −1/S(k) [corresponding to the first peak in S(k)] are indicated by the red solid line. The red dashed lines mark the inflection points of the first peak of S(k), i.e., where . The distance between these dashed lines represents a measure for the width of the peak of S(k). The blue star indicates the point in the (η, k)-plane where the HC liquid becomes marginally stable against kf, i.e., where −1/S(kf) → 0 (see text). The red star marks the closest hexagonal packing at η = 0.906 and . The green lines indicate the locus of the points where −1/S(k) = 0. (b) A plot of the freezing kf as a function of λ to show the effect of the shoulder interaction (in terms of the shoulder width λ) (black and left vertical axis) together with a plot of the freezing packing fraction ηf(λ) (blue and right vertical axis). These are for the scaled shoulder heights ϵ* = ϵλ [see Eq. (24)] of ϵ* = 0.2 (solid), ϵ* = 0.4 (dashed), ϵ* = 0.6 (dotted), and ϵ* = 0.8 (dashed-dotted). The black horizontal line indicates the pure HC case. If we assume that at the freezing transition the particles form a perfect hexagonal lattice (governed by kf), one can derive an ideal packing fraction of such a lattice, given by . The relative error between the predicted ηf and this “ideal” value ηid(kf) is plotted in (c) for the different ϵ* values specified above.
Having considered the pure HC system, we now move on to discuss the ϵ > 0 case. Additional wavevectors are introduced by the shoulder interactions. Those that lie in the peak region of the HC fluid are reinforced (amplified), especially at higher densities, while the smaller wavevectors are increasingly suppressed as the density is increased. For example, in the case where λ = 1.3 discussed above, kc ≈ j2,1/1.3 ≈ 3.95 so that it is in the peak region of the HC fluid up to η ≈ 0.6. This is why this is also roughly the position of the maximum in kc(η), displayed in Fig. 3(c).
In Fig. 4(b), we illustrate the impact of the square-shoulder interaction on the freezing transition. This plot shows how kf and ηf change as λ is varied. We compare the behaviors for different values of the scaled shoulder height ϵ* = λϵ, showing results for ϵ* ∈ {0.2, 0.4, 0.6, 0.8}. Scaling in this way leads to comparable amplitudes in the shoulder contributions to S(k)—see the combination of coefficients in the second term of Eq. (24) and also Fig. 3. For large η (i.e., η > 0.8), only one (freezing) peak is observable in the HCSS S(k), as long as ϵ remains reasonably small. The freezing wavevector kf displayed in Fig. 4(b) is determined by increasing η until the global maximum in −1/S(k) at this value reaches zero, i.e., satisfying −1/SFMT(kf) = 0. Of course, this also corresponds to the peak of the largest maximum of S(k) diverging, indicating that the uniform liquid is unstable, and so in reality, freezing has already occurred at a lower density. Nonetheless, it is instructive to note the value of kf at which S(k) diverges and to examine the behavior with varying λ. We observe that both the predicted freezing kf as well as the freezing packing fraction ηf are functions that oscillate around the HC value (shown as a black, horizontal line) as λ increases. Of all the maxima in −1/S(k), the one that determines kf depends on λ. The role of which ki (i.e., which maxima of S(k)) becomes the freezing wavenumber kf is passed from one to the next roughly at the λ values corresponding to the position of the maxima in the oscillations of kf(λ). For example, at λ = 3.7, it is the wavevector of the fourth maximum, i.e., kf = k4.
Our MC simulations show that, to a fairly good approximation, kf governs the hexagonal lattice of hard cores that emerges when the liquid temperature is decreased at high densities. However, there is a difference between the “ideal” packing fraction of the hexagonal lattice with wavevector kf(λ) (i.e., assuming that every vertex in this hexagonal lattice is occupied by a particle) and the packing fraction ηf(λ) where −1/S(kf) → 0. Moreover, this difference has a consequence. The ideal hexagonal lattice with this wavevector has packing fraction ηid, determined as the ratio of the area occupied by a single particle and the area of the unit cell, giving . We display in Fig. 4(c) the relative difference between ηid(kf) and ηf. We can interpret positive deviations of this quantity as an indicator of the number of vacancy defects in the crystal, which Fig. 4(c) indicates can be rather large. However, this comparison still does not tell us where the freezing transition occurs. This analysis indicates that freezing must occur at a density value lower than that where −1/S(kf) → 0, but the precise density value comes from calculating the grand potential Ω of the liquid and crystalline phases and identifying phase coexistence in the usual way.26
C. The structure of the liquid state
In Secs. IV A and IV B, we showed that the knowledge of the significant wavevectors governing the structure of the HCSS liquid, ki, i = 1, …, 4, provides valuable information about which structures can potentially emerge when the liquid becomes unstable and solid phases form. Furthermore, we showed that the shoulder interactions lead to the introduction of additional wavevectors, beyond those of the pure HC system. The relevance of these is determined by their interplay with the core interactions, i.e., whether the nonlinear couplings in the free energy lead to a reinforcement of these modes or not. Density modes that are likely to be favored can be identified in plots of −1/S(k), such as that in Fig. 3. In this section, we now benchmark the FMT–RPA’s ability to accurately describe the peaks of S(k) on a quantitative level.
Thus, in Fig. 5, we compare S(k) obtained from GCMC simulations [via g(r) and the inverse of Eq. (18)] with the results obtained analytically via the FMT–RPA–DFT [Eqs. (17) and (23)] for three different values of λ (i.e., λ = 2.5, 3.7, and 4.9) and four different values of η (i.e., η = 0.1, 0.2, 0.3, and 0.4). Note the semi-logarithmic axes in Fig. 5. We have chosen liquid states that are located just above the predicted λ-line (see Fig. 3), i.e., at a temperature T = 1.05 × Tmax(λ), where Tmax(λ) is the temperature for which the λ-line attains its maximum. The state-points at which these are calculated are marked in Fig. 3(b).
Static structure factor S(k) as a function of wavenumber k, as obtained in GCMC simulations (top panels) and from the FMT–RPA–DFT (bottom row), for different values of shoulder width, λ, and packing fractions, η [(a) λ = 2.5, (b) λ = 3.7, and (c) λ = 4.9]. Note the semi-logarithmic scales. The temperatures T were chosen to be slightly above the maximum temperature of the corresponding λ-line, Tmax(λ): T = 1.05 × Tmax(λ)—see the points marked in Fig. 3(b). Hence, these states are all above the λ-line, but lie very close for η = 0.2 and η = 0.3. For better visibility, the maxima Si of S(k) are marked by symbols (crosses, circles, and triangles). The two insets show the relative errors between the simulation and DFT results for the peak positions, ki, and peak heights, Si, as functions of λ. In these insets, the same symbols have been used for the different Si as in the plots of S(k). These differences are defined as Δkrel = (kFMT − kSIM)/kSIM and ΔSrel = (SFMT − SSIM)/SSIM, respectively. A diverging second peak erroneously predicted by the FMT–RPA–DFT for λ = 2.5 and η = 0.4 (see dotted vertical lines) has been exempted from these considerations.
Static structure factor S(k) as a function of wavenumber k, as obtained in GCMC simulations (top panels) and from the FMT–RPA–DFT (bottom row), for different values of shoulder width, λ, and packing fractions, η [(a) λ = 2.5, (b) λ = 3.7, and (c) λ = 4.9]. Note the semi-logarithmic scales. The temperatures T were chosen to be slightly above the maximum temperature of the corresponding λ-line, Tmax(λ): T = 1.05 × Tmax(λ)—see the points marked in Fig. 3(b). Hence, these states are all above the λ-line, but lie very close for η = 0.2 and η = 0.3. For better visibility, the maxima Si of S(k) are marked by symbols (crosses, circles, and triangles). The two insets show the relative errors between the simulation and DFT results for the peak positions, ki, and peak heights, Si, as functions of λ. In these insets, the same symbols have been used for the different Si as in the plots of S(k). These differences are defined as Δkrel = (kFMT − kSIM)/kSIM and ΔSrel = (SFMT − SSIM)/SSIM, respectively. A diverging second peak erroneously predicted by the FMT–RPA–DFT for λ = 2.5 and η = 0.4 (see dotted vertical lines) has been exempted from these considerations.
The first thing that we can observe in Fig. 5 is that the number and positions of the individual peaks of S(k) (i.e., ki) are in all cases predicted quite accurately by the FMT–RPA–DFT. A more quantitative measure is displayed in the central inset, where we display the relative differences between the two sets of ki, i.e., we plot Δkrel = (kFMT − kSIM)/kSIM, for the three different values of λ. We find that the values of ki are in general slightly underestimated by the DFT, with the largest error being for small λ-values and higher η-values, notably for η = 0.4 at λ = 2.5.
We also see in Fig. 5 that the heights of the peaks Si are throughout overestimated by the DFT, in some cases only by a little and in other cases by a considerable amount. In the right-hand inset of Fig. 5, we display the relative differences between the two sets of peak heights Si, i.e., we plot ΔSrel = (SFMT − SSIM)/SSIM, for the three different values of λ. We see that the FMT–RPA–DFT greatly overestimates S1 for states in the vicinity of the λ-line, but is relatively accurate for the remaining Si. This is particularly the case for η = 0.2 and η = 0.3, in the vicinity of the λ-line. In contrast to the DFT, the simulations show a tendency to suppress the first peak with increasing η. For λ = 2.5 and η = 0.4, the DFT predicts a very large (in fact diverging) value of S2, indicating that the liquid is linearly unstable against k2. Note that this k2 neither corresponds to kc nor to kf. What is also interesting about this peak is that for intermediate densities and for λ ≳ 2.4, it points to the possibility of the formation of structures with a characteristic wavelength that is associated neither with the core-length scale nor with the shoulder-length scale. This observation tallies with those in Ref. 13, where low-density crystals were identified that are compatible with these observations.
Despite the DFT prediction of structured phases inside the λ-line (e.g., for λ = 2.5 and η = 0.4), all the states investigated within our simulations at the temperatures considered in Fig. 5 are, in fact, in the liquid phase. This can be seen from the corresponding simulation snapshots displayed in Fig. 6. Nonetheless, these snapshots do show significant structuring in the liquid. The characteristic wavelengths of these patterns are predicted accurately by the DFT because the DFT predictions for the peak positions in S(k) are correct. It is at lower temperatures where the system finally freezes into a variety of different solid phases. In other words, the DFT predicts the transition to the structured phases at temperatures higher than they are observed in MC simulations. However, we should emphasize that even within the DFT description, one must expect that coexistence between the liquid and ordered phases to occur before the structure factor diverges. Thus, the DFT is self-consistent, with the structure factor not actually diverging at the transition—except perhaps at the top of the λ-line itself.29,43,58
Typical snapshots from GCMC simulations of all the states whose structure factors S(k) are shown in Fig. 5. All systems are in the liquid phase. Snapshots are shown for λ = 2.5 (top row), λ = 3.7 (central row), and λ = 4.9 (bottom row). The packing fractions are η = 0.1, 0.2, 0.3, and 0.4 (from left to right).
Typical snapshots from GCMC simulations of all the states whose structure factors S(k) are shown in Fig. 5. All systems are in the liquid phase. Snapshots are shown for λ = 2.5 (top row), λ = 3.7 (central row), and λ = 4.9 (bottom row). The packing fractions are η = 0.1, 0.2, 0.3, and 0.4 (from left to right).
In the following section, we discuss the simplest and most commonly observed of these structured phases, namely, the cluster, stripe, and hole phases. A variety of other complex ordered phases can also arise, but a detailed investigation of these will be presented elsewhere. For higher temperature state points further away from the λ-line, the agreement between the DFT and simulations for the values of Si is considerably better.
D. Structured phases
So far, we have largely used the DFT to determine the liquid state structure and to show that it predicts fairly well the key wavenumbers ki. In this section, we now discuss the DFT predictions for properties of structured phases arising in those parts of the phase diagram where both the DFT and the simulations predict structured phases. We focus in particular on comparing wavelengths of density modulations and assessing how the DFT fares in predicting these, in comparison with the simulations. Comparing the equilibrium density profile ρ(r) obtained from full minimization of both the LDA–RPA and the FMT–RPA functionals with data obtained in GCMC simulations, we find that at higher temperatures, both DFTs overestimate the size of the regions in the phase diagram where the periodic phases exist, but overall are in good qualitative agreement with the MC results. The LDA–RPA–DFT predicts a phase diagram that is qualitatively very similar to that of the SALAR system displayed in Ref. 58. However, at lower temperatures, the LDA is unable to predict the freezing of the particles that occurs within the clusters and stripes of the modulated phases, nor can it describe the variety of other crystalline phases that can arise at lower temperatures and/or higher densities. As illustrated in Fig. 7, the FMT–RPA–DFT does fare better for the cluster, stripe, and hole phases. Additionally, at lower temperatures and/or higher densities, in addition to these three periodic phases, a variety of other crystalline phases also occur, depending on the value of λ. However, as mentioned, detailed investigation of these other phases will be presented elsewhere. Figure 7 shows some examples of density profiles obtained at lower densities, where we find good agreement between the FMT–RPA–DFT and the GCMC simulations, although there are some errors made by the DFT that we show are introduced by the underestimation of ki discussed in Sec. IV C.
Comparison of equilibrium density profiles ρ(r) for different HCSS-systems (specified below) obtained via the LDA–RPA–DFT (first row), the FMT–RPA–DFT (second row), and equal-sized snapshots taken from GCMC simulations of much larger sized systems (third row). The bottom panels show the corresponding structure factors S(k) as functions of k = (kx, ky) extracted from GCMC simulations. The white circles correspond to the peak positions of S(k) observed via the FMT–RPA–DFT for a bulk liquid at the same T and μ values. The columns show (a) a cluster phase for μ(η = 0.25) [see Eq. (15)] and T = 0.25; (b) a striped phase for μ(η = 0.35) and T = 0.25; and (c) a holed phase for μ(η = 0.45) and T = 0.38.
Comparison of equilibrium density profiles ρ(r) for different HCSS-systems (specified below) obtained via the LDA–RPA–DFT (first row), the FMT–RPA–DFT (second row), and equal-sized snapshots taken from GCMC simulations of much larger sized systems (third row). The bottom panels show the corresponding structure factors S(k) as functions of k = (kx, ky) extracted from GCMC simulations. The white circles correspond to the peak positions of S(k) observed via the FMT–RPA–DFT for a bulk liquid at the same T and μ values. The columns show (a) a cluster phase for μ(η = 0.25) [see Eq. (15)] and T = 0.25; (b) a striped phase for μ(η = 0.35) and T = 0.25; and (c) a holed phase for μ(η = 0.45) and T = 0.38.
As discussed in the previous subsection, due to the overestimation of the height of the first peak S1, the RPA-based DFT functionals tend to predict the formation of patterned structures at temperatures higher than they actually occur, i.e., within the λ-line, which is at too high temperatures. For instance, the clustering transitions discussed in Sec. IV A are found using GCMC simulations at η ≈ 0.2 to occur at Tc ≈ 0.36 (see Fig. 1), whereas the λ-line for the same η-value is found using the LDA at Tc = 0.6181, and from the FMT, it is at Tc = 0.6979. Similarly, at a higher density, using GEMC simulations, we find that the transition from the liquid to the cluster phase occurs at η ≈ 0.2827 and Tc ≈ 0.3663, whereas for this density, the LDA predicts Tc = 0.5892 and the FMT at Tc = 0.6829. Thus, to compare predictions from the DFTs with simulation results for the structured phases, we must choose state points at relatively low temperatures, where both the LDA–RPA–DFT and FMT–RPA–DFT reliably predict the emerging phases.
As discussed in Sec. IV C, at higher densities and for higher values of k, the FMT–RPA functional tends to be more accurate in its predictions for Si. Hence, in such situations, we can expect this functional to be more reliable for predicting the phase transitions associated with the corresponding ki. Despite the inaccuracy of the DFTs for the values of the peak heights S1, we find that the DFT still accurately predicts the equilibrium structures. This can be seen in Fig. 7, where we compare density profiles for the striped, clustered, and holed phases (other phases can also be observed in some cases; these will be presented elsewhere), together with corresponding simulation snapshots. In the top row, we display equilibrium density profiles ρ(r) obtained from the LDA–RPA–DFT. Note that for better visibility of the profiles in the low-density regions, we plot the logarithm of the density profile. Unsurprisingly, given that these results are from a functional that treats the hard-disk contributions via the simplest LDA, this DFT is only able to resolve the large-scale structuring, characterized by the length-scale ac, induced by the shoulder interactions, which are treated within the RPA.
The second row of Fig. 7 displays ρ(r) obtained from the FMT–RPA–DFT. These are obtained by using the LDA–RPA density profile as the initial guess in the numerical procedure for minimizing the FMT–RPA functional. Note that results calculated in this manner were found to be consistent with those obtained by starting from a random initial guess for the density profile. However, this “educated guess,” based on the LDA–RPA profile, was found out to be computationally more efficient since both density profiles share similar large-scale features, as they use the same RPA treatment of the shoulder. We see in Fig. 7 that the density-profile ρ(r) obtained from the FMT–RPA has the same stripe/cluster/hole interspacing as the LDA–RPA, but additionally resolves the hexagonal ordering formed by the particles being frozen within the sub-structures. At first glance, the FMT–RPA results in the second row appear somewhat “messy.” However, understanding of what these show can be obtained by comparing the FMT–RPA results with equal-sized snapshots extracted from GCMC simulations, which are displayed in the third row. Recall that DFT is a statistical mechanical theory, calculating the ensemble average density profile.25,26 Thus, comparing the second and third rows, one can see that the second row FMT–RPA profiles do indeed resemble ensemble averages of configurations such as those displayed in the third row. Additionally, from comparing the DFT profiles with the respective GCMC snapshots, we can also see that while the local hexagonal order, governed by k4, is very well predicted by the FMT–RPA, both the LDA–RPA and FMT–RPA overestimate the larger length-scale associated with k1 = kc. This can most easily be seen by comparing the spacing between stripes in the three upper panels of Fig. 7(b).
The average packing fractions of the structures displayed in Fig. 7, obtained by the DFT and also from the GCMC simulations using the same chemical potential values μ(ηb) [see Eq. (15)], are given in Table I. We see that the LDA has the tendency to slightly underestimate the densities of the structured phases and to predict their occurrence at even lower densities, as compared to a bulk liquid with the same value of μ. In contrast, the FMT predicts an increase of the density compared to the reference unstable liquid with imposed packing fraction ηb. This is consistent with, but slightly less than the values obtained from the GCMC simulations. Note that for simulations carried out at constant μ, the increase of the density relative to the density of a bulk liquid with equal μ can be used as an indicator for the formation of a structured phase. This can also be seen in the Appendix, where a plot showing a comparison of the bond order parameters and of η for GCMC simulations is displayed.
Phase . | ηb . | LDA–RPA . | FMT–RPA . | GCMC . |
---|---|---|---|---|
Clusters | 0.25 | 0.2499 | 0.2517 | 0.2883 |
Stripes | 0.35 | 0.3466 | 0.3636 | 0.3905 |
Holes | 0.45 | 0.4441 | 0.4640 | 0.4875 |
Phase . | ηb . | LDA–RPA . | FMT–RPA . | GCMC . |
---|---|---|---|---|
Clusters | 0.25 | 0.2499 | 0.2517 | 0.2883 |
Stripes | 0.35 | 0.3466 | 0.3636 | 0.3905 |
Holes | 0.45 | 0.4441 | 0.4640 | 0.4875 |
In the bottom panels of Fig. 7, we display the structure factors S(k) obtained from the GCMC simulations. Similar to in Fig. 2, we have added on top of these heatmap plots the relevant significant wavevectors as predicted by the FMT–RPA–DFT, obtained for a reference bulk liquid with the same T and μ. We find that—even though the structured phases are found well below their freezing temperatures—the wavenumbers extracted from the liquid theory still agree very well with the Bragg-peaks in S(k), especially for k1 and k4 (similar to in Fig. 2), but slightly underestimating the other wavevectors, in accordance with observations made in Sec. IV C. Note that the full simulation box from which these S(k) are obtained contains several distinctly oriented domains, leading to some overlay in the resulting Bragg-peaks. Nonetheless, we observe that the major discrete peaks in the GCMC S(k) are located very closely to the predicted kc = k1 and kf = k4 in SFMT(k).
Regarding both DFTs, the slight underestimation of the smaller wavenumber k1 = kc leads to an overestimation of the interspacing between stripes/clusters/holes. While the error made for k4 has seemingly minimal effects on the predicted equilibrium density profile, the error in k1 results in an overestimation of ac, hence the discrepancy in the spacings of the different structures in Fig. 7 when comparing DFT and simulation. Although the error in k1 is relatively small, due to the reciprocal relation between real- and Fourier-space, any small error δk, made in the determination of a wavevector k, leads to an error in the related length-scale of the real-space equilibrium density profile, which can be much more noticeable. In other words, errors made for smaller k lead to more severe consequences than the ones for bigger wavevectors. Nonetheless, we conclude that we can use the liquid state DFT based theory to successfully predict the wavevectors shaping the various solid states.
V. CONCLUDING REMARKS
In this contribution, we have investigated in detail the pattern formation of a two-dimensional hard-core, square-shoulder (HCSS) system at the point where it becomes unstable with respect to the emerging ordered phases. The interaction potential acting between the particles is characterized by a circular, impenetrable core with an adjacent, repulsive square shoulder interaction; the latter is characterized by the step height (which can be treated as the inverse of the system temperature) and the range of the interaction, λ. Investigations were carried out via intensive Monte Carlo simulations (in the grand-canonical—GCMC—and in the Gibbs ensemble—GEMC) and via classical density functional theory; in the latter case, two different approximations were implemented for the hard core part of the potential (based on the local density approximation—LDA—and on fundamental measure theory—FMT), while the shoulder part of the potential was included in a mean-field-type random phase approximation (RPA).
Our simulation-based investigations have provided evidence that the density modes that dominate fluctuations in the liquid state (in terms of relevant wavevectors) of the HCSS system are the ones that ultimately shape the adjacent (ordered) cluster phase, once the liquid becomes unstable. We have elucidated how the interplay of the characteristic length scales of the hard-core and soft-shoulder interactions gives rise to the dominance of specific density modes. We then have checked the capacities of the two different DFTs for predicting the relevant wavevectors; we find that they agree with the MC simulations reasonably well. Focusing on the liquid phase, we have compared the static structure factor as obtained directly from the simulations with data obtained via DFT: we observe that the positions of the peaks in the static structure factor agree very nicely, while larger deviations between simulation- and DFT-results are observed for the peak heights, in particular in the vicinity of the λ-line.
We have also investigated the structured phases and addressed the question: in what respect do the aforementioned differences between simulation and DFT results affect the density profiles predicted by DFT of the ordered phases? We conclude from our investigations that errors in the small-k regime have a larger impact than deviations occurring for larger wavevectors. However, we consider FMT to be sufficiently accurate to predict the significant wavevectors. Thus, it seems a valid approach to apply DFT to the liquid phase in order to reliably predict those wavevectors that characterize the emerging ordered phases. We conclude that the DFTs developed here can be used in a rather simple manner to navigate (with a relatively small computational cost) through the complex phase diagram of the HCSS system. In return, this feature might also enable us to design interaction potentials that give rise to specific combinations of density modes in emerging, self-assembled ordered phases.
It should be noted that observations similar to the ones here have been made for liquid metals.59 Using effective interatomic pair potentials and reliable liquid state theories60 to calculate the static structure factor of the liquid phase, good agreement to the well-known crystalline structures of the metals was observed.
ACKNOWLEDGMENTS
The computational results presented here were enabled via a generous share of CPU time, offered by the Vienna Scientific Cluster (VSC) under Project No. 71263. The authors thank Ms. Katrin Muck for her guidance related to the use of HPC. A.J.A. gratefully acknowledges support from the EPSRC under Grant No. EP/P015689/1.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Michael Wassermair: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Software (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). Gerhard Kahl: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal). Roland Roth: Conceptualization (equal); Supervision (equal); Writing – review & editing (equal). Andrew J. Archer: Conceptualization (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: FURTHER DETAILS
1. Numerical solution of the DFT—Picard algorithm
The numerical minimization of the DFT functionals used in this contribution via the iterative Picard algorithm31,37 allows for the numerical determination of the equilibrium density profile ρ(r). In the following, we briefly outline the procedure.
The mixing parameter α is chosen to be small enough to ensure convergence of the sequence of density profiles emerging from Eq. (A5). In practice, this value is chosen by trial; in the case of the LDA–RPA functional, a value of α = 10−3–10−2 was sufficient to obtain convergence depending on initial density ρb as well as on the value of the shoulder height ϵ. For the FMT–RPA functional, we implemented an algorithm proposed by Roth in Ref. 31 that evaluates the optimal choice of α in each step.
2. Structure factor
We also compared results for S(k) obtained via Eq. (A6) with the numerical calculation of S(k) via the inverse Fourier transform of Eq. (18), starting from the g(r) obtained in GCMC simulations (comparison not displayed), since the latter can be prone to numerical errors. This is because the two-dimensional Fourier transform requires special care when discretizing k, in view of the resolution in r-space. The results for S(k) obtained via the two routes were found to be in good agreement, except for in the (irrelevant for present purposes) small-k region, for k < k1.
3. Identification of the liquid–cluster transition temperature
Here, we provide evidence of the approximate value of the transition temperature Tc between the liquid and the ordered cluster phase, for the HCSS system with λ = 3.7 and chemical potential μ(η = 0.2), calculated via Eq. (15).
GCMC simulations have been performed, gathering data for the ensemble averaged packing fraction η and the bond orientational order parameter BOO(m = 12). The results are shown in Fig. 8. For temperatures above T ≈ 0.36, we observe a rather flat BOO(m = 12) and only moderately increasing average packing fraction ⟨η⟩ with decreasing temperature. In contrast, both curves show a pronounced and characteristic increase as we pass the T ≈ 0.36 threshold value from above. Thus, we identify T = Tc ≈ 0.36 as a rather reliable estimate for the transition temperature between the two phases for this value of μ. Note that the calculated packing fraction is slightly above the corresponding η-value of the bulk liquid calculated within SPT.
(a) GCMC results for ensemble averaged packing fraction ⟨η⟩ (blue symbols and line; left coordinate) and the BOO(m = 12) (with particles being in shoulder contact) as functions of temperature T; the system is characterized by μ(η = 0.2)—see Eq. (15)—and λ = 3.7. The temperature range considered covers the transition between the liquid and the ordered cluster phase. The black dotted line marks our estimate for the transition temperature Tc ≈ 0.36. (b) Evolution of the energy of the system (normalized by the energy of the initial configuration) as a function of MC moves as obtained in GCMC simulation of the ordered states depicted in Fig. 7.
(a) GCMC results for ensemble averaged packing fraction ⟨η⟩ (blue symbols and line; left coordinate) and the BOO(m = 12) (with particles being in shoulder contact) as functions of temperature T; the system is characterized by μ(η = 0.2)—see Eq. (15)—and λ = 3.7. The temperature range considered covers the transition between the liquid and the ordered cluster phase. The black dotted line marks our estimate for the transition temperature Tc ≈ 0.36. (b) Evolution of the energy of the system (normalized by the energy of the initial configuration) as a function of MC moves as obtained in GCMC simulation of the ordered states depicted in Fig. 7.
In Fig. 8(b), we display the evolution of the system energy, i.e., the number of shoulder overlaps, in the course of the GCMC simulation. The energy is divided by the energy of the initial hexagonal configuration of particles (see Sec. III A for details). After the rapid initial decrease (or increase in the case of the holed phase at ηb = 0.45), we observe equilibration of all states toward constant values. It is hence safe to assume that our solid state simulations are equilibrated after about 1010 moves. For the ensemble averages, states from the end of the simulation are used.
4. Liquid–cluster phase coexistence in the Gibbs-ensemble
The extent of the liquid–cluster coexistence region can be obtained by calculating the average densities in the two simulation boxes used in the GEMC scheme. Thus, we calculate the probability distributions P(η) for each of the boxes to have the average packing fraction η. Having these P(η) form two distinct Gaussian distributions can be seen as a sign of proficient sampling of the Gibbs ensemble. The GEMC simulation of the state at η = 2.827 and T = 0.3663 involving N = 3000 particles was run for 8 × 109 MC-moves. The statistics were obtained using 3000 states from the end of the simulation, from after phase separation in boxes occurred. All sample states are separated by a single simulation block (see Sec. III A for details). P(η) and the corresponding densities of the liquid and clustered phase can be found in Fig. 9. Despite the small extent of the coexistence region, of width Δη = ⟨ηcl⟩ − ⟨ηliq⟩ = 0.0126, the GEMC algorithm is capable of directly identifying phase coexistence.
Probability to find the two individual simulation boxes of the GEMC simulations (see Fig. 2) at a specific average packing fraction η. The solid line denotes the (constant) packing fraction of the full Gibbs ensemble (η = 0.2827) and the dashed line denotes the average packing fraction in the individual boxes, respectively. Here, “Box 1” corresponds to the simulation box hosting the cluster phase (⟨ηcl⟩ = 0.2834), while “Box 2” hosts the liquid phase (⟨ηliq⟩ = 0.2821).
Probability to find the two individual simulation boxes of the GEMC simulations (see Fig. 2) at a specific average packing fraction η. The solid line denotes the (constant) packing fraction of the full Gibbs ensemble (η = 0.2827) and the dashed line denotes the average packing fraction in the individual boxes, respectively. Here, “Box 1” corresponds to the simulation box hosting the cluster phase (⟨ηcl⟩ = 0.2834), while “Box 2” hosts the liquid phase (⟨ηliq⟩ = 0.2821).