Multiple-scale structures: from Faraday waves to soft-matter quasicrystals

Models describing the thermodynamic stability of soft-matter quasicrystals are reviewed and expanded. New analytical methods for treating them are presented, and a number of new stable quasicrystalline structures are reported.


Introduction and outline
The scope of research on quasicrystals 1 has greatly expanded in the last decade, owing mainly to the advent of a host of new experimental systems exhibiting aperiodic structures with long-range order. The ever-growing variety of stable solidstate quasicrystals (Tsai, 2003(Tsai, , 2008Janssen et al., 2007;Steurer & Deloudi, 2009), where quasiperiodic long-range order occurs on the atomic scale, has been joined in recent years by a host of exciting new soft-matter systems that exhibit this quasiperiodicity on a larger mesoscopic scale -typically from a few nanometres to a few micrometres (Zeng et al., 2004;Ungar & Zeng, 2005;Takano et al., 2005;Hayashida et al., 2007;Percec et al., 2009;Talapin et al., 2009;Ungar et al., 2011;Dotera, 2011Dotera, , 2012Fischer et al., 2011;Xiao et al., 2012;Zhang & Bates, 2012;Bodnarchuk et al., 2013;Chanpuriya et al., 2016). In addition to having promising applications, particularly as metamaterials in the optical domain, these substances give us the opportunity to study quasicrystals in ways that were impossible before. The obvious reason for this is indeed the fact that their building blocks -rather than being individual atoms -are composed of large synthesized particles such as macromolecules, nanoparticles and colloids. At these dimensions it may be possible to track the dynamics of individual particles, manipulate their positions or possibly design the interaction between them. If so, an obvious question to ask is how to design this interaction to obtain a particular desired quasicrystal. To answer this question, one clearly requires an understanding of the spontaneous formation and subsequent stability of such materials.
Phenomenological Landau theories based on ad hoc free energies have been widely applied to study the thermodynamics of phase transitions (Alexander & McTague, 1978) and to explain the stability of different phases, including quasicrystals (Mermin & Troian, 1985;Bak, 1985;Kalugin et al., 1985;Jaric, 1985;Gronlund & Mermin, 1988;Narasimhan & Ho, 1988). In such models, one identifies an order parameter field -which is often a simple scalar function (r) that describes the relative deviation ½cðrÞ À c=c of a coarse-grained density c(r) of a material from its average value c -and uses generic symmetry arguments to formulate a free-energy functional F ½ðrÞ, expressed in powers of the order parameter and its gradients. One assumes that the equilibrium phase is the one that minimizes F, and then all that remains is to find out which structures could minimize such a free energy, or conversely, how to tweak F to obtain the desired structures. For a textbook introduction see, for example, chapter 4 of the book by Chaikin & Lubensky (1995). These models are particularly attractive for soft-matter systems (de Gennes & Prost, 1993;Gompper & Schick, 1994) owing to their mesoscopic building blocks, which render the long-wavelength gradient expansion and the truncation at low order a more valid approximation than in the atomic case, especially when the transition is only weakly first-order.
Important insight into the stability question emerged when Edwards & Fauve (1993) discovered that parametrically driven liquid surfaces, exhibiting standing-wave patterns known as Faraday waves, can become quasiperiodic when driven by a superposition of two temporal frequencies [see also Kudrolli et al. (1998), Gollub & Langer (1999) and Arbell & Fineberg (2002)]. The realization that these temporal frequencies impose two spatial length scales on the stable structures that form prompted Lifshitz & Petrich (1997), henceforth LP, to generalize the Swift-Hohenberg equation (Swift & Hohenberg, 1977) and introduce a rather simple Landau free-energy expansion (or a Lyapunov functional) of a scalar field (r) in two dimensions of the form In Section 2 we carefully review the basic features of this free energy. Here, we only wish to point out its two main features: (i) The first integral in F LP , containing the non-local gradient expansion, is responsible for selecting two length scales, whose ratio is given by the parameter q. This is because density modes with wavenumbers differing from unity or q increase its value.
(ii) The second integral, containing the local expansion in powers of the order parameter field, contains an odd, cubic, power that breaks the Z 2 symmetry of ! À. It is this term that is responsible for stabilizing structures with two length scales, depending on the value of q, as it has the ability to lower the free energy if there exist triplets of density modes with wavevectors that add up to zero. These are known in the Faraday wave literature as triad resonances, and amount to effective three-body interaction in the coarse-grained density context.
In particular, by setting the value of the wavenumber ratio q to one can form triplets or 'triangles' containing two unit wavevectors separated by 2/n and a third wavevector of length k n . This may sufficiently lower the free energy of structures with N-fold rotational symmetry (where N is equal to n or 2n for even or odd n, respectively), making them the absolute minimum of F LP . In Section 3 we repeat and extend the calculations of LP of the free energies of candidate structures, setting q = k n for different values of n and assuming LP's limit of c ! 1, which leads to exact length-scale selectivity. In doing so we provide a more complete and definitive calculation, while correcting a few discrepancies in their results that have caused some confusion over the years. Owing to its simplicity and clarity in explaining the stability of the decagonal (10-fold) and dodecagonal (12-fold) quasicrystals that it exhibits, as well as the ease with which one can numerically simulate the dynamical equation that it generates via simple relaxation @ t = ÀF LP =, the LP model has been studied in depth since its original publication and extended in a number of different ways (Wu et al., 2010;Mkhonta et al., 2013;Achim et al., 2014;Jiang & Zhang, 2014;Jiang et al., 2015Jiang et al., , 2016Jiang et al., , 2017Subramanian et al., 2016).
Here we further extend the LP model as follows: (i) Jiang & Zhang (2014) and Jiang et al. (2015) improved on the free-energy calculations of LP by relaxing the exact length selection imposed by the c ! 1 limit, using a highdimensional numerical evaluation scheme which they call the 'projection method'. In Section 4 we introduce an approximation scheme for calculating the LP free energy (1) with finite c that allows competing structures to contain two length scales that are roughly, rather than exactly, equal to unity or q, thus improving their competitiveness and reducing the regions in parameter space where the quasicrystalline structures are stable. This qualitatively captures the importance of lengthscale selectivity, but is quantitatively accurate only in the dodecagonal case. Nevertheless, it is much simpler to evaluate and provides further analytical insight about the model.
(ii) The only quasicrystals which can be stabilized by the original LP model, which allows for two length scales in the structures, are the decagonal and dodecagonal phases. In Section 5, we show that increasing the number of allowed length scales from two to four allows for the stabilization of octagonal (8-fold) and octadecagonal (18-fold) quasicrystals.
One can improve on the Landau expansions by using a density functional mean-field theory (Ramakrishnan & Yussouff, 1979), which is valid for all orders, by rigorously coarse-graining a system of interacting discrete particles. For a textbook introduction see, for example, chapter 5 of the book by Fredrickson (2006). Such theories were also considered in the early studies of quasicrystals (Sachdev & Nelson, 1985). A particularly simple coarse-grained free-energy functional of the form containing the familiar mean-field terms of pair interaction and ideal entropy, was used by Barkan, Diamant & Lifshitz (2011), henceforth BDL, as an extension of the LP model, to study the stability of soft-matter quasicrystals, initially in two dimensions. Again, one assumes that the equilibrium density field is the one that minimizes F CG for the given thermodynamic parameters -such as temperature T and chemical potential -and is then left with the question of how to design the pair potential U ðrÞ to obtain the desired structures, giving us the ability to address our starting question.
To do so, BDL followed an earlier conjecture of Lifshitz & Diamant (2007), who attributed the stability of certain soft quasicrystals to the same mechanism that stabilizes the Faraday wave structures, namely the existence of two length scales in the pair potential, combined with effective manybody interactions. That stable quasicrystals may require the existence of two length scales in their effective interaction potentials U ðrÞ is not a new idea (Olami, 1990;Smith, 1991). Many two-length-scale potentials have been investigated numerically over the years and found to exhibit stable quasiperiodic phases (Dzugutov, 1993;Jagla, 1998;Skibinsky et al., 1999;Quandt & Teter, 1999;Roth & Denton, 2000;Engel & Trebin, 2007;Keys & Glotzer, 2007;Archer et al., 2013;Dotera et al., 2014;Engel et al., 2015;Pattabhiraman & Dijkstra, 2017a;Damasceno et al., 2017). The novelty and emphasis of BDL was in their quantitative understanding of the stabilization mechanism -comparing the nonlocal pairwise interaction and the local entropy terms of F CG to the nonlocal gradient expansion and local power expansion terms of F LP , respectively. This allowed them to pinpoint regions of stability in the parameter spaces of different potentials instead of performing exhaustive searches. Indeed, Barkan et al. (2014) confirmed these predictions by employing molecular dynamics simulations with particles that interact through pair potentials that were designed according to the principles of BDL. By properly setting the two length scales in these potentials, they were able to generate periodic crystals with square and hexagonal symmetry, quasicrystals with decagonal and dodecagonal symmetry, and a lamellar (or striped) phase.
The inclusion of a second length scale in U ðrÞ, imitating the gradient term of F LP , provides greater control over the selfassembly of desired structures than can be achieved with just a single scale, and turns out to be the key to obtaining stable quasicrystals and other novel structures. Yet, calculating the exact value of the coarse-grained free energy F CG turned out to be a challenge. Instead, BDL expanded the logarithmic entropy term in a power series to fourth order in (r) = ½cðrÞ À c=c and mapped the resulting approximate free energy onto the LP free energy, thus obtaining a rough estimate of the physical parameters that stabilize the different structures using the results known for the LP model.
Here we present new insight into the stability of soft-matter quasicrystals, by significantly improving upon the original BDL analysis as follows: (i) In Section 6, we introduce the 'density distribution method' for evaluating the free energy of candidate structures with non-polynomial local free-energy terms.
(ii) Section 7 applies this technique to the coarse-grained free energy F CG and uses it to point out the differences in the stabilities of different structures between the BDL and LP models, and in particular to explain the previously surprising robustness of decagonal structures in the BDL model.
(iii) Finally, equipped with this new understanding of the effect of the local free-energy term, we again use the density distribution method in Section 8 to generate an artificial local free-energy term that can stabilize quasicrystals with 6n-fold symmetry, with arbitrarily large n, using only two length scales.
Note that the vast majority of the stable two-dimensional quasicrystals that have been discovered to date have symmetry orders no greater than 12-fold. Possible explanations for this have been suggested by Levitov (1988) and Mikhael et al. (2010). Exceptions are the octadecagonal quasicrystal discovered by Fischer et al. (2011), those found numerically (Dotera et al., 2014;Engel & Glotzer, 2014;Pattabhiraman & Dijkstra, 2017b) and the one discussed in Section 5, as well as the numerically discovered icositetragonal (24-fold) quasicrystals (Dotera et al., 2014;Engel & Glotzer, 2014).
For completeness, we should note three additional extensions of the LP and BDL models that we do not discuss here. First, in the present work we focus solely on the question of thermodynamic stability (or metastability), searching for the minimum free-energy states, without considering any actual dynamics. LP used purely relaxational dynamics, also known as Model A of Hohenberg & Halperin (1977) starting with random initial conditions to confirm numerically that the steady states were indeed the targeted ones. This also established that quasicrystals were not as difficult to obtain as solutions of simple partial differential equations as one may have thought. Recent authors have been using conserved dynamics of the form @ t = Dr 2 ðF =Þ, also known as Model B of Hohenberg & Halperin (1977), or slight variations of Model B, known as dynamic density functional theory (DDFT) or phase-field crystal (PFC) models (Archer et al., 2013;Achim et al., 2014Achim et al., , 2015Subramanian et al., 2016). 2 Interestingly, the PFC model of Achim et al. (2014) uses F LP without the cubic term, but still genrates the same structures. At first sight, this seems to contradict LP's explanation of the stability of their quasicrystals. However, it turns out that this apparent restoration of the Z 2 symmetry ! À in the local term of the free energy is destroyed by the conservation of mass condition, which constrains the average density to be a positive constant. This, in turn, generates an effective cubic term in the local free energy, with eff / (Barkan, 2015).
Extensions of the LP and BDL models, which we intend to pursue elsewhere, include: (i) The application of these models in three dimensions. This has already been shown by some authors to produce stable icosahedral quasicrystals using two length scales (Subramanian et al., 2016;Jiang et al., 2017). 3 (ii) The generalization to two (or more) interacting densities or order-parameter fields. The use of two coupled fields or two coupled Swift-Hohenberg equations, where each field carries one of the length scales, was considered very early on (Mermin & Troian, 1985;Sachdev & Nelson, 1985;Narasimhan & Ho, 1988;Mü ller, 1994) and has been resumed recently in the context of binary and ternary soft-matter systems (Dotera, 2007;Barkan, 2015;Jiang et al., 2016), with new insight gained from results of the LP and BDL models.
Soft-matter quasicrystals provide rich and versatile platforms for the realization of relatively simple theoretical models as classical particles interacting via pre-designed pair potentials, treated either by molecular dynamics simulations or by coarse-grained mean-field theories and their Landau expansions. Such theoretical tools may be inadequate for treating atomic-scale quasicrystals, yet perfect for the fundamental study of the basic notions of the physics of quasicrystals as they appear in soft condensed matter. Armed with renewed insight from soft-matter systems and the potential to realize them directly in the laboratory, some of the outstanding fundamental questions in the field can be treated afresh, allowing one to get closer than ever to their resolution. Admittedly, our discussion here may apply only to soft condensed matter, but intriguing new analogies between softmatter and solid-state systems continue to emerge (Lee et al., 2014;Lifshitz, 2014), possibly enlarging our scope.

The Lifshitz-Petrich model and its immediate generalizations
Following the original analysis by Lifshitz & Petrich (1997), we define a scalar field (r) on the two-dimensional Cartesian plane. The Swift-Hohenberg free energy of this field (Swift & Hohenberg, 1977) can be written as where f ðÞ is a local contribution to the free energy, which may or may not be symmetric under the operation that replaces by À (Cross & Greenside, 2009). In what follows, we use (r) to refer to the field and for specific scalar values the field can take on at a given point. The LP free energy (1) changes this to and sets explicitly breaking the ! À symmetry, where c is assumed positive and q is the ratio of the two selected length scales, which generally satisfies 1 < q 2 . Note that the coefficient of the cubic term in equation (1) has been scaled to unity by measuring the field amplitude in units of and consequently measuring energy in units of 4 . The parameter c, which sets the length-scale selectivity of the system, and the control parameter are then measured in units of 2 . By substituting the Fourier transform into the first terms of free energies like the ones in equations (5) and (6), they can be written as  V V LP ðkÞ for c = 1 and q = k 5 and k 12 : note thatṼ V LP ðkÞ is positive for all k except k = 1 or q, where it is zero. Also note that the barrier between the minima when q = k 12 is significantly larger than the barrier when q = k 5 . This disparity is important for the study of the finite-c case discussed in Section 4.
which is sketched in Fig. 1. 4 The (horizontally stretched) local quartic free-energy density (7) is plotted as the solid colored lines in Fig. 2 for different values of its single parameter . The reader should note that all the position-and momentum-space integrals are implicitly normalized to give a free energy per unit area. After rescaling, the LP model is left with only three free parameters, q, c and . Given specific values for these, we seek the configurations (r) that minimize F LP . This is easiest in the limit where c is taken to infinity. Because this infinite-c limit is also generally favorable for the formation of quasicrystals, we adopt it throughout this paper, with the exception of Section 4. In this limit,Ṽ V LP ðkÞ = 0 if k belongs to the setṼ V 0 = {1, q} and is otherwise infinite. Thus, we immediately conclude that lim c!1 restricting the support of ðkÞ to lie entirely on two concentric circles of radii unity and q, centered about the origin. Given this restriction, the free energy is simply Upon substituting the Fourier transform (8) of the fieldwhich for quasiperiodic density fields is supported on a countable set of wavevectors, changing the integral into a sum -this becomes where the tilde indicates the restriction that the magnitude of all the wavevectors and their sum must belong to the setṼ V 0 . The products of Fourier coefficients on wavevectors that add up to zero, appearing in this expression for the free energy, are known in crystallography as 'structure invariants'. The sums can easily be evaluated on a computer using symbolic algebra. A wise choice of q can make use of the triplets, or wavevector triangles, on the second line for stabilizing the desired twoscale structures, as mentioned earlier. The benefit of adding such triplets usually comes at the cost of more quadruplets on the third line that generally increase the free energy. In Section 3, essentially by counting triplets and quadruplets, we repeat the calculation of Lifshitz & Petrich (1997) and show that this simple free energy is able to stabilize periodic square and hexagonal crystals, decagonal and dodecagonal quasicrystals, and lamellae, also called stripes. In addition to the obvious generalization to three dimensions, the free energy in equation (9) can immediately be generalized by modifyingṼ VðkÞ, f ðÞ or both, in a number of ways: (i) While remaining in the infinite-c limit,Ṽ V LP ðkÞ can be changed by modifying the setṼ V 0 . We show in Section 5 that doubling the cardinality ofṼ V 0 , i.e. going from two to four concentric circles, allows us to stabilize octagonal and octadecagonal quasicrystals.
(ii) Subramanian et al. (2016) modifyṼ V LP ðkÞ in a particular way in order to gain control of the relative heights of the two minima (see Fig. 1), while leaving c finite.
(iii)Ṽ VðkÞ, as indicated by its notation, can be thought of as the radial Fourier transform (also known as the Hankel transform) of an isotropic interaction potential UðrÞ of a pair of particles in real space, possibly scaled and shifted by a constant. This, along with a replacement of f ðÞ by the local entropy term from F CG in equation (3), forms the basis of the BDL model, which we consider and expand our understanding of in Section 7. A comparison of these two choices for f ðÞ is shown in Fig. 2.
(iv) In Section 8 we show that an artificial yet judicious choice of f ðÞ can actually stabilize 6n-fold quasicrystals for any n ! 2, with just two length scales. Local free-energy densities f ðÞ used in this work: solid colored lines show the scaled ( = 1) quartic local free-energy function (7) used by LP, for values of above, equal to and below the spinodal value of = 0. The full logarithmic local free-energy function f CG in equation (49), as used by BDL, is shown with dashed lines for temperatures above, below and equal to the spinodal temperature, which is scaled to T = 1. The values of and T are related by T = 4=ð3 þ 4Þ, as explained in Section 7. In order to demonstrate visually the resulting fourth-order agreement, the solid lines have been stretched horizontally by 50% while the dashed lines have been compressed vertically by a factor of 27T/16. Note that the LP quartic free-energy density penetrates into the < À 1 region and diverges from the BDL logarithmic free-energy density for > 1. Finally, the solid black line shows the free-energy density used in Section 8 to stabilize 6n-fold quasicrystals with arbitrary n.
These two-scale structures are in thermodynamic competition with the uniform liquid phase (r) = 0, and with singlescale and trivial two-scale periodic structures consisting of two degenerate lamellar phases varying in their spatial scale (set by which circle the two peaks lie on 5 ), four degenerate hexagonal configurations, two of which are regular and two distorted, containing two length scales, and infinitely many oblique, rectangular and square structures consisting of a sum of two cosines with an arbitrary relative orientation, whose wavevectors are taken from the set {1, q}. These competing structures are shown schematically in Figs. 3(a)-3(h). The targeted structures and the competing ones are also listed in Table 1. As is typical for these kinds of stability calculations, one cannot be certain that the list of candidate and competing structures is exhaustive. We believe it is complete, however, not only due to intuitive symmetry considerations, but also based on repeated computational simulations. As one example, the Model A dynamics of equation (4) have been applied to many realizations of random initial conditions, thereby exploring the space of likely minimum free-energy states over a wide range of model parameters.
As noted by LP, because all the candidate structures are centrosymmetric, and because there are no screw rotations in two dimensions, we may always take each of the Fourier coefficients on a given circle to be equal, and their phases may all be chosen such that they are either 0 or , corresponding to positive and negative real values, respectively. 6 The minimization of the free energy (13) is therefore always with respect to no more than two real variables, which we denote in structures with a single scale and 1 and q in the two-scale structures. For example, the larger regular hexagonal phase has a real-space structure of The free energy of equation (13) then becomes a quartic function of two variables, given by for two-scale structures, where all k i = j k i j2 f1; qg, and a similar function of the single variable for single-scale structures. Computer-assisted symbolic algebra is used to evaluate the sums, and then to minimize them with respect tõ or 1 and q . The structure that has the lowest free energy, given the value of , is the thermodynamically stable one, assuming we have not overlooked any additional competing structures.
3.1.2. Calculation of metastability bounds. The majority of work on the LP model has been focused on finding the freeenergy minimizing structure for various choices of the parameters. However, as seen in the work by Barkan et al. (2014), the phase transitions between these structures exhibit significant hysteresis. As a first attempt at evaluating the meta-  Fourier spectra of the candidate structures: structures (a)-(g) use the arbitrary ratio q = k 5 , while all other ratios are specified explicitly. The uniform phase (a) has no Fourier modes. An example of a stabilizing triangle is included in the k 5 decagonal structure (l). stability bounds on , we can imagine writing a structure as a linear combination of multiple components such as aligned lamellar, hexagonal and dodecagonal modes on a circle. The coefficients A i are set by the local minimum of the free energy, where, for a given phase, some of the A i 's will be zero. These represent the potential 'directions' in which the structure can decay. The spinodal decomposition of a phase, where it is no longer metastable, occurs when that point on the free-energy landscape transitions from a local minimum to a saddle point. This occurs when the determinant of the Hessian of the free energy of this structure, becomes zero. In our calculations, we include all of the candidates as potential decay directions, but we have no proof that these are the only options, so the reader should take the metastability bounds reported below as tentative results.

Stability and metastability bounds in the original LP model
3.2.1. Single-scale phases. First, we consider the single-scale lamellar, hexagonal and square phases, along with the uniform liquid state. Recall, also, that this includes the two distorted hexagonal phases that have the same free energy as the two regular ones, even though, strictly speaking, they consist of two length scales. The Fourier spectra of these structures are depicted in Figs. 3(a)-3(e) and 3(h). The stability bounds are summarized in the top section of Table 1 and plotted in Fig. 4.
The uniform phase always has a free energy of and decays spinodally when > 0. For the lamellar phase, the free-energy equation (15) gives Minimizing F LAM over shows that, for < 0, = 0, thus giving a uniform phase. The lamellar phase therefore only exists for positive , wherein The free energy of the hexagonal phase is given by It only exists for ! À1=15, below which the nontrivial minima of equation (22) Table 1 Stability and metastability boundaries in the infinite-c LP model.
Note that the stability regions of the single-scale structures are modified by their competition with the q-tuned two-scale structures. The metastability regions are left unchanged, except in the k 6 and k 12 cases where the lower bound becomes 0. Recall that here is in units of 2 , as in equation (7); without scaling one should multiply the quoted values by 2 . The choices of q = k 8 and k 10 fail to stabilize octagonal and decagonal quasicrystals, respectively. However, they can exhibit regions of metastability, which might be observable in experiments or simulations given proper initial conditions.

Figure 4
Free energies of the single-scale periodic structures in the infinite-c LP model: as is increased, the uniform liquid is the equilibrium phase until it reaches À8/135, at which point a first-order transition to the hexagonal structure occurs. This persists until reaches (6 3/2 + 14)/15 where the lamellar structure becomes the equilibrium phase. The gray region corresponds to the forbidden zone below the lower bound, calculated in Section 6.6.
state is the uniform one with = 0. Above this point -a saddle node in the context of dynamical bifurcation theory -we have the nontrivial minimum For generic q, as is increased, the system first undergoes a first-order transition from the uniform phase to the hexagonal phase at = À8/135 ' À0.05926 and then to the lamellar phase at = (6 3/2 + 14)/15 ' 1.913. At the second transition, F = À[84(6 1/2 ) + 206]/675. This behavior is shown in Fig. 4. The hexagonal phase is metastable when À0.06667 ' À1/15 16/3 ' 5.333, and the lamellar phase is metastable for all ! 4/3 ' 1.333.
The single-scale square structure in Fig. 3(h) and its infinitely many degenerate oblique, rectangular and square structures, consisting of the sum of two cosines with an arbitrary relative orientation, all have a free energy of which leads to a minimized free energy of Because the square structure has additional quadruplets compared with the lamellar phase (21) without any compensating triplets, its free energy is always higher, and it is therefore never in thermodynamic equilibrium. 3.2.2. Two-scale periodic phases. Next, we consider the two-scale square, hexagonal and striped superstructures, for q = k 4 , k 6 and k 1 , respectively. Their Fourier spectra are shown in Figs. 3(i)-3(k). 7 For these structures there are no simple expressions for the minimized 1 and q values, which are generally unequal, nor for their minimized free energies and critical values. Thus, we provide only numerical results for the stability bounds in the middle section of Table 1 and plot these bounds in Fig. 5.
The free energies from which these bounds are obtained are given by where F SQU ð q ; Þ, F LAM ð q ; Þ and F HEX ð q ; Þ are given in equations (24), (22) and (20), respectively. 3.2.3. Two-scale quasiperiodic phases. Finally, we consider the two-scale quasicrystals with q = k 5 , k 12 , k 8 and k 10 , whose Fourier spectra are shown in Figs. 3(l)-3(o), respectively. We find that the k 8 and k 10 structures are never global minima of the free energy and are therefore unstable. We do not give the detailed calculation of their free energies here, and only note that they may exhibit regions of metastability. Thus, one could potentially observe them in experiment or simulation given proper initial conditions. The stability bounds for the k 5 decagonal quasicrystal and the k 12 dodecagonal quasicrystal are included in the latter half of Table 1 and plotted in Fig. 6. The stability bounds reported here should be taken in place of the original bounds reported by Lifshitz & Petrich (1997), as they missed the existence of a stable decagonal quasicrystal with q = k 5 rather than k 10 and miscalculated the stability boundaries of the dodecagonal one. 8 The simplest expressions for the upper stability bounds, for both of these phases, involve roots of quintic polynomials that are provided below for the first time. The metastability bounds for all of the structures studied in this section are reported here for the first time as well. Because all of the results are summarized in Table 1, readers who are not interested in the detailed calculation itself may skip to the next section.
The free energy of the k 5 decagonal phase is given by Free energies of the two-scale periodic structures in the infinite-c LP model: in panel (a), with q = k 4 , the square superstructure dominates when À0.09205 < $ < $ 0.7689. In panel (b), with q = k 6 , the hexagonal superstructure dominates when À0.1143 < $ < $ 2.074. Finally, in panel (c), with q = k 1 , the lamellar superstructure dominates when > $ À0.06904. The two-scale superstructure always has a lower free energy than its single-scale analogue. The gray region is the same as in Fig. 4.
It exists only when ! À3/31 ' À0.09677. For À3/31 < $ 0.5245, 1 = q and Above this approximate upper bound, for which there is no simple expression, the free energy continues to decrease, but 1 no longer equals q . The free energy at the transition is approximately À0.04889.
It is interesting to note that the free energies of both the decagonal and dodecagonal quasicrystals, in the infinite-c limit, have an additional, accidental, Z 2 symmetry associated with the exchange of 1 and q . In both cases, when minimizing the free energy with respect to these amplitudes there is one solution branch that maintains this symmetry with 1 = q and a second branch where the Z 2 symmetry is sponta-neously broken and 1 6 ¼ q . Yet, as it turns out, in the decagonal quasicrystal this symmetry-breaking transition is first order, while in the dodecagonal case it is a continuous phase transition. We compare the free energies of these quasicrystal phases, given the requisite value of q, with those of the uniform and hexagonal phases in Fig. 6, which shows that the decagonal phase is stable for À0.08602 ' À8/93 < $ 0.02290. The upper bound is given by the second real root of the quintic 3 5 Â 25 Â 31 2 Â 43 4 5 À 16 Â 3 5 Â 5 Â 31 Â 43 Â 103 703 4 À 64 Â 81 Â 4 938 418 073 3 À 2 10 Â 27 Â 5 Â 11 798 281 2 À 2 9 Â 3 Â 5 Â 1 046 081 þ 2 12 Â 5 Â 67 303 ¼ 0; ð31aÞ and the free energy at the transition is approximately À3.694 Â 10 À3 . The dodecagonal phase is stable for À0.1009 ' À128/1269 < $ 0.03055. This upper bound is given by the second real root of and the free energy at the transition is approximately À4.180 Â 10 À3 . The decagonal phase with 1 = q is metastable when À0.09677 ' À3/31 763/972 ' 0.7850. The decagonal phase with 1 6 ¼ q appears to be metastable for all > $ À0.01493.
The dodecagonal phase with 1 = q is metastable for À0.1135 ' À16/141 208/867 ' 0.2399. The phase with 1 6 ¼ q appears to be metastable for all ! 208/867. Additionally, when q = k 12 , the lower bound of the hexagonal metastability region is increased to zero.

Relaxing the requirement of exact length-scale selection 4.1. The two-ring approximation
Despite the convenience of taking the infinite-c limit in the analysis of F LP , as given by equations (6) and (7), realistic systems can never fully extinguish all unwanted Fourier modes. It is therefore important to examine the LP model with finite-c values. Evaluating the finite-c LP free energy with quantitative precision requires an approach like the projection method of Jiang & Zhang (2014), which has been applied to the LP model (Jiang et al., 2015). However, one can obtain a fair understanding of the role of length-scale selectivity by employing a simple 'two-ring' approximation.
We restrict ðkÞ to lie within two rings centered about the origin. This is in contrast with allowing ðrÞ to vary freely, as most numerical simulations do (Lifshitz & Petrich, 1997;Barkan et al., 2011), or restricting ðkÞ to some subset of a twodimensional or higher-dimensional lattice, as in the projection method (Jiang et al., 2015). This two-ring approximation compromises the numerical accuracy of our results, but what Free energies of the quasiperiodic structures in the infinite-c LP model: in panel (a), with q = k 5 , the decagonal structure dominates when À8/93 < $ 0.02290. In panel (b), with q = k 12 , the dodecagonal structure dominates when À128/1269 < $ 0.03055. The gray region is the same as in Fig. 4. we lose in quantitative correctness we make up for in the simplicity with which we demonstrate the qualitative importance of length-scale selectivity in stabilizing quasicrystals via two preferred scales.
With this approximation, the target decagonal or dodecagonal quasicrystals have their Fourier amplitudes on exact circles of radii unity and q as before, and so their free energies are unchanged. The competing lamellar and hexagonal phases are rescaled by a factor so as to position both their first and second harmonics to fit, as well as possible, within two finitewidth rings near the minima ofṼ V LP ðkÞ from equation (10), which is plotted with c = 1 in Fig. 1. This lowers the free energies of these phases by adding triplets to the calculation of the local contribution to the free energy in equation (9), as with the superstructures in Section 3.2.2, but comes at a cost in the integral overṼ V LP ðkÞ, where n and k are both 2 for the lamellar phase, and 6 and 3 1/2 , respectively, for the hexagonal one. Altogether, this gives and These are minimized numerically over , k and for each value of c and and compared with the free energies of the decagonal or dodecagonal quasicrystals, so that the minimum free-energy phase can be identified.

c-Dependent phase diagrams for decagonal and dodecagonal quasicrystals
The c-dependent phase diagrams, calculated using the tworing approximation, are shown in Figs. 7 and 8, for the decagonal and dodecagonal quasicrystals, respectively. In both cases, one clearly observes that length selectivity, as parameterized by c, is a key factor contributing to quasicrystal stability. As c is decreased, the competing phases change from a solution where = 1 and k = 0 to one where is shifted and both and k are nonzero, in order to take advantage of both minima ofṼ V LP ðkÞ in an optimal way. This causes the upper bound of for quasicrystal stability to constrict with decreasing c until the quasicrystalline phase vanishes. This vanishing occurs at a uniform-hexagonal-quasicrystal triple point. Below the triple point, the critical for the uniformhexagonal transition continues to decrease towards the value  Phase diagram of the LP model with q = k 5 in the two-ring approximation: the thick dark-green hexagonal-decagonal boundary, which was calculated using the projection method by Jiang et al. (2015), indicates that this phase diagram should be considered only qualitatively. Nevertheless, note the expected uniform-hexagonal-decagonal triple point, and the possibility that the lamellar-hexagonal coexistence line has a maximum for intermediate c before it reaches its expected value in the limit of infinite c.

Figure 8
Phase diagram of the LP model with q = k 12 in the two-ring approximation: in the dodecagonal case, the thick dark-green hexagonal-dodecagonal boundary, which was calculated using the projection method of Jiang et al. (2015), shows very good agreement with the results of the two-ring approximation. Note the uniformhexagonal-dodecagonal triple point and the lamellar-hexagonal coexistence curve exhibiting a turning point at intermediate c with a minimum value of . $ À0.1143 at zero c, as calculated earlier for the two-scale hexagonal phase in Section 3.2.2.
The portion of the hexagonal to quasicrystal phase boundary, calculated by Jiang et al. (2015) using the more accurate projection method, 9 is shown on both phase diagrams using thick dark-green lines. These lines indicate that, while both approximate phase diagrams qualitatively agree with the projection-method calculation, only the dodecagonal phase diagram agrees with it quantitatively. This is because, at the relevant c values, the free-energy barrier between the two minima ofṼ V LP ðkÞ in the decagonal (q = k 5 ) case, shown in Fig. 1, is sufficiently low that many higher-harmonic peaks appear in this region and stabilize the decagonal phase relative to the hexagonal one, which has no additional Fourier peaks there. This enlarges the stability region of the decagonal phase compared with what we calculated by restricting its Fourier coefficients to two exact circles. On the other hand, the freeenergy barrier between the two minima ofṼ V LP ðkÞ in the dodecagonal (q = k 12 ) case is much steeper, preventing additional rings from forming. This leads to the two-ring approximation and its predictions for the position of the triple point and the precise shapes of the phase-boundary curves, being quantitatively reasonable when q = k 12 but only qualitatively valid when q = k 5 . Lifshitz & Petrich (1997) speculated that, with more than just two length scales, the LP model could stabilize quasicrystals with higher orders of symmetry than the dodecagonal structures they obtained, such as 18-or 24-fold. In the meantime, such structures have been observed experimentally (Fischer et al., 2011) and in simulations (Engel & Glotzer, 2014;Dotera et al., 2014;Pattabhiraman & Dijkstra, 2017b). In addition, Arbell & Fineberg (2002) discovered patterns with 8-fold symmetry in Faraday wave experiments using three driving frequencies. Here, we study the infinite-c LP model with four length scales, by modifyingṼ V 0 in equation (11) from {1, q} to {q 1 , 1, q 2 , q 3 }. We consider two cases: (i) octagonal quasicrystals, with q 1 = k 8/3 = 2cos(3/8) = (2 À 2 1/2 ) 1/2 ' 0.7654, q 2 = k 4 and q 3 = k 8 = (2 + 2 1/2 ) 1/2 ' 1.848; and (ii) octadecagonal quasicrystals, with q 1 = k 18/7 ' 0.6840, q 2 = k 18/5 ' 1.286 and q 3 = k 18 ' 1.970. The anticipated diffraction spectra of these two structures are shown in Fig. 9.

Four-scale octagonal and octadecagonal quasicrystals
With more than two length scales, one must carefully check for competing structures, additional to the single-scale phases in Section 3.2.1. For the octagonal case, we must consider the two-scale square superstructure shown in Fig. 3(i), allowed by the fact that q 2 = k 4 .
For the octadecagonal case, additional competing structures stem from the fact that q 1 + q 2 = q 3 . This allows for the 'modified' lamellar and hexagonal candidates shown in Fig. 10. These have free energies of and where n q AE = P 3 i¼1 n q i and q Å = Q 3 i¼1 q i . Minimizing these equations in the relevant range shows that the coincidental lamellar phase has q 1 = q 2 = q 3 and a free energy degenerate with the single-scale hexagonal phase. The coincidental hexagonal phase has only one ring of active modes when its free energy is minimized and so does not have its free energy reduced relative to the single-scale hexagonal structure. Thus, we can continue treating the usual single-scale hexagonal phase as the only candidate competing with the octadecagonal quasicrystal for the relevant values of .
Applying equation (15) to the octagonal structure in Fig. 9(a) gives a free energy of Fourier spectra of the four-scale quasicrystals|: (a) for the octagonal quasicrystal, the radii of the circles from inside to outside are k 8/3 , 1, k 4 and k 8 ; (b) for the octadecagonal quasicrystal, the radii are k 18/7 , 1, k 18/5 and k 18 .

Figure 10
Fourier spectra of additional structures competing with the four-scale octadecagonal quasicrystal: the radii of the circles are the ones listed for Fig. 9(b). 9 We thank the authors for graciously sharing their raw data with us.
where n AE = n 1 + n q AE . While it is lengthy, it is not difficult for a computer to minimize this quartic numerically over the four 's for each value of . Interestingly, despite the quartic lacking q 1 $ q 2 and q 2 $ q 3 symmetry, the minima in the relevant small range all satisfy q 1 = q 2 = q 3 .
As shown in Fig. 11(a), the octagonal quasicrystal is expected to be thermodynamically stable when À0.1119 < $ < $ À0.06227. In this range, the optimized 1 is larger than the equal q 's by a factor varying between roughly 2.1 and 1.7. The free energy at the transition to the two-scale supersquare phase is approximately À1.094 Â 10 À3 . A finite section of this quasicrystal is shown in Fig. 12.
The free energy of the octadecagonal structure is As shown in Fig. 11(b), the octadecagonal quasicrystal is expected to be stable when À0.06882 < $ < $ À0.05140. In this range, the optimized 1 is larger than the q 's by a factor varying between roughly 2.2 and 2.6. The q 's are almost, but not exactly, equal, varying by about 1%. The free energy at the transition to the hexagonal phase is approximately À2.085 Â 10 À4 . A finite section of this quasicrystal is shown in Fig. 13. 6. The density distribution method 6.1. The notion of a density distribution and its quantum density of states analogy When the local free-energy function f ðÞ is a finite polynomial such as equation (7)  Free energies of four-scale quasicrystals in the infinite-c LP model: (a) the octagonal structure is the equilibrium phase when À0.1119 < $ < $ À0.06227; (b) the octadecagonal structure is the equilibrium phase when À0.06882 < $ < $ À0.05140.
given candidate configuration can be evaluated using the approach of equation (15). However, this is not the case when the local free-energy function is non-polynomial. We describe here an alternative technique, which we call the 'density distribution method', that not only allows us to evaluate such free energies, but also provides a novel way of understanding the stability of the various periodic and quasiperiodic phases. This understanding is applied in Section 7 to calculate the free energy (3) of the candidates under the BDL model and explain the surprising stability of certain decagonal quasicrystals, and in Section 8 to aid in the artificial design of a new local-energy function f ðÞ which stabilizes arbitrarily highorder quasicrystals with only two length scales. Rather than evaluating the integral (12) -used to calculate the contribution of the local term to the free energy -over space, this term can be summed differently by integrating over the set of possible values that (r) may take on, where is normalized such that Z PðÞ d ¼ 1; and is the Dirac delta function. PðÞ is the probability density function of (r) being . Intuitively, it is essentially a histogram of the position-space (r) values obtained when r is selected by blindly throwing a dart at the entire Cartesian plane on which the structure is defined. This reformulation of the free energy is conceptually reminiscent of Lebesgue integration. It can also be thought of as taking a uniformweight inner product of f and P over the vector space of real functions. An intriguing analogy exists between the density distribution and the density of energy eigenstates of a quantum Hamiltonian for a single particle in a periodic potential. There, it is the energy dispersion, or band structure, E(k), that plays the role of our density field (r). When the Hamiltonian is that of a nearest-neighbor tight-binding model for a particle hopping on a lattice, corresponding to one of our candidate structures, the band structure E(k) assumes the same form taken by our density field (r) and the analogy becomes exact. 10 Consequently, the formal expression for the calculation of the density distribution in one dimension is similar to that of the density of states 11 where in d dimensions the sum is replaced by an integral over the d À 1-dimensional equal-surface. We can take the analogy one step further if we noticeusing standard complex analysis -that the density distribution PðÞ ¼ ÀImfGðÞg=, where It turns out that this function GðÞ is the on-site lattice Green's function, obtained for the nearest-neighbor tightbinding Hamiltonian of a particle hopping on the corresponding lattice, with (r) replaced by E(k). The real and imaginary parts of the lattice Green's function are related to each other by the Kramers-Kronig relations and so encode equivalent information about the crystal structure. The real part can be useful for evaluating the density distribution-based free-energy equation (36) for certain local free-energy density functions f ðÞ, such as the one in Section 7. Many advanced mathematical approaches have been developed for evaluating these lattice Green's functions, including contour integrals (Ray, 2014), hypergeometric functions and Calabi-Yau differential equations (Guttmann, 2010), holonomic functions (Koutschan, 2013;Zenine et al., 2015;Hassani et al., 2016), and Chebyschev polynomials (Loh, 2017).

Rescaling and skewness of the density distribution
As for any normalized probability distribution (38), rescaling the field strength (r), and with it the width of the density distribution, merely rescales the distribution itself by the reciprocal factor, namely P ðÞ = P ðÞ=jj. In the case of single-scale structures, whose overall field strength is determined by a single Fourier amplitude , it is therefore sufficient to calculate the density distribution once for P ¼1 ðÞ, and rescale later if necessary.
For two-scale structures, characterized by two Fourier amplitudes 1 and q , a rescaling of (r) affects both amplitudes together, giving P 1 ; q ðÞ = jjP 1 ; q ðÞ. Thus, the density distributions differ for fields with different ratios q = 1 of the amplitudes, but are otherwise independent of the overall scale of the field.
For all the structures relevant to us, PðÞ has compact support ½ min ; max between the extreme values of (r), which are both finite, because the fields are all finite sums of harmonic functions. The value = À max = min is a measure of the 'skewness' or 'lopsidedness' of the density distribution. It characterizes the imbalance between the ground state and the highest excited state of the corresponding tight-binding model. It is a useful measure that will serve us in what follows.
Because of the freedom to rescale the density distribution, for single-scale structures can take on at most only two distinct values: for positive and its inverse 1/ for negativẽ . We need only consider positive . On the other hand, with this assumption, for two-scale structures varies continuously as a function ð q = 1 Þ of the ratio of the two Fourier amplitudes. One may need to resort to numerical sampling of the field (r) in order to generate the density distribution PðÞ when analytical methods for calculating equations (37) or (39) prove difficult. For periodic crystals this is readily achieved by uniformly sampling the unit cell of the crystal in both spatial directions. Quasicrystals lack periodicity, so this approach would, in principle, require a uniform sampling of the entire two-dimensional plane.
An alternative approach for the periodic case, which is easier to generalize to quasicrystals, is to remain at the origin of the two-dimensional plane and shift the field itself, until a full unit cell is sampled. This procedure samples the origin of the degenerate minimum free-energy states, which in the periodic case merely differ by a rigid translation within the unit cell.
One can sample the minimum free-energy states in terms of the Fourier coefficients of the fields by ensuring that the value of the free energy -like the one in equation (13) -does not change. This implies that one may generally shift the phases of the (complex) Fourier coefficients as long as the sum of these phases is zero for all possible structure invariants. This amounts to performing a Rokhsar-Wright-Mermin gauge transformation (Rokhsar et al., 1988), as explained elsewhere (Lifshitz, 2011). Thus, one may freely shift the phases of the Fourier coefficients on wavevectors that are linearly independent over the integers. All the phase shifts of the remaining Fourier coefficients are then determined by the structure invariants. For periodic crystals in two dimensions there are two such independent phases. For the decagonal and dodecagonal quasicrystals of interest here there are four independent phases. Shifting these phases uniformly from 0 to 2 yields the uniform sampling that we seek. 12 6.4. Density distributions for the candidate phases 6.4.1. Single-scale structures. Trivially, P UNIF ðÞ = ðÞ. We therefore begin with the single-scale lamellar field which, after scaling to unity, is given by (r) = 2cos x. Because jj never exceeds two, P LAM ðÞ vanishes when jj > 2. We express it analytically between these bounds using equation (39), by substituting À2 sinðxÞx x for rðrÞ and cos À1 ð=2Þ for x and normalizing according to equation (38). This gives This calculation is demonstrated schematically in Fig. 14. While calculating P LAM is straightforward, doing so for P HEX is quite difficult. The mathematics necessary to do so was worked out by Ramanujan (1914) using one of his theories of elliptic functions to alternative bases. The connection to lattice Green's functions was introduced by Horiguchi (1972). We simply provide the result, where K is the complete elliptic integral of the first kind. These density distributions are plotted in Fig. 15, showing that the lamellar phase is unskewed with LAM = 1 as expected, while the hexagonal phase is skewed, with HEX = 2.
6.4.2. Two-scale structures. The density distributions of the decagonal and dodecagonal fields are calculated numerically by measuring them at the origin, as explained earlier, while sampling the set of all degenerate minimum free-energy states via appropriate phase shifts of the harmonic functions. The phase-shifted fields are given by Graphical evaluation of the lamellar density distribution: the left-hand panel shows a half-period of the sinusoidal lamellar spatial structure. The horizontal lines coming off it are evenly spaced in the horizontal direction. Their vertical density determines the distribution in the righthand panel. Note how the stationary regions of the structure, where the gradient vanishes at 0 and , lead to the inverse square root Van Hove singularities in the density distribution at jj = 2.
We call 1 the value of (r = 0; { i }) when 1 = 1 and q = 0 and likewise call q the value when 1 = 0 and q = 1, so that = 1 1 + q q . We numerically sample the ordered pairs ð 1 ; q Þ uniformly over i 2 ½0; 2Þ, with 96 sampling points along each phase, to give a total of roughly 85 million samples. If we write this distribution function as P ð 1 ; q Þ, then Upon numerically maximizing the skewness parameter over the ratio of amplitudes q = 1 for the decagonal and dodecagonal phases, we find that the optimal ratio is unity, where 1 = q > 0, in both cases. The density distributions obtained for this ratio are shown in Fig. 16, where it can be seen that the decagonal and dodecagonal phases have maximal values of 4 and 3, respectively.
The density distributions in Figs. 15 and 16, along with equation (36), allow us to understand the stability of the candidate phases more generally than before. The uniform phase simply has F UNIF = f(0). The lamellar phase can potentially dominate this by heavily sampling the local freeenergy density f ðÞ far from = 0. For example, a local freeenergy density with a strongly negative quadratic component would likely favor the lamellar phase. Indeed, this is what we observe in Section 3.2.1 when exceeds (6 3/2 + 14)/15 and the free energy of the lamellar phase is lower than that of the hexagonal phase. On the other hand, the lopsided nature of the hexagonal structure allows it to take advantage of the odd components of the local free-energy density.
Similarly, quasicrystalline structures also attain stability through the skewness of their density distributions. In particular, their extremes can be more lopsided than those of the hexagonal phase. In other words, they have > HEX = 2. In Section 8, we use this feature of quasicrystals to stabilize arbitrarily high-order structures using only two length scales.

Van Hove singularities
As shown in Figs. 15 and 16, the density distributions exhibit a variety of Van Hove singularities (Van Hove, 1953). Note the inverse square root Van Hove singularities exhibited by the lamellar structure and the zeroth-order step discontinuities and logarithmic singularities in the hexagonal density histogram. The logarithmic singularity at = À2 is due to the corresponding stationary point in the structure function being a saddle, to leading order, unlike the extrema which lead to the step discontinuities.
By numerically minimizing the magnitude of the gradients of the effective four-dimensional periodic functions [equations (42a) and (42b)], we can identify Van Hove singularities at min and max , at AE4/5 in the decagonal structure, and at À1/2, À3/8 and 0 in the dodecagonal structure.
With the exception of the discontinuity in the decagonal distribution at min = À1, all of the quasicrystal Van Hove singularities are first order. This zeroth-order Van Hove singularity turns out to be crucial for stabilizing the decagonal structure under the BDL model in Section 7 and can be understood in terms of the spatial Hessian of the effective four-dimensional function. One of the phase-shifted structures Density distributions of the uniform, lamellar and hexagonal structures: note the Van Hove singularities (i) the inverse square root singularities of the lamellar distribution at min = À2 and max = 2, (ii) the logarithmic singularity of the hexagonal distribution at = À2, and (iii) the discontinuous jumps from 1/(3 1/2 ) and 1/[4(3 1/2 )] to zero at min = À3 and max = 6, respectively. The skewness of the latter two distributions is given by LAM = 1 and HEX = 2.

Figure 16
Density distributions for decagonal and dodecagonal quasicrystals: for these quasicrystals, is maximized when 1 = q = 1/5 or 1/8, respectively. Observe that DEC = 4 and DOD = 3, and note the interior Van Hove singularities at AE4/5 and À1/2, À3/8 and 0, and the zeroth-order discontinuity in the decagonal distribution at min = À1. This final Van Hove singularity is analyzed in Section 6.5 and is a key factor leading to the decagonal structure's stability under the BDL model. that has this minimum of À1 at its origin is given by { i } min = {sec À1 (À4), 0, 2 sec À1 (4), 0}. At this point, the Hessian is 8 6 4 2 6 12 18 9 4 18 32 16 2 9 16 8 The nullity, or the rank of the kernel, of this matrix is two. In this case, this indicates that the manifold of the minima is twodimensional. This generates a quadratic minimum of effective dimension two, the rank of the Hessian. A quadratic d-dimensional minimum generates a singularity of order d=2 À 1, so the decagonal density distribution has a zerothorder discontinuity at its minimum. As explained by Van Hove (1953), topological considerations in Morse theory require the existence of a certain number of stationary points with each quadratic signature, although degenerate stationary points such as the one analyzed in the previous paragraph complicate the situation.

A lower bound on the LP free energy
Using the language of density distributions, we can calculate a lower bound for the free energy in the LP model. Clearly, if it were not for the infinite-c penalty at k = 0, which requires the average density to be zero, the best possible density distribution would have a single delta-function peak, corresponding to a uniform field (r) = 0 that minimizes the local freeenergy density (7) everywhere. To satisfy the zero-average requirement we must add a compensating delta function peak at some negative value ' < 0, and possibly allow the positive value r > 0 to shift away from 0 . This yields a field with sharp boundaries between two allowed values and a density distribution of the form PðÞ = P ' ð À ' Þ + P r ð À r Þ.
Given sufficient harmonics, structures of arbitrary symmetry can attain this sharp form. Of course, with too many allowed length scales, the physical requirements on the lengthscale selectivity c are stricter, and even so the set of competing candidate phases can increase, so the results below provide only an extreme theoretical lower bound on the free energy.
All that is left is to minimize the free energy (36) under the constraints of zero-averaging, P ' ' + P r r = 0, and the normalization P ' þ P r = 1 of the density distribution. Solving for the left-hand side variables gives P ' = 1 À P r and ' = ðP r r Þ=ðP r À 1Þ. Substituting them into PðÞ and minimizing the free energy (36) over P r and r gives Essentially, we have fitted the right peak into the wells of the solid colored lines in Fig. 2, while remembering that it must be balanced out by a corresponding peak at negative . This gives a lower bound on the free energy of F ! Àð9 þ 2Þ 2 =324.
Note that this implies that must be greater than À2/9 to allow for the possibility of structures with negative free energy. This 'forbidden zone' is shown as the grayed-out regions in Figs. 4-6. Furthermore, the lowest possible metastability bound is = À1/3.

7.
Mean-field theory for soft interacting particles 7.1. The Barkan-Diamant-Lifshitz model Equipped with the density distribution method and the ability to calculate free energies with non-polynomial local terms, we can perform a more detailed and informed analysis of the coarse-grained free energy (3) used in the BDL model (Barkan et al., 2011), which contains a local logarithmic entropy term. Assuming a sufficiently dense system of soft particles that interact via a Fourier transformable isotropic pair potential UðrÞ -implying that it does not diverge at a higher order than the usual 1/r electrostatic potential as the particles get closer together -one can express the BDL coarse-grained free energy in the form of equation (9) with V VðkÞ ¼Ũ UðkÞ ÀŨ U min 8 2 jŨ U min j ; and whereŨ UðkÞ is the Hankel transform of UðrÞ andŨ U min is its minimum value. The temperature T is measured here in units of the spinodal temperature T sp = ÀcŨ U min =k B , where c is the average number density of the particles and k B is the Boltzmann constant. Recall that T sp is the lower metastability boundary of the uniform liquid phase, below which the system must become ordered, and note that the minimumŨ U min of U UðkÞ must be negative for this temperature to be positive. 13 Finally, the value of the field (r) is here constrained to > À 1 by the fact that the number density c(r) = c ½ðrÞ þ 1 of the particles cannot be negative. This 'vacuum constraint' ensures that the logarithm in equation (49) does not diverge.
BDL proceeded to take the fourth-order Taylor expansion of f CG to obtain and mapped this resulting quartic free energy onto the LP free energy, giving them a rough estimate of the physical parameters that might stabilize the different targeted structures, based on the LP results. By rescaling and F , we can make f 4 equivalent to the LP form in equation (7). The correspondence between T and necessary to do so is then given by T = 4=ð3 þ 4Þ. Note that the range À4/3 < < 0 corresponds to scaled temperature values of T > 1 above the spinodal decomposition, and that positive values of correspond to values of T < 1 below it. The cases of À4/3 and T < 0 are unphysical for the coarse-grained free-energy model. The success of the estimates obtained by BDL through this mapping were somewhat fortuitous, as a comparison between even properly rescaled plots of f LP and f CG in Fig. 2 reveals that they are very different outside of the radius of convergence jj > 1. As the transition from the uniform liquid to the ordered state is first order, the field (r) generally contains regions with large values, making f 4 a poor approximation even at the transition. It is therefore important that we can now evaluate the exact local free energy f CG . As for the evaluation of the nonlocal term of the free energy, in order to make the analytical calculations more tractable, we again work in the limit of exact length-scale selection -analogous to the infinite-c limit of the LP model -corresponding here to the smallŨ U min and therefore small T sp limit. As seen below, we indeed find important differences between the behaviors of the LP and BDL models.

Single-scale structures
For the single-scale lamellar and hexagonal phases, substituting the density distributions of equations (41a) and (41b) and the logarithmic local free-energy density of equation (49) into the free-energy equation (36) gives and where 2 F 1 is a hypergeometric function. 14 Now, minimizing F LAM over yields which is shown in Fig. 17. At absolute zero, the free energy is exactly À1/4. In the first temperature range, where F LAM is linear in T, = 1/2, which is the maximum value it can obtain without violating the vacuum constraint. Increasing it further would violate the non-negative density condition. At the transition to the second temperature range, the free energy is ð1 À 2 ln 2Þ=4. From here, as T increases to unity, decreases to zero. At T = 1, a second-order phase transition to the uniform phase occurs. The result of minimizing F HEX , which is obtained numerically, is also shown in Fig. 17, displaying similar behavior. At absolute zero, the free energy is exactly À1/3. For 0 T < $ 0.8514, has its maximum allowed value, which is 1/3, and the free energy is linear in T. When T ' 0.8514, F ' À0.05278, and in the next temperature range decreases monotonically to $0.1923 at T ' 1.063, at which point the system undergoes a first-order transition to the uniform phase.
Note that the hexagonal phase always has a lower free energy than the lamellar phase. Therefore, we do not need to consider the single-scale lamellar candidate when evaluating the stability of quasicrystals in the BDL model, as it is never the equilibrium phase. This is qualitatively different from the behavior with the quartic local free energy f 4 , analogous to that of the LP model, where the lamellar phase would be expected to take over at T (34 À 6 3/2 )/47 ' 0.4107.

Two-scale structures
In the original LP model, one can set q to k 4 , k 6 and k 1 to stabilize two-scale square, hexagonal and lamellar superstructures, respectively. Unsurprisingly, the BDL model can also do this, as demonstrated in Fig. 18, and it stabilizes these two-scale periodic phases all the way down to absolute zero. Their free energies are calculated using the same techniques as the quasicrystalline structures treated below, but we omit a detailed analysis of their behavior.
Again, as in the LP model, setting q to k 8 and k 10 fails to stabilize octagonal and decagonal quasicrystals, as their free energies, shown in Fig. 18, are always greater than that of the single-scale hexagonal phase. Only decagonal and dodecagonal structures with q = k 5 and k 12 occur as stable quasicrystalline states in the BDL model. Their free energies are calculated using the sampled distribution of ordered pairs Free energies of the single-scale periodic structures in the BDL model: in this model, the free energy of the hexagonal structure is always lower than that of the lamellar structure, so we need not consider it further. The hexagonal structure is the equilibrium phase up to T ' 1.063, where it undergoes a first-order transition to the uniform liquid phase. Colored dots correspond to temperatures below which the Fourier amplitude reaches its maximum value and the free energy becomes a linear function of the temperature.
Pð 1 ; q Þ as explained in Section 6.4.2, and plotted in Fig. 19 as functions of T.
We would like to emphasize that the minimization of the free energy with respect to 1 and q is performed subject to the vacuum constraint, which is more difficult to take into account for the two-scale structures. Each ð 1 ; q Þ pair from Section 6.4.2 not only gives a small free-energy contribution, but also imposes the linear constraint 1 1 + q q ! À1 on the values of 1 and q . The problem of finding the intersection of these half-planes is dual to that of finding the convex hull of the ð 1 ; q Þ points (de Berg et al., 2008). If ð 1A ; q A Þ and ð 1B ; q B Þ are adjacent extremal points on the convex hull, then the point is a vertex of the polygonal boundary of the set of allowed ð 1 ; q Þ values that do not violate the vacuum constraint. The feasible sets for the decagonal and dodecagonal quasicrystals are displayed in Fig. 20. As shown there, we are able to find exact values for all the polygonal vertices bounding the regions. A simple constrained descent algorithm is used to minimize the free energy over these convex sets.
For the decagonal phase with q = k 5 , a portion of which is shown in Fig. 21, as the temperature is lowered the system undergoes a first-order phase transition from the uniform liquid to the quasicrystal at T ' 1.125, at which point 1 = q ' 0.1352. These 's increase together until T ' 0.8977, where they reach their maximal allowed value of 1/5. At this point, the free energy of the decagonal phase is approximately À0.06746. This continues to be the equilibrium phase all the way down to absolute zero, where the free energy becomes exactly À2/5.
For the dodecagonal phase with q = k 12 , the first-order transition from the uniform liquid occurs at T ' 1.159, wherẽ 1 = q ' 0.1220.

Figure 18
Free energies of the two-scale periodic and unstable quasiperiodic structures in the BDL model: plotted for comparison are the free energies of the uniform and single-scale hexagonal phases. q takes the values k 4 , k 6 , k 8 , k 10 and k 1 , as labeled. As in the LP model, the free energies of the two-scale square (k 4 ), hexagonal (k 6 ), lamellar (q = k 1 ), decagonal (k 5 ) and dodecagonal (k 12 ) structures are lower than that of single-scale hexagonal structure, but unlike the LP model the decagonal structure remains the equilibrium phase down to zero temperature, without undergoing a second phase transition. On the other hand, still in line with the LP model, the free energies of the octagonal (k 8 ) and decagonal (k 10 ) quasicrystals are always higher than that of the single-scale hexagonal structure, and are therefore never the equilibrium phase.

Figure 19
Free energies of the stable two-scale quasiperiodic structures in the BDL model: the decagonal structure (k 5 ) is the equilibrium phase for T < $ 1.125 and the dodecagonal structure (k 12 ) is the equilibrium phase for 0.8697 < $ T < $ 1.159. Black dots mark the phase transitions, while colored dots mark transitions between linear and nonlinear regimes of the free energy, where the Fourier amplitudes 1 and q reach a stationary pair of values.

Figure 20
Feasible ð 1 ; q Þ sets for decagonal and dodecagonal quasicrystals in the BDL model: these are the convex regions over which the Fourier coefficients can vary without violating the vacuum constraint. There is a mirror line given by 1 = q . Solid dark lines in the first quadrant show the path of the minimizing values of ð 1 ; q Þ as the temperature is changed. While the polygonal vertices are exact, we believe that there is no simple expression for the two smooth curves in the dodecagonal feasibility boundary. Note that the peakedness of the decagonal structure jutting far into the first quadrant essentially explains its surprising stability in the BDL model. Similar shapes exist for the additional phases considered in Fig. 18, but are omitted here.
10 À3 and the 's reach their maximum allowed value of 1/8, and the free energy as a function of temperature enters a linear regime. The 's remain at 1/8 until T ' 0.8697, at which point the hexagonal phase takes over at a free energy of approximately À0.04697. Regardless, if we continue to restrict our attention to the dodecagonal structure, it undergoes a second-order phase transition-like event where the Z 2 symmetry 1 $ q is broken at T ' 0.8446 and F ' À0.05086. After this point, the F 's move along their maximum allowed sum line 1 + q = 1/4 until T ' 0.7022, where they land in either of the two degenerate states ð 1 ; q Þ = (7/36, 1/18) or (1/18, 7/36) with a free energy of approximately À0.07973. The structure remains in one of these two minima until absolute zero, where the free energy becomes exactly À53/216.

Skewness and the vacuum constraint
While the dodecagonal quasicrystal shows qualitatively similar stability under the LP and BDL models, the decagonal quasicrystal exhibits significantly different behavior, showing surprisingly robust stability in the BDL model. Until now, it was understood that soft quasicrystals are stabilized by threebody (or more generally odd-body) interactions that break the Z 2 symmetry of the free energy F, namely, the ! À symmetry of f ðÞ. However, the fact that the decagonal quasicrystal dominates even at T = 0, where there are no three-body interactions, demonstrates that this cannot be the whole story.
The primary quasicrystal stabilizing factor in the BDL model is the important ! À1 vacuum constraint, which is a very effective alternative way of breaking the Z 2 symmetry of F . The decagonal structure with 1 = q possesses the most lopsided of any of the density distributions examined, giving it a skewness of four, and allows the decagonal phase to take maximal advantage of the the highly negative f CG ðÞ values at large 's without violating the vacuum constraint. The important Van Hove singularity at the vacuum minimum, which allows this high skewness to occur, is analyzed in Section 6.5. Note that the structure itself, shown in Fig. 21, contains well separated very highly peaked positive red spots in a shallow blue sea of negative values. This argument also explains the result of the molecular dynamics simulations performed by Barkan et al. (2014) with q = k 5 , showing that the decagonal quasicrystal remains stable as T is lowered to absolute zero, without undergoing a transition to the hexagonal phase as would be predicted by a naïve correspondence to the LP model. 8. Higher-order quasicrystals 8.1. Skewness of the 6n-fold two-scale structures Finally, we examine the skewness of the density distributions of quasicrystals of order 6n for arbitrary n ! 2. Then, using the information obtained, we judiciously design an artificial local free-energy function f ðÞ that allows for the stabilization of quasicrystals of these orders. While the local free energies previously used by Lifshitz & Petrich (1997) and Barkan et al. (2011) were physically justified by entropic considerations and their truncated polynomial expansions, the one we construct below is engineered with the sole goal of stabilizing higher-order quasicrystal phases. However, as we argue below, it might not be impossible to design a physical system with sufficiently similar behavior to stabilize some of these higher-order quasicrystals, particularly when n is not too large.
We first demonstrate that, for all n, two-scale 6n-fold crystals, like the ones shown in Figs. 3(j) and 3(m) for n = 1 and 2, all have their skewness from Section 6.2 greater than two, when the two amplitudes 1 and q are equal. We restrict our attention to this case and scale both 1 and q to unity.
By inverse Fourier transforming its momentum-space representation according to equation (8), the position-space field representing this quasicrystal can be written as Predicted decagonal quasicrystal with 1 = q = 1/5: red shades correspond to positive field values 4, whereas blue shades represent negative values which just barely scrape against the vacuum when = À1. Note the abundance of blue or white areas which are interspersed with bright red spots. This provides the skewness which makes this structure so stable in the BDL model.
have length |k j | = 1, and the sum of two consecutive vectors has a length |k j + k j+1 | = q = k 6n . 15 In principle, the additional phases j and j,j+1 are free to vary so as to minimize the free energy (12), yet for similar arguments mentioned in Section 3.1.1 there is always a representative structure, within the set of all degenerate minimum free-energy states, for which the phases within each circle are all equal and can be taken to be 0 or . We limit our attention to structures where the phases are the same on both circles, taking them all to be zero without any further loss of generality, owing to our freedom to change the sign of the cubic term in f ðÞ accordingly. With this choice, with all the 's set to zero, the field obtains its maximum value max = 12n at the origin. We now show that min > À6n so that > 2.
As explained in Section 6.3, we sample the field by staying at the origin r = 0 and shifting the phases, while keeping the structure invariants constant, thereby performing Rokhsar-Wright-Mermin (Rokhsar et al., 1988) gauge transformations. This requires the sum of the phase shifts at wavevectors that add to zero to vanish, and immediately implies that j,j+1 j + j+1 , where '' stands for equality modulo 2. Thus, the phase shifts on the outer circle are all fixed by the choice of shifts on the inner circle.
In addition, within the inner circle, each vector and its negative impose the constraint j + j+3n 0, and each triplet of wavevectors adding to zero imposes the constraint j + j+2n j+n . This leaves at most two independent phases on each of the n sextets forming the inner circle that still need to satisfy additional constraints imposed by each additional prime divisor of n other than 2 or 3. The resulting number of independent phase shifts is given by Èð6nÞ, where È is the Euler totient function.
This can be used to rewrite the shifted field at the origin as Each of the triplets of cosines in the two sums of equation (56) can be written in the form cos x + cos y + cos (x + y). This function has a minimum of À3/2 when both x y AE2/3. However, this bound is unattainable for all 2n cosine triplets: If we set m m+2n AE2/3 for m = 1 . . . n, so as to obtain the À3/2 minimum for all the triplets of the first sum, then no matter what the sign choices are, at least one of the triplets in the second sum must have phases of zero and therefore does not achieve the À3/2 minimum. This implies that min > À6n, and so > 2. Indeed, numerical sampling for 1 = q shows that = 3, 3, 36/13 ' 2.769 and $2.634, for n = 1 . . . 4, respectively. Asymptotically, appears to decrease no faster than 2 + 2/n.

A local free-energy density that stabilizes 6n-fold quasicrystals
We set our local free-energy density to be which is shown as the black lines in Figs. 1 and 22. Using the density distributions calculated in Section 6.4, it is not difficult to show that the uniform and single-scale lamellar and hexagonal phases must have non-negative free energies with this local free-energy function. In particular, because HEX = 2, the hexagonal phase cannot probe the negative f ð > 2Þ ¼ À1 region without suffering greater free-energy penalties from the positive f ð < À 1Þ ¼ 2 region, as can be inferred graphically from Fig. 22. Furthermore, because the freeenergy penalty is twice the free-energy decrease, the hexagonal phase is unable to reach negative free energies through negative values of either. For the 6n-fold two-scale structures considered above, we scale 1 = q until min reaches À1. Then, because > 2, max is also greater than two. This, together with the density distribution (36) and the local free-energy function (57), implies that the free energy F is negative for this structure. Thus, a 6n-fold quasicrystal is the minimum free-energy state Density distribution of the two-scale octadecagonal quasicrystal: the black line is the local free-energy density from equation (57). Here q = k 18 ' 1.970 and 1 = q = 1/13. Note that the purple-colored octadecagonal density distribution extends into the positive values of > 2, where the free-energy density is negative, making its overall free energy F ' À1.223 Â 10 À3 < 0. On the other hand, the density distribution of the single-scale hexagonal structure, which is plotted here for reference, cannot extend beyond = 2 without running into the barrier at < À1, which would force its free energy to become positive. This approach for forcing the quasicrystal structure to be the minimum free-energy state succeeds theoretically for all 6n-fold quasicrystals, where n ! 2, although they become increasingly fragile.
of the system. Indeed, the calculated free energies of the octadecagonal and icositetragonal quasicrystals are approximately À1.223 Â 10 À3 and À6.093 Â 10 À5 , respectively. These free energies are simply the fraction of the density distribution above a density of two when the 's are scaled such that min = À1, as can be seen graphically in Fig. 22.
Only some of the features of the contrived local free-energy density (57) are necessary for this stabilization to occur, and lower orders will be much more tolerant of imprecision than higher orders. For arbitrarily high orders, the flat region which includes = 0 in the local free-energy functional is essential to this argument, as is some degree of favorable free energy for high values of and a somewhat greater free-energy penalty for sufficiently negative ones. These deviations from zero do not have to be sudden jumps. If the quasicrystalline can only be proven to be greater than two, as we have done here for 6nfold structures, the ratio of the onset of these latter two effects must be exactly two, but if it can be shown to be even higher there will be some room for error.
Reproducing these results in the laboratory is likely to be challenging for at least three reasons: (i) High length-scale selectivity will be required.
(ii) The thermodynamic stability of a stable state does not necessarily imply that it is kinetically accessible within a reasonable time frame.
(iii) Engineering an effective local free-energy density function like the one in equation (57) may be difficult. We suggest using a system similar to the interacting particles in the BDL model, where the vacuum constraint (r) ! À1 implements the required barrier for negative concentrations relative to the average. A drop in the local free-energy density for sufficiently high concentrations could potentially be implemented through a kind of local phase change which sets in at a critical density, or some other highly nonlinear effect, such as the formation of oscillons (Umbanhowar et al., 1996;Arbell & Fineberg, 2000).

Closing remarks
In closing, we wish to emphasize the ease with which one can stabilize quasicrystals in rather simple isotropic models of interacting particles or their mean-field descriptions. It was appreciated from the outset that one needs to introduce multiple length scales into the interaction potentials of the constituent particles. Yet the ability to do so in a quantitatively predictive and controlled manner has only emerged in the last two decades, based on the understanding of how the multiple scales 'work together' to produce the targeted quasicrystalline structures.
The Faraday wave experiments of Edwards & Fauve (1993) led to the understanding of Lifshitz & Petrich (1997) that one needs to break the ! À symmetry of the Landau freeenergy expansion in order to allow the two length scales to couple via triad resonances or effective three-body interactions. This understanding was then generalized by Barkan et al. (2011) with their symmetry-breaking logarithmic entropy term.
Here, enabled by the density distribution method for calculating such non-polynomial free energies, we have come to an even deeper understanding that the breaking of ! À symmetry favors the formation of structures with skewness. Indeed, quasicrystalline structures attain stability through the large skewness of their density distributions. Importantly, their extremes can be more lopsided than those of the hexagonal phase, which also takes advantage of its skewness to compete with the lamellar and uniform states. We have taken this idea to the extreme in Section 8 to design a local free energy which allows arbitrarily high-order quasicrystals to be stabilized.
Quasicrystals are stabilized by local free energies which take advantage of the unique skewed shape of quasicrystalline density distributions. Three-body interactions are responsible for this in the LP model, but any symmetry-breaking term may do the job.