Quasicrystal formation in binary soft matter mixtures

Using a strategy that may be applied in theory or in experiments, we identify the regime in which a model binary soft matter mixture forms quasicrystals. The system is described using classical density functional theory combined with integral equation theory. Quasicrystal formation requires particle ordering with two characteristic lengthscales in certain particular ratios. How the lengthscales are related to the form of the pair interactions is reasonably well understood for one component systems, but less is known for mixtures. In our model mixture of big and small colloids confined to an interface, the two lengthscales stem from the range of the interactions between pairs of big particles and from the cross big-small interactions, respectively. The small-small lengthscale is not significant. Our strategy for finding quasicrystals involves tuning locations of maxima in the dispersion relation, or equivalently in the liquid state partial static structure factors.

For systems of soft particles of a single type, there is growing understanding of the ingredients required for the self-assembly of quasicrystals (QCs). The necessary features are rather special, which explains why QCs are rare in nature. These are best seen by considering the particle pair interaction potentials and pair correlation functions in Fourier space, where one observes that there are two characteristic peaks at wave numbers k 1 and k 2 , with the ratio k 1 /k 2 taking certain special values which are geometric in origin [1][2][3][4][5][6][7][8][9][10][11]; e.g. for two dimensional (2D) dodecagonal QCs, k 1 /k 2 = 2 cos(π/12) ≈ 1.93. These features also manifest in the dispersion relation ω(k), which characterises the growth or decay of density modulations with wave number k in the liquid state.
For one-component systems in the liquid state with number density ρ 0 we may express a perturbation in the density profile δρ(r, t) ≡ ρ(r, t) − ρ 0 as a Fourier sum of modes with wave vector k of form ∼ exp(ik · r + ωt). The equation for the time evolution of the density ρ can be written ∂ t δρ = Lδρ + O(δρ 2 ), where ∂ t is a partial time derivative and L is a spatial operator [12][13][14]. Linearising, we see that ω is the eigenvalue of L acting on the eigenfunction exp(ik · r), so if ω < 0 for all k = |k|, then all modes decay and the uniform state is linearly stable. However, if ω > 0 for some k, then those modes grow over time. How ω(k) is related to the form of the soft pair potentials is well established [12][13][14][15][16].
To obtain QCs, the dispersion relation should exhibit maxima at k 1 and k 2 , with roughly equal peak values that are as close to zero as possible. This is equivalent to requiring that the static structure factor S(k) [17] should exhibit two prominent peaks at k 1 and k 2 , since for a bulk fluid ω(k) ∝ −k 2 /S(k) [6,8,9,16,18]. Additionally, ω(k) must be sufficiently negative at the reciprocal lattice vectors of k 1 and k 2 that are involved in stabilising competing periodic crystal structures (e.g. with wavenumbers √ 3k 1 , √ 3k 2 , etc), so that these are suppressed [11].
In addition to characterising how a uniform bulk liq-uid evolves after being perturbed, the dispersion relation is important because it gives crucial understanding of what wave number density modulations are favourable and which are likely be present in any incipient nonuniform crystalline or QC states [11]. In systems that are near to or beyond freezing, the characteristic modes that form the crystal or QC either grow or decay slowly. In one-component QC forming systems, this occurs in the vicinity of the point in the phase diagram where the system is marginally unstable at both k 1 and k 2 [1][2][3][4][5][6][7][8][9][10][11].
On the face of it, QCs should occur more widely in two component systems, since these intrinsically have at least two lengthscales, originating from the different particle sizes. Indeed, the vast majority of QCs discovered so far are metallic alloys with at least two components, e.g. Al-Mn or Ni-Cr [19,20]. For mixtures where the particles have a well defined (hard) core, requiring certain geometrical motifs as minimal energy structures in local particle arrangements can be a fruitful way to find QCs [21][22][23][24]. However, there is not an established 'recipe' for finding them, at least in soft matter. The three-step strategy we follow and advocate here, which works for the colloidal mixture model considered below and which we expect to be more generally applicable, is as follows: (i) Obtain the liquid state partial static structure factors. (ii) Select the parameters or state point such that they exhibit two peaks at wave numbers k 1 and k 2 with the specific ratio corresponding to the desired QC, whilst also checking that there is no peak at k = 0, as this is a signature of demixing, which can overtake the desired QC formation. (iii) Tune the parameters or state point so that the maxima in ω(k) are similar in height and as close to zero as possible. We successfully apply this strategy for finding QCs in a binary mixture modelled using a simple classical density functional theory (DFT) with direct correlations functions obtained from the hypernetted chain (HNC) Ornstein-Zernike integral equation theory [17]. This is only qualitatively correct for the system to which it is applied (see below), but the simplicity makes it an ideal test system on which to develop our approach. Our strategy has similarities to the approach used previously to find QCs in one-component systems, but for binary mixtures there are several additional complexities to overcome. Mixtures of big (b) and small (s) particles generically have at least three (not two) lengthscales present. These are the characteristic ranges of the b-b and s-s interactions and also the b-s cross interaction. In the beautiful work [25], results for several model binary mixtures of soft particles are presented, treated using a phase field crystal type theory [26]. This predicts the mixtures to form a variety of QCs, with the b-b and s-s interaction potentials providing the two lengthscales needed for the QC formation. However, treating these mixtures with a more accurate DFT, which retains the logarithmic ideal gas free energy instead of approximating via a Taylor expansion [26,27], we find that these systems actually just phase separate and do not form QCs (results not displayed). In the mixture considered here, the two lengthscales required for QC formation originate in the b-b and b-s interaction lengthscales. The s-s lengthscale seems to be irrelevant.
To determine the growth or decay rate of density perturbations in a uniform binary fluid mixture (i.e. the dispersion relation), where the bulk densities of the two species are ρ 0,b and ρ 0,s , one must consider the time evolution of δρ = (δρ b , δρ s ) = (ρ b − ρ 0,b , ρ s − ρ 0,s ). Fourier transforming the coupled dynamical equations for the two density profiles we obtain ∂ tρ (k, t) = Lρ(k, t) + O(ρ 2 ), whereρ is the Fourier transform of δρ and L is a 2 × 2 matrix [13]. In the Appendix we give an explicit expression for L, which depends onĉ ij (k), the Fourier transforms of the fluid pair direct correlation functions c ij (r), where i, j = b, s, for the case when the particles have Brownian equations of motion, i.e. where we can use dynamical DFT to describe the dynamics [16,[28][29][30][31]. For binary mixtures, the dispersion relation has two distinct branches, ω + (k) and ω − (k). The three partial structure factors S ij (k) are closely linked toĉ ij (k) [17] (see the Appendix) and therefore ω + and ω − depend crucially on the form of S ij (k). In Fig. 1 we display examples of S ij (k) and also ω + (k) and ω − (k) for the model defined below. Note that the peak locations in S ij (k) are where the peaks in ω + occur, i.e. these are the wave numbers of the slowest decaying density modes (or growing, if ω + (k) > 0). The lower panel of Fig. 1 shows how ω + and ω − vary as the inverse-temperature-like parameter J (defined below) is varied. We compare these with the ideal gas case, where ω(k) ∝ −k 2 , which highlights how the particle interactions are responsible for the shape of ω + and ω − . Since ω + ≥ ω − , the most important branch is ω + .
The system we consider is a 2D binary mixture of charged colloids adsorbed on a flat oil-water interface [32][33][34][35]. The interactions between particles can be modelled by the pair potentials where r is the distance between particles and β = 1/k B T , where k B is Boltzmann's constant and T is the temper- particles, Γ is the dimensionless interaction strength between b-particles and m i is the dipole moment ratio of species i = s, b relative to that of species b (i.e. m b = 1 and m s < 1). This system exhibits a rich variety of 2D crystal structures [33][34][35][36][37][38][39]. We describe the system using DFT with the Ramakrishnan-Yussouff (RY) approximation [40]. The RY DFT was used in Ref. [35] to obtain density profiles for system (1) at state points where periodic crystals occur, with the c ij (r) that are inputs to the DFT obtained from an accurate but computationally intensive theory. Here instead the c ij (r) are obtained from HNC theory [17], which is less accurate but simpler and so much faster [35]. The RY DFT and HNC theory are described in the Appendix. The speed of this is important because tuning the pair potential parameters following the three steps of the recipe above for finding QCs requires numerous calculations. So, although the theory we use is at best qualitatively accurate [35], it is used because of its speed and the fact that our aim is to test the efficacy of the strategy for finding QCs, not the accuracy of the theory. If more accuracy were required, one could use the DFT in Refs. [41,42].
Density profiles for a QC state are displayed in Fig. 2, with parameters ρ 0,b l 2 = 1, ρ 0,s l 2 = 2, Γ = 42 and m s = 0.025, identified by following the three-step strategy given above. The partial structure factors and dispersion relations displayed in Fig. 1 also correspond to  FIG. 2. The logarithm of the two densities and the total density: in (a) we display ln(ρ b l 2 ), in (b) ln(ρsl 2 ) and in (c) ln(ρ b l 2 + ρsl 2 ). In (d) is the Fourier transform of the b-particle density. The highest peaks in (a) and (c) correspond to the positions of the b-particles, while the s-particles in (b) are much more delocalised and almost fluid-like. The state point is the same as in Fig. 1, with J = 1. In (d), the inner circle has radius k1 and the outer has radius k2 = 1.93k1. this system. The branch ω + (k) has a peak at k = k 1 with ω + (k 1 ) approaching zero from below, but the second peak at k 2 > k 1 is not as high. The pair potential parameters were tuned so that k 1 /k 2 ≈ 1.93, in order to observe dodecagonal QCs. The aim was of course to have both peaks at the same height and as close to zero as possible, but the physical constraints stemming from the form of the potentials in Eq. (1) prevents this. In particular, φ bs (r) cannot be varied independently of φ bb (r) and φ ss (r). Of course, the HNC theory also fails before To obtain the QC density profiles after calculating the pair direct correlation functions, in the DFT we initially replace c ij (r) → Jc ij (r), where J is a constant scaling parameter. Recall that c ij (r) ∼ −βφ ij (r) for large r [17], so increasing J is much like decreasing the temperature. This makes freezing easier and the uniform density state to be linearly unstable, so using an initial guess for the two density profiles consisting of the desired average density value plus a small amplitude random field is sufficient to observe the QC formation (or periodic crystals, at other state points). See Ref. [35] for further details on this approach involving J for calculating solid state density profiles. To obtain the results in Fig. 2, we initially calculate for J = 1.5 and then take the resulting profiles as our initial guess at the physical value J = 1. The density profile for the b-particles exhibits sharp peaks, with the QC structure clearly visible. The Fourier transform in Fig. 2(d) shows the characteristic 12-fold symmetry. In contrast, the s-particles are much more delocalised and fluid-like, acting as a 'stabiliser' for the structure. For these particles, the density peaks represent preferred locations where particles might be found some of the time, as there are more such locations than particles. Similar behaviour was observed at other state-points, where periodic crystals are the equilibria [35].
The c ij (r) obtained from the HNC theory are displayed in Fig. 3. Inspecting these, one can roughly identify a typical lengthscale (effective diameter) as the range r beyond which c ij (r)/c ij (0) becomes small. The effective diameter obtained from c bb (r) is ≈ l and from c bs (r) is ≈ 0.5l. The ratio of these is ≈ 2, but one needs to go to Fourier space ( Fig. 1) to see much more precisely the ratio k 1 /k 2 ≈ 1.93, characteristic of QC formation. Note too that the effective diameter from c ss (r) is ≈ 0.25l. This corresponds to a wave vector k s l ≈ 2π/0.25 ≈ 25. Fig. 1 shows that neither ω(k) nor S ij (k) have significant features near k s , indicating that this lengthscale is irrelevant to the QC formation. A deeper understanding of the observed QC formation can be obtained by considering the phase diagram in the concentration χ ≡ ρ 0,b /(ρ 0,b +ρ 0,s ) versus J plane, calculated using only the (scaled with J) pair direct correlation functions displayed in Fig. 3 from the state point (χ, J) = (1/3, 1) as input to the RY DFT. The result is displayed in Fig. 4. We should emphasise that because the phase diagram is calculated by rescaling the c ij (r) from the state point (χ, J) = (1/3, 1) to all other state points, in a strict sense, this is the only physically relevant state point in Fig. 4. However, by exploring this theoretical model phase diagram, we obtain important insight into the observed QC formation that we would not obtain otherwise. At small J the uniform density liquid state is found; recall that decreasing J is like increasing T . At higher J the system freezes to form one of three different solid phases. For small χ, i.e. where the s-particles dominate, the system forms a hexagonal crystal with lattice spacing ≈ 2π/k s , which we refer to as s-Hex; see Fig. 4. A portion of a typical example is displayed in Fig. 5(a). The defects originate from the the random initial conditions. Increasing χ, the system forms a hexagonal crystal with much larger lattice spacing ≈ 2π/k 1 , determined by the range of c bb (r) (b-Hex in Fig. 4). These crystal structures are discussed in detail in [33]. See also the DFT results in [35]. In these, the s-particles can either be fluid-like, leading to a honeycomb like density distribution surrounding the peaks of the b-particles -see Fig. 5(c). Alternatively, they can be more localised, so that the b-particle density peaks are surrounded by density peaks from the s-particlessee Fig. 5(b). Moving to even higher χ, we find QCs.
An example is displayed in Fig. 2. Note that in this model [obtained by rescaling the c ij (r) from the state point (χ, J) = (1/3, 1)] the QCs extends right up to χ = 1 (where ρ 0,s = 0). We believe this is because the influence of the s-particles is still present in the rescaled c bb (r) that is calculated at (χ, J) = (1/3, 1), although it could also be because −c bb (r) is somewhat akin to the soft effective pair potential of the monodisperse system in Ref. [43], which forms QCs.
This simplified model enables us to easily calculate the linear stability threshold for the uniform liquid, i.e. the locus in the phase diagram where either ω + (k 1 ) = 0 or ω + (k 2 ) = 0, or ω + (k s ) = 0. These are the lines in Fig. 4. Those satisfying the first two of these conditions meet at (χ, J) = (0.85, 2.715), where the system is marginally unstable at both k 1 and k 2 (right hand cusp on the solid line in Fig. 4). Recall that for monodisperse systems points of this type are intimately connected with QC formation [1][2][3][4][5][6][7][8][9][10][11]. Thus, finding this point explains much of why we observe QCs in the present binary mixture. Importantly, notice that this point exists in the theoreticallyconstructed (χ, J) plane, rather than in the physical parameter space of the original system. So, although in the physical parameter space the two peaks at k 1 and k 2 in ω(k) are not at the same height and nor is the second peak in S ij (k) at k 2 as prominent as the first at k 1 , nonetheless there is still the influence of the interaction between density modes with wavenumber k 1 and k 2 to stabilise the QC state. Much insight on such two-mode interactions is in the pattern formation literature related to Faraday waves [1,10,[44][45][46][47][48][49][50][51][52][53][54][55][56][57] and though binary mixtures have the added complication of consisting of two coupled fields, much of this insight still applies.
In Fig. 4 the boundaries between regions of the different phases are only guides for the eye. For J = 2 we have determined the states at coexistence between the b-Hex and QC phases and found the width of the coexistence region to be ∆χ ≈ 0.04 (not displayed). Since this is small, it justifies our approximate approach for identifying the locations of the phase boundaries (see the Appendix).
To summarise, we have proposed a 'recipe' for finding QCs in soft matter mixtures. The key quantities for inspection are the partial static structure factors and the dispersion relation. In the model system studied here the two lengthscales required for QC formation arise from the b-b and b-s particle interactions. In principle, these could instead arise from the b-b and s-s interactions, but from our studies of soft-particle models (not shown), phase separation occurs much more readily than QC formation in this case. Using the cross-interaction b-s lengthscale as one of the key QC lengthscales helps to avoid this. We thank Ron Lifshitz, Daniel Ratliff, Alastair Rucklidge, and Priya Subramanian for fruitful discussions at various points over the course of this work. We also thank Ron Lifshitz for sending us Ref. [25].

APPENDIX The Ornstein-Zernike equation and liquid state structure
The Ornstein-Zernike (OZ) equations for the total correlation functions h ij (r) of a binary fluid mixture are: (2) where c ij (r) are the pair direct correlation functions and ρ 0,i for i = b, s are the bulk fluid densities of the two species [17]. The radial distribution functions are related to the total correlation functions via g ij (r) = 1 + h ij (r). These coupled equations must be solved in conjunction with the following (exact) closure relations (3) where B ij (r) are the so-called bridge-functions, φ ij (r) are the pair potentials and β = 1/k B T [17]. The hypernetted chain (HNC) approximation consists of setting B ij (r) = 0 for all r. Due to the convolutions in (2), on Fourier transforming we obtain the following set of algebraic equationŝ whereĥ ij (k) andĉ ij (k) are the Fourier transforms of h ij (r) and c ij (r), respectively. The partial static structure factors are related to these as follows [13,17]: S bb (k) = 1 + ρ 0,bĥbb (k), S ss (k) = 1 + ρ 0,sĥss (k), From (4) we obtainĥ with the numerators given by and the common denominator (8) For the stable liquid, D(k) > 0 for all k. However, if this is not the case, then the liquid state is unstable. Thus, we can determine the stability threshold for the uniform liquid from solving for the locus in the phase diagram where a solution to the equation D(k) = 0 appears.

Density functional theory for binary mixtures
The density profiles ρ i (r) are obtained using classical density functional theory (DFT). The grand potential of the system is [15,17] where F is the intrinsic Helmholtz free energy functional, V ext i (r) is the one-body external potential acting on species i (here we set V ext i (r) ≡ 0 for i = b, s, in order to study bulk phases) and µ i are the chemical potentials. The intrinsic Helmholtz free energy can be split into two terms where the first term is the ideal gas contribution, where Λ i is the (irrelevant) thermal de Broglie wavelength and d is the dimensionality of the system. The second term in Eq. (10) is the excess Helmholtz free energy, arising from the interactions between the particles. Following Ramakrishnan and Yussouff [40], the approximation we use here is to expand this functional around the homogeneous fluid state in a functional Taylor expansion and truncate at second order, giving i=b,s j=b,s drδρ i (r)c ij (| r − r |)δρ j (r ), (12) where δρ i (r) = ρ i (r)−ρ 0,i and µ ex i = µ i −k B T ln ρ 0,i Λ d i are the excess chemical potentials. We further approximate the pair direct correlation functions c ij (r) via those obtained from the HNC theory. The equilibrium density profiles are those which minimise the grand potential Ω and which therefore satisfy the following pair of coupled Euler-Lagrange equations for i = b, s.
Dynamics: the growth or decay of small amplitude density perturbations When the equations of motion of the particles can be approximated by stochastic Brownian equations of motion, then dynamical density functional theory (DDFT) shows that the non-equilibrium density distributions for the two species of particles ρ i (r, t) is described by [16,[28][29][30]: where the mobility coefficient γ i = βD i and where D i is the diffusion coefficient of species i. Note that if instead the particles evolve according to Newton's equations of motion, then the equations for the time evolution of the density profiles are more complicated, but in dense systems one can argue that Eq. (14) still governs the long time (on diffusive timescales) behaviour [31]. If we consider the growth or decay of small amplitude density perturbations around the bulk value of the form δρ i (r, t) = ρ i (r, t) − ρ 0,i , then we can expand Eqs. (14) to obtain [12][13][14][15][16]: Linearising this equation and then Fourier transforming, we obtain whereρ i (k, t) is the Fourier transform of δρ i (r, t) and k = |k|. Assumingρ i (k, t) ∝ exp(ω(k)t), then Eq. (16) becomes [13]: whereρ = (ρ b ,ρ s ) and the matrix L = ME, where the two matrices M and E are defined as and Solving Eq. (17) for the dispersion relation ω(k), one obtains two branches of solutions, ω ± (k). These are given by ω ± (k) = 1 2 Tr(ME) ± 1 4 Tr(ME) 2 − det(ME). (20) Further details of this derivation can be found in Ref. [13]. Note that the equation det(E) = 0 is entirely equivalent to solving D(k) = 0, from Eq. (8).
It is worth recalling that the values of the diffusion coefficients D b and D s do not ever determine which structure is the thermodynamic equilibrium state, i.e. the minimum of the free energy. Therefore, the values of D b and D s are not involved in determining the phase diagram in Fig. 4 of the main text. Nor do the values of D b and D s determine the locations of the linear stability threshold lines in the phase diagram, i.e. the lines in Fig. 4 where either ω + (k 1 ) = 0 or ω + (k 2 ) = 0 or ω + (k s ) = 0. This is because these lines come from solving the equation det(E) = 0, whilst the values of the diffusion coefficients only enter the mobility matrix M in Eq. (18). That said, the precise value of the ratio D b /D s does influence the dispersion relation curves, but does not affect where the peaks occur (i.e. does not change k 1 or k 2 ). Thus, the value of the ratio D b /D s is only relevant to the non-equilibrium dynamics of the system. However, since here we are solely ultimately interested in the equilibrium phase behaviour of the system, which does not depend on D b /D s , we therefore set this ratio equal to 1, i.e. we set D b = D s = D.
Note on the width of the coexistence region between the QC and the b-Hex phase In the main text we comment briefly on the fact that in the phase diagram in Fig. 4 the coexistence region between the QC and the b-Hex phase is fairly small. It is worth expanding on those comments here. That the coexistence region is narrow is important, because it implies that in large portions of the phase diagram (as displayed in Fig 4), the QC is the thermodynamic equilibrium. In the main text we give the width of the coexistence region ∆χ ≈ 0.04 for J = 2. For lower values of J the coexistence region becomes a little broader (e.g. at J = 1.5 the width of the coexistence region ∆χ ≈ 0.06) and for higher J it is narrower. Other model systems where the coexistence gap between the QC and hexagonal phases is very narrow include the systems described in Refs. [6,11], so based on our experience with those systems, the narrowness in the present system is perhaps not too surprising.
Another observation on this issue worth noting is the following: If one initiates the system in the QC state and then decreases χ in small steps, following the QC branch of solutions, one eventually falls off that branch onto the b-Hex phase branch of solutions. For example, for J = 1.5 this occurs at χ ≈ 0.3. Some authors would refer to this as the "spinodal" point for the QC phase. In other words, for J = 1.5 and χ < 0.3 the QC state is no longer a stable solution to the model equations. In a similar way, if one initiates the system in the b-Hex state and then increases χ in small steps, following the b-Hex branch of solutions, one eventually falls off that branch onto a state that is a periodic approximant for the QC state. For J = 1.5 this b-Hex spinodal point occurs at χ ≈ 0.37. In other words, for J = 1.5 and χ > 0.37 the b-Hex state is no longer a stable solution to the model equations. This fact that the system falls from b-Hex branch of solutions onto a branch related to the QC state is a very strong indicator that the QC is the thermodynamic equilibrium state. Moreover, the dis-tance in the phase diagram between these two spinodal points 0.37 − 0.30 = 0.07, is an upper bound for the coexistence region width ∆χ.