- 1. Introduction and outline
- 2. The Lifshitz–Petrich model and its immediate generalizations
- 3. Stable periodic and quasiperiodic crystals in the original LP model with exact length-scale selection
- 4. Relaxing the requirement of exact length-scale selection
- 5. Four-scale octagonal and octadecagonal quasicrystals
- 6. The density distribution method
- 7. Mean-field theory for soft interacting particles
- 8. Higher-order quasicrystals
- 9. Closing remarks
- References

- 1. Introduction and outline
- 2. The Lifshitz–Petrich model and its immediate generalizations
- 3. Stable periodic and quasiperiodic crystals in the original LP model with exact length-scale selection
- 4. Relaxing the requirement of exact length-scale selection
- 5. Four-scale octagonal and octadecagonal quasicrystals
- 6. The density distribution method
- 7. Mean-field theory for soft interacting particles
- 8. Higher-order quasicrystals
- 9. Closing remarks
- References

## topical reviews

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

^{a}Condensed Matter Physics, California Institute of Technology, Pasadena, CA 91125, USA, ^{b}Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA, and ^{c}Raymond and Beverly Sackler School of Physics and Astronomy, Tel Aviv University, Tel Aviv 69978, Israel^{*}Correspondence e-mail: sam@savitz.org, ronlif@tau.ac.il

For many years, quasicrystals were observed only as solid-state metallic alloys, yet current research is now actively exploring their formation in a variety of soft materials, including systems of macromolecules, nanoparticles and colloids. Much effort is being invested in understanding the thermodynamic properties of these soft-matter quasicrystals in order to predict and possibly control the structures that form, and hopefully to shed light on the broader yet unresolved general questions of *via* multiple-scale pair potentials. The article reviews, and substantially expands, the quantitative predictions of these models, while correcting a few discrepancies in earlier calculations, and presents new analytical methods for treating the models. In so doing, a number of new stable quasicrystalline structures are found with octagonal, octadecagonal and higher-order symmetries, some of which may, it is hoped, be observed in future experiments.

Keywords: quasicrystals; soft matter; pattern formation.

### 1. 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 *solid-state* quasicrystals (Tsai, 2003, 2008; Janssen *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 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, 2011, 2012; Fischer *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 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 field – which is often a simple scalar function ρ(**r**) that describes the relative deviation of a coarse-grained density *c*(**r**) of a material from its average value – and uses generic symmetry arguments to formulate a free-energy functional , expressed in powers of the and its gradients. One assumes that the equilibrium phase is the one that minimizes , and then all that remains is to find out which structures could minimize such a free energy, or conversely, how to tweak 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 , 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 *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 . In Section 3 we repeat and extend the calculations of LP of the free energies of candidate structures, setting

*q*=

*k*for different values of

_{n}*n*and assuming LP's limit of , 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 = , 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.*, 2015, 2016, 2017; Subramanian *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 limit, using a high-dimensional 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 length-scale 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 ), 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 for the given thermodynamic parameters – such as temperature *T* and μ – and is then left with the question of how to design the pair potential 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 many-body interactions. That stable quasicrystals may require the existence of two length scales in their effective interaction potentials 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, 2017*a*; 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 terms of to the nonlocal gradient expansion and local power expansion terms of , 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 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 , imitating the gradient term of , provides greater control over the self-assembly 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 turned out to be a challenge. Instead, BDL expanded the logarithmic ρ(**r**) = 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 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 6*n*-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 discovered by Fischer *et al.* (2011), those found numerically (Dotera *et al.*, 2014; Engel & Glotzer, 2014; Pattabhiraman & Dijkstra, 2017*b*) 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 = , 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.*, 2014, 2015; Subramanian *et al.*, 2016).^{2}

Interestingly, the PFC model of Achim *et al.* (2014) uses 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 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 (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 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 soft-matter and solid-state systems continue to emerge (Lee *et al.*, 2014; Lifshitz, 2014), possibly enlarging our scope.

### 2. 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 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 . 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

where in the function is given by the octic polynomial,

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 . 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, = 0 if *k* belongs to the set = {1, *q*} and is otherwise infinite. Thus, we immediately conclude that

restricting the support of 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 field – which 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 . 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 two-scale 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 , or both, in a number of ways:

(i) While remaining in the infinite-*c* limit, can be changed by modifying the set . We show in Section 5 that doubling the cardinality of , *i.e.* going from two to four concentric circles, allows us to stabilize octagonal and octadecagonal quasicrystals.

(ii) Subramanian *et al.* (2016) modify 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) , 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 of a pair of particles in real space, possibly scaled and shifted by a constant. This, along with a replacement of by the local , 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 is shown in Fig. 2.

term from in equation (3)(iv) In Section 8 we show that an artificial yet judicious choice of can actually stabilize 6*n*-fold quasicrystals for any , with just two length scales.

### 3. Stable periodic and quasiperiodic crystals in the original LP model with exact length-scale selection

#### 3.1. Notation and method of calculation

##### 3.1.1. Calculation of stability bounds

We set *q* equal to *k _{n}* from equation (2) with

*n*3, so that 1

*q*2, and the upper limit of 2 is obtained for . Our goal is to stabilize

*N*-fold symmetric structures whose Fourier coefficients are confined to two circles of radii unity and

*k*. Each circle is expected to contain

_{n}*N*equally separated Bragg peaks, with

*N*=

*n*or 2

*n*when

*n*is even or odd, respectively. These targeted structures are shown schematically in Figs. 3(

*i*)–3(

*o*) for

*k*

_{4}= 2

^{1/2},

*k*

_{6}= 3

^{1/2}, = 2,

*k*

_{5}= (1 + 5

^{1/2})/2,

*k*

_{12}= (2 + 3

^{1/2})

^{1/2},

*k*

_{8}= (2 + 2

^{1/2})

^{1/2}and

*k*

_{10}= [(5 + 5

^{1/2})/2]

^{1/2}, respectively.

These two-scale structures are in thermodynamic competition with the uniform liquid phase ρ(**r**) = 0, and with single-scale 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 and 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}* = , 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 to or and . 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 free-energy 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 metastability 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 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.

#### 3.2. 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 .

For the lamellar phase, the free-energy equation (15) gives

Minimizing over shows that, for , = 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 , below which the nontrivial minima of equation (22) are complex, and so the only possible 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, = −[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 , respectively. Their Fourier spectra are shown in Figs. 3(*i*)–3(*k*).^{7} For these structures there are no simple expressions for the minimized and 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 , and 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 and the *k*_{12} dodecagonal 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 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

It exists only when −3/31 ≃ −0.09677. For −3/31 0.5245, = and

Above this approximate upper bound, for which there is no simple expression, the free energy continues to decrease, but no longer equals . The free energy at the transition is approximately −0.04889.

The free energy of the dodecagonal phase is given by

It exists only for −16/141 ≃ −0.1135, where

for , and

for . In the first case, , but in the second, . The free energy at the transition between these two free-energy minima is -2^{9} ×73/(3^{3} ×17^{4}).

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, symmetry associated with the exchange of and . In both cases, when minimizing the free energy with respect to these amplitudes there is one solution branch that maintains this symmetry with = and a second branch where the symmetry is spontaneously broken and . Yet, as it turns out, in the decagonal this is first order, while in the dodecagonal case it is a continuous phase transition.

We compare the free energies of these *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

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 = is metastable when −0.09677 ≃ −3/31 763/972 ≃ 0.7850. The decagonal phase with appears to be metastable for all −0.01493.

The dodecagonal phase with = is metastable for −0.1135 ≃ −16/141 208/867 ≃ 0.2399. The phase with 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.

### 4. 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 , 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 to lie within two rings centered about the origin. This is in contrast with allowing to vary freely, as most numerical simulations do (Lifshitz & Petrich, 1997; Barkan *et al.*, 2011), or restricting to some subset of a two-dimensional 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 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 finite-width rings near the minima of 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 ,

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 , 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.

#### 4.2. *c*-Dependent phase diagrams for decagonal and dodecagonal quasicrystals

The *c*-dependent phase diagrams, calculated using the two-ring 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 stability. As *c* is decreased, the competing phases change from a solution where ν = 1 and = 0 to one where ν is shifted and both and are nonzero, in order to take advantage of both minima of in an optimal way. This causes the upper bound of for stability to constrict with decreasing *c* until the quasicrystalline phase vanishes. This vanishing occurs at a uniform–hexagonal–quasicrystal Below the the critical for the uniform–hexagonal transition continues to decrease towards the value ∼ −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 *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 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 free-energy barrier between the two minima of 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 and the precise shapes of the phase-boundary curves, being quantitatively reasonable when *q* = *k*_{12} but only qualitatively valid when *q* = *k*_{5}.

### 5. Four-scale octagonal and octadecagonal quasicrystals

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, 2017*b*). 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 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.

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 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 = and = . Minimizing these equations in the relevant range shows that the coincidental lamellar phase has = = 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

for the relevant values of .Applying equation (15) to the octagonal structure in Fig. 9(*a*) gives a free energy of

where = + . 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 and symmetry, the minima in the relevant small range all satisfy = = .

As shown in Fig. 11(*a*), the octagonal is expected to be thermodynamically stable when −0.1119 −0.06227. In this range, the optimized is larger than the equal '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 is shown in Fig. 12.

| Figure 12 has |

The free energy of the octadecagonal structure is

As shown in Fig. 11(*b*), the octadecagonal is expected to be stable when −0.06882 −0.05140. In this range, the optimized is larger than the 's by a factor varying between roughly 2.2 and 2.6. The '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 is shown in Fig. 13.

| Figure 13 has |

### 6. The density distribution method

#### 6.1. The notion of a density distribution and its quantum analogy

When the local free-energy function is a finite polynomial such as equation (7) and the structure ρ(**r**) consists of a finite number of harmonic components, the free energy of a 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 which stabilizes arbitrarily high-order 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

and δ is the Dirac delta function. 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 uniform-weight inner product of *f* and *P* over the 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 notice – using standard complex analysis – that the density distribution , where

It turns out that this function is the on-site lattice Green's function, obtained for the nearest-neighbor tight-binding 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 The real part can be useful for evaluating the density distribution-based free-energy equation (36) for certain local free-energy density functions , 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).

#### 6.2. 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 = . 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 , and rescale later if necessary.

For two-scale structures, characterized by two Fourier amplitudes and , a rescaling of ρ(**r**) affects both amplitudes together, giving = . Thus, the density distributions differ for fields with different ratios of the amplitudes, but are otherwise independent of the overall scale of the field.

For all the structures relevant to us, has compact support between the extreme values of ρ(**r**), which are both finite, because the fields are all finite sums of harmonic functions. The value γ = is a measure of the `skewness' or `lopsidedness' of the density distribution. It characterizes the imbalance between the ground state and the highest 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 negative . We need only consider positive . On the other hand, with this assumption, γ for two-scale structures varies continuously as a function of the ratio of the two Fourier amplitudes.

#### 6.3. Numerical sampling of the field

One may need to resort to numerical sampling of the field ρ(**r**) in order to generate the density distribution when analytical methods for calculating equations (37) or (39) prove difficult. For periodic crystals this is readily achieved by uniformly sampling the 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

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, = .

We therefore begin with the single-scale lamellar field which, after scaling to unity, is given by ρ(**r**) = 2cos*x*. Because never exceeds two, vanishes when 2. We express it analytically between these bounds using equation (39), by substituting for and 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

and

We call the value of ρ(**r** = 0; {χ_{i}}) when = 1 and = 0 and likewise call the value when = 0 and = 1, so that = + . We numerically sample the ordered pairs uniformly over , with 96 sampling points along each phase, to give a total of roughly 85 million samples. If we write this distribution function as , then

Upon numerically maximizing the skewness parameter γ over the ratio of amplitudes for the decagonal and dodecagonal phases, we find that the optimal ratio is unity, where = 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*(0). The lamellar phase can potentially dominate this by heavily sampling the local free-energy density far from = 0. For example, a local free-energy 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.

#### 6.5. 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 (42*a*) and (42*b*)],

we can identify Van Hove singularities at and , at ±4/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 = −1, all of the and can be understood in terms of the spatial Hessian of the effective four-dimensional function. One of the phase-shifted structures 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

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 two-dimensional. 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 zeroth-order 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.

#### 6.6. 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**) = that minimizes the local free-energy density (7) everywhere. To satisfy the zero-average requirement we must add a compensating delta function peak at some negative value , and possibly allow the positive value to shift away from . This yields a field with sharp boundaries between two allowed values and a density distribution of the form = + .

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 length-scale 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, + = 0, and the normalization *P*_{l} + *P*_{r} = 1 of the density distribution. Solving for the left-hand side variables gives *P*_{l} = 1 - *P*_{r} and = . Substituting them into and minimizing the free energy (36) over *P*_{r} and gives

and

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 . 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 term. Assuming a sufficiently dense system of soft particles that interact *via* a Fourier transformable isotropic pair potential – 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

and

where is the Hankel transform of and is its minimum value. The temperature *T* is measured here in units of the spinodal temperature *T*_{sp} = , where is the average of the particles and *k*_{B} is the 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 of must be negative for this temperature to be positive.^{13} Finally, the value of the field ρ(**r**) is here constrained to by the fact that the *c*(**r**) = 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 , 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* = . Note that the range −4/3 0 corresponds to scaled temperature values of *T* 1 above the 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 . 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 and therefore small *T*_{sp} limit. As seen below, we indeed find important differences between the behaviors of the LP and BDL models.

#### 7.2. Single-scale structures

For the single-scale lamellar and hexagonal phases, substituting the density distributions of equations (41*a*) and (41*b*) and the logarithmic local free-energy density of equation (49) into the free-energy equation (36) gives

_{2}

*F*

_{1}is a hypergeometric function.

^{14}

Now, minimizing over yields

which is shown in Fig. 17. At absolute zero, the free energy is exactly −1/4. In the first temperature range, where 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 ln2)/4. From here, as *T* increases to unity, decreases to zero. At *T* = 1, a second-order to the uniform phase occurs.

The result of minimizing , 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, ≃ −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.

#### 7.3. Two-scale structures

In the original LP model, one can set *q* to *k*_{4}, *k*_{6} and 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 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 and is performed subject to the vacuum constraint, which is more difficult to take into account for the two-scale structures. Each pair from Section 6.4.2 not only gives a small free-energy contribution, but also imposes the linear constraint + −1 on the values of and . The problem of finding the intersection of these half-planes is dual to that of finding the convex hull of the points (de Berg *et al.*, 2008). If and are adjacent extremal points on the convex hull, then the point

is a vertex of the polygonal boundary of the set of allowed 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 from the uniform liquid to the at *T* ≃ 1.125, at which point = ≃ 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, where = ≃ 0.1220. At *T* ≃ 1.152, the free energy ≃ −1.123 × 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 symmetry is broken at *T* ≃ 0.8446 and ≃ −0.05086. After this point, the 's move along their maximum allowed sum line + = 1/4 until *T* ≃ 0.7022, where they land in either of the two degenerate states = (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.

#### 7.4. Skewness and the vacuum constraint

While the dodecagonal *T* = 0, where there are no three-body interactions, demonstrates that this cannot be the whole story.

The primary γ skewness of four, and allows the decagonal phase to take maximal advantage of the the highly negative 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 simulations performed by Barkan *et al.* (2014) with *q* = *k*_{5}, showing that the decagonal 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 6*n*-fold two-scale structures

Finally, we examine the γ skewness of the density distributions of quasicrystals of order 6*n* for arbitrary *n* 2. Then, using the information obtained, we judiciously design an artificial local free-energy function 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 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 6*n*-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 and are equal. We restrict our attention to this case and scale both and to unity.

By inverse Fourier transforming its momentum-space representation according to equation (8), the position-space field representing this can be written as

where the wavevectors

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 accordingly. With this choice, with all the χ's set to zero, the field obtains its maximum value = 12*n* at the origin. We now show that −6*n* 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 , 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* ±2π/3. However, this bound is unattainable for all 2*n* cosine triplets: If we set χ_{m} χ_{m+2n} ±2π/3 for *m* = , 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 −6*n*, and so γ 2.

Indeed, numerical sampling for = shows that γ = 3, 3, 36/13 ≃ 2.769 and ∼2.634, for *n* = , respectively. Asymptotically, γ appears to decrease no faster than 2 + 2/*n*.

#### 8.2. A local free-energy density that stabilizes 6*n*-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 region without suffering greater free-energy penalties from the positive region, as can be inferred graphically from Fig. 22. Furthermore, because the free-energy penalty is twice the free-energy decrease, the hexagonal phase is unable to reach negative free energies through negative values of either.

For the 6*n*-fold two-scale structures considered above, we scale = until reaches −1. Then, because γ 2, is also greater than two. This, together with the density distribution (36) and the local free-energy function (57), implies that the free energy is negative for this structure. Thus, a 6*n*-fold is the minimum free-energy state 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 = −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 6*n*-fold 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).

### 9. 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 free-energy 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 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.

### Footnotes

^{1}For definitions, basic terminology and some background on the from a liquid phase to a see, for example, Lifshitz (2003, 2007, 2011).

^{2}For a comparison of some of the different approaches see, for example, van Teeffelen *et al.* (2009).

^{3}Additionally, we have found that, for the single-length-scale case in arbitrarily high dimensions, the set of potential minimal free-energy candidate phases corresponds to the ADE Lie algebra root lattices. The minimal free-energy states are precisely the laminated lattices in no more than eight dimensions.

^{4}The particular form of the function , with its octic ultraviolet divergence, comes from duplicating the gradient term of the Swift–Hohenberg equation (5), which possesses a quartic divergence, and whose purpose is to favor modes with wavenumbers approximately equal to unity or *q*. Conveniently, after taking the infinite-*c* limit, as we do everywhere outside of Section 4, it diverges almost everywhere, and its only remaining features are its minima at *k* = 1 and *k* = *q* where = 0. On the other hand, when viewed as the Fourier transform of a real-space interaction potential, it is incompatible with the kind of soft-matter particles we have in mind, as it introduces a strongly diverging interaction at the origin. One may, technically, tame this interaction by multiplying it with a Gaussian, as was done by Barkan *et al.* (2014).

^{5}Jiang *et al.* (2015) termed these `sibling periodic crystals'.

^{6}More technically, there is always a Rokhsar–Wright–Mermin gauge transformation (Rokhsar *et al.*, 1988) that can be applied to the free-energy minimizing structure to obtain these phase values, without altering the free energy. For specific information regarding the required gauge transformation, and a definition of a screw operation in the case of quasicrystals, see, for example, Rabson *et al.* (1991).

^{7}The single-scale lamellar, square and hexagonal structures, and the two-scale square and hexagonal superstructures, correspond to the *A*_{1}, *A*_{1} × *A*_{1}, *A*_{2}, *B*_{2} and *G*_{2} Lie algebra root systems, respectively.

^{8}These discrepancies led to some confusion in the subsequent literature [see, *e.g.*, p. 3 and footnote 32 of Barkan *et al.* (2011), p. 2 and footnote 26 of Barkan *et al.* (2014), and pp. 7 and 9 of Jiang *et al.* (2015)].

^{9}We thank the authors for graciously sharing their raw data with us.

^{10}This analogy is reminiscent of that between the shape of micelles in real space and electronic Fermi surfaces in momentum space, suggested by Lee *et al.* (2014) and discussed by Lifshitz (2014).

^{11}See, for example, equation (8.63) of Ashcroft & Mermin (1976).

^{12}Equivalently, one can think of this process as uniformly sampling a of the higher-dimensional periodic structure through which the can be constructed as a slice at an irrational slope (Senechal, 1995).

^{13}It must also be noted that BDL denote the spinodal temperature as *T*_{c}, while Barkan *et al.* (2014), who confirmed the BDL results using simulations, denote it as *T*_{sp} but leave the temperature *T* unscaled.

^{14}Alternatively, one could evaluate these expressions using a double integral of an expression involving the analytically continued real part of the lattice Green's function (40) and no logarithms, instead of a single integral with a logarithm.

^{15}We note that, for 6*n* 46, there are additional inequivalent arrangements of wavevectors to consider, corresponding to distinct Bravais lattices, which we ignore here. See Mermin *et al.* (1987) for additional information.

### Acknowledgements

S. Savitz thanks Gil Refael, the Institute of Quantum Information and Matter, the Caltech Student–Faculty Programs office and Marcella Bonsall for their support. R. Lifshitz thanks Dean Petrich, Gilad Barak, Kobi Barkan, Yoni Mayzel, Michael Engel, Haim Diamant and Mike Cross for their fruitful collaboration on these problems over the years. We again extend our gratitude to Jiang *et al.* (2015) for sharing their data with us. Thanks to Marcus Bintz for pointing out the connections to the Kramers–Kronig relation and Morse theory in Section 6. The calculations in Sections 6 and 7 utilized the open-source computational geometry software library CGAL (The CGAL Project; https://www.cgal.org/).

### Funding information

Funding for this research was provided by: Israel Science Foundation (grant No. 1667/16).

### References

Achim, C. V., Schmiedeberg, M. & Löwen, H. (2014). *Phys. Rev. Lett.* **112**, 255501. Web of Science CrossRef Google Scholar

Alexander, S. & McTague, J. (1978). *Phys. Rev. Lett.* **41**, 702–705. CrossRef CAS Web of Science Google Scholar

Arbell, H. & Fineberg, J. (2000). *Phys. Rev. Lett.* **85**, 756–759. Web of Science CrossRef CAS Google Scholar

Arbell, H. & Fineberg, J. (2002). *Phys. Rev. E*, **65**, 036224. Web of Science CrossRef Google Scholar

Archer, A. J., Rucklidge, A. M. & Knobloch, E. (2013). *Phys. Rev. Lett.* **111**, 165501. Web of Science CrossRef Google Scholar

Archer, A. J., Rucklidge, A. M. & Knobloch, E. (2015). *Phys. Rev. E*, **92**, 012324. Web of Science CrossRef Google Scholar

Ashcroft, N. W. & Mermin, N. D. (1976). *Solid State Physics.* New York: Holt, Rinehart and Winston. Google Scholar

Bak, P. (1985). *Phys. Rev. Lett.* **54**, 1517–1519. CrossRef CAS Web of Science Google Scholar

Barkan, K. (2015). *Theory and Simulation of the Self Assembly of Soft Quasicrystals.* PhD thesis, Tel Aviv University, Israel. Google Scholar

Barkan, K., Diamant, H. & Lifshitz, R. (2011). *Phys. Rev. B*, **83**, 172201. Web of Science CrossRef Google Scholar

Barkan, K., Engel, M. & Lifshitz, R. (2014). *Phys. Rev. Lett.* **113**, 098304. Web of Science CrossRef PubMed Google Scholar

Berg, M. de, Cheong, O., van Kreveld, M. & Overmars, M. (2008). *Comput. Geom.* 3rd ed., ch. 4, pp. 63–93. Heidelberg: Springer. Google Scholar

Bodnarchuk, M., Erni, R., Krumeich, F. & Kovalenko, M. (2013). *Nano Lett.* **13**, 1699–1705. Web of Science CrossRef CAS Google Scholar

Chaikin, P. M. & Lubensky, T. C. (1995). *Principles of Condensed Matter Physics.* Cambridge University Press. Google Scholar

Chanpuriya, S., Kim, K., Zhang, J., Lee, S., Arora, A., Dorfman, K. D., Delaney, K. T., Fredrickson, G. H. & Bates, F. S. (2016). *ACS Nano*, **10**, 4961–4972. Web of Science CrossRef CAS Google Scholar

Cross, M. & Greenside, H. (2009). *Pattern Formation and Dynamics in Nonequilibrium Systems*, ch. 5. Cambridge University Press. Google Scholar

Damasceno, P. F., Glotzer, S. C. & Engel, M. (2017). *J. Phys. Condens. Matter*, **29**, 234005. Web of Science CrossRef Google Scholar

Dotera, T. (2007). *Philos. Mag.* **87**, 3011–3019. Web of Science CrossRef CAS Google Scholar

Dotera, T. (2011). *Isr. J. Chem.* **51**, 1197–1205. Web of Science CrossRef CAS Google Scholar

Dotera, T. (2012). *J. Polym. Sci. B Polym. Phys.* **50**, 155–167. Web of Science CrossRef CAS Google Scholar

Dotera, T., Oshiro, T. & Ziherl, P. (2014). *Nature*, **506**, 208–211. Web of Science CrossRef CAS Google Scholar

Dzugutov, M. (1993). *Phys. Rev. Lett.* **70**, 2924–2927. CrossRef CAS Web of Science Google Scholar

Edwards, W. & Fauve, S. (1993). *Phys. Rev. E*, **47**, R788–R791. CrossRef CAS Web of Science Google Scholar

Engel, M., Damasceno, P. F., Phillips, C. L. & Glotzer, S. C. (2015). *Nat. Mater.* **14**, 109–116. Web of Science CrossRef CAS Google Scholar

Engel, M. & Glotzer, S. (2014). *Nat. Phys.* **10**, 185–186. Web of Science CrossRef CAS Google Scholar

Engel, M. & Trebin, H.-R. (2007). *Phys. Rev. Lett.* **98**, 225505. Web of Science CrossRef Google Scholar

Fischer, S., Exner, A., Zielske, K., Perlich, J., Deloudi, S., Steurer, W., Lindner, P. & Förster, S. (2011). *Proc. Natl Acad. Sci. USA*, **108**, 1810–1814. Web of Science CrossRef CAS PubMed Google Scholar

Fredrickson, G. H. (2006). *The Equilibrium Theory of Inhomogeneous Polymers.* Oxford University Press. Google Scholar

Gennes, P. G. de & Prost, J. (1993). *The Physics of Liquid Crystals*, 2nd ed. New York: Oxford University Press. Google Scholar

Gollub, J. P. & Langer, J. S. (1999). *Rev. Mod. Phys.* **71**, S396–S403. Web of Science CrossRef Google Scholar

Gompper, G. & Schick, M. (1994). *Self-Assembling Amphiphilic Systems.* Vol. 16 of *Phase Transitions and Critical Phenomena.* London: Academic Press. Google Scholar

Gronlund, L. & Mermin, N. D. (1988). *Phys. Rev. B*, **38**, 3699–3710. CrossRef CAS Web of Science Google Scholar

Guttmann, A. J. (2010). *J. Phys. A Math. Theor.* **43**, 305205. Web of Science CrossRef Google Scholar

Hassani, S., Koutschan, C., Maillard, J.-M. & Zenine, N. (2016). *J. Phys. A Math. Theor.* **49**, 164003. Web of Science CrossRef Google Scholar

Hayashida, K., Dotera, T., Takano, A. & Matsushita, Y. (2007). *Phys. Rev. Lett.* **98**, 195502. Web of Science CrossRef PubMed Google Scholar

Hohenberg, P. C. & Halperin, B. I. (1977). *Rev. Mod. Phys.* **49**, 435–479. CrossRef CAS Web of Science Google Scholar

Horiguchi, T. (1972). *J. Math. Phys.* **13**, 1411–1419. CrossRef Web of Science Google Scholar

Jagla, E. A. (1998). *Phys. Rev. E*, **58**, 1478–1486. Web of Science CrossRef CAS Google Scholar

Janssen, T., Chapuis, G. & de Boissieu, M. (2007). *Aperiodic Crystals.* Oxford University Press. Google Scholar

Jarić, M. V. (1985). *Phys. Rev. Lett.* **55**, 607–610. CrossRef PubMed CAS Web of Science Google Scholar

Jiang, K., Tong, J. & Zhang, P. (2016). *Commun. Comput. Phys.* **19**, 559–581. Web of Science CrossRef Google Scholar

Jiang, K., Tong, J., Zhang, P. & Shi, A.-C. (2015). *Phys. Rev. E*, **92**, 042159. Web of Science CrossRef Google Scholar

Jiang, K. & Zhang, P. (2014). *J. Comput. Phys.* **256**, 428–440. Web of Science CrossRef Google Scholar

Jiang, K., Zhang, P. & Shi, A.-C. (2017). *J. Phys. Condens. Matter*, **29**, 124003. Web of Science CrossRef Google Scholar

Kalugin, P. A., Kitaev, A. Y. & Levitov, L. C. (1985). *JETP Lett.* **41**, 145. Google Scholar

Keys, A. S. & Glotzer, S. C. (2007). *Phys. Rev. Lett.* **99**, 235503. Web of Science CrossRef Google Scholar

Koutschan, C. (2013). *J. Phys. A Math. Theor.* **46**, 125005. Web of Science CrossRef Google Scholar

Kudrolli, A., Pier, B. & Gollub, J. (1998). *Physica D*, **123**, 99–111. Web of Science CrossRef CAS Google Scholar

Lee, S., Leighton, C. & Bates, F. S. (2014). *Proc. Natl Acad. Sci. USA*, **111**, 17723–17731. Web of Science CrossRef CAS Google Scholar

Levitov, L. S. (1988). *Europhys. Lett.* **6**, 517–522. CrossRef CAS Web of Science Google Scholar

Lifshitz, R. (2003). *Found. Phys.* **33**, 1703–1711. Web of Science CrossRef Google Scholar

Lifshitz, R. (2007). *Z. Kristallogr.* **222**, 313–317. Web of Science CrossRef CAS Google Scholar

Lifshitz, R. (2011). *Isr. J. Chem.* **51**, 1156–1167. Web of Science CrossRef CAS Google Scholar

Lifshitz, R. (2014). *Proc. Natl Acad. Sci. USA*, **111**, 17698–17699. Web of Science CrossRef CAS Google Scholar

Lifshitz, R. & Diamant, H. (2007). *Philos. Mag.* **87**, 3021–3030. Web of Science CrossRef CAS Google Scholar

Lifshitz, R. & Petrich, D. (1997). *Phys. Rev. Lett.* **79**, 1261–1264. CrossRef CAS Web of Science Google Scholar

Loh, Y. L. (2017). *J. Phys. A Math. Theor.* **50**, 405203. Web of Science CrossRef Google Scholar

Mermin, D., Rokhsar, D. & Wright, D. (1987). *Phys. Rev. Lett.* **58**, 2099–2101. CrossRef CAS Web of Science Google Scholar

Mermin, N. D. & Troian, S. M. (1985). *Phys. Rev. Lett.* **54**, 1524–1527. CrossRef CAS Web of Science Google Scholar

Mikhael, J., Schmiedeberg, M., Rausch, S., Roth, J., Stark, H. & Bechinger, C. (2010). *Proc. Natl Acad. Sci. USA*, **107**, 7214–7218. Web of Science CrossRef CAS Google Scholar

Mkhonta, S. K., Elder, K. R. & Huang, Z.-F. (2013). *Phys. Rev. Lett.* **111**, 035501. Web of Science CrossRef Google Scholar

Müller, H. W. (1994). *Phys. Rev. E*, **49**, 1273–1277. Google Scholar

Narasimhan, S. & Ho, T. L. (1988). *Phys. Rev. B*, **37**, 800–809. CrossRef CAS Web of Science Google Scholar

Olami, Z. (1990). *Phys. Rev. Lett.* **65**, 2559–2562. CrossRef CAS Web of Science Google Scholar

Pattabhiraman, H. & Dijkstra, M. (2017*a*). *J. Phys. Condens. Matter*, **29**, 094003. Web of Science CrossRef Google Scholar

Pattabhiraman, H. & Dijkstra, M. (2017*b*). *J. Chem. Phys.* **146**, 114901. Web of Science CrossRef Google Scholar

Percec, V., Imam, M. R., Peterca, M., Wilson, D. A., Graf, R., Spiess, H. W., Balagurusamy, V. S. K. & Heiney, P. A. (2009). *J. Am. Chem. Soc.* **131**, 7662–7677. Web of Science CrossRef PubMed CAS Google Scholar

Quandt, A. & Teter, M. P. (1999). *Phys. Rev. B*, **59**, 8586–8592. Web of Science CrossRef CAS Google Scholar

Rabson, D. A., Mermin, N. D., Rokhsar, D. S. & Wright, D. C. (1991). *Rev. Mod. Phys.* **63**, 699–733. CrossRef CAS Web of Science Google Scholar

Ramakrishnan, T. & Yussouff, M. (1979). *Phys. Rev. B*, **19**, 2775–2794. CrossRef CAS Web of Science Google Scholar

Ramanujan, S. (1914). *Q. J. Math.* **45**, 350–372. Google Scholar

Ray, K. (2014). arXiv: 1409.7806. Google Scholar

Rokhsar, D. S., Wright, D. C. & Mermin, N. D. (1988). *Acta Cryst.* A**44**, 197–211. CrossRef CAS Web of Science IUCr Journals Google Scholar

Roth, J. & Denton, A. R. (2000). *Phys. Rev. E*, **61**, 6845–6857. Web of Science CrossRef CAS Google Scholar

Sachdev, S. & Nelson, D. (1985). *Phys. Rev. B*, **32**, 4592–4606. CrossRef CAS Web of Science Google Scholar

Senechal, M. (1995). *Quasicrystals and Geometry.* Cambridge University Press. Google Scholar

Skibinsky, A., Buldyrev, S. V., Scala, A., Havlin, S. & Stanley, H. E. (1999). *Phys. Rev. E*, **60**, 2664–2669. Web of Science CrossRef CAS Google Scholar

Smith, A. (1991). *Phys. Rev. B*, **43**, 11635–11641. CrossRef CAS Web of Science Google Scholar

Steurer, W. & Deloudi, S. (2009). *Crystallography of Quasicrystals: Concepts, Methods and Structures.* Heidelberg: Springer. Google Scholar

Subramanian, P., Archer, A. J., Knobloch, E. & Rucklidge, A. M. (2016). *Phys. Rev. Lett.* **117**, 075501. Web of Science CrossRef Google Scholar

Swift, J. & Hohenberg, P. (1977). *Phys. Rev. A*, **15**, 319–328. CrossRef Web of Science Google Scholar

Takano, A., Kawashima, W., Noro, A., Isono, Y., Tanaka, N., Dotera, T. & Matsushita, Y. (2005). *J. Polym. Sci. B Polym. Phys.* **43**, 2427–2432. Web of Science CrossRef CAS Google Scholar

Talapin, D. V., Shevchenko, E. V., Bodnarchuk, M. I., Ye, X., Chen, J. & Murray, C. B. (2009). *Nature*, **461**, 964–967. Web of Science CrossRef PubMed CAS Google Scholar

Teeffelen, S. van, Backofen, R., Voigt, A. & Löwen, H. (2009). *Phys. Rev. E*, **79**, 051404. Google Scholar

Tsai, A. P. (2003). *Acc. Chem. Res.* **36**, 31–38. Web of Science CrossRef PubMed CAS Google Scholar

Tsai, A. P. (2008). *Sci. Technol. Adv. Mater.* **9**, 013008. Web of Science CrossRef Google Scholar

Umbanhowar, P. B., Melo, F. & Swinney, H. L. (1996). *Nature*, **382**, 793–796. CrossRef CAS Web of Science Google Scholar

Ungar, G., Percec, V., Zeng, X. & Leowanawat, P. (2011). *Isr. J. Chem.* **51**, 1206–1215. Web of Science CrossRef CAS Google Scholar

Ungar, G. & Zeng, X. (2005). *Soft Matter*, **1**, 95–106. Web of Science CrossRef CAS Google Scholar

Van Hove, L. (1953). *Phys. Rev.* **89**, 1189–1193. CrossRef CAS Web of Science Google Scholar

Wu, K.-A., Adland, A. & Karma, A. (2010). *Phys. Rev. E*, **81**, 061601. Web of Science CrossRef Google Scholar

Xiao, C., Fujita, N., Miyasaka, K., Sakamoto, Y. & Terasaki, O. (2012). *Nature*, **487**, 349–353. Web of Science CrossRef CAS Google Scholar

Zeng, X., Ungar, G., Liu, Y., Percec, V., Dulcey, A. E. & Hobbs, J. (2004). *Nature*, **428**, 157–160. Web of Science CrossRef CAS Google Scholar

Zenine, N., Hassani, S. & Maillard, J. M. (2015). *J. Phys. A Math. Theor.* **48**, 035205. Web of Science CrossRef Google Scholar

Zhang, J. & Bates, F. (2012). *J. Am. Chem. Soc.* **134**, 7636–7639. Web of Science CrossRef CAS Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.