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

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 quasicrystal formation and stability. Moreover, the ability to control the self-assembly of soft quasicrystals may contribute to the development of novel photonics or other applications based on self-assembled metamaterials. Here a path is followed, leading to quantitative stability predictions, that starts with a model developed two decades ago to treat the formation of multiple-scale quasiperiodic Faraday waves (standing wave patterns in vibrating fluid surfaces) and which was later mapped onto systems of soft particles, interacting 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.

1. Introduction and outline

The scope of research on quasicrystals¹ 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 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, 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 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 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 order parameter 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 order parameter field, contains an odd, cubic, power that breaks the 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 . 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 , 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) quasi­crystals.

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 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 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, 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 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 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 do­decagonal 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 entropy term in a power series to fourth order in ρ(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 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 = , 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).²

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).³

(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 quasi­crystals 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 α⁴. The parameter c, which sets the length-scale selectivity of the system, and the control parameter are then measured in units of α².

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

Figure 1 for c = 1 and q = k5 and k12: note that is positive for all k except k = 1 or q, where it is zero. Also note that the barrier between the minima when q = k12 is significantly larger than the barrier when q = k5. This disparity is important for the study of the finite-c case discussed in Section 4.

Figure 2 Local free-energy densities 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 fCG 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 = , 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 region and diverges from the BDL logarithmic free-energy density for . Finally, the solid black line shows the free-energy density used in Section 8 to stabilize 6n-fold quasicrystals with arbitrary n.

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 entropy term from 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 is shown in Fig. 2.

(iv) In Section 8 we show that an artificial yet judicious choice of can actually stabilize 6n-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_n. Each circle is expected to contain N equally separated Bragg peaks, with N = n or 2n when n is even or odd, respectively. These targeted structures are shown schematically in Figs. 3(i)–3(o) for k₄ = 2^1/2, k₆ = 3^1/2, = 2, k₅ = (1 + 5^1/2)/2, k₁₂ = (2 + 3^1/2)^1/2, k₈ = (2 + 2^1/2)^1/2 and k₁₀ = [(5 + 5^1/2)/2]^1/2, respectively.

Figure 3 Fourier spectra of the candidate structures: structures (a)–(g) use the arbitrary ratio q = k5, 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 k5 decagonal structure (l).

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

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₆ and k₁₂ cases where the lower bound becomes 0. Recall that ∊ here is in units of α², as in equation (7); without scaling one should multiply the quoted values by α². The choices of q = k₈ and k₁₀ 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.

	Structure	q	Stable	Metastable	Figures
	Uniform				3(a)
One-scale periodic	Hexagonal				3(d)–3(g), 4, 15, 22
	Lamellar				3(b), 3(c), 4, 14, 15
	Square		Unstable	Unstable	3(h)
Two-scale periodic	Square	k4 = 2^1/2 ≃ 1.414			3(i), 5(a)
	Hexagonal	k6 = 3^1/2 ≃ 1.732			3(j), 5(b)
	Lamellar	= 2			3(k), 5(c)
Two-scale quasicrystal	Decagonal	k5 = (5^1/2 + 1)/2 ≃ 1.618			1, 3(l), 6(a), 16, 21
	Dodecagonal	k12 = (2 + 3^1/2)^1/2 ≃ 1.932			1, 3(m), 6(b), 16
		k8, k10 etc.	Unstable	See caption	3(n), 3(o)

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.⁶ 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 meta­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.

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.

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.

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₄, k₆ and , respectively. Their Fourier spectra are shown in Figs. 3(i)–3(k).⁷ 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.

Figure 5 Free energies of the two-scale periodic structures in the infinite-c LP model: in panel (a), with q = k4, the square superstructure dominates when −0.09205 0.7689. In panel (b), with q = k6, the hexagonal superstructure dominates when −0.1143 2.074. Finally, in panel (c), with q = , 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.

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₅, k₁₂, k₈ and k₁₀, whose Fourier spectra are shown in Figs. 3(l)–3(o), respectively. We find that the k₈ and k₁₀ 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₅ decagonal quasicrystal and the k₁₂ 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₅ rather than k₁₀ and miscalculated the stability boundaries of the dodecagonal one.⁸

Figure 6 Free energies of the quasiperiodic structures in the infinite-c LP model: in panel (a), with q = k5, the decagonal structure dominates when −8/93 0.02290. In panel (b), with q = k12, the dodecagonal structure dominates when −128/1269 0.03055. The gray region is the same as in Fig. 4.

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₅ 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⁹ ×73/(3³ ×17⁴).

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

and the free energy at the transition is approximately −3.694 × 10⁻³.

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⁻³.

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₁₂, 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 do­decagonal 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 quasicrystal 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 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 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.

Figure 7 Phase diagram of the LP model with q = k5 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 = k12 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 uniform–hexagonal–dodecagonal triple point and the lamellar–hexagonal co­existence curve exhibiting a turning point at intermediate c with a minimum value of ∊.

The portion of the hexagonal to quasicrystal phase boundary, calculated by Jiang et al. (2015) using the more accurate projection method,⁹ 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₅) 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₁₂) 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₁₂ but only qualitatively valid when q = k₅.

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, 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 in equation (11) from {1, q} to {q₁, 1, q₂, q₃}. We consider two cases: (i) octagonal quasicrystals, with q₁ = k_8/3 = 2cos(3π/8) = (2 − 2^1/2)^1/2 ≃ 0.7654, q₂ = k₄ and q₃ = k₈ = (2 + 2^1/2)^1/2 ≃ 1.848; and (ii) octadecagonal quasicrystals, with q₁ = k_18/7 ≃ 0.6840, q₂ = k_18/5 ≃ 1.286 and q₃ = k₁₈ ≃ 1.970. The anticipated diffraction spectra of these two structures are shown in Fig. 9.

Figure 9 Fourier spectra of the four-scale quasicrystals|: (a) for the octagonal quasicrystal, the radii of the circles from inside to outside are k8/3, 1, k4 and k8; (b) for the octadecagonal quasicrystal, the radii are k18/7, 1, k18/5 and k18.

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₂ = k₄.

For the octadecagonal case, additional competing structures stem from the fact that q₁ + q₂ = q₃. 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 quasicrystal for the relevant values of .

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

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 quasicrystal 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⁻³. A finite section of this quasicrystal is shown in Fig. 12.

Figure 11 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.

Figure 12 Predicted four-scale octagonal quasicrystal: blue and red shades correspond to negative field values −0.2306 and positive values 1.117, respectively. This quasicrystal has ∊ = −0.09. At this ∊, the minimization of the quartic energy in equation (35a) gives ≃ 0.05427 and = = ≃ 0.02856.

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 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⁻⁴. A finite section of this quasicrystal is shown in Fig. 13.

Figure 13 Predicted four-scale octadecagonal quasicrystal: blue and red shades correspond to negative field values −0.1615 and positive values 0.6008, respectively. This quasicrystal has ∊ = −0.06. At this ∊, the minimization of the quartic energy in equation (35b) gives ≃ 0.02960, ≃ 0.01234, ≃ 0.01243 and ≃ 0.01246.

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 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 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.¹⁰ Consequently, the formal expression for the calculation of the density distribution in one dimension is similar to that of the density of states¹¹

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

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) = 2cosx. 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.

Figure 14 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 right-hand 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 = 2.

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.

Figure 15 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 = −2 and = 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 = −3 and = 6, respectively. The skewness of the latter two distributions is given by γLAM = 1 and γ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.

Figure 16 Density distributions for decagonal and dodecagonal quasicrystals: for these quasicrystals, γ is maximized when = = 1/5 or 1/8, respectively. Observe that γDEC = 4 and γDOD = 3, and note the interior Van Hove singularities at ±4/5 and −1/2, −3/8 and 0, and the zeroth-order discontinuity in the decagonal distribution at = −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.

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

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 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 that has this minimum of −1 at its origin is given by {χ_i}_min = {sec⁻¹(−4), 0, 2 sec⁻¹(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 entropy 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 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 of must be negative for this temperature to be positive.¹³ Finally, the value of the field ρ(r) is here constrained to by the fact that the number density 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₄ 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 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 . 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₄ 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 (41a) and (41b) and the logarithmic local free-energy density of equation (49) into the free-energy equation (36) gives

and

where ₂F₁ is a hypergeometric function.¹⁴

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 phase transition to the uniform phase occurs.

Figure 17 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.

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₄, 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₄, k₆ 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.

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 k4, k6, k8, k10 and , as labeled. As in the LP model, the free energies of the two-scale square (k4), hexagonal (k6), lamellar (q = ), decagonal (k5) and dodecagonal (k12) 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 (k8) and decagonal (k10) quasicrystals are always higher than that of the single-scale hexagonal structure, and are therefore never the equilibrium phase.

Again, as in the LP model, setting q to k₈ and k₁₀ 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₅ and k₁₂ 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.

Figure 19 Free energies of the stable two-scale quasiperiodic structures in the BDL model: the decagonal structure (k5) is the equilibrium phase for T 1.125 and the dodecagonal structure (k12) 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 and reach a stationary pair of values.

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 inter­section 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.

Figure 20 Feasible 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 = . Solid dark lines in the first quadrant show the path of the minimizing values of 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.

For the decagonal phase with q = k₅, 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 = ≃ 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.

Figure 21 Predicted decagonal quasicrystal with = = 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.

For the dodecagonal phase with q = k₁₂, 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⁻³ 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 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 three-body (or more generally odd-body) interactions that break the symmetry of the free energy , namely, the symmetry of . 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 symmetry of . The decagonal structure with = 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 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₅, 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 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 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 quasicrystal 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.¹⁵

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 = 12n at the origin. We now show that −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 , 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 cosx + cosy + cos(x + y). This function has a minimum of −3/2 when both x y ±2π/3. However, this bound is unattainable for all 2n 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 −6n, 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 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 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.

Figure 22 Density distribution of the two-scale octadecagonal quasicrystal: the black line is the local free-energy density from equation (57). Here q = k18 ≃ 1.970 and = = 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 ≃ −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.

For the 6n-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 6n-fold quasicrystal is the minimum free-energy state of the system. Indeed, the calculated free energies of the octadecagonal and icositetragonal quasi­crystals are approximately −1.223 × 10⁻³ and −6.093 × 10⁻⁵, 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 6n-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 inter­actions. 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.

Footnotes

¹For definitions, basic terminology and some background on the symmetry-breaking transition from a liquid phase to a quasicrystal see, for example, Lifshitz (2003, 2007, 2011).

²For a comparison of some of the different approaches see, for example, van Teeffelen et al. (2009).

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

⁴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).

⁵Jiang et al. (2015) termed these `sibling periodic crystals'.

⁶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).

⁷The single-scale lamellar, square and hexagonal structures, and the two-scale square and hexagonal superstructures, correspond to the A₁, A₁ × A₁, A₂, B₂ and G₂ Lie algebra root systems, respectively.

⁸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)].

⁹We thank the authors for graciously sharing their raw data with us.

¹⁰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).

¹¹See, for example, equation (8.63) of Ashcroft & Mermin (1976).

¹²Equivalently, one can think of this process as uniformly sampling a unit cell of the higher-dimensional periodic structure through which the quasicrystal can be constructed as a slice at an irrational slope (Senechal, 1995).

¹³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 molecular dynamics simulations, denote it as T_sp but leave the temperature T unscaled.

¹⁴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.

¹⁵We note that, for 6n 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.