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.

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.

In our two-dimensional system, the particles interact via the spherically symmetric hard-core, square-shoulder (HCSS) potential, which can be split up into the hard disk interaction (with diameter σ) and an adjacent soft, repulsive, square shoulder interaction corona of diameter λσ and an interaction strength ϵ > 0. The pair potential Φ(r) thus reads as follows:
(1)
Throughout this work, we use the hard-core diameter σ as the unit length of the system (i.e., σ ≡ 1). Furthermore, we introduce the temperature T via β = 1/(kBT), with kB being the Boltzmann constant. The properties of the system solely depend on the dimensionless combination βϵ, so for simplicity, we henceforth set β ≡ 1 and refer to the reciprocal shoulder height as the temperature, i.e., T = 1/ϵ. We introduce the average bulk density ρb with its dimensionless counter-part ρb=ρbσ2. Furthermore, we define the packing fraction η = πρbσ2/4.

We analyze our system with both extensive MC based simulations23,24 and classical DFT based calculations.25,26 The essential conceptual elements of both approaches are summarized in this section.

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 AIV C, whereas we use rectangular boxes with aspect ratio 1:3 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.

1. DFT formalism and thermodynamics

The thermodynamic and structural properties of our system can alternatively be calculated within the framework of classical DFT. This approach is based on the grand potential functional Ω[ρ], which in the absence of any external potentials reads as follows:25,26
(2)
Ω[ρ] is a functional of the one-body density profile, ρ(r), and is formed from a Legendre transform of the Helmholtz free energy functional F[ρ], where μ is the chemical potential of the system.

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.

The Helmholtz free energy can be split into two contributions as follows:
(3)
The first term is the (exact) ideal gas contribution, which can be written as25,26
(4)
where Λ is the thermal de Broglie wavelength. The second term in Eq. (3) is the excess Helmholtz free energy functional Fex[ρ], which incorporates all contributions beyond that of the ideal gas, due to the particle interactions. For the overwhelming majority of systems, an expression for Fex[ρ] can only be given in an approximate manner. To do this, we split the particle pair interaction potential in Eq. (1) into a hard core (index “hc”) and a square shoulder (index “ss”) part as follows:
(5)
(6)
Thus, the full potential is Φ(r) = Φhc(r) + Φss(r). Correspondingly, we can split Fex[ρ] into a contribution due to the hard core and a remainder, associated with the shoulder as follows:
(7)
Furthermore, we assume that we can treat (in the sense of the random phase approximation—RPA26) the contribution due to the square shoulder potential as a perturbation to the hard core. The resulting RPA for Fexss[ρ] reads26 
(8)
Note that via Eqs. (6) and (8), we have followed a standard practice of extending the shoulder of the potential “inside” the hard core. This mean-field approximation is justified as long as Φss(r) is fairly long ranged and the shoulder height ϵ in Eq. (6) is reasonably small. However, in practice, the RPA can often do better than one might expect, even when these conditions are only loosely met.26,29,30

For the contribution to the free energy Fexhc[ρ], 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.

The LDA functional assumes that the density profile varies sufficiently slowly and that at every point r, where the density is ρ(r), we can approximate the local contribution to the free energy by that of a bulk liquid with corresponding bulk density ρ. This requires the corresponding bulk hard particle fluid equation-of-state (EOS); here, we use the one obtained within scaled particle theory (SPT).26,33,34 For the two-dimensional hard disk fluid, the LDA functional together with the SPT EOS is
(9)
where the local packing fraction η(r) = πρ(r)/4. Experience shows that the above functional fails when the density profile ρ(r) varies on length scales close to that of the core diameter, e.g., for fluids in the vicinity of hard walls or for crystalline states,35 but for clusters and aggregations of multiple particles, we can expect it to be qualitatively reliable.10,29,36 The results presented below show well when the LDA is sufficient and when it is not, in which case one must use the more sophisticated FMT discussed next.
The FMT functional for hard disks32 has been shown to successfully describe the entropic freezing of pure hard disks into a hexagonal crystal, with densely packed arrangement of these particles, having density profiles ρ(r) with sharp Gaussian-like peaks at the appropriate particle positions in the hexagonal lattice. In sharp contrast, unlike the FMT, the LDA is unable to predict the freezing of the pure hard-disk (or hard-sphere) fluid. The FMT formalism uses so-called weight functions ωα(r) and ω(m)(r), which by convolution with the one-particle density ρ(r) lead to weighted densities, that are either scalar
(10)
or tensorial in nature,
(11)
In the above relations, the symbol “⊗” represents spatial convolutions in two dimensions; the weight functions ωα(r) and ω(m)(r) are specified below.
One can then use these weighted densities to generate an approximation for the hard core contribution to the free energy Fexhc[ρ] as an integral over a function of these weighted densities,32 
(12)
The integrand of the above relation, i.e., the excess free energy density kBTΨ(r), is obtained from an ansatz, which is motivated by dimensional analysis, assuming that it must have dimension length−2 and so only has combinations of weight functions consistent with this. This ansatz is then combined with a set of differential equations for the bulk pressure, arising from SPT.31,32 The free coefficients in the solutions to these differential equations are determined by requiring that the resulting Helmholtz free energy yields the SPT EOS of for hard disks with bulk density ρ and is well behaved in certain limiting cases of extreme confinement. The resulting equation [Eq. (12)] is valid over a broad range of densities up to and beyond freezing and recovers the exact low-density limit.
The scalar weight functions introduced above are32 
(13)
where R = σ/2 = 1/2, δ(r) is the Dirac delta-function, and Θ(r) is the Heaviside step-function. The tensorial weight functions are given by
(14)
where the tensor-rank m of the weight function arises from (m − 1) tensor products of the unit vector r̂ with itself.
To calculate the equilibrium density profile ρ(r), we use numerical minimization of the DFT functionals via the iterative Picard algorithm described, e.g., in Refs. 31 and 37. Further details are also given in the  Appendix. Since both the LDA–RPA [Eqs. (9) and (8)] and the FMT–RPA [Eqs. (8) and (12)] are constructed so that they reproduce the SPT equation of state in the uniform density limit, for both DFTs, the chemical potential for a uniform bulk liquid is
(15)

2. Structure

Classical DFT not only provides thermodynamic properties but also gives structural (pair correlation) information about the system. One route to this is via the density profiles, typically using the test-particle method, where a single particle is fixed, treating it as an external potential. The DFT is then minimized to obtain the density profile of the surrounding fluid in the presence of this fixed particle. The resulting profile then directly yields the radial distribution function—see, e.g., Ref. 38—for further background on this approach. Here, we follow another approach, going via the direct correlation functions. To be more specific, the two-particle direct correlation function c(2)(r, r′) for an inhomogeneous fluid can be readily obtained as the following functional derivative of the excess Helmholtz free energy:26 
(16)
This function characterizes the two-point density correlations in the system due to the interactions between the particles. In the homogeneous bulk liquid where ρ(r) = ρb, this correlation function becomes isotropic, i.e., c(2)(r, r′) = c(2)(r = |rr′|) = c(r). The static structure factor, S(k), of such a homogeneous liquid is related to the Fourier transform (FT, indicated by a hat) of the direct pair correlation function as follows:
(17)
The static structure factor S(k) is in turn related to the radial distribution function g(r) as follows:26 
(18)
Due to the dimensionality of the problem at hand, the transformation from k- to r-space is realized via the inverse Hankel transform, with J0(x) being the Bessel function of the first kind of order 0.
Since the Helmholtz free energy functional has been split in an additive manner into the hard-core and the square-shoulder contribution—see Eq. (7)—a related splitting of the direct correlation function (in both r- and k-space) can be done by virtue of Eq. (16) as follows:
(19)
The simplicity of the last term arises due to the quadratic nature of FexRPA[ρ]—see Eq. (8). The required Fourier transform of the shoulder potential is
(20)
where J1(x) is the Bessel function of first kind of order 1.

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 ĉhc(k).

Using Eqs. (9) and (16), one finds in the case of the LDA the constant value
(21)
In the case of the FMT functional, the situation is more involved, where the direct correlation function is obtained in terms of derivatives of the excess Helmholtz free energy density Ψ(r) with respect to the weighted densities—see Eq. (12). The subsequent Fourier transform of the weight functions leads to39 
(22)
by which, together with Eq. (12), one eventually obtains
(23)

Since the static structure factor S(k) can easily be calculated from ĉ(k) via Eq. (17) and since in both approximations ĉ(k) is given by analytic expressions, one can—in principle—write down the expression for S(k).

In the LDA, this is straightforward, giving the following closed expression:
(24)
A related expression is found for SFMT(k), which is considerably more lengthy [in view of the complexity of ĉhcFMT(k), Eq. (23)] and so is not reproduced here.

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., ddkS(k)1k=kc=0, and (ii) the divergence condition, i.e., S(kc)10.

By virtue of the simplicity of the LDA expression for the static structure factor [see Eq. (24)], the extremum condition can be formulated in a closed expression,
(25)
where ji,n is the nth zero of the Bessel function of the first kind of ith order. Since kc is evaluated via an extremum condition, it is worth mentioning that for the LDA–RPA not only is the position of the main peak of S(k) given by the zeros of this specific Bessel function but also all the other maxima and minima. Consequently, the number of critical wavevectors that have to be considered in a stability analysis strongly depends on the interaction range λ since for larger λ, more local minima of S(k)−1 arise in the relevant interval 0 < k ≲ 2π. Here, we focus on the instability of the bulk liquid against the single wavenumber kc, corresponding to the global minimum of S(k)−1. However, we should mention that one can tune the system so that there are two (or more) critical wavenumbers, a topic we will elaborate elsewhere. Using the divergence conditions and inserting the critical wavenumber kc into Eq. (24), we obtain the following parametric curve for the λ-line:
(26)
The corresponding calculation for the FMT–RPA is not analytically tractable, but using numerical methods, the critical wavenumber kc and the λ-line can be obtained for the FMT–RPA, by following the general approach outlined above. Discussion and comparison of the LDA–RPA with the FMT–RPA results are made in Sec. III B 3.

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 khex=4π/37.26 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.

FIG. 1.

(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λ.

FIG. 1.

(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λ.

Close modal
FIG. 2.

(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 khex=4π/3 (see text). The small green circles in (d) highlight the Bragg-peaks associated with the hexagonal patterns observed between and within clusters (see text).

FIG. 2.

(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 khex=4π/3 (see text). The small green circles in (d) highlight the Bragg-peaks associated with the hexagonal patterns observed between and within clusters (see text).

Close modal

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 < kikhex (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.

In Fig. 1(c), we display the bond orientational parameters BOO(m) for a selected range of m-values. For 2D systems, these are defined as follows:51,52
(27)
where NN is the number of neighbors for a tagged particle at position r and θr is the angle between an arbitrary but fixed axis and the line connecting the particle at hand with its neighbors within a specified range. In evaluating BOO(m), we consider for a tagged particle only those neighbors that are (i) separated by distances r ≤ 1.5 (a distance which corresponds roughly to the position of the first peak in g(r), so these are the nearest-neighbors) or (ii) all particles with distances rλ from the tagged particle (i.e., up to shoulder contact). These two distances are marked in Fig. 1(c) by dashed and solid vertical lines, respectively. The average ⟨⋯ ⟩ in Eq. (27) indicates a statistical average performed over all particles in the system. In addition to the four temperatures considered in Figs. 1(a) and 1(b), in Fig. 1(c), we have also added the results for an even lower temperature state (i.e., for T = 0.250) for which the orientations of the particles within each cluster are fully aligned. In the case of nearest neighbors, the two peaks of the BOO(m) (for m = 6 and m = 12) become more pronounced as the temperature drops below Tc. In contrast, the transition from the liquid to the cluster phase can be associated with an increase in BOO(12) at shoulder contact. Note that due to the range of the shoulder interaction, the averaging procedure in the evaluation of BOO(m) for shoulder contact involves in general 10 to 20 particles.

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 khex=4π/3 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 BIV 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.

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.

FIG. 3.

(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 khex=4π/3 (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.

FIG. 3.

(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 khex=4π/3 (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.

Close modal

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 khex=4π/3, representing the wavenumber that is responsible for a closest hexagonal packing, occurring at ηhex=π/(23)0.907. 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 kkf. However, at smaller k, the smooth curve of SHC(k) is superposed by modulations with maxima occurring at 2πn/λ, 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 (ϵλ2)1 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 kcj2,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 kf=khex=4π/3 (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).

FIG. 4.

(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 d2S(k)dk2=0. 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 khex=4π/3. 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 ηid(kf)=3kf2/32π. The relative error between the predicted ηf and this “ideal” value ηid(kf) is plotted in (c) for the different ϵ* values specified above.

FIG. 4.

(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 d2S(k)dk2=0. 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 khex=4π/3. 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 ηid(kf)=3kf2/32π. The relative error between the predicted ηf and this “ideal” value ηid(kf) is plotted in (c) for the different ϵ* values specified above.

Close modal

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, kcj2,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 ηid(kf)=3kf2/(32π). 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 

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).

FIG. 5.

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 = (kFMTkSIM)/kSIM and ΔSrel = (SFMTSSIM)/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.

FIG. 5.

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 = (kFMTkSIM)/kSIM and ΔSrel = (SFMTSSIM)/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.

Close modal

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 = (kFMTkSIM)/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 = (SFMTSSIM)/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

FIG. 6.

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).

FIG. 6.

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).

Close modal

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.

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.

FIG. 7.

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.

FIG. 7.

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.

Close modal

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.

TABLE I.

Average packing fractions of the three investigated solid phases displayed in Fig. 7, obtained via DFT from GCMC simulations. For both the DFT and MC calculations, a fixed chemical potential μ(ηb) determined via Eq. (15) was imposed.

PhaseηbLDA–RPAFMT–RPAGCMC
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ηbLDA–RPAFMT–RPAGCMC
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 δa2πk+δk2πk 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.

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.

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.

The authors have no conflicts to disclose.

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).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

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.

Minimization of the grand-potential functional (2),
(A1)
leads (in the absence of an external potential) to the Euler–Lagrange equation,
(A2)
with c(1)(r)δF[ρ]ex/δρ(r) being the one-body direct correlation function.26 Evaluating the latter for a uniform liquid allows us to express this reference bulk density ρb in terms of the chemical potential μ,
(A3)
With this relation at hand, we can rearrange Eq. (A2) to obtain the following expression for the density profile:
(A4)
The index “P” indicates that this expression can be used in the iterative algorithm to obtain an update for the density profile at a given step, when the density profile of the preceding step is inserted into the right-hand side of Eq. (A4). It should be emphasized that ρb is only equal to the average density of the system ρ̄=ρ(r) for uniform (fluid) states. More generally, in the case of the non-uniform phases, the average density ρ̄ρb, although differences can be small. Thus, Eq. (A4) represents the basis for the iterative Picard algorithm as described, e.g., in Refs. 31 and 61.
To be more specific, we start from an initial guess for ρ(r), and then, a new density profile is constructed by mixing the density profile of the preceding step via a mixing parameter α with the density profile ρP(r), given Eq. (A4). For reasons of numerical stability, one must introduce a small value for the mixing parameter α. Thus, the density profile at step (i + 1) of the Picard iteration is given by
(A5)
This algorithm is iterated until the squared difference between the density profiles of two consecutive steps integrated over the domain, i.e., if the dimensionless quantity σ2 |ρ(i)(r) − ρ(i−1)(r)|2dr, drops below a predefined threshold value. The actual value depends on the system size, but typically amounts to (10−10–10−14). Simultaneously, the grand potential Ω[ρ] converges toward a constant value.

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

To calculate the static structure factor, in addition to using the method described in Sec. IV A, we also used the more standard definition,26 
(A6)
This relation requires as input the positions of all the particles {ri}, obtained in our MC simulations. The brackets denote an ensemble average in the related ensemble (here the grand-canonical or Gibbs ensemble).

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.

FIG. 8.

(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.

FIG. 8.

(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.

Close modal

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.

FIG. 9.

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).

FIG. 9.

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).

Close modal
1.
N.
Vogel
,
C.
Fernández-López
,
J.
Pérez-Juste
,
L. M.
Liz-Marzán
,
K.
Landfester
, and
C. K.
Weiss
,
Langmuir
28
,
8985
(
2012
).
2.
K.
Geisel
,
L.
Isa
, and
W.
Richtering
,
Angew. Chem.
126
,
5005
(
2014
).
3.
M.
Rey
,
M. Á.
Fernández-Rodríguez
,
M.
Steinacher
,
L.
Scheidegger
,
K.
Geisel
,
W.
Richtering
,
T. M.
Squires
, and
L.
Isa
,
Soft Matter
12
,
3545
(
2016
).
4.
K.
Volk
,
F.
Deißenbeck
,
S.
Mandal
,
H.
Löwen
, and
M.
Karg
,
Phys. Chem. Chem. Phys.
21
,
19153
(
2019
).
5.
M.
Rey
,
A. D.
Law
,
D. M. A.
Buzza
, and
N.
Vogel
,
J. Am. Chem. Soc.
139
,
17464
(
2017
).
6.
J.
Menath
,
J.
Eatson
,
R.
Brilmayer
,
A.
Andrieu-Brunsen
,
D. M. A.
Buzza
, and
N.
Vogel
,
Proc. Natl. Acad. Sci. U. S. A.
118
,
e2113394118
(
2021
).
7.
M.
Ickler
,
J.
Menath
,
L.
Holstein
,
M.
Rey
,
D. M. A.
Buzza
, and
N.
Vogel
,
Soft Matter
18
,
5585
(
2022
).
8.
G.
Malescio
and
G.
Pellicane
,
Nat. Mater.
2
,
97
(
2003
).
10.
M. A.
Glaser
,
G. M.
Grason
,
R. D.
Kamien
,
A.
Košmrlj
,
C. D.
Santangelo
, and
P.
Ziherl
,
Europhys. Lett.
78
,
46004
(
2007
).
11.
J.
Dobnikar
,
J.
Fornleitner
, and
G.
Kahl
,
J. Phys.: Condens. Matter
20
,
494220
(
2008
).
12.
J.
Fornleitner
and
G.
Kahl
,
Europhys. Lett.
82
,
18001
(
2008
).
13.
J.
Fornleitner
and
G.
Kahl
,
J. Phys.: Condens. Matter
22
,
104118
(
2010
).
14.
H.
Pattabhiraman
and
M.
Dijkstra
,
Soft Matter
13
,
4418
(
2017
).
15.
T.
Dotera
,
T.
Oshiro
, and
P.
Ziherl
,
Nature
506
,
208
(
2014
).
16.
P.
Ziherl
and
T.
Dotera
,
Proc. Int. Sch. Phys. “Enrico Fermi”
193
,
307
(
2016
).
17.
T.
Dotera
,
S.
Bekku
, and
P.
Ziherl
,
Nat. Mater.
16
,
987
(
2017
).
18.
P. C.
Hemmer
and
G.
Stell
,
Phys. Rev. Lett.
24
,
1284
(
1970
).
20.
E. A.
Jagla
,
J. Chem. Phys.
110
,
451
(
1999
).
21.
E. G.
Strekalova
,
D.
Corradini
,
M. G.
Mazza
,
S. V.
Buldyrev
,
P.
Gallo
,
G.
Franzese
, and
H. E.
Stanley
,
J. Biol. Phys.
38
,
97
(
2012
).
22.
W. R.
Somerville
,
A. D.
Law
,
M.
Rey
,
N.
Vogel
,
A. J.
Archer
, and
D. M. A.
Buzza
,
Soft Matter
16
,
3564
(
2020
).
23.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
2017
).
24.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
(
Elsevier
,
2023
).
26.
J. P.
Hansen
and
I. R.
McDonald
,
Theory of Simple Liquids
(
Academic Press
,
Oxford
,
2013
).
27.
R.
Evans
,
M.
Oettel
,
R.
Roth
, and
G.
Kahl
,
J. Phys.: Condens. Matter
28
,
240401
(
2016
).
28.
A. Z.
Panagiotopoulos
,
Mol. Phys.
61
,
813
(
1987
).
29.
A. J.
Archer
,
C.
Ionescu
,
D.
Pini
, and
L.
Reatto
,
J. Phys.: Condens. Matter
20
,
415106
(
2008
).
30.
A. J.
Archer
,
B.
Chacko
, and
R.
Evans
,
J. Chem. Phys.
147
,
034501
(
2017
).
31.
R.
Roth
,
J. Phys.: Condens. Matter
22
,
063102
(
2010
).
32.
R.
Roth
,
K.
Mecke
, and
M.
Oettel
,
J. Chem. Phys.
136
,
081101
(
2012
).
33.
E.
Helfand
,
H. L.
Frisch
, and
J. L.
Lebowitz
,
J. Chem. Phys.
34
,
1037
(
1961
).
34.
H.
Reiss
,
H. L.
Frisch
, and
J. L.
Lebowitz
,
J. Chem. Phys.
31
,
369
(
1959
).
35.
R.
Evans
, in
Fundamentals of Inhomogeneous Fluids
, edited by
D.
Henderson
(
Taylor & Francis
,
1992
), Chap. 3, pp.
85
175
.
36.
A. J.
Archer
,
A. M.
Rucklidge
, and
E.
Knobloch
,
Phys. Rev. Lett.
111
,
165501
(
2013
).
37.
A. P.
Hughes
,
U.
Thiele
, and
A. J.
Archer
,
Am. J. Phys.
82
,
1119
(
2014
).
38.
A. J.
Archer
and
R.
Evans
,
J. Chem. Phys.
118
,
9726
(
2003
).
39.
A. L.
Thorneywork
,
S. K.
Schnyder
,
D. G.
Aarts
,
J.
Horbach
,
R.
Roth
, and
R. P.
Dullens
,
Mol. Phys.
116
,
3245
(
2018
).
40.
U. M. B.
Marconi
and
P.
Tarazona
,
J. Chem. Phys.
110
,
8032
(
1999
).
41.
A. J.
Archer
and
R.
Evans
,
J. Chem. Phys.
121
,
4246
(
2004
).
42.
M.
te Vrugt
,
H.
Löwen
, and
R.
Wittkowski
,
Adv. Phys.
69
,
121
(
2020
).
43.
A. J.
Archer
,
D.
Pini
,
R.
Evans
, and
L.
Reatto
,
J. Chem. Phys.
126
,
014104
(
2007
).
44.
A. J.
Archer
,
M. J.
Robbins
,
U.
Thiele
, and
E.
Knobloch
,
Phys. Rev. E
86
,
031603
(
2012
).
45.
A.
Imperio
and
L.
Reatto
,
J. Phys.: Condens. Matter
16
,
S3769
(
2004
).
46.
J.-P.
Hansen
and
L.
Verlet
,
Phys. Rev.
184
,
151
(
1969
).
47.
L.
Costigliola
,
T. B.
Schrøder
, and
J. C.
Dyre
,
Phys. Chem. Chem. Phys.
18
,
14678
(
2016
).
48.
J. Q.
Broughton
,
G. H.
Gilmer
, and
J. D.
Weeks
,
Phys. Rev. B
25
,
4651
(
1982
).
49.
Z.
Wang
,
W.
Qi
,
Y.
Peng
,
A. M.
Alsayed
,
Y.
Chen
,
P.
Tong
, and
Y.
Han
,
J. Chem. Phys.
134
,
034506
(
2011
).
50.
G.
Jiménez-Serratos
,
C.
Vega
, and
A.
Gil-Villegas
,
J. Chem. Phys.
137
,
204104
(
2012
).
51.
K. J.
Strandburg
,
J. A.
Zollweg
, and
G. V.
Chester
,
Phys. Rev. B
30
,
2755
(
1984
).
52.
K. J.
Strandburg
,
Rev. Mod. Phys.
60
,
161
(
1988
).
53.
A.
Imperio
and
L.
Reatto
,
J. Chem. Phys.
124
,
164712
(
2006
).
54.
D.
Pini
,
A.
Parola
, and
L.
Reatto
,
J. Phys.: Condens. Matter
18
,
S2305
(
2006
).
55.
A. J.
Archer
and
N. B.
Wilding
,
Phys. Rev. E
76
,
031501
(
2007
).
56.
B.
Chacko
,
C.
Chalmers
, and
A. J.
Archer
,
J. Chem. Phys.
143
,
244904
(
2015
).
57.
S.-C.
Lin
and
M.
Oettel
,
Phys. Rev. E
98
,
012608
(
2018
).
58.
59.
J.
Hafner
and
G.
Kahl
,
J. Phys. F: Met. Phys.
14
,
2259
(
1984
).
60.
G.
Kahl
and
J.
Hafner
,
Phys. Rev. A
29
,
3310
(
1984
).
61.
A. J.
Archer
,
A. M.
Rucklidge
, and
E.
Knobloch
,
Phys. Rev. E
92
,
012324
(
2015
).
Published open access through an agreement with JISC Collections