Spontaneous formation of thermodynamically stable Al-Cu-Fe icosahedral quasicrystal from realistic atomistic simulations

Icosahedral quasicrystals spontaneously form from the melt in simulations of Al-Cu-Fe alloys. We model the interatomic interactions using oscillating pair potentials tuned to the specific alloy system based on a database of density functional theory (DFT)-derived energies and forces. Favored interatomic separations align with the geometry of icosahedral motifs that overlap to create face-centered icosahedral order on a hierarchy of length scales. Molecular dynamics simulations, supplemented with Monte Carlo steps to swap chemical species, efficiently sample the configuration space of our models, which reach up to 9846 atoms. Exchanging temperatures of independent trajectories (replica exchange) allows us to achieve thermal equilibrium at low temperatures. By optimizing structure and composition we create structures whose DFT energies reach to within ∼2 meV/atom of the energies of competing crystal phases. Free energies obtained by adding contributions due to harmonic and anharmonic vibrations, chemical substitution disorder, phasons, and electronic excitations, show that the quasicrystal becomes stable against competing phases at temperatures above 600 K. The average structure can be described succinctly as a cut through atomic surfaces in six-dimensional space that reveal specific patterns of preferred chemical occupancy. Atomic surface regions of mixed chemical occupation demonstrate the proliferation of phason fluctuations, which can be observed in real space through the formation, dissolution and reformation of large-scale icosahedral motifs—a picture that is hidden from diffraction refinements due to averaging over the disorder and consequent loss of information concerning occupancy correlations.


I. INTRODUCTION
Since the discovery of quasicrystals as a distinct phase of matter [1], and recognition of their quasiperiodicity [2], two fundamental questions remain to be definitively answered: Where are the atoms [3]? What stabilizes their quasiperiodic order? Excellent descriptions of their average structures are possible in terms of cuts through higher-dimensional periodic lattices obtained by single-crystal diffraction refinements [4]. However, quasicrystalline structures can only be reliably equilibrated at high temperatures; consequently, these models contain ambiguous atomic positions with uncertain occupation and chemistry. They omit important correlations in the case of mixed or partial occupation, and they omit atomic vibrations and diffusion. As regards their thermodynamic stability, local icosahedral motifs are clearly favored energetically. Local preference need not force long-range quasiperiodicity, as is illustrated by the prevalence of periodic "approximants," which mimic quasiperiodic order locally within an ordinary crystalline unit cell that in turn repeats periodically [5,6].
An intriguing puzzle that has eluded researchers for three decades is identifying the mechanism that selects an ordered yet nonperiodic state. One possible explanation is that quasicrystals are energetic minima, whose structure is forced by specific interatomic interactions [7], by maximizing the density of some favorable motif [8,9], or by creating a deep pseudogap in the electronic density of states [10][11][12]. Another possibility is that the structural ambiguity is an intrinsic characteristic of quasicrystals [13][14][15]. In this view, the entropy to be gained from chemical or positional fluctuations serves to reduce the free energy relative to competing phases whose energies (without entropy) are lower. Quasiperiodicity might arise spontaneously because it allows icosahedral symmetry, and this high symmetry maximizes the entropy.
The icosahedral phase of Al-Cu-Fe is an excellent place to seek theoretical insight. Experimentally, the i-phase of Al-Cu-Fe was the first well-ordered thermodynamically stable quasicrystal to be discovered [16]. It exhibits a particular symmetry classified as face-centered-icosahedral, and exhibits strong pseudogap in electronic states density near Fermi energy [17]. Recently, these quasicrystals attracted renewed attention following their discovery in a meteorite [18].
In this paper, we report computer simulations leading to spontaneous formation of icosahedral quasicrystals from the melt. Such simulations have been previously reported [19][20][21], but only for artificial models that do not describe actual chemical species and hence yield no direct insight into specific experimentally studied compounds. Here we model Al-Cu-Fe ternary alloys using a combination of chemically accurate density functional theory (DFT) calculations and accurate interatomic pair potentials. Based on structure models and thermodynamic data provided by our simulation, we calculate a temperature-dependent phase diagram for the Al-rich region of the alloy system showing that energy and entropy conspire in the emergence of the quasicrystal phase as thermodynamically stable at elevated temperatures. Figure 1 illustrates the diffraction patterns of our simulated structures and verifies that we indeed obtain a quasicrystal with the expected face-centered icosahedral symmetry.

II. SIMULATION METHODS
Three important ingredients enable the success of our atomistic simulations: realistic DFT-derived interatomic interactions; appropriately sized simulation cells with periodic boundary conditions; efficient hybrid Monte Carlo/molecular dynamics augmented by replica exchange. We exploit the data generated by our simulation to present optimized structure models revealing icosahedral order on a hierarchy of length scales. Combining DFT-calculated formation enthalpies with entropies derived from fluctuations, we calculate the absolute free energies of the quasicrystal phase and competing crystalline phases. Then, from the convex hull of the set of free energies we predict a temperature-dependent phase diagram that shares characteristics with the experimentally assessed behavior, including the emergence of the quasicrystal as a high temperature stable phase.
Oscillating interatomic pair potentials accurately describe elemental metals and alloys characterized by weakly-bound s electrons, even in the presence of s-pd hybridization. While these can be derived analytically within electronic density functional theory [22][23][24], we instead employ a parametrized empirical form known as EOPP [25] that we fit to a database of DFT energies and forces (Appendix A). This form has found success modeling many binary and ternary alloys [26][27][28]. EOPP also lead to spontaneous formation of singlecomponent icosahedral quasicrystals [21]. Henley [29] connected the oscillating potentials with pseudogap near E F via the Hume-Rothery scenario, and argued that their second minima participate in formation of fundamental clusters. Indeed, under some circumstances the second minima can even create local matching rules favoring quasiperiodicity [30]. Figure 2 illustrates our fitted potentials and the comparison with interatomic separations in our simulated structures. Details of our fitting procedure are provided in Appendix B.
Because the quasicrystal is aperiodic it cannot be precisely represented in a finite-size system. Fortunately, a series of "rational approximants" are known that capture the local quasicrystal structure and minimize the deviation caused by application of periodic boundary conditions. Reflecting the hierarchical nature of quasiperiodic order, these special sizes grow geometrically as a cub = 2a q τ n / √ τ + 2, where τ = (1 + √ 5)/2 ≈ 1.61803... is the golden mean and a q is "quasilattice parameter" or Penrose rhombohedron edge length [5]. The volume grows rapidly as a 3 cub ∼ τ 3n . Strictly, face-centered icosahedral symmetry implies τ 3 scaling for FIG. 2. Empirical oscillating pair potentials V αβ (r) (thick solid lines, energy axis on the left) and partial pair correlation functions g αβ (r) (dashed lines, vertical axis on the right) for Al-Cu-Fe quasicrystals. The EOPP are fit to a database of DFT energies and forces (see Appendix B). The pair correlation functions are obtained from single snapshot of the 9846-atom 8/5 approximant EOPP MD simulation. self similarity (volume ∼τ 9 ) along fivefold and threefold directions, but τ 1 scaling along twofold icosahedral directions. Here, we label the approximants with a ratio of successive Fibonacci numbers F n+1 /F n , with 128 atoms for "2/1" up to 9846 for "8/5." Our naming convention is drawn from the definition of Henley's canonical cells designed for packing icosahedral clusters [31]. Since the fundamental clusters in AlCuFe are τ -times smaller than the proper Mackay clusters decorating the B.C.C. lattice in α-AlMnSi [5,32,33], our "2/1" approximant has about the same edge length a cub ∼ 12.3 Å as α-AlMnSi.
For small cell sizes (2/1 and 3/2-2/1-2/1), imposing this special length encourages nucleation of the quasicrystal from the melt [34]. In larger cells (3/2, 5/3, and 8/5), the entropic barrier to nucleation is hard to overcome; instead we seed the structure using the previous approximant size. Because of the discrete cell size scaling, a single unit cell of the seed occupies less than 24% of the cell volume. Thus we take a supercell of the smaller approximant and enforce the periodic boundary of the bigger approximant. Near the large approximant boundary we let the small approximant overlap itself in a region that is 10% of the large approximant cell size. This introduces a 25% excess of atoms in the large approximant placed at unphysically short separations. We remove the excess atoms through a fixed-site lattice gas annealing [35] to reach the desired atomic density.
To efficiently anneal both chemical and positional order we utilize a hybrid Monte Carlo/molecular dynamics (MC/MD [36]) method that samples continuous evolution of atomic positions through molecular dynamics while enabling interchange of chemical species through Monte Carlo swaps. Species swaps are accepted with the Boltzmann probability exp (− E /k B T ) with E the change in total energy for the swap. Our simulations are performed in the canonical ensemble with constant temperature, volume, and numbers of atoms of each species. To achieve equilibration at low temperatures and enhance sampling of the configurational ensemble at all temperatures, we supplement our MC/MD simulation with replica exchange [37]. MC/MD simulations are performed in parallel at many temperatures T i . At fixed intervals we suspend the simulations and consider swapping configurations at adjacent temperatures T i and T i+1 . The swap is accepted with a Boltzmann-like probability based on the energy difference between the configurations. Although the temperature of a configuration jumps during the swap, the configuration remains a properly weighted member of the equilibrium ensemble at its instantaneous temperature. Further details are in Appendix C.

III. RESULTS
We analyze the simulated structures to demonstrate their quasiperiodicity. The clearest demonstration is their diffraction pattern (Fig. 1) which shows characteristic 2×, 3×, and 5× patterns with sharp peaks near the characteristic positions for face-centered icosahedral symmetry (minor deviations occur due to the finite-size approximant). The face-centering causes certain diffraction peak positions to occur in ratios of τ 3 , while the remainder show τ 1 scaling.
Examining the structure in real space, we observe that small approximants solidify into well-ordered structures that can be conveniently described as packings of two cluster types. One is a small 13-atom (Al 12−x Cu x )Cu icosahedron that we denote as I [see Fig. 3(a)]. The other is a larger pseudo-Mackay icosahedron (pMI) [ Fig. 3(b)] with a partially occupied Al 12−x Fe inner shell (x ∼ 1-3 due to Al-Al repulsion), and a second shell made up of two subshells: a (Cu,Fe) 12 "unit-icosahedron" on the fivefold axes at 4.45 Å from the center, and an Al-rich (Al,Cu) 30 icosidodecahedron on the twofold axes. In the example shown in Fig. 3(b), the Cu and Fe atoms on the unit-icosahedron segregate to break the symmetry from icosahedral to fivefold. I clusters connect along twofold icosahedral directions with a spacing of b = 7.55 Å, and alternate with pMI along threefold directions (c = 6.54 Å spacing). This even/odd alternation implements the face-centered icosahedral order. For the 2/1 approximant this packing is a unique structure producing A-, B-, and C-type canonical cells [31] (CCT), while the 3/2-2/1-2/1 approximant also contains a symmetry-broken D-cell.
Larger approximants avoid the bulkiest canonical cell D by introducing a new three-shell cluster extending to ∼7.7 Å [twofold radius; see Fig. 3(c)]. This cluster is entirely bounded by CCT Y -faces, hence it extends, rather than violates, the CCT concept. Its innermost shell is Al 12−x Fe as in usual pMIs, but the second shell, albeit topologically similar to pMI, is Cu-rich with only 25% Al and no Fe atoms. Finally, the third shell is made up from three icosahedral subshells: Al 60 (at 6.4 Å), a Fe 30−x Cu x (x ∼ 6) τ -icosidodecahedron (at twofold radius 7.7 Å) and a Cu 12 τ 2 -icosahedron (at fivefold radius 7.2 Å). These 12 Cu atoms are in fact all centers of the small I clusters; upon including them the whole cluster has 282 atoms.
The three clusters provide a simple, highly economical zero-order description of the structure, since they cover practically all atoms in the structure (99% in 3/2, 98% in 5/3, and 97% in 8/5 approximant). A typical example of these clusters as they appear in our simulations is shown in Fig. 3(d), taken from the equilibrium ensemble at T = 1200 K, followed by relaxation. Mixed chemical occupation breaks the cluster symmetry and can serve as a stabilizing source of entropy, while, together with the cluster covering, it is potentially a means of forcing quasiperiodicity [8,9].
The ensemble of structures can be represented using the six-dimensional (6D) cut and project scheme [3,4]. This is an efficient way to represent the average structure in a manner that automatically enforces perfect quasiperiodicity. Briefly, we define a set of three-dimensional (3D) volumes, termed "atomic surfaces," placed on a 6D hypercubic lattice. A 3D "parallel space" cuts through the 6D space, oriented at an angle such that the 6D lattice edges project onto the six fivefold axes of a 3D icosahedron in the parallel space. Then, atoms are placed in parallel space at the points of intersection with the atomic surfaces and then averaged over the ensemble. (c) Four-shell τ -pMI cluster with icosahedral Al 60 Cu 12 (Fe,Cu) 30 outer shell encapsulating an Al 12−x Fe inner shell and a Cu-rich (Al,Cu) 42 pMI-like second and third shell. (d) Snapshot of 5/3 approximant from a simulation at T=1200 K that we then relaxed. Small icosahedra I are outlined by green circles, pMI by red, and an azimuth of the τ -pMI by black. Bonds between small icosahedra shown in blue lie along 2-fold directions, bonds between pMIs shown in red lie along fivefold directions, and bonds between Is and pMIs shown in orange lie along threefold directions. Tick marks are placed at 1 Å intervals. than the wells of the Al-TM potential. Three atomic surfaces emerge (see Fig. 4), two large ones (AS 1 and AS 2 ) sit at hypercubic lattice sites ("nodes"), one even and the other odd. The remaining small one (B 1 ) sits at the hypercubic body center. Each atomic surface has a particular pattern of chemical occupation. AS 1 is primarily Fe, concentrated at the center, with Cu surrounding and finally traces of Al. AS 2 is primarily Al, with a small concentration of Fe at the center and Cu separating the Al from the Fe. The remaining surface B 1 , at the body center, is primarily Cu. The contrast between AS 1 and AS 2 , along with absence of B 2 , reflects the strong symmetry breaking from simple icosahedral to face-centered lattice. There is a unique connection between the three fundamental clusters (Fig. 3) and the three atomic surfaces. The B 1 surface Cu atoms are centers of the I clusters, the central Cu/Fe part of the AS 2 surfaces occupy pMI clusters, and the Fe core atoms of the AS 1 surface are centers of the large τ -pMI clusters.
Notice the smooth variations in color (i.e., chemical occupancy). Cu separates Al from Fe on the atomic surfaces while blending continuously into each. Curiously, the location of Cu at the boundary of Fe and Al on the atomic surface is consistent with the position of Cu in the periodic table between transition (d-band) metals and nearly free electron (sp-band) metals. Mixed occupation implies swaps in chemical occupation in real space, and low atomic surface densities correspond to partial site occupation. Because these lead to spreading of the occupation domains in perpendicular space, we may identify these chemical species swaps and fractional occupation as types of phason fluctuations.
Our simulated atomic surface occupation largely agrees with the popular Katz-Gratias model [38]. Specifically, the node vertex surfaces (AS 1 and AS 2 ) both transition from Fe at the center, through Cu, to Al around the edges. Only a single body center (B 1 ) is occupied, solely by Cu. Shortdistance constraints determined specific shapes of the KG atomic surfaces. In our model these constraints are obeyed through positional correlations so our surfaces have outer shapes that differ from the KG model. Note that the outer shapes are defined by regions of low occupation probability.
After annealing down to low temperatures using the interatomic potentials, we apply DFT to relax the structures to T = 0 K and compute their enthalpies of formation, H relative TABLE I. Chemistry, instability E (letter "S" for stable), and formation enthalpy H for quasicrystal approximants and competing ternary phases at T = 0 K. The two variants of the B2-type β-phase, resulting from annealing simulations, are a 2 × 2 × 2 Fericher supercell with Pearson symbol tP16 and a 3 × 3 × 3 supercell with ordered vacancies and Pearson symbol cP50). to pure elements, and energetic instabilities, E relative to the tie-plane of competing crystal structure enthalpies. Table I summarizes the compositions and formation enthalpies of several competing phases. The small approximants exhibit an unusual electronic density of states (eDOS), with a wide and deep pseudogap as is usual in Al-based quasicrystals, and in addition a deeper and very narrow second pseudogap inside the broad pseudogap (see Appendix F). For optimal density and composition, which we achieve by replacing Cu with a combination of Al and Fe, the Fermi level lies inside this second pseudogap. Specifically, we find the effective valence rules of Al = +3, Cu = +1, and Fe = −2 apply, so we can raise or lower E F by 1 electron without altering the number of atoms through targeted chemical substitutions such as 2Cu ↔ Al + Fe. Composition can be shifted at constant effective valence through substitutions such as 3Al + 2Fe ↔ 5Cu. This rule matches the slope of the quasicrystal and approximant phase fields in the ternary composition space [39]. We discovered that neighboring Fe-Fe pairs lead to states in the pseudogap that can be removed by avoiding these pairs. These optimizations can substantially lower the total energy. However, these structures remain unstable by 2-4 meV/atom relative to competing ordinary crystal phases in the triangle η 2 (AlCu)-λ(Al 3 Fe)ω(Al 7 Cu 2 Fe) as described in Table I. Simulated larger approximants, and supercells of the 2/1 approximant, exhibit only the broad pseudogap. Apparently the higher entropy available in larger simulation cells introduces disorder that washes out the detailed structure leading to the narrow second pseudogap, trading a gain in energy for a compensating gain in entropy. It is conceivable that the EOPP potentials are not sufficiently accurate to capture the interactions responsible for the narrow pseudogap feature.
Notice the sequences of enthalpies of formation, H, that decreases monotonically with increasing approximant size. This suggests a possible energetic mechanism favoring quasiperiodicity. However, this is not yet clear, as we have not demonstrated that the energetically optimized structures are more perfect in their quasiperiodicity than representative high temperature structures. Indeed, the decreasing enthalpy is primarily a reflection of increasing Fe content. From the energies relative to the convex hull, E , which are positive and not systematically decreasing, it is likely that any quasicrystal model will be unstable relative to competing crystal phases at low temperatures. Hence, to explain the formation of the quasicrystals from the melt at high temperatures we must seek either a kinetic or a thermodynamic argument.
We consider the larger approximants, 3/2 (552 atoms) and 5/3 (2324 atoms), as representative of the true quasicrystal. They exhibit interesting structures and dynamics at elevated temperatures. Because clusters are difficult to identify in individual snapshots due to chemical disorder and atomic displacements, it is best to examine time averages of the structure. In Fig. 5 the inner shells of the pMIs are clearly visible as smeared circles due to the high mobility of the Al atoms, whose positions are frustrated by the incompatible length scale of the icosahedral potential produced by the outer shells with the short-range repulsion of the Al-Al potential. An azimuth of the three-shell τ -pMI clusters with Cu-rich interior is indicated by large black circles. Cu-centered small icosahedra are also clearly visible. As time evolves, the identity of these clusters shifts, with some becoming more distinct while others dissolve. Occasionally, the structures take pleasing hexagon-boat-star-decagon (HBSD [35,40]) tiling forms as as outlined in yellow in Fig. 5, but these structures, in turn, further evolve. We provide a video illustrating the evolving structure in our Supplemental Material [41].
We find a curious behavior at very high temperatures as the quasicrystals melt. One aspect is that the melting point grows as the approximant size increases. The 2/1 approximant (in 2 × 2 × 2 supercell) melts at T = 1669 K, the 3/2 at T = 1683 K, the 5/3 at T = 1723 K, and the 8/5 at T = 1788 K. Additionally, the 2/1 and 3/2 melt in a single step, while the 5/3 melts in two steps (the first broad heat-capacity maximum at 1590 K), and the 8/5 have three additional smaller but sharp heat capacity peaks at 1266 K, 1471 K, and 1680 K that precede melting. Melting of the larger approximants begins around 1266 K with small Al-Cu-rich regions [see the Curich sides of the pMI and τ -pMI clusters in Figs. 3(b) and 3(c) and the center of the τ -pMI in Fig. 3(d)] that leave behind an Fe-richer, and hence more mechanically stable, solid quasicrystal phase. Because we perform our simulations at fixed volumes, the actual melting occurs at high pressure-we estimate P = 7 GPa at the melting point of 8/5. Note that this is consistent with recent experiments of Bindi [42]. A video showing liquid-quasicrystal phase coexistence is provided in our Supplemental Material [41]. In this video the quasicrystal melts and resolidifies as the liquid interface advances into the solid and then recedes.
Finally, we seek to explain the thermodynamic stability. As shown in Table I our lowest energy quasicrystal models remain above the convex hull of energy by 2-5 meV/atom. Thus at low temperatures we anticipate phase separation into the competing phases, η 2 (AlCu)-λ(Al 3 Fe)-ω(Al 7 Cu 2 Fe). Indeed, this is what is shown in the standard phase diagram for the Al-Cu-Fe system [43]. Finite temperature stability is given by the convex hull of free energy. As discussed in Appendix G the free energy includes harmonic vibrational free energy F h derived from phonons, electronic free energy  Comparing free energies of the competing ordinary crystals, small quasicrystalline approximants and the quasicrystal (considered as a large approximant) we predict the stable phases at various temperatures by evaluating the convex hull of free energies over the composition space (see Fig. 6). Details are presented in Appendices E and F. Notably, we observe that the 2/1 approximant gains stability at T = 0 K when quantum zero point vibrational energy is included, and the 5/3 approximant, which we take as a proxy for the quasicrystal emerges as a stable phase at temperatures above T = 600 K owing to the excitations contained in F a .

IV. DISCUSSION
We have provided four pieces of evidence demonstrating the appearance of a thermodynamically stable quasicrystal state in our model of Al-Cu-Fe. First is the spontaneous formation of quasicrystal approximants, directly from the melt in the case of small approximants, and with the assistance of smaller approximant seeds in the case of large approximants. Second is the appearance of progressively larger clusters and superclusters with increasing approximant size. Third is the decreasing total energy of highly optimized structures with increasing approximant size, suggesting that energy could favor quasiperiodicity. Finally, we have calculated free energies for quasicrystal approximants and competing ordinary crystalline phases and foud that the 5/3 approximant, taken as a proxy for the true quasicrystal, acquires thermodynamic stability above T = 600 K.
Our findings shed light on the underlying reasons for quasicrystal formation and suggest there is not a single explanation but rather a coincidence of favorable conditions. Although the formation enthalpies decrease (i.e., become more negative; see Table I) with increasing approximant size, they actually lose stability relative to competing crystal states owing to the slope of the convex hull with increasing Fe content. The largest cluster for which we have DFT energies, namely, 5/3, is the smallest approximant that can accomodate the τ -pMI supercluster including its complete surrounding I clusters. If we postulate that this is an energetically favorable motif, then it may be that DFT energies are needed for yet larger approximants before we can claim existence of an energetic preference for the quasicrystal as compared with finite approximants. Further, the appearance of superclusters is a consequence of quasiperiodicity, not a cause. Indeed, both matching rule and random tiling models share this feature provided that these are viewed in a time-averaged sense as in the present case. While our explanation for thermodynamic stability at high temperatures is an example of entropic stabilization, the entropy includes phonon anharmonicity in addition to phason-related chemical species swaps and tile flips.
Our simulation goes beyond usual structure models of quasicrystals in that it captures the full dynamic evolution, including correlations among positions of partial and mixed occupation that represent phason fluctuations. At the same time it surpasses previous quasicrystal simulation models in the degree of fidelty to the specific chemical interactions of the constituent elements that we capture at or near the DFT level of accuracy. This allows us to evaluate the free energy of the real quasicrystal with accuracy sufficient to explain its formation.
Limited experimental information concerning the quasicrystal structure exists at present. The most important is the 6D Patterson function (a type of positional correlation function) measured by neutron diffraction on a single crystal sample [44]. Our results are consistent with that data, which was the first to identify the three atomic surfaces similar to those that we show in Fig. 4. A clear need exists for additional diffraction studies whose refinements could be compared more precisely with our average structure. Additional predictions could conceivably be tested experimentally, for example free energies could be measured, and our predicted low temperature stable 128-atom structure could possibly be produced by experiments that circumvent the diffusion barriers.
In conclusion, although our explanation for quasicrystal stability is not simple, it illuminates the complex interplay of multiple factors. These include the composition-dependence of competing phase energies, as well as multiple sources of entropy. Chemical preferences for site classes in cluster motifs and on atomic surfaces are favored by our EOPP interatomic potentials. We remark that cluster overlap together with cluster symmetry breaking provides a possible mechanism to force quasiperiodicity, but this appears insufficient to stabilize the Al-Cu-Fe quasicrystal against competing ordinary crystalline phases at low temperature. In our model, the quasicrystal is a high temperature phase.

APPENDIX A: DFT
First principles calculations within the approximations of electronic density functional theory (DFT) lie at the foundation of our structural and thermodynamic models. We employ projector augmented wave potentials [45] in the PW91 generalized gradient approximation [46] as implemented in the plane-wave code VASP [47]. Our k-point meshes are increased to achieve convergence to better than 1 meV/atom with tetrahedron integration. We employ the default energy cutoffs. For T = 0 K enthalpies, all internal coordinates and lattice parameters are fully relaxed. Finite difference methods are applied to calculate interatomic force constants that we need for low temperature vibrational free energies. Ab initio molecular dynamics (AIMD) was performed to generate additional energy and force data for fitting interatomic pair potentials. AIMD calculations contained 544 atoms in cubic cell and utilized only a single k-point. For accurate cohesive energies, approximant k-meshes were converged to 10 3 kpoints/BZ for 128-atom 2/1, 6 3 for 552-atom 3/2, and 2 3 for 2324-atom 5/3.

APPENDIX B: EMPIRICAL OSCILLATING PAIR POTENTIALS
We choose a six-parameter empirical oscillating pair potentials (EOPP [25]) of the form to fit a DFT-derived database of force components and energies. The database contains binary compounds Al 2 Cu (both tetragonal and cubic), Al 3 Fe, Al 5 Fe 2 , ternary tetragonal Al 7 Cu 2 Fe, orthorhombic Al 23 CuFe 4 , as well as the ternary  Table II. extension of Al 3 Fe, and a number of approximant structures. The database contains a significant portion of AIMD data at elevated temperatures, ranging from T = 300 K up to 2000 K, covering both solid and liquid configurations. In total we use 13176 force-component data points and 63 energy differences (see Fig. 7).
We initialized the fit from parameter values that fit GPT potentials [35] for a similar system (Al-Co-Ni). The fit quickly converged, with RMS deviation 0.16 eV/Å for forces, and 9.4 meV/atom for energy differences. Since Fe-Fe and Fe-Cu potentials are prone to softening at nearest-neighbor distances due to the lack of data in Al-rich systems, we increased repulsion term coefficients C 1 manually for Fe-Fe and Fe-Cu. Final parameters of our potentials are listed in Table II. Our cutoff radius is taken as 7 Å.
As a demonstration of the accuracy of our pair potentials, we computed the vibrational densities of states (VDOS) for a 208-atom "3/2-2/1-2/1" approximant (see Fig. 8). The three partials computed by EOPP semiquantitatively match the DFT result, with accuracy comparable to the Sc-Zn case [27].
It should be noted that the potentials are valid only for a particular density of the free-electron sea, and they should be used exclusively at constant volume; or in constant-pressure simulations, with additional external pressure set to a value leading to the same equilibrium volume.

APPENDIX C: REPLICA EXCHANGE SIMULATIONS
To enhance the efficiency with which our simulations explore the configurational ensemble, we employ a replica exchange mechanism [37], also known as parallel tempering, in which we perform many runs simultaneously at different temperatures. The probability for a configuration i of energy E i to occur at temperature T i is where (E ) is the configurational density of states, Z i is the partition function at temperature T i , and β i = 1/k B T i . The joint probability for configuration i at T i and configuration j at T j is Now consider swapping temperature between configurations i and j. The joint probability for configuration i to occur in the equilibrium ensemble at temperature T j and configuration j to occur at temperature T i is Hence, if the swap is performed with probability P /P = e β E , then equilibrium is preserved following the swap. This process works most efficiently if energy fluctuations are sufficiently large that the energy distributions H (E ) at adjacent temperatures overlap, so that swaps occur frequently. In this case, a given run (sequence of consecutive configurations) will diffuse between low and high temperatures. Rapid structural evolution at high temperatures thus provides an ongoing source of independent configurations for low temperatures where structural change is intrinsically slow.
We perform isochoric replica-exchange atomistic simulations in the temperature range of 200-1800 K, for several sizes where α is a multiplicative coefficient and N a number of atoms; such spacing guarantees uniform replica exchange acceptance rates assuming constant heat capacity. Consequently, for large systems an added load is due to the increasingly finer temperature grid. The parameters of most important simulations are summarized in Table III.

APPENDIX D: HYPERSPACE RECONSTRUCTION
Lifting a raw atomic structure to hyperspace amounts to associating the position of each atom to an ideal point in a 6D hypercubic crystal. We require that the projection of the ideal position into physical space match the actual atomic position to within a tolerance r core . Additionally, we require that the ideal position lie close in 6D space to the physical 3D space within "atomic surfaces" of definite positions, shapes and sizes. Specifically, we choose τ 2 -triacontahedra at 6D nodes, and unit triacontahedra at 6D body-centers, of radius r 6D =11.7 Å in the fivefold direction. In real space, the projected sites are separated by at least 1.05 Å. We obtain unambiguous registration with r core = 0.45 Å.
We register each structure by: (i) identifying all Al 12 Cu icosahedra in the actual atomic structure (centers of perfectly icosahedral clusters will have the smallest positional deviation from ideal sites); (ii) mapping these to the even-only bodycenter sublattice (hence resolving the even/odd ambiguity); (iii) Given the relative shift obtained from step (ii), we attempt to map the entire set of atoms positions to the ideal 6D sites. In practice about 80% of sites map, with the exceptions being primarily the disordered inner shells of the pMI .
Working with periodically bounded boxes in real space results in finite resolution of the perp-space, namely d ⊥ res = a q / √ F 2 n + F 2 n−1 , F n are Fibonacci numbers and a q = 4.462 Å. For our 8/5 approximant, d ⊥ res = a q / The atomic surfaces that result from the registration for our largest 8/5 approximant, after averaging over 3000 configurations, are shown in Fig. 4.

APPENDIX E: TERNARY PHASE DIAGRAMS
We model phase stability by exploring the full composition space. In addition to the quasicrystal phase, we include the pure species in their favored structures, all known binary Al-Cu and Al-Fe phases, and all known ternary phases: λ-Al 3 Fe.mC102, β-AlCuFe, Al 6 Mn structure type τ 1 -Al 23 CuFe 4 .oC28, ω-Al 7 Cu 2 Fe, φ-Al 10 Cu 10 Fe. All of the phases have known structures with the exception of β and φ, exhibiting vacancy and chemical disorder; the φ phase was claimed to be a superstructure of Al 3 Ni 2 structure [48] but belongs to the same B2-type family.
To represent the β phase family, we evaluated DFT cohesive energies for every member of the 2 × 2 × 2-supercell ensemble of the cubic B2 structure type, under the constraint of fixing the cube-vertex occupation as Al. For a given Cu content, we created a list of all symmetry-independent Cu/Fe orderings on the body-center sublattice. Ground states with 1-5 Cu atoms per 16-atom supercell revealed that 2-4 Cu atom range yields stable structures, while 1/16 or 5/16 Cu atom compositions (with 7/16 and 3/16 Fe atom content, respectively) are unstable. By examining the electronic DOS we concluded that the ternary β phase is electronically stabilized at low T by pushing E F beyond the steep Fe-d-band shoulder, optimally by substituting 3Fe →3Cu (per 16-atoms). This ternary β-phase is predicted to be stable at T = 0 K in the composition range 12.5-25% of Cu.
At the Cu-rich composition, the β phase takes a vacancyordered form with experimental composition Al 10 Cu 10 Fe. Since we could not find any promising T = 0 K structure in the 2 × 2 × 2 supercell, and larger supercells ensembles are inaccessible to our direct DFT method, we proceeded with EOPP potentials and fixed-site lattice-gas annealing [35], in which atoms are constrained to occupy fixed lattice of sites, but pairs of species with different chemistry are allowed to swap their positions. Using this method we discovered a T = 0 K stable state in a 3 × 3 × 3 supercell, whose composition can be described by a single parameter x = 4/(2×3 2 ) ∼0.074 and composition (Al 0.5−x Cu x )(Cu 0.5−2x Fe x Vac x ), where the parentheses separate cube-vertex/body-center sublattices, respectively.
The quasicrystal family of structures is represented by 2/1 and 5/3 approximants. The former turns out to be the most stable T = 0 K structure within the family (δE = +1.8 meV/atom), while the 5/3 model at (δE = +4.0 meV/atom is our best representation of the quasicrystal phase. To predict the phase diagram at finite temperature, T > 0, we add the free energy F Tot in Eq. (G1) to the DFT-calculated enthalpy for each structure considered. Because full phonon calculations for large QC approximants are prohibitive, we assume that they all share the same F h as the 2/1 approximant. We then compute the convex hull of the set of free energies (see Fig. 6). Vertices of the convex hull are predicted to be stable. Line segments and enclosed triangles are predicted tie-lines and tie-planes connecting coexisting phases. We also include a few structures whose free energies lie slightly above the convex hull by up to 4 meV/atom. Structures are labeled using their phase name followed by their Pearson symbol.
All known structures lie on or near the convex hull, both at T = 0 K and at T = 600 K, with the exception of Al 2 Fe in its observed structure of Pearson type aP19. Instead we find the hypothetical Al 2 Fe structure of Pearson type tI6 that is believed to be more energetically stable [49]. Some binaries extend into the ternary composition space, for example Al(Cu,Fe).cP2. The two quasicrystal approximants, iQC-2/1 and iQC-5/3 (reference numbers 21 and 20, respectively), swap stability between low and high temperature leading to changes in the vertices, edges, and facets of the convex hull as can be seen by comparing the two diagrams in Fig. 6. This occurs because the smaller 2/1 approximant has an optimized structure and composition at which a deep pseudogap appears and the energy reaches the convex hull, whereas the larger 5/3 approximant lacks a unique optimized structure and composition but instead enjoys a large anharmonic contribution to its entropy. The 2/1 approximant appears on the convex hull at T = 0 despite having E > 0 according to DFT because of the quantum zero point vibrational (or competing phases) that is contained in F h . The 5/3 approximant is the largest for which we can reliably obtain its low temperature enthalpy by DFT, and hence we take it as a proxy for the true quasicrystal state.
Thus we predict that the 2/1 approximant should be stable at low temperatures but transform to the icosahedral quasicrystal at elevated temperatures of 600 K and above. Correspondingly, the quasicrystal loses stability at low temperature and transforms into the 2/1 approximant. Many actual or implicit transformations have been reported for Al-Cu-Fe quasicrystals [50][51][52][53]. Fine detail of the phase diagram at 600 C [39] revealed that except for small window of single-phase stable quasicrystal composition around Al 62 Cu 25.5 Fe 12.5 , the quasicrystal phase transforms into a rhombohedral approximant (R-5/3 in our notation), and this transformation is reversible [54]. The quasicrystal phase field shrinks with decreasing temperature [43] but remains finite at T = 560 C, a temperature above our predicted transformation. At low temperatures the kinetics becomes slow and the transformation will be inhibited, so the quasicrystal remains metastable at low temperatures.
According to the experimental phase diagram, our energy minimizing, electronically optimized composition (Al 65 Cu 22.5 Fe 12.5 ) with E F in the center of the pseudogap lies in a coexistence region of three phases: λ, ω, and R, and the latter should have composition Al 63.5 Cu 24 Fe 12.5 in agreement with Ref. [54]. However, the center of the R-approximant phase field is around x Cu = 0.26 and x Fe = 0.12, not far from the composition of our low-temperature winner, the 2/1 approximant Al 63.3 Cu 25.8 Fe 10.9 , which lies just outside the R-phase stability range.
We simulated the R-phase structure in a cell of 718 atoms at our optimized composition of Al 65 Cu 22.5 Fe 12.5 . Although the optimal structure places E F at the center of a pseudogap, the energy remains higher than the cubic 5/3 approximant by ∼ 3 meV/atom. Due to kinetic barriers at low temperatures, existence of the 2/1 approximant cannot be ruled out. Adequate evaluation of all competing phases around 800-1000 K would require systematic variation of composition and density of all competing phases.

APPENDIX F: ENERGETIC OPTIMIZATION OF QUASICRYSTAL APPROXIMANTS
The most direct comparison between approximants is achieved under the constraint of equal density/composition revealing the impact of structure alone. Also, these are the conditions under which the EOPP energies are most meaningful (see Appendix B). Expressing atomic density per b 3 volume, where b ∼ 12.2 Å is the side of the 2/1 cubic cell, we chose a density of 129.51 atoms/b 3 and composition Al 65 Cu 22.5 Fe 12.5 . To satisfy this constraint accurately for all approximants, we worked with the 2/1 approximant in a 2 × 2 × 2 supercell (1032 atoms), the 3/2 approximant in a √ 2 × √ 2 × 1 supercell with 1096 atoms, and the 5/3 in its unit cell with 2324 atoms. Supercells were also required to counter the size effect when measuring anharmonic heat capacity, F a . The resulting optimized energies are presented in Table IV. The sequence of formation enthalpies, both for EOPP and full DFT calculations, appears to favor larger approximants that minimize the phason strain and accommodate larger superclusters.
Since the optimal density and composition could vary between approximants, we varied these individually for each approximant within its conventional unit cell, with the results presented previously in Table I. Again the enthalpy is found to be a decreasing function of approximant size, both for EOPP and for DFT. However, the energy relative to the convex hull is minimized for the smallest approximant. This is because the larger approximants favor greater Fe content, and the strong bonding of Fe (see Fig. 2) causes a strong slope of the convex hull facets in the direction of increasing Fe.
The 2/1 approximant system size (128-129 atoms) allowed for full exploration of all degrees of freedom. Under EOPP the icosahedral structure forms easily from the melt, and we can apply full DFT refinement to simultaneously explore compositional and density variation for EOPP preoptimized models. The structure with lowest E occured at increased Cu content, and density of 128 atoms/cell (identical with the Katz-Gratias model prediction). Starting from this model, we then performed AIMD for 5000 steps (5-fs time step) of MD annealing at 1100, 900, and 700 K, and quenched several snapshots from each temperature. The lowest energy snapshot was from the 900 K annealing batch and yielded the FIG. 9. Electronic density of states for 2/1, 3/2, and 5/3 approximants at fixed composition Al 65.0 Cu 22.5 Fe 12.5 . 2/1 and 3/2 approximants represented by supercells (see Table IV) contain amount of frozen disorder comparable to the 5/3 approximant. Resolution with Gaussian σ = 0.02 eV (a), 0.01 eV (b), and 0.006 eV (c). best atomic structure. This optimal structure exhibits a deep pseudogap (0.015 states/eV/atom according to tetrahedron method calculation; see Fig. 10) centered on the Fermi energy. The pseudogap becomes shallower and broader for structures in the equlibrium ensemble at higher temperatures.
In the 3/2 approximant cell, we explored several densities: 544 (KG model density), 552, and 448 atoms/unit cell. The composition was refined by seeking deepening of the pseudogap and controlling the Fermi energy assuming a rigid-band picture and simple valence rules. The best structure was a 552-atom model, greater by 1.5% than the KG-model density. This structure was then annealed under AIMD for 5000 steps at 1250 K; we observed strong atomic diffusion and some Al atoms moved as far as 6 Å. Subsequently, we cooled gradually from 1250 K to 800 K in another 5000 steps, and finally from 800 K down to 300 K in 1000 steps. At the end, we found that maximal displacement was 5.1 Å for Al, 2.9 Å for Cu, and 0.5 Å for Fe atoms; 20 Cu and 100 Al atoms displaced by more than 0.5 Å. Despite that, the energy of the final annealed configuration was nearly identical to the initial configuration. Our best 3/2 approximant has a narrow true gap (according to the tetrahedron method) at E F . The 5/3 approximant system size did not allow for systematic variations, but we did explore several densities (2304, 2324, 2338 atoms) and compositions, using pseudogap depth and Fermi level position as guides. The final structure was optimized under ab initio relaxation in ∼70 ionic steps. AIMD annealing was not feasible. In contrast to the smaller approximants, all 5/3 models considered had broader and less deep pseudogaps. Nonetheless, the 5/3 approximant achieves the lowest formation enthalpy of the approximants, most likely as a result of its greatest Fe content (as seen in Fig. 2 the two   FIG. 10. Electronic density of states for 2/1, 3/2, and 5/3 approximants at their optimal compositions (see Table I), eigenenergies density smeared using Gaussian σ = 0.02 eV (a). Zoom into fineaccuracy tetrahedron method DOS calculation using 10×10 × 10 and 6 × 6 × 6 meshes, respectively, for 2/1 and 3/2 approximants in (b) and (c). 5/3 approximant (2×2 × 2 calculation) with resolutions σ = 0.006 eV (b) and 0.002 eV (c). strongest bonds are near-neighbor Al-Fe and next-nearestneighbor Fe-Fe).

APPENDIX G: THERMODYNAMICS
The Helmholtz free energy F (N, V, T ) can be approximately decomposed into a relaxed T = 0 K energy E 0 , plus corrections due to harmonic vibrational free energy, F h , anharmonic positional disorder F a , and electronic excitations F elect . Thus, we write The harmonic vibrational free energy of a single phonon mode of frequency ω is Notice that as k B T hω, f h (ω) →hω/2, which is the zero point vibrational energy. As k B T hω, f h (ω) → k B T ln (hω/k B T ), which is the classical limit of the free energy. The full harmonic free energy The anharmonic contribution includes corrections due to shifts in phonon frequency with large amplitude of oscillation and additional discrete degrees of freedom connected to chemical substitution and possible additional tiling flips. At low temperature these contributions can be neglected, so we only include them beyond a temperature, T 0 , which we set at 200 K. At these elevated temperatures positional degrees of freedom behave nearly classically, so we will evaluate F a from our classical MC/MD simulations using EOPP.
In the canonical NV T ensemble, the Helmholtz free energy F (N, V, T ) = U − T S has the differential dF (N, V, T ) = −SdT − pdV + μdN. (G4) In particular, S = −∂F/∂T . The entropy is also related to the heat capacity C = ∂U/∂T through C/T = ∂S/∂T . Hence, Conveniently, C can be obtained at high temperatures from classical MC/MD simulation through fluctuations of the energy, (G6) C, thus obtained, includes contributions from both harmonic and anharmonic atomic vibrations, and potentially also from chemical substitution and tile flipping. Because we seek the anharmonic component of the free energy, we define C a ≡ C − 3Nk B for T > T 0 , and C a ≡ 0 for T < T 0 . We now integrate C a once to obtain and then integrate once more to obtain Summing over electronic states, we obtain free energy (G10) The overlapping energy distributions at different temperatures created by replica exchange provide an opportunity for accurate calculation of heat capacity, entropy, and free energy through the method of histogram analysis. At a single temperature T , the frequency distribution of simulated energies E is proportional to the Boltzmann probability This equation can be inverted to obtain the density of states (E ) ∼ H T (E )e +E /k B T , where H T (E ) is a normalized histogram of energies obtained from a simulation at fixed temperature T . Given the density of states, the partition function may be calculated (up to an undetermined multiplicative factor) by integrating The free energy is determined (up to an additive linear function of T) from and all other thermodynamic functions can be obtained by differentiation. Notice that Z (T ) and F (T ) are obtained as continuous functions of temperature T over a range of temperatures surrounding the original simulation temperature. The same approach interpolates between the fixed simulation temperatures by consistently merging densities of states obtained from each temperature [55]. Up to an unknown multiplicative constant, we have where the free energies F T must be obtained self-consistently with (E ) from Eqs. (G14) and (G13). By setting the value of F and its derivative S = −∂F/∂T to values determined from first principles methods at the lowest simulation temperature, we obtain absolute free energy across the entire simulated temperature range.