Novel Ground-State Crystals with Controlled Vacancy Concentrations: From Kagomé to Honeycomb to Stripes
Robert D. Batten, David A. Huse, Frank H. Stillinger, Salvatore Torquato
aa r X i v : . [ c ond - m a t . m t r l - s c i ] O c t Novel Ground-State Crystals with Controlled VacancyConcentrations: From Kagom´e to Honeycomb to Stripes
Robert D. Batten, David A. Huse, Frank H. Stillinger, and Salvatore Torquato
2, 3, 4, 5, 6, 7, ∗ Department of Chemical Engineering,Princeton University, Princeton, NJ 08544, USA Department of Physics, Princeton University, Princeton, NJ, 08544, USA Department of Chemistry, Princeton University, Princeton, NJ, 08544, USA Princeton Institute for the Science and Technology of Materials,Princeton University, Princeton, NJ, 08540, USA Program in Applied and Computational Mathematics,Princeton University, Princeton, NJ, 08544, USA Princeton Center for Theoretical Science,Princeton University, Princeton, NJ, 08644, USA School of Natural Sciences, Institute for Advanced Study, Princeton, NJ, 08544 USA (Dated: June 20, 2018) bstract We introduce a one-parameter family, 0 ≤ H ≤
1, of pair potential functions with a single relativeenergy minimum that stabilize a range of vacancy-riddled crystals as ground states. The “quinticpotential” is a short-ranged, nonnegative pair potential with a single local minimum of height H at unit distance and vanishes cubically at a distance of √
3. We have developed this potentialto produce ground states with the symmetry of the triangular lattice while favoring the presenceof vacancies. After an exhaustive search using various optimization and simulation methods, webelieve that we have determined the ground states for all pressures, densities, and 0 ≤ H ≤
1. Forspecific areas below 3 √ /
2, the ground states of the “quintic potential” include high-density andlow-density triangular lattices, kagom´e and honeycomb crystals, and stripes. We find that theseground states are mechanically stable but are difficult to self-assemble in computer simulationswithout defects. For specific areas above 3 √ /
2, these systems have a ground-state phase diagramthat corresponds to hard disks with radius √
3. For the special case of H = 0, a broad rangeof ground states is available. Analysis of this case suggests that among many ground states, ahigh-density triangular lattice, low-density triangular lattice, and striped phases have the highestentropy for certain densities. The simplicity of this potential makes it an attractive candidate forexperimental realization with application to the development of novel colloidal crystals or photonicmaterials. PACS numbers: . 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 us-ing DNA linkers, to construct photonic crystals using binary systems of colloids, and toprototype a self-organzing colloidal battery. Recently, there have been significant efforts to design pair potentials that yield self-assembly of targeted ground-state (potential energy minimizing) structures using “inversemethods.”
With an inverse approach, one chooses a targeted structure or property anduses optimization methods to develop a robust pair potential function. We envision thatcolloids will be the experimental testbed for the optimized interactions because colloidsoffer a significant range of repulsive and attractive interactions that can be tailored in thelaboratory. For these models, the many-body interactions of a system are reduced to pairinteractions. For a pair potential v ( r ), where r is the distance between two particles, thepotential energy per particle u of a many-body configuration r N of N particles reduces to u ( r N ) = N P i FIG. 1: (Color online) Quintic potential, Eq. 6 for various H values. The scales of energy andlength are dimensionless. riddled and low-coordinated lattices are energetically degenerate or possibly favorable to thetriangular lattice for certain densities. Each quintic potential function curve, illustrated inFig. 1, is steeply repulsive for small r , has a local minimum at r = 1 so that v (1) = H , isnonnegative, and is smoothly truncated to be zero at r = √ a = √ / r = 1 and second neighbors at r = √ v ( r ) is at positive energy, the removal of a particle will lowerthe total potential energy of the system. However, to maintain the constraint on a , thelattice spacing must decrease, which increases the total potential energy. This family ofpotentials is designed so that the creation of a vacancy and the subsequent reduction inlattice spacing yields a net decrease in the potential energy per particle.4owever, 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 contribu-tion to the free energy. At T = 0, they can be “frozen in” during cooling, though theyare not present in the equilibrium ground states of these systems. For typical potentialswhose ground states are ordered, some pair potentials, such as the Lennard-Jones interac-tion, vacancies cost energy due to missing “bonds” of negative energy, while with repulsivepotentials 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, strainswill arise so that the system can reduce its potential energy and achieve mechanical stabil-ity. In typical materials such as Lennard-Jones systems, these relaxations cannot reduce theenergy below that of the ground state. When more than one vacancy is present, the strainfields interact and subsequently mediate long-ranged vacancy-vacancy interactions. Theseinteractions are highly nontrivial even in the case of two vacancies in a crystal. In twodimensions, several studies have examined these effective interaction energies between twovacancies as a function of separation distance for a number of potentials including the elec-tron crystal (Coulombic potential with positive background charge), Lennard-Jones, Gaussian core, and a repulsive r − potential. In the two-dimensional electron crystal, vacancy-vacancy interactions are attractive andfall off as r − v , where r v is the separation distance between two vacancies, which agreeswith continuum elasticity theory. For the Lennard-Jones potential, vacancy-vacancy inter-actions are also attractive at least for short distances. These effective vacancy interactionpotentials are not necessarily monotonic as a function of lattice spacing. The displace-ment fields can be highly anisotropic near a vacancy. In the vicinity of the vacancy, theatomistic details account for the discrepancy between numerical results for the displacementfields and those predicted by continuum elasticity theory. Beyond fifteen lattice spacings,this continuum theory accurately reproduced the results of numerical studies under someboundary conditions. The arrangement of low concentrations of vacancies may not simply be controllable be-cause of the nontrivial coupling between the pair potential function, the local deformationsthat it produces when a vacancy is introduced into a perfect crystal, and the effectivevacancy-vacancy interactions. However, a number of pair potential functions have been5ound 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 crys-tal as a robust ground-state structure. At positive temperatures, the honeycomb struc-ture remained stable for a large temperature and pressure region of the phase diagram. This double-well potential was generalized as the Lennard-Jones-Gaussian potential, whichshowed a number of complex crystal structures including pentagonal, hexagonal, nonago-nal, and decagonal unit cells as the depths and locations of the energy wells were varied. Colloidal particles with three “patches,” or attractive interaction sites on the surface of a par-ticle, have recently been shown to yield a honeycomb structure for some pressures at T = 0in addition to other low-coordinated crystal structures. In addition, the square-shoulderpotential, a hard-core potential with a soft corona of variable length, shares some of the sameground-state features as the quintic potential. In particular, low-coordinated structures suchas a honeycomb-like trivalent configuration and lanes are ground-state features that areshared with the quintic potential. Such a potential can also form lamellar and micellarphases for relatively long-ranged coronas as was shown theoretically. Lastly, we note that several potentials qualitatively similar to the quintic potential havebeen studied in other contexts. For example, a potential that gives rise to quasiperiodicorder and the square lattice in two dimensions has one local minimum and local maximumand finite cutoff, as does the Dzugutov potential which gives rise to a number of com-plex local clusters in three dimensions. However, both of these potentials have a localattractive minimum ( i.e. negative potential energy). A piecewise linear potential with ahard core, single local minimum, and single local maximum showed the formation of chainsand labyrinths at positive temperature, though the ground state is not fully characterized. However, in contrast to the above potentials, the quintic potential is short-ranged, positive,isotropic, continuous and differentiable, which may be beneficial to experimentalists hopingto realize this potential in a laboratory setting. We find that with this potential, a numberof low-coordinated structures arise as ground states including the honeycomb and kagom´ecrystals as well as a “striped” phase. This represents a significant achievement and a firststep toward controlling vacancies in ground-state structures.The remainder of this paper is as follows. The quintic potential is defined and detailed inSec. II while our methodology is detailed in Sec. III. In Sec. IV, we define the phase diagramin the specific area- H plane and in the pressure- H plane. We discuss the characteristics of6he phases, the results of simulation, and their mechanical stability. In Sec. V, we focus onthe special case where H = 0. This case yields anomalous behavior because the potentialenergy vanishes inside and outside the energy barrier. Lastly, we discuss the implications ofthis 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 truncatedbeyond √ f ( r ) = (cid:20) L ( m ) (cid:16) r − √ (cid:17) + K ( m ) (cid:16) r − √ (cid:17) − (cid:16) r − √ (cid:17) (cid:21) , r ≤ √ m is an unscaled height of the local minimum. The coefficients L ( m ) and K ( m ) are chosen so that r = 1 is a local minimum and f (1) = m , and, as afunction of m , are given respectively as L ( m ) = 4 m (cid:16) √ − (cid:17) − − (cid:16) √ − (cid:17) − , (2) K ( m ) = 5 m (cid:16) √ − (cid:17) − − (cid:16) √ − (cid:17) − . (3)The condition that f ′′ (1) > m must be constrainted to m < (cid:0) √ − (cid:1) ≈ . . (4)The position of the energy barrier (the relative maximum) occurs at r b = √ − (cid:0) √ − (cid:1) h − m (cid:0) √ − (cid:1) − i . (5)The generalized pair potential is rescaled so that v ( r b ) = 1 to set a uniform energy scale, v ( r ) = f ( r ) /f ( r b ) . (6)Upon rescaling, the height of the relative minimum H is defined so that v (1) = H andvaries from zero to unity. When H = 1, the potential is monotonically decreasing and is flatat r = 1. In this case, there is no local minimum or maximum. These soft potentials arerepulsive for small r , and the first and second derivatives vanish at r = √ 3. We construct7 m H FIG. 2: The relationship between m and H for the quintic potential. the potential such that the first and second derivatives vanish at r = √ e.g. phonons) arewell-defined. Examples from the family of quintic potentials are shown in Fig. 1. Becausedeveloping a simple expression relating m and H is difficult, a table relating m and H wasused. The relation between m and H is shown in Fig. 2. III. METHODS At T = 0, the free energy per particle g is identical to the enthalpy per particle andis related to the average potential energy per particle u via g = u + pa . To map theground-state phase diagram as a function of specific area a , pressure p , and H , we utilizedthe double-tangent construction method to find the lowest free-energy structure amongcandidate ground-state configurations. The double-tangent construction is demonstrated inFig. 3.We considered several crystal structures as candidate ground states including the trian-gular lattice (TRI), square lattice (SQ), honeycomb crystal (HC), kagom´e crystal (KAG),and its close relative, the “anti-kagom´e” crystal (AKG). Whereas the kagom´e crystal hasvacancies located on a triangular lattice separated by two nearest neighbor spacings, theanti-kagom´e crystal has vacancies located on a rectangular lattice. For all densities, theanti-kagom´e has a potential energy greater than or equal to the kagom´e lattice due to thediffering local coordination numbers. We considered these crystal structures because they8re relative simple crystal structures with varying degrees of local coordination. The con-struction of the potential, with a well at r = 1 and vanishing energy at r = √ n -particle crystals, peri-odic configurations containing n particles per unit cell, while allowing for shape deformationof the unit cell. For example, consider a two-particle crystal with lattice vectors [ d , 0] and[ d , d ] and a basis where particle is at (0 , 0) and particle is located at ( x , y ). For thissystem, the potential energy per particle is given as u = v ( r ) + 12 X i [ v ( r ′ ) + v ( r ′ ) + v ( r ′ ) + v ( r ′ )] , (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 u for a fixedarea per particle a and eliminating d = na/d casts this as an unconstrained minimizationthat is function of d , d , x and y . The function is minimized using the conjugate gradientalgorithm. It is straightforward to generalize this to a larger number of particles in the unitcell.To ensure that we obtain globally optimal solutions, we use a large number of initialconditions slowly and systematically varying d , d , x and y from zero to n √ 3. After eachminimization, we ensure that the unit cell is not significantly sheared. A highly shearedunit cell requires summation over a large number of unit cells to ensure that all interactionsare accounted for and therefore we discard those solutions with angles less than 10 ◦ . Thisprocedure is repeated for nearly two thousand values of a on the range 0 . ≤ a ≤ √ / u versus a curve is generated.These minimizations are a function of 2 n variables, and therefore many minimizationsat a fixed area per particle can be performed easily. However, as the number of particlesper unit cell grows, the number of initial conditions required to ensure global optimalitygrows as 2 n . We performed these minimizations systematically for up to four particles perunit cell, beyond which it becomes too intense to systematically explore thousands of initialconditions for thousands of specific areas.We also used slow cooling via molecular dynamics (MD) to obtain candidate ground-state9tructures for larger and possibly disordered systems. With a system of 800 particles in a pe-riodic box, the system was simulated using the Verlet algorithm and Andersen thermostat. The system was initialized as a high-temperature liquid and slowly cooled according to a lin-ear temperature schedule from T = 0 . × time steps. Afterthe simulation was terminated, the system was quenched to a potential energy minimumusing the conjugate gradient algorithm.Lastly, we performed lattice Monte Carlo (MC) optimizations using several hundredlattice sites. In these optimizations, the specific area was fixed and the lattice sites wereeither occupied or unoccupied by particles. The optimization variables were the occupationstate of the lattice sites, the angle between the lattice vectors and the length of one latticevector. The length of the other lattice vector was constrained by the area, number ofparticles, angle between the lattice vectors, and the length of the first lattice vector. PossibleMonte Carlo moves included inserting a particle to a vacant site, removing a particle froman occupied site, swapping a particle from an occupied site to an unoccupied site, perturbingthe angle of the lattice, and perturbing the length of one of the lattice vectors. The systemwas then optimized via simulated annealing. It was given an initial “temperature,” or energyscale, and moves were accepted or reject based upon the Metropolis algorithm. Ten thousandMC trials were attempted in each cycle and several hundred cycles were performed for eachoptimization.It is important to note that although many interesting structural characteristics arosefrom the results of molecular dynamics simulation and Monte Carlo optimization, no fi-nal configurations were identified as ground-state structures. The lattice sums and crystaloptimization methods always provided the lowest free-energy structures.The double-tangent construction requires care and precision since other phases can comevery close to the coexistence free energy. The p - g and u - a curves are discretized over severalthousand data points and linearized between between data points. The identification ofcoexisting phases can be done by using the lower, concave envelope of the p - g curve. Thecoexistence points and ranges of stability were then calculated using a tabulation of u , a , p and g and linear interpolation. 10 .5 1 1.5 2 2.5 3 a u Opt 3 PartOpt 2 PartOpt BravHCTriKagMD + QuenchMC FIG. 3: (Color online) Potential energy per particle u as a function of specific area a for selectedcrystals and the optimal 2- and 3-particle crystals for H = 0 . 5. The dense triangular, honeycomb,striped, and open triangular crystals are stable ground states for H = 0 . 5. The square and anti-kagom´e structures are omitted for clarity. Dashed lines represent the double tangent construction. IV. RESULTSA. Phase Diagram The double-tangent construction revealed two different types of phase behavior dependingon the value of H . In no case did the configurations resulting from MD or MC yield a groundstate. For 0 < H < . u - a curves for H = 0 . 5. In the figure, the optimal n -particle crystals are not visible when theyare identical to the triangular, honeycomb, or kagom´e crystals. As shown in in Fig. 3, athighest pressure, a dense triangular lattice (TRI), where neighboring particles lie inside thepotential energy barrier of other neighboring particles, is the ground state. At a reductionof pressure to p = 1 . p against g reveals an intersection between the curves for theTRI and HC structures. The kagom´e lattice has a higher free energy than both structuresin this density range. 11urther reduction of the pressure to 0.75263 yields a “striped” (ST) or lane phase coex-isting with the honeycomb phase. This phase is an affinely-stretched triangular lattice. Thefirst few coordination shells contain two, four, and two particles. The curves representingthe lowest-energy crystals with a two- and three-particle basis either the triangular lattice,honeycomb crystal, or kagom´e crystals generally for a < . a , these curves appear to come close to the coexistence line between the HC andST phases or the ST and open TRI phases. The p - g curves reveal that the free energy isalways above that of the coexisting phases. These near-ground states are generally crystalsthat nearly mimic the coexistence. For example, the optimal three-particle crystal betweenthe HC and ST phases is an alternating sequence of lines with honeycomb-like coordinationand striped coordination. These are suboptimal structures because the specific neighborlengths do not match exactly. At sufficiently low pressure, for H = 0 . 50 and p = 0 . √ a ≥ √ / 2. Inthis region, the ground states behave like hard disks (HD) of radius √ 3. There would bea hard-disk crystal and liquid regime with specific area scaled to the hard-disk equation ofstate, which has been well characterized. For 0 . ≥ H ≥ 1, the kagom´e crystal emerges as a ground state structure for a nar-row range of stability. This is an especially important result. In previous work, researchersused an inverse methodology to design pair interactions for targeted ground states. Al-though they were successful in engineering potential for the honeycomb and square crystals,they were unable to do so for the kagom´e. The area and pressure ranges of stability increasewith H . The emergence of the kagom´e crystal as a ground state is shown in the p - a curvefor H = 1 . 00 in Fig. 4. As the pressure is reduced, the sequence of stable phases is the denseTRI, KAG, HC, ST, open TRI, and HD phases.The full phase diagrams in the p - H plane and the a - h plane are shown in Figs. 5 and 6,respectively. Figure 5 shows that the coexistence lines are nearly linear for small H . Thepressure range of stability for the open TRI, ST, and KAG phases widens as H increases.When H is sufficiently large to include the KAG phase, the pressure range of stability forthe HC phase narrows. As H increases, the area range of stability increases for most phases,12 .5 1 1.5 2 2.5 a u Opt 3 PartOpt 2 PartOpt BravHCTriKagMD + QuenchMC FIG. 4: (Color online) Potential energy per particle u as a function of specific area a for selectedcrystals and the optimal 2- and 3-particle crystals for H = 1 . 00. The dense triangular, kagom´e,honeycomb, striped, and open triangular crystals are stable ground states for H = 0 . H = 1 . 00, the kagom´e also emerges as a stable ground state. The square and anti-kagom´e structuresare omitted for clarity. Dashed lines represent the double tangent construction. excluding the dense TRI phase, which is evident in Fig. 6.Near H = 1 . 00, the slopes of the curves change dramatically. This is due to the softeningof the pair potential function as H approaches unity and the fact that the relative minimumand maximum come together much more rapidly as H approaches unity. Fig. 1 shows thatthe H = 1 . 00 potential is significantly softer near r = 1 than other potentials.The softening of the potential and the rise of the relative minimum are important featuresfor the inclusion of the KAG phase. In comparing Fig. 3 to Fig. 4, one can see that the curvefor the kagom´e crystal initially begins above the double tangent connecting the TRI and HCphases. As the potential softens and the relative minimum increases, the first minimum ineach curve increases proportionally to H . However, the softening of the potential bends theleft-side of the curves in such a way that the KAG curve lies at lower free energy than thatTRI-HC coexistence line. The combination of the height of the relative minimum and thesoftening are directly related to the potential energy and pressure of the system, the twocomponents of the T = 0 free energy.The KAG phase emerges at a triple point at H = 0 . p = 2 . p H HD TRIKAGHCTRI ST FIG. 5: (Color online) Phase diagram in the p - H plane as calculated by double tangent method.The hard-disk-like state (HD) occurs for p = 0 and H > a H TRI TRI Hard DiskST-TRITRI-HC HC-STTRI-KAGKAG KAG-HCHC VRLST FIG. 6: Phase diagram in the a - H plane as calculated by the double tangent method. .8 1 1.2 1.4 1.6 1.8 2 2.2 r -1-0.500.511.52 v (r) ; - v ’ (r) / v(r)-v’(r) /2Tri 6,6HC 3,6Kag 4,4 FIG. 7: (Color online) Comparison of the location and coordination number of the triangularlattice, honeycomb crystal and kagom´e crystal for a = 0 . H = 0 . H and p = 2 . TRI, HON, and KAG phases form a triple point. We have examined the subtle interplaybetween the shape of the pair potential function and its derivative at this coexistence pointin Fig. 7. The figure shows v ( r ) and − v ′ ( r ) / H = 0 . p ∗ and free energy g ∗ , the system of equations becomes g ∗ = 3 " v ( d t ) + v ( d t √ − d t v ′ ( d t ) − d t √ v ′ ( d t √ (8) p ∗ = − √ d t h v ′ ( d t ) + √ v ′ ( d t √ i (9) g ∗ = 3 " v ( d h ) + v ( d h √ − d h v ′ ( d h ) − d h √ v ′ ( d h √ (10) p ∗ = − d h √ (cid:20) v ′ ( d h ) + √ v ′ ( d h √ (cid:21) (11) g ∗ = 2 v ( d k ) + 2 v ( d k √ − d k v ′ ( d k ) − d k √ v ′ ( d k √ 3) (12) p ∗ = − √ d k h v ′ ( d k ) + √ v ′ ( d k √ i , (13)where d i is the nearest-neighbor distance for each structure and is necessarily less than unityfor the quintic potential. Evidently, the quintic potential for H = 0 . B. Stability The positivity of the squared frequency of all phonons ensures mechanical stability of acrystal. We have examined the phonon spectra for the KAG, HC, and ST phases. Thesephases are confirmed to be mechanically stable within the relevant density and H ranges.Figure 8 illustrates the phonon spectra for certain points in the reduced first Brillouin zonefor the KAG, HC, and ST structures for H = 0 . ω max / ω min at M point a simple metric for the extent of mechanical stability, this ratio stays relatively fixedfor the KAG lattice as H increases. This marks an achievement since previous attemptsto stabilize the kagom´e crystal had been unsuccessful. For the HC structure, this ratio ishigher for H = 1 . 00 than for that which is plotted in Fig. 8, indicating that the crystal ismore stable as H increases. In general, the quintic potential and the honeycomb potential yield similar ratios. The striped phase whose lattice vectors for a = 1 . 51 and H = 0 . 875 are[1 . , 0] and [0 . , . ω K Γ M* (a) ω K Γ M* (b) ω Γ Γ Γ k y k x (c) FIG. 8: (Color online) Phonon spectra for (a) the kagom´e crystal at a = 1 . 035 and (b) thehoneycomb crystal at a = 1 . 2, and (c) the striped phase at a = 1 . 51 for wave vectors in thereduced first Brillouin zone for H = 0 . other “quadrants” are identical to the quadrant displayed in Fig. 8(c) due to the symmetryof the Brillouin zone. C. Low-Energy Structures Although those structure obtained by MD and MC methods were suboptimal, Figures 3and 4 show that the free-energy differences between these structures and the ground statesare small. Several structural motifs emerge for this potential and are shown in Fig. 9 for H = 0 . 5. These metastable states may have technological value if the freezing kinetics aresufficiently slow. Preliminary simulations suggest the freezing behavior varies with density.For example, with a < . 2, systems typically freeze into a triangular lattice with vacanciesrandomly distributed as in Fig. 9(a). Although this is not a ground-state structure, it yieldsa vacancy-riddled lattice. The vacancies appear to have no particular order. The freezingtransition, the temperature at which the system changes from a high-density liquid to anordered phase, appears to be first order. The drop in potential energy as the temperatureis reduced is sharp in this density range.Systems at densities where coexistence between the HC and ST phases are the groundstates tend to exhibit a freezing from the liquid phase to a rigid, structured phase. However,17he structured phase, which we believe to be metastable, is characterized by rings and stringsas in Fig. 9(b). The six-particle ring is a characteristic of the HC structure and the stringsare characteristics of the ST phase. These metastable states are disordered due to the fluidnature of the strings. The drop in potential with temperature is not as sharp for thesedensities compared to those at higher densities.Lastly, at lower densities, as with those associated with the ST and open TRI phases, themetastable, low-energy states adopt labyrinthine characteristics, Fig. 9(c), and eventuallycolloidal polymers and monomers, 9(d). In this density range, the drop in potential energyassociated with freezing is weak. The phenomena and characteristics of low-energy statesare common to all H . However, as H increases so does the propensity to form labyrinthinecharacteristics. This is due to the lower energy barrier and the decrease in r b , the locationof the energy barrier. A qualitatively similar, but piecewise linear, potential also yieldsa labyrinth phase at positive temperature. Although the positive temperature behaviorhas not been fully explored in this paper, we anticipate that the equation of state for thisfamily of potentials will have interesting behavior. Manipulation of the kinetics to achievesuch unusual structures represents one path to achieving kinetically stable materials withcontrolled vacancy concentrations. V. SPECIAL CASES: H = 0 The H = 0 pair potential has interactions that vanish at a pair distance of unity andbeyond √ 3. At positive pressure, the ground state is the dense triangular lattice. However,at zero pressure and a ≥ √ / 2, ground-state configurations will have a vanishing poten-tial energy and pressure (enthalpy vanishes). Thus, the type of available ground states isdependent on the area.A number of interesting structures arise as ground states. For a > √ / 2, dilutions of thetriangular lattice with unit neighbor spacing are ground states, including the honeycomband kagom´e crystals and other lattice gases. For a = √ / 2, the striped phase, a Bravaislattice with lattice vectors of [1 , 0] and [1 / , √ / √ a > √ / 2, dilutions of this lattice are also ground states. In addition, any smallexpansion in the direction normal to the stripes, will allow each stripe to gain some fluidity.18 a) a = 0 . a = 1 . a = 1 . a = 2 . FIG. 9: (Color online) Metastable configurations obtained via molecular dynamics for H = 0 . Complex liquids of strings have found as ground states using molecular dynamics for at least a ≥ . 15. However, the geometric problem is highly nontrivial since for several runs with a > . 15, ground states could not always be obtained. For a = 3 √ / 2, an open triangularlattice with neighbor spacing of √ e.g. lattice gas) and continuous en-tropy ( e.g. hard-disk crystal). In the classical case, continuous entropy is uncountable whilediscrete entropy is not. Making the analogy to hard disks and hard spheres, we believe thatthe highest entropy phases for a specified area are those that have the availability for contin-uous entropy. Configurations with the highest-dimensional configuration space dominate the T = 0 entropy. We estimate that, for a certain a , configurations with the fewest constraineddegrees of freedom will have the highest dimensional configuration space, and therefore thehighest entropy. For a ≤ √ / 2, all degrees of freedom per particle are constrained. For a ≥ √ / 2, no degrees of freedom per particle are constrained.In order to gain continuous entropy, the system must have the ability to expand the √ √ √ 3. These are depicted in Table Ialong with their shorthand notations.Because there are no constraints beyond distances of √ 3, the edges with lengths of √ • 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 thevertices also match. • The (1,1,r) and (1,r,r) tiles cannot share a √ i.e. two vertices are separated by a length between unity and √ ABLE I: Available triangular tiles for the H = 0 tiling problem.Name Shape Angles ( ◦ ) Area(1,1,1) 60, 60, 60 √ / √ / √ / √ / in Fig. 10(a) where each vertex is decorated with a particle. These tilings are simply di-lutions of the dense triangular lattice, or coexistence between a dense and open triangularlattice. The (1,1,r) tiles have two roles. They can act as intermediaries between domainsof the dense triangular lattice and the open triangular lattice, and they can dimerize toform a rhombus making a domain of the dense triangular lattice. This coexistence betweenthe dense and open triangular lattices is represented as a double-tangent line connectingthe fully constrained dense triangular lattice to the unconstrained open triangular lattice asshown in Fig. 11.Next we consider those tilings that include (1,r,r) tiles. Two periodic tilings with aspecific area of a = √ / ≈ . √ √ a) (b) (c)(d) (e) FIG. 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 thestriped phase, and (e) coexistence between the open triangular lattice and the striped phase. the the striped phase and the open triangular lattice. The double tangent construction ofthese coexistences is displayed as the solid red line in Fig. 11. These coexistences between thestriped phase and the dense triangular lattice and the striped phase and the open triangularlattice maximize the number of unconstrained degrees of freedom and likely maximize theentropy of the ground state. In addition, these system can take on discrete entropy by wayof stacking variants of striped phase and the lattice phase. Using “S” and “T” to denote astripe and a triangular lattice layer, a STTTTSSS system is degenerate with a STSTSTSTsystem.The zig-zag phase can only form local coexistence with the open triangular lattice. To doso, the (1,r,r) tiles must form a trimer which exposes the √ .5 1 1.5 2 2.5 3 a = A/N C on f i g . S p ace D i m . / N Dense/Open Tri CoexDense/ST, ST/Open Coex FIG. 11: (Color online) Unconstrained degrees of freedom per particle for candidate ground-statestructures for H = 0 potential. Configurations with the maximum number of unconstrained degreesof 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 potentialthat 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 fullphase 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 honeycombpotential, there was no gas-liquid coexistence in the low-temperature, low-density part ofthe phase diagram, nor was there a liquid-liquid phase coexistence We expect most ofthe solid-solid transitions that we find at zero temperature will remain at small nonzerotemperature. Our preliminary calculations suggest that strings, or polymers, may arise inequilibrium at low densities.In addition, we have interest in the mobility of vacancies, particularly in the cases wherevacancy concentrations are dilute, ( e.g. H = 0 and a ≈ √ / For a just above √ / H = 0, we expect there to be a smallnumber of vacancies in the system at low-temperature. We have estimated the potential23nergy required for a particle to “hop” from a lattice site to a vacant site. By initializing asystem with one vacancy with a “hopping” particle along the transition to the vacant site,in an otherwise undistorted lattice, we minimized F ( r N ) = |∇ u ( r N ) | , where r N representsparticle coordinates, using conjugate gradient minimization. The resulting configurationrepresents a saddle point in the energy landscape. The difference in potential energy of thesaddle point and the ground state is the energy barrier required for the particle and vacancyto swap positions. This barrier sets an activation energy for classical thermal motion. Forthe quantum case, it would also enter in the tunneling rate. For H = 0 potential, there is asingle saddle point in the vicinity of the numerous initial conditions in the energy landscapewith a total energy barrier height of 5.569015. The diffusing particle is midway between theorigin and destination sites while the bracing particles are displaced off the lattice line toaccommodate the jump. Using molecular dynamics, we expect to relate this saddle pointenergy to a vacancy diffusion coefficient.The quintic potential can further be generalized by varying the distance at which thefunction is truncated. We set this cutoff distance to be √ 3. However, allowing this cutoffdistance to vary introduces a larger class of potentials. A systematic study of the cutoffradius on the robustness of the ground states is necessary for experimental realization. Thehard-core plus square shoulder potential has ground states that vary significantly dependingon the relative lengthscale of the hard-core distance and the square-shoulder distance It is expected that the ground states of the generalized quintic potential would be sensitiveto the location of the minimum and the cutoff distance, though it is currently unknownhow sensitive the ground-state phase diagram is to these parameters. Understanding thissensitivity is important to experimentalists who want a simple, robust potential. Developingan 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 develop-ing potentials for targeted material properties. We intend to develop a general and broadinverse optimization technique to target specific vacancy arrangements by accounting forand/or manipulating long-ranged vacancy-vacancy interactions. For example, one might de-velop an objective function whose variables include the strength, sign (attractive/repulsive),and angular dependence of vacancy-vacancy interactions. Alternatively, one might considera two-component system system of a heavy particle and a light particle and apply an in-24erse optimization technique to this system. Using a broad family of potentials, one couldthen optimize over the available parameters to achieve a dilute concentration of effectivelyrepulsive vacancies. VII. ACKNOWLEDGMENTS We thank Lawrence Cheuk for initial work on a related model. S.T. thanks the Institutefor Advanced Study for its hospitality during his stay there. This work was supported by theOffice of Basic Energy Sciences, U.S. Department of Energy, Grant DE-FG02-04-ER46108.. ∗ Corresponding author; Electronic address: [email protected] J. Sharma, R. Chhabra, A. Cheng, J. Brownell, Y. Liu and H. Yan, Science , 2009, , 112–116. A. P. Hynninen, J. H. J. Thijssen, E. C. M. Vermolen, M. Dijkstra and A. van Blaaderen, Nature Materials , 2007, , 202–205. Y. K. Cho, R. Wartena, S. M. Tobias and Y. M. Chiang, Adv. Func. Mat. , 2007, , 379–389. S. Torquato, Soft Matter , 2009, , 1157–1173. H. Cohn and A. Kumar, Proc. Nat. Acad. Sci. , 2009, , 9570–9575. W. B. Russel, D. A. Saville and W. R. Schowalter, Colloidal dispersions , Cambridge UniversityPress: Cambridge, 1989. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett. , 2005, , 228301. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E , 2006, , 011406. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. E , 2006, , 021404. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev, E , 2007, , 031403. M. C. Rechtsman, F. H. Stillinger and S. Torquato, Phys. Rev. Lett. , 2008, , 85501. M. C. Rechtsman, F. H. Stillinger and S. Torquato, J. Phys. Chem. A , 2007, , 12816–12821. R. D. Batten, F. H. Stillinger and S. Torquato, J. App. Phys. , 2008, , 033504. J. D. Joannopoulos, S. G. Johnson, J. N. Winn and R. D. Meade, Photonic crystals: moldingthe flow of light , Princeton University Press: Princeton, N.J., 2008. R. C. Agrawal and R. K. Gupta, J. Mat. Sci. , 1999, , 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, , 1107. E. Kim and M. H. W. Chan, Nature , 2004, , 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, , 4692–4712. E. Cockayne and V. Elser, Phys. Rev. B , 1991, , 623–629. L. Cˆandido, P. Phillips and D. M. Ceperley, Phys. Rev. Lett. , 2001, , 492–495. L. Modesto, D. L. S. J´unior, J. N. Teixeira Rabelo and L. Cˆandido, Solid State Communications ,2008, , 355–358. W. Lechner and C. Dellago, Soft Matter , 2009, , 646–659. W. Lechner and C. Dellago, Soft Matter , 2009, , 2752–2758. W. Lechner, E. Sch¨oll-Paschinger and C. Dellago, J. Phys. Condens. Matter , 2008, , 404202. A. P. Hynninen, A. Z. Panagiotopoulos, M. C. Rechtsman, F. H. Stillinger and S. Torquato, J.Chem. Phys. , 2006, , 024505. M. Engel and H. R. Trebin, Phys. Rev. Lett. , 2007, , 225505. M. Engel and H. R. Trebin, Zeitschrift f¨ur Kristallographie , 2008, , 721–725. G. Doppelbauer, E. Bianchi and G. Kahl, J. Phys.: Condensed Matter , 2010, , 104105. J. Fornleitner and G. Kahl, Eur. Phys. Lett. , 2008, , 18001. J. Fornleitner and G. Kahl, J. Phys: Condes. Matter , 2010, , 104118. M. A. Glaser, G. M. Grason, R. D. Kamien, A. Koˇsmrlj, C. D. Santangelo and P. Ziherl, Eur.Phys. Lett. , 2007, , 46004:1–5. A. Quandt and M. P. Teter, Phys. Rev. B , 1999, , 8586–8592. M. Dzugutov, Phys. Rev. A , 1992, , 2984–2987. J. P. K. Doye, D. J. Wales and S. I. Simdyankin, Faraday Discussions , 2001, , 159–170. M. D. Haw, Phys. Rev. E , 2010, , 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, , 359–361., 359–361.