Phase diagrams of simple models of colloidal nanocrystals in two dimensions

The self-assembly of colloidal nanocrystals at interfaces provides a bottom-up approach to create functional materials for developing next-generation flexible optoelectronic devices and sensors. In this work, we report phase diagrams of simple models of colloidal nanocrystals confined at a flat interface. By performing extensive computer simulations we elucidate the mesoscale organization that takes place as different parameters are varied. Our simulation results uncover rich phase diagrams where hexagonal, rhomboid, honeycomb and stripe phases as well as hierarchical self-assembly are found. Our results could serve as a guideline for experimentalists to design colloidal nanocrystal arrangements to target specific applications.


Introduction
Colloidal nanocrystals are a fascinating class of materials with many possible applications, from catalysis [1] and oil-recovery [2] to flexible optoelectronic devices [3,4]. Each colloidal particle is composed by an inorganic solid core coordinated at its surface by organic ligand shells, typically small polymer chains [5]. Once these colloids are synthesized, they can be dispersed into solutions that can be used to coat flexible or rigid substrates to form very thin films [6,7]. These materials were originally referred to as 'artificial atoms' because their physical properties can be tuned by adjusting the nanocrystal's composition, shape and size [8], thus, they are promising candidates for designer materials [9]. The ultimate goal is to achieve the hierarchical organization of these colloidal nanocrystals into complex arrangements by a bottom-up approach. The primary premise of this goal is that the complexity associated to the nanocrystals self-organization will be reflected in creating materials with remarkable new functionalities. For example, there is a lot of interest on creating hollow superstructures for applications in electrochemistry and energy storage [10][11][12], as well as vehicles in drug-delivery applications [13]. The potential of these hollow superstructures relies on the nanoparticle assembly into two-dimensional manifolds, and their organization in these subspaces. To fabricate these superstructures, researchers are using a wide range of strategies, going from the use of hard and soft templates, template free and self-templated methods [14]. However, a deep understanding of the important factors governing the assembly at the different length scales is lacking [14]. Unveiling these factors will lead to controlled design of hollow nanoparticle assemblies for advanced applications.
An important insight was the realization of the fundamental role played by the organic ligands on the collective organization of colloidal nanocrystals [15], as well as on the mechanical properties of the nanocrystal aggregates [16,17]. Thus, studies have been performed to elucidate the effect of different kinds of organic ligands, from DNA [18] to block copolymers [19], on the colloidal nanocrystals self-assembly. The versatility of organic molecular components that can be used by polymer chemists offers a rich platform to design these ligand shells. By carefully choosing these organic components, it has been possible to create nanocrystal assemblies, in the low concentration regime, with different shapes and geometries [19][20][21].
Indeed, solution self-assembly of colloidal nanocrystals has become a strong methodology to create ultra-thin nanocrystal superstructures by harnessing the effect of solvents on the organic ligands solubility behavior [22]. In fact, this has motivated the idea of creating stimuli-responsive nanocrystals assemblies for biomedical applications [23]. The self-assembly behavior and response of these materials is encoded into the surface ligands on the colloidal nanocrystals. These colloidal assemblies can be designed to respond to a variety of stimuli (temperature, pH, light, etc) [23]. It should be highlighted that not only the organic ligands play an important role on the self-organization of colloidal nanocrystals, but also the shape is a critical parameter. The associated symmetry of the nanocrystal can also be translated into the global symmetry of the lattices. Thus, for example, it has been possible to create lattices with tetragonal symmetry by using cuboidal nanocrystals [24][25][26]. Also, experimental work have found that it is possible to control the self-assembly of nanoplates into ordered superlattices by using dendrimer ligands [27]. Moreover, this work showed the importance of the grafting distribution on the ordering of these nanocrystals.
To achieve the full potential of these materials, the mechanisms governing the self-assembly process and response need to be better understood. A route to achieve this fundamental understanding is by performing sophisticated experiments complemented by modeling studies. For example, in Ref [28]. both experiments and simulations have been performed to provide physical insights on the effect of the shape of gold nanocrystals on the colloidal large-scale organization in the bulk. In optoelectronic applications, the type of ordering of these colloidal nanocrystals at the substrates is important, thus, tools that allow the prediction of the self-assembly at interfaces of these colloidal materials would be of great benefit to experimentalists. However, given the complexity of these materials, a first-principles description is incapable to address the time and length scales associated to the self-assembly process. Coarse-grained approaches, on the other hand, have the potential of providing physical insights about the role of structural/chemical parameters on the self-organization of these colloidal materials. In this regard, models that included polymer chains explicitly have been used to explore the self-assembly and mechanical properties of colloidal nanocrystals lattices [29,30]. Also, effective theories have been used to predict the self-assembly behavior of nanocrystal superlattices [31]. On the other hand, simple coarse-grained models rely on a few effective parameters representing the physical system of interests. These coarse-grained approaches are useful as the computational cost of exploring full phase diagrams is low, and could serve as zeroth-order approximations to the experimentally-relevant phase diagrams. These kinds of models have been successfully used to study the self-assembly behavior of colloidal nanocrystals both in bulk [32][33][34][35][36] as well as under confinement [37].
Historically, computer simulations started by performing calculations of simple model systems such as hard spheres (HS) in three dimensions (3D) or hard disks (HD) in two dimensions (2D) [38][39][40]. The reason behind this choice is the idea that microscopically-different systems might share same macroscopic phenomenology if they belong to the same universality class. Thus, by choosing the simplest models, it would be possible to provide predictions about the macroscopic behavior of many different systems. Indeed, by using models with simple interactions (purely attractive, repulsive, or mixed potentials) computer simulations have found a variety of different mesophases [41][42][43][44][45][46][47][48], that are also found by using more elaborated models [49][50][51][52][53]. Also, it was assumed that an important ingredient to obtain complex lattices is the introduction of anisotropic interactions between constituent particles [54][55][56][57]. However, experimental evidence showed the formation of complex lattices in spherical colloidal particles (including colloidal nanocrystals) grafted with polymers and/or DNA strands [5,[58][59][60]. These results suggest that simple spherical models could still provide useful information regarding the collective behavior of self-assembling colloids.
In this work, we use a simple coarse-grained representation of colloidal nanocrystals that captures the salient features of these materials: a solid core and a soft shell representing the grafted organic ligands. We have used this model to study the nanocrystals phase behavior for a relatively small soft shell [61]. As shown in our previous work, that instance displays a rich set of morphologies, including hexagonal and square lattices, as well as quasicrystalline order [61]. In this work, we have performed extensive computer simulations to elucidate phase diagrams of three instances of this model, representing different physical systems in which the organic ligand shells are varied in size. These phase diagrams were constructed by analyzing different thermodynamic and structural descriptors, namely: pressures, internal energies, compressibility factors, orientational order parameters, structure factors, radial distribution functions and diffusion coefficients. As reported below, a rich variety of lattices were found that could be realized experimentally.

Model and methods
As mentioned above, we consider a simple coarse-grained model to represent a colloidal nanocrystal. The solid core is represented by an impenetrable sphere of diameter σ, whereas the associated organic ligands' degrees of freedom are integrated out and represented by a soft corona of thickness ∆ (see figure 1(a)). Therefore, the effective interaction between two colloidal nanocrystals is defined by a characteristic distance, λ = σ + 2∆. In the case of polymer-grafted colloids, self-consistent polymer field theory calculations have shown that this simple approach represents quite well the effective interactions between colloidal particles [62]. Thus, the size of the organic shell and the associated chemistry are encoded into the characteristic distance of the soft interaction, and the corresponding interaction energies. Specifically, the mathematical expression of the effective pair interaction potential between colloidal nanocrystals, considered in this work, is defined as [1]: where r is the distance between two particles, ε is the height of the soft repulsive shoulder, and σ is the hard core diameter. Parameter k accounts for the steepness of the repulsive shoulder, and λ = σ + 2∆ is the characteristic length scale of soft repulsion (see figure 1). In this work, we set kσ = 10, and explore the phase behavior of three different instances of the family of models represented by equation (1), namely, λ = 1.7σ, 2.0σ, and 3.0σ. We use σ and ε as length and energy units, respectively. In a previous work, we used this model to explore the phase behavior of nanocrystals with λ = 1.35σ [61], in here we are extending the exploration to larger values of λ to have a full picture of the effect of this parameter on the phase behavior of this family of models.
To mimic the self-organization of these colloidal materials at interfaces, all particles are constrained to move on the XY plane. We have performed computer simulations in the NAT ensemble with a Nosé-Hoover thermostat [64,65], where N is the total number of colloidal nanocrystals, A is the total area of the simulation box, and T is the temperature. All calculations were performed using HOOMD-blue [66,67], and periodic boundary conditions were imposed. The systems are composed by N = 2 15 colloids at reduced densities ρ * ≡ Nσ 2 /A, and reduced temperatures T * ≡ k B T/ϵ. For each different instance of the model, we explored the phase behavior in the parameter space (ρ * , T * ), specifically ρ * ∈[0.1, 1.0] and T * ∈1.2, 1.0, 0.8, 0.6, 0.5, 0.4, 0.3, 0.2, 0.18, 0.16, 0.14, 0.12, 0.1, where the density was varied in steps of δρ * = 0.05, unless otherwise noted. To avoid quenching effects at low temperatures, an annealing procedure was used: systems were prepared at a disordered state at the highest temperature, then equilibrated at that condition, the final configuration is used as initial configuration for the next lower temperature, and so on. By following this protocol, we expect that quenching effects are minimized and the configurations we get are equilibrium states. For each simulation, the system is equilibrated by using 5 × 10 7 time steps, then properties are sampled during the next 8 × 10 7 time steps. Post-processing of data is performed by using 10 sub-blocks of the whole trajectory. The reduced time step was set at δt * = δt ϵ/mσ 2 = 0.001, and the thermostat coupling constant to τ = Q/gk B T = 1.0, where Q is the Nosé mass, g is the thermostat's number of degrees of freedom, k B is the Boltzmann constant, and T is the temperature.
The thermodynamic behavior of the system is characterized by using the reduced internal energy U * ≡U/Nε, reduced pressure P * ≡ Pσ 2 /ϵ, compressibility factor Z ≡ P * /(ρ * T * ), and reduced specific heat where the heat capacity, C v , is defined by: To characterize the symmetry of the equilibrated structures, we use the bond orientational order (BOO) parameter, and the structure factor (SF). The global bond orientational order parameter is computed as: [68] where m is an integer characterizing the symmetry, θ j is the angle formed between a fixed axis and the imaginary bond connecting the jth pair of particles separated by no more than a distance r ⋆ ; N B is the number of neighbor bonds, and the brackets indicate ensemble average. The value of r ⋆ is determined as the first minimum after the global maximum of the radial distribution function, at a given temperature and density.
The structure factor is computed by using the following expression: [69] S where ⃗ k is a wave vector, and⃗ r j is the vector position of the jth particle. We also computed the radial distribution function (RDF) as: [70] where A is the area of the simulation box, and r ij is the distance between the ith and jth colloids. Finally, the dynamical properties are characterized by computing the mean squared displacement (MSD): the associated self-diffusion coefficient D was also computed. To obtain D, a linear fit was performed over the MSD vs t curves at long times. Then, we used the following relationship where d is the dimensionality of the system.

Results and discussion
As we mentioned in the Introduction, we are interested in exploring the phase behavior of colloidal nanocrystals as different structural and thermodynamic parameters are varied. In this work, we focus on the following parameters: temperature, density and thickness of the soft corona. We selected three different corona thicknesses: ∆ = 0.35σ, 0.5σ, and 1.0σ, where σ is the diameter of the hard core. For each of these instances, we explore the previously defined parameter space (ρ * , T * ). At every single point of this parameter space, the internal energy U * , pressure P * , compressibility factor Z * , and specific heat c * v were computed. Besides these thermodynamic variables, we also computed the following structural descriptors: bond orientational order (BOO) parameter, structure factor (SF), and the radial distribution function (RDF). To characterize the dynamical properties, the self-diffusion coefficient, D, was computed. All these quantities were used to categorize the different colloidal arrays obtained in our simulations, and to build the phase diagrams reported below.

Small corona, ∆ = 0.35σ
Based on our simulation results, the phase behavior of this instance can be divided into low and high temperature regimes. In the high-temperature regime, T * ∈[0.2, 1.2], the behavior of U * and P * is relatively simple for densities in the range 0.1 ≤ ρ * ≤ 0.7 (see figures S1a-c in Supplementary Information (available online at stacks.iop.org/JPMATER/4/015006/mmedia)). The simple monotonic behavior displayed by these thermodynamic variables, the absence of any global order as deduced from the BOO and SF descriptors (data not shown), and the finite value of the diffusion coefficient indicate that a single isotropic fluid phase dominates the whole thermodynamic region. As the density increases to values larger than 0.7, then all U * and P * isotherms display signatures of a phase transition taking place at different temperatures as the density is varied. According to the BOO parameter, the dominant symmetry after the transition is m = 6 (see figure 2(a)), thus a solid with a 6-fold symmetry is formed. This symmetry-breaking transition from an isotropic phase to a phase with a 6-fold symmetry is analogous to that observed in a hard disks (HD) system. The observed temperature-dependence can be rationalized by considering that the model used in this work has a weak repulsion before reaching the algebraic functional dependence of the hard repulsion (see figure 1(b)). At high temperatures, the colloids have enough energy to overcome that weak repulsion and feel only the hard core, thus, their behavior approaches that one displayed by HD systems. As the temperature is lowered, the weak repulsion becomes more pronounced, therefore, leading to a larger effective diameter, thus decreasing the density where the transition to the solid phase occurs.
In the low-temperature regime, a richer scenario is found: U * (ρ * ) and P * (ρ * ) isotherms indicate the existence of additional phase transitions at intermediate densities (see figures S1d-f in SI). Indeed, the specific heat c * v displays a noticeable increase around those points where the isotherms change in a non-monotonic way (with kinks or apparent discontinuities, see figure S2a in SI). According to the BOO and SF, for densities around ρ * = 0.3, the system exhibits a phase with a 6-fold symmetry (see figure 2(b)), the stability region of this phase is surrounded by a disordered phase. Representative particle configurations of these two phases are shown in figure 3(a) and 3(b). Notice the peculiar colloidal organization in the disordered phase: both the structure factor (inset of figure 3(b)), and the global BOO parameter (see figure 2(b)) indicate the absence of global order in these states. However, as can be seen in the particle configurations, there is a strong local organization where colloids organize into short strings along which particles are close to each other by a distance smaller than the soft shoulder range. These strings, in turn, are separated by a second larger scale, with no orientational correlations. Thus, they behave like a system of fully-flexible polymers. Diffusion coefficient calculations suggest that both ordered and disordered phases are solids. At higher densities (ρ *0 .56), the colloids undergo a phase transformation from the disordered phase to a honeycomb lattice ( figure 3(c)). By increasing the density even further, it is found that the added particles start to fill the empty regions enclosed by the honeycomb framework (see figure 3(d)), and eventually, a high-density hexagonal phase is formed (figures 3(e) and (f)). We used all information provided by the isotherms, BOO, SF and particle configurations to build the phase diagram displayed at figure 4(a). Visual inspection was used to help classify the different geometric structures found in our calculations. In figure 4(b) we present a color map of the self-diffusion coefficient computed for the whole parameter space explored in this work. As can be seen, the diffusion coefficient vanishes both at the ordered phases, and in some isotropic regions at low temperatures and intermediate densities. The two hexagonal phases arise as a consequence of the two length scales associated to the colloids: for intermediate densities and low temperatures, each colloidal nanocrystal behaves as a colloid with an effective particle size equal to the range of the soft interaction, however, at high densities each colloidal particle just sees each other hard core. Therefore, at the respective regimes, they will tend to form the triangular lattice with the appropriate lattice constant. In the supplementary information we have added two representative plots of the RDF (see figures S2b and S2c) that illustrate the change in the correlation length before and after the formation of the low-density hexagonal phase.
It should be noted here that in the phase diagram presented in figure 4(a) we have included several colloidal arrays as a single 'honeycomb' phase. As mentioned above, once the honeycomb lattice is formed, if the density is increased a little, then the global honeycomb structure remains, and the extra colloids that do not find a lattice site to occupy start to fill the empty space enclosed by the honeycomb framework (see sequence of configurations shown in figure 5). Notice that, as the density is increased, the SF evolves towards that expected for the high-density hexagonal phase. Note, however, that even though the SF looks very similar to that of a high-density hexagonal phase, there are still many vacancies in the structure (see figure 5(d)). To define the phase 'boundary' between the honeycomb phases and high-density hexagonal phases (dashed line in figure 4(a)), we computed the coordination number, Υ, defined as where ρ * and g(r * ) are the density and RDF defined in the Methods section, whereas r * min is defined as the distance to the minimum after the first peak of the RDF. Thus, Υ is an estimate of the average number of neighbors contained inside the first shell of radius r * min . Note, however, that the location of the first peak of the RDF (and therefore r * min ) depends on the density and temperature. Thus, to obtain the coordination number as a function of the density, Υ(ρ * ), we need to extract r * min for every RDF at a given temperature and density, we illustrate this on Figure 6(a). In this figure we display the radial distribution functions, g(r * ), for three different densities, namely ρ * = 0.3, 0.4 and 0.6, at T * = 0.1. As can be seen in figure 6(a), for densities ρ * < 0.4, the first maximum of g(r * ) is located at distances larger or equal to the range of the full interaction potential (equation (1)). However, when ρ * = 0.4, the first maximum appears at shorter distances, but its height is smaller that the second peak. At ρ * = 0.6 the location of the first peak does not change much, however its relative height increases whereas the second peak becomes smaller. Figure 6(b) shows the coordination number as a function of the density, Υ(ρ * ), for two different temperatures. As can be seen, at the highest temperature, T * = 1.2, the coordination number increases steadily until it reaches a value of Υ≈ 6 at ρ * ≈ 0.9, and remains in that value as the density increases even further. This density zone corresponds to the high-density hexagonal phase observed at high temperatures (see phase diagram in figure 4(a). On the other hand, at the lowest temperature, the behavior of Υ(ρ * ) changes drastically ( figure 6(b)), in this case, a non-monotonic behavior of Υ(ρ * ) is obtained. The first feature to notice is that Υ(ρ * ) increases rapidly at small densities, and reaches Υ≈ 6 at ρ * ≈ 0.3, which corresponds to the stability region of the low-density hexagonal phase. For ρ * > 0.35 there is a sharp drop in the coordination number, reaching values of Υ≈ 1. This sudden change occurs because the first maximum of g(r * ) is now located at shorter distances and its relative height is smaller than the second peak (see figure 6(a)). Further increase in the density leads to a continuous growth of the coordination number, and it reaches values close to Υ≈ 3 at approximately ρ * ≈ 0.55, the density at which the honeycomb phase appears for the first time. The fact that Υ < 3 for ρ * < 0.6 is due to the existence of vacancies in the honeycomb lattice. On the other hand, Υ > 3 for ρ * > 0.6, and this occurs because the extra colloids can not be placed into the honeycomb lattice, and start filling the empty spaces, therefore some colloids will have more than 3 neighbors, in average. Note that by increasing the density even more, a kink appears at ρ * ≈ 0.65 and then Υ starts to grow at a faster rate, reaching and saturating at Υ≈ 6 around a density of ρ * ≈ 0.85, which we designate as the transition to the high-density hexagonal lattice. We use this criterion to draw the dashed phase boundary displayed in figure 4(a).

Medium size corona, ∆ = 0.5σ
By increasing the thickness of the corona to a larger value, ∆ = 0.5σ, a richer phase diagram is found (see figure 7(a)). As before, we analyze the phase behavior in terms of two temperature regimes. The high temperature regime is the region defined by T * ∈[0.4, 1.2], whereas the low temperature regime includes T * ∈[0.1, 0.4). Isotherms of U * (ρ * ) and P * (ρ * ) in both temperature regimes are shown in figures S3a-f of the supplementary information. For temperatures above 0.7, the scenario is similar to that of the small corona system at high temperatures: at low and intermediate densities the colloidal nanocrystals are in a fluid thermodynamic state, however there is a critical density value where the fluid lost stability, and a new high-density hexagonal solid phase appears. As before, the transition density depends on temperature, and it is associated to the softness of the potential close to the hard core algebraic functional dependence. Close to the lower limit of the high temperature regime, T * ≈ 0.42, the last scenario changes by the occurrence of honeycomb lattices at intermediate densities, which appear for the first time, when ρ * ≈ 0.7 at this temperature. By increasing the density, these honeycomb lattices transform into high-density hexagonal  phases following the same mechanism as that displayed by the small corona nanocrystals. As in the previous instance, to define the transition line between these two phases (dashed line in figure 7(a)), we computed the coordination number as defined above. Two representative Υ(ρ * ) at the lowest and highest temperatures are displayed in figure 8(c). Once again, we classify the states as high-density hexagonal phases if Υ≈ 6, and the corresponding BOO and SF also indicate that m = 6 is the leading symmetry. Interestingly, in the interval of temperatures T * ∈[0.45, 0.65] and densities ρ * ∈[0.8, 0.95], there is a region (light blue on figure 7(a)), where the thermodynamic states are characterized by BOO parameters and structure factors that indicate that these states correspond to globally-isotropic phases. However, at local level, there are strong orientational correlations, with m = 6 as the dominant symmetry. Importantly, the associated diffusion coefficients, although small, are different from zero, we call these states highly-ordered liquids.
As can be deduced by observing the behavior of the thermodynamic descriptors U * , P * , and Z * (see figures S3a-f in SI), as the temperature decreases below 0.4, a series of phase transformations take place. The fact that strong structural rearrangements are taking place in the low temperature regime can be already deduced by the behavior of the specific heat, and RDFs (see figures S4a-b in SI). At T * = 0.1, the system transforms from an isotropic state to a low-density hexagonal phase at ρ * = 0.225 with the characteristic length scale dictated by the soft repulsion range. By increasing the density a bit more, once again disordered states are obtained, however, at ρ * ≈ 0.35 the symmetry is broken again, and colloids collectively arrange themselves into ordered linear arrays, for which m = 2 is the leading global symmetry (see figure 8(a)). This  stripe phase remains stable up to ρ * ≈ 0.425, where the colloidal system rearrange into honeycomb phases, which change as the density is increased even further, as described for the small corona systems, until high-density hexagonal phases are formed. Interestingly, at ρ * ≈ 0.9 a new structural transformation takes place and the nanocrystals form a rhomboid solid (see a representative configuration in figure 8(b), which remains stable up to ρ * ≈ 0.975. At higher densities, we observed the coexistence between the rhomboid solid and the high-density hexagonal phase. The complete phase diagram is presented in figure 7(a), along with the associated diffusion coefficient map, figure 7(b).

Large size corona, ∆ = 1.0σ
The last instance to be considered in this work are colloidal nanocrystals with large coronas, specifically, we set ∆ = 1.0σ. For this case, at temperatures lower than T * = 0.16, the phase behavior becomes too complex that a detailed description at these low temperature states would require δρ * and δT * to be much smaller than the current values being used in this work, thus, increasing the computational cost. Due to this, we confine our analysis to temperatures T * ≥ 0.16. We should mention that a low-density hexagonal phase was also found for this instance, but it appears at temperatures below T * = 0.12 and densities ρ * ≈ 0.1. The phase diagram associated to this system is displayed in figure 9(a).
The first important feature to notice is that the stability region of high-density hexagonal phases is smaller. In fact, by comparing the three phase diagrams elucidated in this work (figures 4(a), 7(a) and 9(a)), one can notice that the stability region of such phase becomes smaller as the interaction range between colloidal nanocrystals increases. Another interesting feature of these colloidal nanocrystals is that hierarchical organization is very common. For example, at T * = 0.16 and ρ * = 0.225, nanocrystals self-organize into a complex lattice where two characteristic length scales are present. At the shortest length scale, three colloids come together to form the basic structural motif that makes up the large scale organization. These triads, in turn, self-assemble into a lattice where the center of mass of each triad is placed on the sites of a hexagonal super-lattice (see figure 10(a)). Note also that each triad forms a triangle, and that strong orientational correlations are apparent in the orientation of these triangles. The corresponding structure factor reflects this complexity by displaying a hierarchical arrangement of scattering peaks ordered as hexagonal patterns embedded in hexagonal patterns at larger scales (see Inset of figure 10(a). At the same low temperature, but higher densities (0.275 ≤ ρ * ≤ 0.425) a new hierarchical structure becomes the most stable configuration, a double-strand stripe phase ( figure 10(b)). As in the previous case, this structure is also characterized by a two-length scale organization of the colloidal nanocrystals: at the shorter scale, the nanocrystals form dumbbells, which in turn self-assemble into stripes (see figure 10(b)). The scattering pattern associated to this particle configuration reflects clearly this two-length scale organization and the associated linear arrangement (Inset of figure 10(b)). By increasing the density to ρ * = 0.45, the system rearranges the stripes into a structure similar to a honeycomb lattice, but now the framework surrounding the empty spaces is formed by more than one particle (compare figures 3(c) and 10(d)). By increasing the density even further, those empty spaces fill up, but instead of self-organizing into a high-density hexagonal phase, colloids arrange into a rhomboid lattice. It was also found that at low temperatures a coexistence region between hexagonal and rhomboid phases is present. Finally, at the largest densities explored in this work, a high-density double stripe phase is obtained (see figure 10(e)). By taking a close inspection to this configuration, one can see that each colloid is in contact with five other colloids in average, a number that is corroborated by our coordination number calculations presented in figure 10(f). Figure 9(b) displays a color map corresponding to the self-diffusion coefficient computed for each state in the thermodynamic parameter space explored in this work.
An important point to notice is the fact that, for the corona sizes studied in here, we did not find quasicrystalline structures on the explored thermodynamic parameter space. In fact, Dijkstra et al [71] and Gruhn et al [72] have shown that both the range and type of the interaction potential have a strong effect on the stability region of quasicrystals (QCs) formation (see, for example, figure 11 in reference [71]). Note, in particular, that soft interaction potentials tend to shift the temperature region where QCs are stable to lower temperatures. In reference [71], the authors have also shown that for the hard, discontinuous model, the interaction range alters the stability of QCs (see figure 9 in reference [71]). For a corona size ∆ = 0.175, they reported the coexistence of square lattice+QC, and pure QC phases as the density increases, which is exactly what our simulation results found in our previous study using the soft model with a similar corona size [61], however our calculations showed that those phases appeared at lower temperatures, in agreement with the observations by Dijkstra et al regarding the effect of the softness of the interaction. Also, in Ref [73], Dijkstra et al used Monte Carlo simulations to study the phase behavior of a hard, discontinuous model with a corona size ∆≈ 0.5, and found that a quasicrystalline approximant, the so-called sigma phase (a periodic lattice), was stable instead of a QC phase. Our simulations also did not uncover a QC phase in the explored thermodynamic parameter space. The main difference is that our simulations predict the existence of a rhomboid phase instead of the sigma phase. This discrepancy can be attributed to: (a) the soft potential used in this work; (b) the small systems considered in reference [73], in there, they used ≈ 256 particles in their Monte Carlo calculations when starting from a crystal structure, and 1600 particles when starting from a disordered state; or (c) the detailed procedure they used to equilibrate the samples. In our calculations we have used substantially bigger systems (2 15 ≈ 32, 000 particles), and performed a careful annealing procedure in which a total of 1.3 × 10 8 time steps were used at each thermodynamic condition. Thus, we are convinced that finite size or finite time effects are not present in our calculations. The existence of large free energy barriers, on the other hand, is a difficult matter to address, and require advanced sampling methods that allow the evaluation of kinetic pathways and metastable states [74,75]. It should also be highlighted that in all other similar works, QCs appear at low temperatures, typically reported to be T * = 0.1 − 0.01 (see below). Our work focuses on higher temperatures, thus it is not unexpected that QCs are not present for the conditions considered in this work.
Before ending this section, we would like to highlight the richness of colloidal nanocrystal arrays that can be generated by changing few structural and thermodynamic parameters. As reported above, colloidal nanocrystals can self-assemble into simple and complex lattices, but also, can display fluid-like states with strong local order. Also, re-entrant transitions are present in the different instances studied in this work. This is a common feature of models with soft interactions and two-length scale potentials [61,[76][77][78][79][80][81][82][83][84]. Note that similar models have also been used to study the formation of quasicrystalline states, and other complex lattices, both with particle-based simulations [68,[71][72][73][85][86][87][88] and field-based modeling [89][90][91][92][93][94][95]. Most of these studies have focused on very low temperatures and high densities where these solid states become stable. The aim of this work was to explore the other region (low and intermediate densities and intermediate temperatures) that has not been studied in detail, thus providing a global picture of the phase behavior of a model that could be a good approximation to colloidal nanocrystals.
We would like to comment about the nanocrystal-nanocrystal interaction potential. The interaction model used in this work is quite simple, that is why we consider our numerical predictions to be a zeroth-order approximation to experimental phase diagrams. A more realistic interaction potential should include a contribution, coming from ligand-to-ligand interactions, that increases monotonically as two nanocrystals approach to each other [96,97]. Ideally, an experimentally-or atomistically-informed interaction potential would be the best option. In this regard, Clancy et al [98] have performed molecular dynamics simulations using united atom force fields, and obtained effective interaction potentials between nanocrystals in the presence of solvents, which are experimentally-relevant conditions. It would be really interesting to map out that effective interaction potential, into the coarse-grained description used in our modeling approach, and explore how much these phase diagrams change. In two recent works [72,99], an interaction potential that increases as two particles approach to each other, was used. Those potentials are still 'toy models' in the sense that the functional forms are guided by intuition rather than from atomistic/experimental data. However, in both cases, the reported simulation results showed the presence of complex particle arrangements, such as lines, honeycomb lattices and quasicrystals. Thus, suggesting that some of these nanocrystal patterns could be realized experimentally. Also, Bonnecaze et al [100] have used an interaction potential that describes the interactions between star polymers, as an effective interaction potential between polymer-decorated nanoparticles. The results in the latter work predict complex arrangements of nanoparticles at various packing fractions, very similar to those reported in here. It should be highlighted that in a recent experimental report by Rey et al [101], spherical colloidal particles were assembled at the air/water interface in the presence of block polymers. It was reported the presence of hexagonal and square lattices as well as stripe phases, all of these phases already obtained in the results described above. A quantitative comparison will still be needed, however, these experimental results already suggest that simple models could be useful to guide experimentalists.
Finally, we would like to conclude by highlighting the richness of possible mesophases that could be obtained by extending the dimensionality of the systems studied in here to three-dimensions. Zero-temperature phase diagrams obtained for the discontinuous square-shoulder model display interesting colloidal arrangements as the size of the corona is increased [102,103]. On the other hand, Pattabhiraman et al [104] performed Monte Carlo simulations at finite temperature for colloids with a corona size of ∆ = 0.5σ in 3D. In the latter, a rich phase behavior displayed by the system was found, including the existence of a pyrochlorelike phase (a photonic crystal). Thus, a study of the model presented in this work in 3D could provide with interesting new phases, phase diagrams at finite temperature and dynamics would be two important aspects to explore.

Conclusions
We have performed extensive computer simulations of simple models of colloidal nanocrystals, and explored the self-assembly behavior as several thermodynamic and structural parameters were varied. Our simulation results suggest that the thickness of the organic soft shell of these versatile materials has a profound impact on the collective organization of these colloids. In particular, we found that hierarchical lattices are common when the size of the corona is comparable to the diameter of the hard core. Our results suggest that complex lattices could be designed by tailoring the structural and thermodynamic parameters. Among the possible configurations we found (low-and high-density) hexagonal lattices, honeycomb, rhomboid and stripe phases. In future work, we aim to study more realistic models of colloidal nanocrystals and explore binary mixtures of these interesting materials.
As mentioned above, most of previous studies have focused on low-temperature, high-density states, because it is in there where some complex patterns have been found: quasicrystals, sigma phases etc. To achieve such a low-temperature states, an experimentalist would have to change temperature to really low values (compared to thermal energy), or to increase the effective repulsion between nanocrystals to high values, which can be achieved by increasing ligand grafted density to really high numbers, which can be difficult to do from a synthesis point of view. Therefore, it could be that those previously predicted states can not be realized experimentally, and high-temperature states (as those described in this work) are the experimentally-relevant configurations. Thus, our work completes the big picture of the organization of these materials by providing this missing information. Our work also provides important insights regarding the behavior at the-so far neglected-low-density regime, where we found hierarchical self-assembly, which is of much interest in different areas of materials science.