Novel GroundState Crystals with Controlled Vacancy Concentrations: From Kagomé to Honeycomb to Stripes
Abstract
We introduce a oneparameter family, , of pair potential functions with a single relative energy minimum that stabilize a range of vacancyriddled crystals as ground states. The “quintic potential” is a shortranged, 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 highdensity and lowdensity triangular lattices, kagomé and honeycomb crystals, and stripes. We find that these ground states are mechanically stable but are difficult to selfassemble in computer simulations without defects. For specific areas above , these systems have a groundstate 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 highdensity triangular lattice, lowdensity 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 selfassembly provides numerous examples of using atoms, molecules, colloids, and polymers as building blocks in the development of selforganizing 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 selforganzing colloidal battery.Cho et al. (2007)
Recently, there have been significant efforts to design pair potentials that yield selfassembly of targeted groundstate (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 manybody 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 manybody configuration of particles reduces to .
Inverse approaches have been applied to produce lowcoordinated 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 manybody 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 manyparticle systems.Batten et al. (2008) The ability to control vacancy concentrations and arrangements in groundstate structures is one property yet to be achieved via inverse methods. We expect the design of vacancyriddled 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 vacancyriddled lattices and lowcoordinated crystals as ground states at high density and striped patterns as ground states at lowdensity. Although not a result of a true inverse method, this family of potential functions arose through a selection process so that vacancyriddled and lowcoordinated 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 groundstate structures.
We find that ground states can vary from a highdensity triangular lattice, the kagomé and honeycomb crystals, a striped phase, a lowdensity triangular lattice, and a harddisklike 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 longranged vacancyvacancy interactions. For the LennardJones and many shortranged 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 LennardJones 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 LennardJones 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 longranged vacancyvacancy 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) LennardJones,Modesto et al. (2008) Gaussian core,Lechner and Dellago (2009) and a repulsive potential.Lechner and Dellago (2009)
In the twodimensional electron crystal, vacancyvacancy 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 LennardJones potential, vacancyvacancy 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 vacancyvacancy 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 groundstate 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 doublewell potential was generalized as the LennardJonesGaussian 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 lowcoordinated crystal structures.Doppelbauer et al. (2010) In addition, the squareshoulder potential, a hardcore potential with a soft corona of variable length, shares some of the same groundstate features as the quintic potential. In particular, lowcoordinated structures such as a honeycomblike trivalent configuration and lanes are groundstate featuresFornleitner and Kahl (2008, 2010) that are shared with the quintic potential. Such a potential can also form lamellar and micellar phases for relatively longranged 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 shortranged, 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 lowcoordinated 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 groundstate 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 fifthorder polynomial that is truncated beyond to be zero. It has the form
(1) 
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
(2)  
(3) 
The condition that ensures that the stationary point is a relative minimum, and therefore the parameter must be constrainted to
(4) 
The position of the energy barrier (the relative maximum) occurs at
(5) 
The generalized pair potential is rescaled so that to set a uniform energy scale,
(6) 
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 secondorder expansions of the potential energy required for stability calculations (e.g. phonons) are welldefined. 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.
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 groundstate phase diagram as a function of specific area , pressure , and , we utilized the doubletangent construction method to find the lowest freeenergy structure among candidate groundstate configurations. The doubletangent 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 “antikagomé” crystal (AKG). Whereas the kagomé crystal has vacancies located on a triangular lattice separated by two nearest neighbor spacings, the antikagomé crystal has vacancies located on a rectangular lattice. For all densities, the antikagomé 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 energyminimizing particle crystals, periodic configurations containing particles per unit cell, while allowing for shape deformation of the unit cell. For example, consider a twoparticle 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
(7) 
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 groundstate 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 hightemperature 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 groundstate structures. The lattice sums and crystal optimization methods always provided the lowest freeenergy structures.
The doubletangent 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
The doubletangent 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 affinelystretched triangular lattice. The first few coordination shells contain two, four, and two particles. The curves representing the lowestenergy crystals with a two and threeparticle 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 nearground states are generally crystals that nearly mimic the coexistence. For example, the optimal threeparticle crystal between the HC and ST phases is an alternating sequence of lines with honeycomblike 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 freeenergy vanishes. This occurs for . In this region, the ground states behave like hard disks (HD) of radius . There would be a harddisk crystal and liquid regime with specific area scaled to the harddisk 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 leftside of the curves in such a way that the KAG curve lies at lower free energy than that TRIHC 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.
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 nextnearest 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
(8)  
(9)  
(10)  
(11)  
(12)  
(13) 
where is the nearestneighbor 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.
iv.3 LowEnergy Structures
Although those structure obtained by MD and MC methods were suboptimal, Figures 3 and 4 show that the freeenergy 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 groundstate structure, it yields a vacancyriddled lattice. The vacancies appear to have no particular order. The freezing transition, the temperature at which the system changes from a highdensity 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 sixparticle 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, lowenergy 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 lowenergy 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.
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 , groundstate 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. harddisk 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 highestdimensional 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 firstneighbor or secondneighbor 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 
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 doubletangent 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 “zigzag” tilings respectively. In these phases, any small expansion perpendicular to the stripe or zigzag 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 zigzag 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 zigzag phase.
Vi Discussion and Concluding Remarks
In this paper, we have developed the groundstate phase diagram of a new potential that gives rise to a number of novel phases that include lowcoordinated 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 LennardJonesGauss, or honeycomb potential, there was no gasliquid coexistence in the lowtemperature, lowdensity part of the phase diagram, nor was there a liquidliquid phase coexistenceHynninen et al. (2006) We expect most of the solidsolid 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 lowtemperature. 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 hardcore plus square shoulder potential has ground states that vary significantly depending on the relative lengthscale of the hardcore distance and the squareshoulder 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 groundstate 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 selfassembly more robust would be particularly useful. In addition, an extension to threedimensions 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 longranged vacancyvacancy interactions. For example, one might develop an objective function whose variables include the strength, sign (attractive/repulsive), and angular dependence of vacancyvacancy interactions. Alternatively, one might consider a twocomponent 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 DEFG0204ER46108..
References
 J. Sharma, R. Chhabra, A. Cheng, J. Brownell, Y. Liu and H. Yan, Science, 2009, 323, 112–116.
 A. P. Hynninen, J. H. J. Thijssen, E. C. M. Vermolen, M. Dijkstra and A. van Blaaderen, Nature Materials, 2007, 6, 202–205.
 Y. K. Cho, R. Wartena, S. M. Tobias and Y. M. Chiang, Adv. Func. Mat., 2007, 17, 379–389.
 S. Torquato, Soft Matter, 2009, 5, 1157–1173.
 H. Cohn and A. Kumar, Proc. Nat. Acad. Sci., 2009, 106, 9570–9575.
 W. B. Russel, D. A. Saville and W. R. Schowalter, Colloidal dispersions, Cambridge University Press: Cambridge, 1989.
 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett., 2005, 95, 228301.
 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E, 2006, 73, 011406.
 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E, 2006, 74, 021404.
 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev, E, 2007, 75, 031403.
 M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett., 2008, 101, 85501.
 M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. Phys. Chem. A, 2007, 111, 12816–12821.
 R. D. Batten, F. H. Stillinger and S. Torquato, J. App. Phys., 2008, 104, 033504.
 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.
 R. C. Agrawal and R. K. Gupta, J. Mat. Sci., 1999, 34, 1131–1162.
 S. Torquato, Random Heterogeneous Materials: Microstucture and Macroscopic Properties, SpringerVerlag: New York, 2002.
 A. F. Andreev and L. M. Lifshitz, Sov. Phys. JETP, 1969, 29, 1107.
 E. Kim and M. H. W. Chan, Nature, 2004, 427, 225–227.
 N. W. Ashcroft and N. D. Mermin, Solid State Physics, Saunders College: Philadelphia, Pa, 1976.
 D. S. Fisher, B. I. Halperin and R. Morf, Phys. Rev. B, 1979, 20, 4692–4712.
 E. Cockayne and V. Elser, Phys. Rev. B, 1991, 43, 623–629.
 L. Cândido, P. Phillips and D. M. Ceperley, Phys. Rev. Lett., 2001, 86, 492–495.
 L. Modesto, D. L. S. Júnior, J. N. Teixeira Rabelo and L. Cândido, Solid State Communications, 2008, 145, 355–358.
 W. Lechner and C. Dellago, Soft Matter, 2009, 5, 646–659.
 W. Lechner and C. Dellago, Soft Matter, 2009, 5, 2752–2758.
 W. Lechner, E. SchöllPaschinger and C. Dellago, J. Phys. Condens. Matter, 2008, 20, 404202.
 A. P. Hynninen, A. Z. Panagiotopoulos, M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. Chem. Phys., 2006, 125, 024505.
 M. Engel and H. R. Trebin, Phys. Rev. Lett., 2007, 98, 225505.
 M. Engel and H. R. Trebin, Zeitschrift für Kristallographie, 2008, 223, 721–725.
 G. Doppelbauer, E. Bianchi and G. Kahl, J. Phys.: Condensed Matter, 2010, 22, 104105.
 J. Fornleitner and G. Kahl, Eur. Phys. Lett., 2008, 82, 18001.
 J. Fornleitner and G. Kahl, J. Phys: Condes. Matter, 2010, 22, 104118.
 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.
 A. Quandt and M. P. Teter, Phys. Rev. B, 1999, 59, 8586–8592.
 M. Dzugutov, Phys. Rev. A, 1992, 46, 2984–2987.
 J. P. K. Doye, D. J. Wales and S. I. Simdyankin, Faraday Discussions, 2001, 118, 159–170.
 M. D. Haw, Phys. Rev. E, 2010, 81, 031402.
 D. Frenkel and B. Smit, Understanding molecular simulation, Academic Press, Inc. Orlando, FL, USA, 2001.
 B. J. Alder and T. E. Wainwright, Phys. Rev., 1962, 127, 359–361.