Novel Ground-State Crystals with Controlled Vacancy Concentrations: From Kagomé to Honeycomb to Stripes

Novel Ground-State Crystals with Controlled Vacancy Concentrations: From Kagomé to Honeycomb to Stripes


We introduce a one-parameter family, , of pair potential functions with a single relative energy minimum that stabilize a range of vacancy-riddled crystals as ground states. The “quintic potential” is a short-ranged, nonnegative pair potential with a single local minimum of height at unit distance and vanishes cubically at a distance of . We have developed this potential to produce ground states with the symmetry of the triangular lattice while favoring the presence of vacancies. After an exhaustive search using various optimization and simulation methods, we believe that we have determined the ground states for all pressures, densities, and . For specific areas below , the ground states of the “quintic potential” include high-density and low-density triangular lattices, kagomé and honeycomb crystals, and stripes. We find that these ground states are mechanically stable but are difficult to self-assemble in computer simulations without defects. For specific areas above , these systems have a ground-state phase diagram that corresponds to hard disks with radius . For the special case of , a broad range of ground states is available. Analysis of this case suggests that among many ground states, a high-density triangular lattice, low-density triangular lattice, and striped phases have the highest entropy for certain densities. The simplicity of this potential makes it an attractive candidate for experimental realization with application to the development of novel colloidal crystals or photonic materials.

I Introduction

The field of self-assembly provides numerous examples of using atoms, molecules, colloids, and polymers as building blocks in the development of self-organizing functional materials. For example, researchers have used nanoparticles to create stacked rings and spirals using DNA linkers,Sharma et al. (2009) to construct photonic crystals using binary systems of colloids,Hynninen et al. (2007) and to prototype a self-organzing colloidal battery.Cho et al. (2007)

Recently, there have been significant efforts to design pair potentials that yield self-assembly of targeted ground-state (potential energy minimizing) structures using “inverse methods.”Torquato (2009); Cohn and Kumar (2009) With an inverse approach, one chooses a targeted structure or property and uses optimization methods to develop a robust pair potential function. We envision that colloids will be the experimental testbed for the optimized interactions because colloids offer a significant range of repulsive and attractive interactions that can be tailored in the laboratory.Russel et al. (1989) For these models, the many-body interactions of a system are reduced to pair interactions. For a pair potential , where is the distance between two particles, the potential energy per particle of a many-body configuration of particles reduces to .

Inverse approaches have been applied to produce low-coordinated ground states such as honeycomb,Rechtsman et al. (2005) square,Rechtsman et al. (2006) simple cubic,Rechtsman et al. (2006) diamond and wurtzite crystalsRechtsman et al. (2007) as well as many-body configurations on a sphere.Cohn and Kumar (2009) In addition to targeting material structures, the inverse approach has been applied to material properties including negative Poisson’s ratio,Rechtsman et al. (2008) negative thermal expansion,Rechtsman et al. (2007) and scattering characteristics of many-particle systems.Batten et al. (2008) The ability to control vacancy concentrations and arrangements in ground-state structures is one property yet to be achieved via inverse methods. We expect the design of vacancy-riddled crystals to have important technological applications and may provide fundamental insight into certain physical phenomena. The scattering of light,Joannopoulos et al. (2008) ionic conductivity,Agrawal and Gupta (1999) transport processes in heterogeneous materials,Torquato (2002) and possibly supersolid behavior in quantum systemsAndreev and Lifshitz (1969); Kim and Chan (2004) are directly related to vacancies in a material structure.

In this paper, we introduce the “quintic potential” as a pair interaction function that yields vacancy-riddled lattices and low-coordinated crystals as ground states at high density and striped patterns as ground states at low-density. Although not a result of a true inverse method, this family of potential functions arose through a selection process so that vacancy-riddled and low-coordinated lattices are energetically degenerate or possibly favorable to the triangular lattice for certain densities. Each quintic potential function curve, illustrated in Fig. 1, is steeply repulsive for small , has a local minimum at so that , is nonnegative, and is smoothly truncated to be zero at and beyond. The algebraic form of the potential is given in Sec. II. The understanding of the phase behavior associated with the quintic potential is the first step toward developing a more comprehensive inverse approach to controlling vacancies in ground-state structures.

Figure 1: (Color online) Quintic potential, Eq. 6 for various values. The scales of energy and length are dimensionless.

We find that ground states can vary from a high-density triangular lattice, the kagomé and honeycomb crystals, a striped phase, a low-density triangular lattice, and a hard-disk-like fluid as the pressure decreases. The positivity of the local relative minimum is a key feature that allows vacancies to be stabilized in the ground state. Consider a system of particles arranged in a triangular lattice at a specific area (area per particle) so that each of the particle’s first neighbors lie at and second neighbors at . Because the local minimum in is at positive energy, the removal of a particle will lower the total potential energy of the system. However, to maintain the constraint on , the lattice spacing must decrease, which increases the total potential energy. This family of potentials is designed so that the creation of a vacancy and the subsequent reduction in lattice spacing yields a net decrease in the potential energy per particle.

However, the general control of vacancies is a challenging task due to long-ranged vacancy-vacancy interactions. For the Lennard-Jones and many short-ranged repulsive potentials, vacancies arise in equilibrium solids at positive temperature due to their entropic contribution to the free energy.Ashcroft and Mermin (1976) At , they can be “frozen in” during cooling, though they are not present in the equilibrium ground states of these systems. For typical potentials whose ground states are ordered, some pair potentials, such as the Lennard-Jones interaction, vacancies cost energy due to missing “bonds” of negative energy, while with repulsive potentials and the electron crystal, vacancies cost energy due to strain energy in the crystal.

If a vacancy is introduced into an otherwise perfect crystal at positive pressure, strains will arise so that the system can reduce its potential energy and achieve mechanical stability. In typical materials such as Lennard-Jones systems, these relaxations cannot reduce the energy below that of the ground state. When more than one vacancy is present, the strain fields interact and subsequently mediate long-ranged vacancy-vacancy interactions. These interactions are highly nontrivial even in the case of two vacancies in a crystal. In two dimensions, several studies have examined these effective interaction energies between two vacancies as a function of separation distance for a number of potentials including the electron crystal (Coulombic potential with positive background charge),Fisher et al. (1979); Cockayne and Elser (1991); Cândido et al. (2001) Lennard-Jones,Modesto et al. (2008) Gaussian core,Lechner and Dellago (2009) and a repulsive potential.Lechner and Dellago (2009)

In the two-dimensional electron crystal, vacancy-vacancy interactions are attractive and fall off as ,Cockayne and Elser (1991) where is the separation distance between two vacancies, which agrees with continuum elasticity theory.Fisher et al. (1979) For the Lennard-Jones potential, vacancy-vacancy interactions are also attractive at least for short distances.Modesto et al. (2008) These effective vacancy interaction potentials are not necessarily monotonic as a function of lattice spacing.Cândido et al. (2001); Modesto et al. (2008) The displacement fields can be highly anisotropic near a vacancy.Lechner et al. (2008); Lechner and Dellago (2009) In the vicinity of the vacancy, the atomistic details account for the discrepancy between numerical results for the displacement fields and those predicted by continuum elasticity theory. Beyond fifteen lattice spacings, this continuum theory accurately reproduced the results of numerical studies under some boundary conditions.Lechner and Dellago (2009)

The arrangement of low concentrations of vacancies may not simply be controllable because of the nontrivial coupling between the pair potential function, the local deformations that it produces when a vacancy is introduced into a perfect crystal, and the effective vacancy-vacancy interactions. However, a number of pair potential functions have been found to give rise to ground states with high concentrations of vacancies. For example, the “honeycomb potential” was developed via inverse methods to yield the honeycomb crystal as a robust ground-state structure.Rechtsman et al. (2005, 2006) At positive temperatures, the honeycomb structure remained stable for a large temperature and pressure region of the phase diagram.Hynninen et al. (2006) This double-well potential was generalized as the Lennard-Jones-Gaussian potential, which showed a number of complex crystal structures including pentagonal, hexagonal, nonagonal, and decagonal unit cells as the depths and locations of the energy wells were varied.Engel and Trebin (2007, 2008) Colloidal particles with three “patches,” or attractive interaction sites on the surface of a particle, have recently been shown to yield a honeycomb structure for some pressures at in addition to other low-coordinated crystal structures.Doppelbauer et al. (2010) In addition, the square-shoulder potential, a hard-core potential with a soft corona of variable length, shares some of the same ground-state features as the quintic potential. In particular, low-coordinated structures such as a honeycomb-like trivalent configuration and lanes are ground-state featuresFornleitner and Kahl (2008, 2010) that are shared with the quintic potential. Such a potential can also form lamellar and micellar phases for relatively long-ranged coronas as was shown theoretically.Glaser et al. (2007)

Lastly, we note that several potentials qualitatively similar to the quintic potential have been studied in other contexts. For example, a potential that gives rise to quasiperiodic order and the square lattice in two dimensions has one local minimum and local maximum and finite cutoff,Quandt and Teter (1999) as does the Dzugutov potentialDzugutov (1992) which gives rise to a number of complex local clusters in three dimensions.Doye et al. (2001) However, both of these potentials have a local attractive minimum (i.e. negative potential energy). A piecewise linear potential with a hard core, single local minimum, and single local maximum showed the formation of chains and labyrinths at positive temperature, though the ground state is not fully characterized.Haw (2010) However, in contrast to the above potentials, the quintic potential is short-ranged, positive, isotropic, continuous and differentiable, which may be beneficial to experimentalists hoping to realize this potential in a laboratory setting. We find that with this potential, a number of low-coordinated structures arise as ground states including the honeycomb and kagomé crystals as well as a “striped” phase. This represents a significant achievement and a first step toward controlling vacancies in ground-state structures.

The remainder of this paper is as follows. The quintic potential is defined and detailed in Sec. II while our methodology is detailed in Sec. III. In Sec. IV, we define the phase diagram in the specific area- plane and in the pressure- plane. We discuss the characteristics of the phases, the results of simulation, and their mechanical stability. In Sec. V, we focus on the special case where . This case yields anomalous behavior because the potential energy vanishes inside and outside the energy barrier. Lastly, we discuss the implications of this potential and identify extensions to this work in Sec. VI.

Ii Quintic Potential

The generalized quintic potential consists of a fifth-order polynomial that is truncated beyond to be zero. It has the form


and zero otherwise, where is an unscaled height of the local minimum. The coefficients and are chosen so that is a local minimum and , and, as a function of , are given respectively as


The condition that ensures that the stationary point is a relative minimum, and therefore the parameter must be constrainted to


The position of the energy barrier (the relative maximum) occurs at


The generalized pair potential is rescaled so that to set a uniform energy scale,


Upon rescaling, the height of the relative minimum is defined so that and varies from zero to unity. When , the potential is monotonically decreasing and is flat at . In this case, there is no local minimum or maximum. These soft potentials are repulsive for small , and the first and second derivatives vanish at . We construct the potential such that the first and second derivatives vanish at so that the second-order expansions of the potential energy required for stability calculations (e.g. phonons) are well-defined. Examples from the family of quintic potentials are shown in Fig. 1. Because developing a simple expression relating and is difficult, a table relating and was used. The relation between and is shown in Fig. 2.

Figure 2: The relationship between and for the quintic potential.

Iii Methods

At , the free energy per particle is identical to the enthalpy per particle and is related to the average potential energy per particle via . To map the ground-state phase diagram as a function of specific area , pressure , and , we utilized the double-tangent construction method to find the lowest free-energy structure among candidate ground-state configurations. The double-tangent construction is demonstrated in Fig. 3.

We considered several crystal structures as candidate ground states including the triangular lattice (TRI), square lattice (SQ), honeycomb crystal (HC), kagomé crystal (KAG), and its close relative, the “anti-kagomé” crystal (AKG). Whereas the kagomé crystal has vacancies located on a triangular lattice separated by two nearest neighbor spacings, the anti-kagomé crystal has vacancies located on a rectangular lattice. For all densities, the anti-kagomé has a potential energy greater than or equal to the kagomé lattice due to the differing local coordination numbers. We considered these crystal structures because they are relative simple crystal structures with varying degrees of local coordination. The construction of the potential, with a well at and vanishing energy at would intuitively favor such geometries. We use other techniques to systematically explore other crystals with a small number of particles in the unit cell. With these crystal structures, we performed lattice sums over the relevant range of specific area.

We also performed an exhaustive search for energy-minimizing -particle crystals, periodic configurations containing particles per unit cell, while allowing for shape deformation of the unit cell. For example, consider a two-particle crystal with lattice vectors and and a basis where particle is at and particle is located at . For this system, the potential energy per particle is given as


where the summation is over all lattice sites except the origin site and represents particle in a cell other than the origin cell (and likewise for particle ). Minimizing for a fixed area per particle and eliminating casts this as an unconstrained minimization that is function of , , and . The function is minimized using the conjugate gradient algorithm. It is straightforward to generalize this to a larger number of particles in the unit cell.

To ensure that we obtain globally optimal solutions, we use a large number of initial conditions slowly and systematically varying , , and from zero to . After each minimization, we ensure that the unit cell is not significantly sheared. A highly sheared unit cell requires summation over a large number of unit cells to ensure that all interactions are accounted for and therefore we discard those solutions with angles less than 10. This procedure is repeated for nearly two thousand values of on the range so that a smooth versus curve is generated.

These minimizations are a function of variables, and therefore many minimizations at a fixed area per particle can be performed easily. However, as the number of particles per unit cell grows, the number of initial conditions required to ensure global optimality grows as . We performed these minimizations systematically for up to four particles per unit cell, beyond which it becomes too intense to systematically explore thousands of initial conditions for thousands of specific areas.

We also used slow cooling via molecular dynamics (MD) to obtain candidate ground-state structures for larger and possibly disordered systems. With a system of 800 particles in a periodic box, the system was simulated using the Verlet algorithm and Andersen thermostat.Frenkel and Smit (2001) The system was initialized as a high-temperature liquid and slowly cooled according to a linear temperature schedule from to 0.025 over the course of 1.6 time steps. After the simulation was terminated, the system was quenched to a potential energy minimum using the conjugate gradient algorithm.

Lastly, we performed lattice Monte Carlo (MC) optimizations using several hundred lattice sites. In these optimizations, the specific area was fixed and the lattice sites were either occupied or unoccupied by particles. The optimization variables were the occupation state of the lattice sites, the angle between the lattice vectors and the length of one lattice vector. The length of the other lattice vector was constrained by the area, number of particles, angle between the lattice vectors, and the length of the first lattice vector. Possible Monte Carlo moves included inserting a particle to a vacant site, removing a particle from an occupied site, swapping a particle from an occupied site to an unoccupied site, perturbing the angle of the lattice, and perturbing the length of one of the lattice vectors. The system was then optimized via simulated annealing. It was given an initial “temperature,” or energy scale, and moves were accepted or reject based upon the Metropolis algorithm. Ten thousand MC trials were attempted in each cycle and several hundred cycles were performed for each optimization.

It is important to note that although many interesting structural characteristics arose from the results of molecular dynamics simulation and Monte Carlo optimization, no final configurations were identified as ground-state structures. The lattice sums and crystal optimization methods always provided the lowest free-energy structures.

The double-tangent construction requires care and precision since other phases can come very close to the coexistence free energy. The - and - curves are discretized over several thousand data points and linearized between between data points. The identification of coexisting phases can be done by using the lower, concave envelope of the - curve. The coexistence points and ranges of stability were then calculated using a tabulation of , , and and linear interpolation.

Iv Results

iv.1 Phase Diagram

Figure 3: (Color online) Potential energy per particle as a function of specific area for selected crystals and the optimal 2- and 3-particle crystals for . The dense triangular, honeycomb, striped, and open triangular crystals are stable ground states for . The square and anti-kagomé structures are omitted for clarity. Dashed lines represent the double tangent construction.

The double-tangent construction revealed two different types of phase behavior depending on the value of . In no case did the configurations resulting from MD or MC yield a ground state. For , there exist five distinct structures. Figure 3 illustrates the - curves for . In the figure, the optimal -particle crystals are not visible when they are identical to the triangular, honeycomb, or kagomé crystals. As shown in in Fig. 3, at highest pressure, a dense triangular lattice (TRI), where neighboring particles lie inside the potential energy barrier of other neighboring particles, is the ground state. At a reduction of pressure to , the dense triangular lattice is in coexistence with the honeycomb crystal (HC), which is illustrated with the dashed tangent lines. In Fig. 3, the curve for the kagomé lattice appears possibly to touch the coexistence line due to the finite thickness of the line. However, plotting against reveals an intersection between the curves for the TRI and HC structures. The kagomé lattice has a higher free energy than both structures in this density range.

Further reduction of the pressure to 0.75263 yields a “striped” (ST) or lane phase coexisting with the honeycomb phase. This phase is an affinely-stretched triangular lattice. The first few coordination shells contain two, four, and two particles. The curves representing the lowest-energy crystals with a two- and three-particle basis either the triangular lattice, honeycomb crystal, or kagomé crystals generally for with a corresponding basis. For larger , these curves appear to come close to the coexistence line between the HC and ST phases or the ST and open TRI phases. The - curves reveal that the free energy is always above that of the coexisting phases. These near-ground states are generally crystals that nearly mimic the coexistence. For example, the optimal three-particle crystal between the HC and ST phases is an alternating sequence of lines with honeycomb-like coordination and striped coordination. These are suboptimal structures because the specific neighbor lengths do not match exactly. At sufficiently low pressure, for and , the ST phase coexists with an “open” triangular lattice where the neighbor particles lie on the outside of the potential energy barrier.

Lastly, at vanishing pressure, any configuration in which particles are separated by at least is a ground state since the free-energy vanishes. This occurs for . In this region, the ground states behave like hard disks (HD) of radius . There would be a hard-disk crystal and liquid regime with specific area scaled to the hard-disk equation of state, which has been well characterized.Alder and Wainwright (1962)

For , the kagomé crystal emerges as a ground state structure for a narrow range of stability. This is an especially important result. In previous work, researchers used an inverse methodology to design pair interactions for targeted ground states.Rechtsman et al. (2006) Although they were successful in engineering potential for the honeycomb and square crystals, they were unable to do so for the kagomé. The area and pressure ranges of stability increase with . The emergence of the kagomé crystal as a ground state is shown in the - curve for in Fig. 4. As the pressure is reduced, the sequence of stable phases is the dense TRI, KAG, HC, ST, open TRI, and HD phases.

The full phase diagrams in the - plane and the - plane are shown in Figs. 5 and 6, respectively. Figure 5 shows that the coexistence lines are nearly linear for small . The pressure range of stability for the open TRI, ST, and KAG phases widens as increases. When is sufficiently large to include the KAG phase, the pressure range of stability for the HC phase narrows. As increases, the area range of stability increases for most phases, excluding the dense TRI phase, which is evident in Fig. 6.

Near , the slopes of the curves change dramatically. This is due to the softening of the pair potential function as approaches unity and the fact that the relative minimum and maximum come together much more rapidly as approaches unity. Fig. 1 shows that the potential is significantly softer near than other potentials.

The softening of the potential and the rise of the relative minimum are important features for the inclusion of the KAG phase. In comparing Fig. 3 to Fig. 4, one can see that the curve for the kagomé crystal initially begins above the double tangent connecting the TRI and HC phases. As the potential softens and the relative minimum increases, the first minimum in each curve increases proportionally to . However, the softening of the potential bends the left-side of the curves in such a way that the KAG curve lies at lower free energy than that TRI-HC coexistence line. The combination of the height of the relative minimum and the softening are directly related to the potential energy and pressure of the system, the two components of the free energy.

Figure 4: (Color online) Potential energy per particle as a function of specific area for selected crystals and the optimal 2- and 3-particle crystals for . The dense triangular, kagomé, honeycomb, striped, and open triangular crystals are stable ground states for while for , the kagomé also emerges as a stable ground state. The square and anti-kagomé structures are omitted for clarity. Dashed lines represent the double tangent construction.
Figure 5: (Color online) Phase diagram in the - plane as calculated by double tangent method. The hard-disk-like state (HD) occurs for and and the filled triangle on the lower curve represents the triple point.
Figure 6: Phase diagram in the - plane as calculated by the double tangent method.
Figure 7: (Color online) Comparison of the location and coordination number of the triangular lattice, honeycomb crystal and kagomé crystal for , 1.18745, and 1.06445, respectively, for . At this and , the TRI, HC, an KAG phases form a triple point and are in coexistence as ground-state structures. The delicate balance between potential energy and forces allow for coexistence.

The KAG phase emerges at a triple point at and , where the TRI, HON, and KAG phases form a triple point. We have examined the subtle interplay between the shape of the pair potential function and its derivative at this coexistence point in Fig. 7. The figure shows and for . The plot also shows the location of the nearest and next-nearest neighbors for the crystals. The HC phase has the closest neighbors while the TRI phase has the smallest forces.

The coexistence of these structures requires that a complex system of equations be satisfied. For the coexistence pressure and free energy , the system of equations becomes


where is the nearest-neighbor distance for each structure and is necessarily less than unity for the quintic potential. Evidently, the quintic potential for satisfies these equations. The complexity of the coexistence equations is clear in that each structure weighs certain parts differently. The design of other potentials that stabilize these structures may utilize this system of equations within an appropriate optimization framework.

iv.2 Stability

The positivity of the squared frequency of all phonons ensures mechanical stability of a crystal. We have examined the phonon spectra for the KAG, HC, and ST phases. These phases are confirmed to be mechanically stable within the relevant density and ranges. Figure 8 illustrates the phonon spectra for certain points in the reduced first Brillouin zone for the KAG, HC, and ST structures for . THE HC and KAG structures have one particularly soft acoustic branch, labeled in the figures. Using the ratio / at point a simple metric for the extent of mechanical stability, this ratio stays relatively fixed for the KAG lattice as increases. This marks an achievement since previous attempts to stabilize the kagomé crystal had been unsuccessful.Rechtsman et al. (2006) For the HC structure, this ratio is higher for than for that which is plotted in Fig. 8, indicating that the crystal is more stable as increases. In general, the quintic potential and the honeycomb potentialRechtsman et al. (2006) yield similar ratios. The striped phase whose lattice vectors for and are and is also mechanically stable. The phonon spectra from the other “quadrants” are identical to the quadrant displayed in Fig. c due to the symmetry of the Brillouin zone.

Figure 8: (Color online) Phonon spectra for (a) the kagomé crystal at and (b) the honeycomb crystal at , and (c) the striped phase at for wave vectors in the reduced first Brillouin zone for . The first Brillouin zone for the striped phase is shown in inset in (c). These crystals are mechanically stable at these densities with the quintic potential. The branches containing soft acoustic modes are denoted with a *.

iv.3 Low-Energy Structures

Although those structure obtained by MD and MC methods were suboptimal, Figures 3 and 4 show that the free-energy differences between these structures and the ground states are small. Several structural motifs emerge for this potential and are shown in Fig. 9 for . These metastable states may have technological value if the freezing kinetics are sufficiently slow. Preliminary simulations suggest the freezing behavior varies with density.

For example, with , systems typically freeze into a triangular lattice with vacancies randomly distributed as in Fig. a. Although this is not a ground-state structure, it yields a vacancy-riddled lattice. The vacancies appear to have no particular order. The freezing transition, the temperature at which the system changes from a high-density liquid to an ordered phase, appears to be first order. The drop in potential energy as the temperature is reduced is sharp in this density range.

Systems at densities where coexistence between the HC and ST phases are the ground states tend to exhibit a freezing from the liquid phase to a rigid, structured phase. However, the structured phase, which we believe to be metastable, is characterized by rings and strings as in Fig. b. The six-particle ring is a characteristic of the HC structure and the strings are characteristics of the ST phase. These metastable states are disordered due to the fluid nature of the strings. The drop in potential with temperature is not as sharp for these densities compared to those at higher densities.

Lastly, at lower densities, as with those associated with the ST and open TRI phases, the metastable, low-energy states adopt labyrinthine characteristics, Fig. c, and eventually colloidal polymers and monomers, d. In this density range, the drop in potential energy associated with freezing is weak. The phenomena and characteristics of low-energy states are common to all . However, as increases so does the propensity to form labyrinthine characteristics. This is due to the lower energy barrier and the decrease in , the location of the energy barrier. A qualitatively similar, but piecewise linear, potential also yields a labyrinth phase at positive temperature.Haw (2010) Although the positive temperature behavior has not been fully explored in this paper, we anticipate that the equation of state for this family of potentials will have interesting behavior. Manipulation of the kinetics to achieve such unusual structures represents one path to achieving kinetically stable materials with controlled vacancy concentrations.

Figure 9: (Color online) Metastable configurations obtained via molecular dynamics for . In an attempt to cool to a liquid to the ground state, the system has difficulty crystallizing and/or phase separating into the appropriate ground-state structures. As density is decreased, the metastable states go from ordered structures with vacancies to a labyrinthine network to a monomeric liquid.

V Special Cases:

The pair potential has interactions that vanish at a pair distance of unity and beyond . At positive pressure, the ground state is the dense triangular lattice. However, at zero pressure and , ground-state configurations will have a vanishing potential energy and pressure (enthalpy vanishes). Thus, the type of available ground states is dependent on the area.

A number of interesting structures arise as ground states. For , dilutions of the triangular lattice with unit neighbor spacing are ground states, including the honeycomb and kagomé crystals and other lattice gases. For , the striped phase, a Bravais lattice with lattice vectors of and , becomes available as a ground state. Each particle has two neighbors at unit separation and four neighbors at separation of . For , dilutions of this lattice are also ground states. In addition, any small expansion in the direction normal to the stripes, will allow each stripe to gain some fluidity. Complex liquids of strings have found as ground states using molecular dynamics for at least . However, the geometric problem is highly nontrivial since for several runs with , ground states could not always be obtained. For , an open triangular lattice with neighbor spacing of is a ground state structure.

The question remains as to which type of configurations are entropically (thermodynamically) favorable, or rather which type of system makes up the largest fraction of configuration space. We first make the distinction between discrete (e.g. lattice gas) and continuous entropy (e.g. hard-disk crystal). In the classical case, continuous entropy is uncountable while discrete entropy is not. Making the analogy to hard disks and hard spheres, we believe that the highest entropy phases for a specified area are those that have the availability for continuous entropy. Configurations with the highest-dimensional configuration space dominate the entropy. We estimate that, for a certain , configurations with the fewest constrained degrees of freedom will have the highest dimensional configuration space, and therefore the highest entropy. For , all degrees of freedom per particle are constrained. For , no degrees of freedom per particle are constrained.

In order to gain continuous entropy, the system must have the ability to expand the “bonds,” which can be either first-neighbor or second-neighbor connections between particles. Any slight expansion of the dense triangular lattice is not a ground state because each neighbor must be constrained to unit separation. However, an expansion of the open triangular lattice is still a ground state because the bonds are not constrained from above. The problem can be cast as a tiling problem of four triangular tiles. The possible triangular tiles are those with side lengths of unity or . These are depicted in Table 1 along with their shorthand notations.

Because there are no constraints beyond distances of , the edges with lengths of can be considered “elastic.” For simplicity, we first consider the tiling problem where these edge lengths are fixed. The rules of the tiling problem are as follows:

  • The plane must be tiled with no gaps.

  • For pairs of tiles that share an edge, the edges must be of the same length so that the vertices also match.

  • The (1,1,r) and (1,r,r) tiles cannot share a edge since this violates a second neighbor constraint (i.e. two vertices are separated by a length between unity and ).

Name Shape Angles () Area
(1,1,1) 60, 60, 60
(r,r,r) 60, 60, 60
(1,1,r) 30, 30, 120
(1,r,r) 33.56, 73.22, 73.22
Table 1: Available triangular tiles for the tiling problem.

The internal angle of the (1,r,r) tiles are such that they cannot integrate well with the other tiles. Therefore, in any tiling that includes (1,r,r) tiles, these tiles must “phase separate” from the others. We consider two types of tilings - those with (1,r,r) tiles and those without. First, we consider those without (1,r,r) tiles. An example of such a tiling is shown in Fig. 10 where each vertex is decorated with a particle. These tilings are simply dilutions of the dense triangular lattice, or coexistence between a dense and open triangular lattice. The (1,1,r) tiles have two roles. They can act as intermediaries between domains of the dense triangular lattice and the open triangular lattice, and they can dimerize to form a rhombus making a domain of the dense triangular lattice. This coexistence between the dense and open triangular lattices is represented as a double-tangent line connecting the fully constrained dense triangular lattice to the unconstrained open triangular lattice as shown in Fig. 11.

Next we consider those tilings that include (1,r,r) tiles. Two periodic tilings with a specific area of exist and are shown in Figs. 10 and 10. We deem these the stripe and “zig-zag” tilings respectively. In these phases, any small expansion perpendicular to the stripe or zig-zag will allow the stripes to have fluidity. On average, these configurations will constrain one degree of freedom per particle, as plotted in Fig. 11.

For any tiling consisting of all tiles, the (1,r,r) tiles must segregate so that the the other tiles can fully tile the plane. Dimerizing the (1,r,r) tiles along the edges creates a wide stripe. Building off the wide stripe requires (1,1,1) tiles or (1,1,r) tiles and forms a local coexistence with the dense triangular lattice as shown in Fig. 10. Dimerizing the (1,r,r) tiles along the unit length edge creates a narrow stripe. The exposed edges have length , requiring (r,r,r) tiles to be directly adjacent to the stripe. This forms a coexistence between the the striped phase and the open triangular lattice. The double tangent construction of these coexistences is displayed as the solid red line in Fig. 11. These coexistences between the striped phase and the dense triangular lattice and the striped phase and the open triangular lattice maximize the number of unconstrained degrees of freedom and likely maximize the entropy of the ground state. In addition, these system can take on discrete entropy by way of stacking variants of striped phase and the lattice phase. Using “S” and “T” to denote a stripe and a triangular lattice layer, a STTTTSSS system is degenerate with a STSTSTST system.

The zig-zag phase can only form local coexistence with the open triangular lattice. To do so, the (1,r,r) tiles must form a trimer which exposes the edges on either side of the trimer. Since they are in trimers, their stacking entropy would be less that the stacking entropy available to the coexistence between the stripe and open triangular lattice. Therefore, we believe the stripe phase has a higher entropy than the zig-zag phase.

Figure 10: (Color online) (a) Tiling excluding (1,r,r) tiles, (b) striped phase using only (1,r,r) tiles, (c) zig-zag phase using only (1,r,r) tiles (d) coexistence between dense triangular lattice and the striped phase, and (e) coexistence between the open triangular lattice and the striped phase.
Figure 11: (Color online) Unconstrained degrees of freedom per particle for candidate ground-state structures for potential. Configurations with the maximum number of unconstrained degrees of freedom are presumed to be the ground state (red circles).

Vi Discussion and Concluding Remarks

In this paper, we have developed the ground-state phase diagram of a new potential that gives rise to a number of novel phases that include low-coordinated crystal structures. Given the unusual nature of the ground state, we expect that the equation of state and full phase diagram for these systems to exhibit other unusual behaviors at positive temperature. For example, in the global phase diagram of the the Lennard-Jones-Gauss, or honeycomb potential, there was no gas-liquid coexistence in the low-temperature, low-density part of the phase diagram, nor was there a liquid-liquid phase coexistenceHynninen et al. (2006) We expect most of the solid-solid transitions that we find at zero temperature will remain at small nonzero temperature. Our preliminary calculations suggest that strings, or polymers, may arise in equilibrium at low densities.

In addition, we have interest in the mobility of vacancies, particularly in the cases where vacancy concentrations are dilute, (e.g. and ), due to the possible relation to supersolid behavior.Andreev and Lifshitz (1969); Kim and Chan (2004) For just above and , we expect there to be a small number of vacancies in the system at low-temperature. We have estimated the potential energy required for a particle to “hop” from a lattice site to a vacant site. By initializing a system with one vacancy with a “hopping” particle along the transition to the vacant site, in an otherwise undistorted lattice, we minimized , where represents particle coordinates, using conjugate gradient minimization. The resulting configuration represents a saddle point in the energy landscape. The difference in potential energy of the saddle point and the ground state is the energy barrier required for the particle and vacancy to swap positions. This barrier sets an activation energy for classical thermal motion. For the quantum case, it would also enter in the tunneling rate. For potential, there is a single saddle point in the vicinity of the numerous initial conditions in the energy landscape with a total energy barrier height of 5.569015. The diffusing particle is midway between the origin and destination sites while the bracing particles are displaced off the lattice line to accommodate the jump. Using molecular dynamics, we expect to relate this saddle point energy to a vacancy diffusion coefficient.

The quintic potential can further be generalized by varying the distance at which the function is truncated. We set this cutoff distance to be . However, allowing this cutoff distance to vary introduces a larger class of potentials. A systematic study of the cutoff radius on the robustness of the ground states is necessary for experimental realization. The hard-core plus square shoulder potential has ground states that vary significantly depending on the relative lengthscale of the hard-core distance and the square-shoulder distanceFornleitner and Kahl (2008, 2010) It is expected that the ground states of the generalized quintic potential would be sensitive to the location of the minimum and the cutoff distance, though it is currently unknown how sensitive the ground-state phase diagram is to these parameters. Understanding this sensitivity is important to experimentalists who want a simple, robust potential. Developing an optimization procedure to make self-assembly more robust would be particularly useful. In addition, an extension to three-dimensions may provide additional fundamental insight.

As mentioned earlier, inverse optimization techniques have been effective in developing potentials for targeted material properties. We intend to develop a general and broad inverse optimization technique to target specific vacancy arrangements by accounting for and/or manipulating long-ranged vacancy-vacancy interactions. For example, one might develop an objective function whose variables include the strength, sign (attractive/repulsive), and angular dependence of vacancy-vacancy interactions. Alternatively, one might consider a two-component system system of a heavy particle and a light particle and apply an inverse optimization technique to this system. Using a broad family of potentials, one could then optimize over the available parameters to achieve a dilute concentration of effectively repulsive vacancies.

Vii Acknowledgments

We thank Lawrence Cheuk for initial work on a related model. S.T. thanks the Institute for Advanced Study for its hospitality during his stay there. This work was supported by the Office of Basic Energy Sciences, U.S. Department of Energy, Grant DE-FG02-04-ER46108..


  1. J. Sharma, R. Chhabra, A. Cheng, J. Brownell, Y. Liu and H. Yan, Science, 2009, 323, 112–116.
  2. A. P. Hynninen, J. H. J. Thijssen, E. C. M. Vermolen, M. Dijkstra and A. van Blaaderen, Nature Materials, 2007, 6, 202–205.
  3. Y. K. Cho, R. Wartena, S. M. Tobias and Y. M. Chiang, Adv. Func. Mat., 2007, 17, 379–389.
  4. S. Torquato, Soft Matter, 2009, 5, 1157–1173.
  5. H. Cohn and A. Kumar, Proc. Nat. Acad. Sci., 2009, 106, 9570–9575.
  6. W. B. Russel, D. A. Saville and W. R. Schowalter, Colloidal dispersions, Cambridge University Press: Cambridge, 1989.
  7. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett., 2005, 95, 228301.
  8. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E, 2006, 73, 011406.
  9. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E, 2006, 74, 021404.
  10. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev, E, 2007, 75, 031403.
  11. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett., 2008, 101, 85501.
  12. M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. Phys. Chem. A, 2007, 111, 12816–12821.
  13. R. D. Batten, F. H. Stillinger and S. Torquato, J. App. Phys., 2008, 104, 033504.
  14. J. D. Joannopoulos, S. G. Johnson, J. N. Winn and R. D. Meade, Photonic crystals: molding the flow of light, Princeton University Press: Princeton, N.J., 2008.
  15. R. C. Agrawal and R. K. Gupta, J. Mat. Sci., 1999, 34, 1131–1162.
  16. S. Torquato, Random Heterogeneous Materials: Microstucture and Macroscopic Properties, SpringerVerlag: New York, 2002.
  17. A. F. Andreev and L. M. Lifshitz, Sov. Phys. JETP, 1969, 29, 1107.
  18. E. Kim and M. H. W. Chan, Nature, 2004, 427, 225–227.
  19. N. W. Ashcroft and N. D. Mermin, Solid State Physics, Saunders College: Philadelphia, Pa, 1976.
  20. D. S. Fisher, B. I. Halperin and R. Morf, Phys. Rev. B, 1979, 20, 4692–4712.
  21. E. Cockayne and V. Elser, Phys. Rev. B, 1991, 43, 623–629.
  22. L. Cândido, P. Phillips and D. M. Ceperley, Phys. Rev. Lett., 2001, 86, 492–495.
  23. L. Modesto, D. L. S. Júnior, J. N. Teixeira Rabelo and L. Cândido, Solid State Communications, 2008, 145, 355–358.
  24. W. Lechner and C. Dellago, Soft Matter, 2009, 5, 646–659.
  25. W. Lechner and C. Dellago, Soft Matter, 2009, 5, 2752–2758.
  26. W. Lechner, E. Schöll-Paschinger and C. Dellago, J. Phys. Condens. Matter, 2008, 20, 404202.
  27. A. P. Hynninen, A. Z. Panagiotopoulos, M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. Chem. Phys., 2006, 125, 024505.
  28. M. Engel and H. R. Trebin, Phys. Rev. Lett., 2007, 98, 225505.
  29. M. Engel and H. R. Trebin, Zeitschrift für Kristallographie, 2008, 223, 721–725.
  30. G. Doppelbauer, E. Bianchi and G. Kahl, J. Phys.: Condensed Matter, 2010, 22, 104105.
  31. J. Fornleitner and G. Kahl, Eur. Phys. Lett., 2008, 82, 18001.
  32. J. Fornleitner and G. Kahl, J. Phys: Condes. Matter, 2010, 22, 104118.
  33. M. A. Glaser, G. M. Grason, R. D. Kamien, A. Košmrlj, C. D. Santangelo and P. Ziherl, Eur. Phys. Lett., 2007, 78, 46004:1–5.
  34. A. Quandt and M. P. Teter, Phys. Rev. B, 1999, 59, 8586–8592.
  35. M. Dzugutov, Phys. Rev. A, 1992, 46, 2984–2987.
  36. J. P. K. Doye, D. J. Wales and S. I. Simdyankin, Faraday Discussions, 2001, 118, 159–170.
  37. M. D. Haw, Phys. Rev. E, 2010, 81, 031402.
  38. D. Frenkel and B. Smit, Understanding molecular simulation, Academic Press, Inc. Orlando, FL, USA, 2001.
  39. B. J. Alder and T. E. Wainwright, Phys. Rev., 1962, 127, 359–361.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description