This item was submitted to Loughborough’s Institutional Repository (https://dspace.lboro.ac.uk/) by the author and is made available under the following Creative Commons Licence conditions. For the full text of this licence, please go to: http://creativecommons.org/licenses/by-nc-nd/2.5/ PHYSICAL REVIEW E 85, 061408 (2012)Modeling the structure of liquids and crystals using one- and two-component modifiedphase-field crystal modelsM. J. Robbins,1 A. J. Archer,1 U. Thiele,1 and E. Knobloch21Department of Mathematical Sciences, Loughborough University, Leicestershire LE11 3TU, United Kingdom2Department of Physics, University of California, Berkeley, California 94720, USA(Received 28 March 2012; published 21 June 2012)A modified phase-field crystal model in which the free energy may be minimized by an order parameterprofile having isolated bumps is investigated. The phase diagram is calculated in one and two dimensions and welocate the regions where modulated and uniform phases are formed and also regions where localized states areformed. We investigate the effectiveness of the phase-field crystal model for describing fluids and crystals withdefects. We further consider a two-component model and elucidate how the structure transforms from hexagonalcrystalline ordering to square ordering as the concentration changes. Our conclusion contains a discussion ofpossible interpretations of the order parameter field.DOI: 10.1103/PhysRevE.85.061408 PACS number(s): 82.70.Dd, 61.50.Ah, 05.70.Fh, 02.30.OzI. INTRODUCTIONModeling materials at the atomic scale is a task which,for example, may be performed using molecular dynamicssimulations. This involves solving coupled equations ofmotion to calculate the position of each particle at every timestep. The resulting calculations can be very computationallyexpensive, especially when one seeks to consider phenomenawhich involve a large number of particles. Only short atomictime scales can be feasibly accessed with this or other suchapproaches. However, there are some instances where it isimportant to consider materials on the atomic length scalefor much longer diffusive time scales, for example, wheninvestigating freezing or glass transitions. One approach tosuch problems that may be adopted is to develop a phase-fieldmodel capable of describing the structure of materials onthe scale of the individual particles. In contrast to traditionalphase-field models, the recently developed phase-field crystal(PFC) models are capable of just such a description and arenow widely used in the literature to model crystalline struc-tures [1–5]. The PFC model consists of a Swift-Hohenberg-like equation [6], but with conserved dynamics rather thanthe nonconserved dynamics of the regular Swift-Hohenbergequation. Similar models with conserved dynamics arise inquite different contexts as well [7,8]. The regular PFC modelis governed by the equation∂φ(x,t)∂t= α∇2 δF [φ]δφ(x,t) , (1)where the free energy functionalF [φ] =∫dxf (φ), (2)withf (φ) = φ2[r + (q2 + ∇2)2]φ + φ44, (3)where α is the mobility coefficient, r is the undercoolingparameter that decreases with decreasing temperature, q is aconstantwhich determines the typicalmicroscopic length scalein the system, and φ(x,t) is the order parameter. For certainparameter values this free energy functional is minimized byan order parameter profile consisting of a periodic array ofbumps which somewhat resembles the density distributionof particles in a crystalline material. This interpretation isbolstered by the fact that it has been shown that the PFCmodel [Eqs. (1) and (2)] may be derived from the densityfunctional theory of freezing [9] and the dynamical densityfunctional theory for colloidal particles [10,11] with certainapproximations. The free energy is minimized by eitherperiodic structures or by a homogeneous flat profile, dependingon the values of q, r , and ¯φ = 1Ld∫dxφ(x,t), where Ldis the size of the system. In two dimensions (d = 2), fourdifferent phases are observed: a homogeneous phase, twohexagonal phases (one consisting of hexagonally orderedbumps and the other of hexagonally ordered holes), and astripe phase [1–3,12]. The literature largely focuses on theregion of the two-dimensional (2D) phase diagram whichcontains the hexagonally arranged bumps and their transitionto the homogeneous state [1–5,9,10,12]. The uniform profileφ(x,t) = ¯φ represents the order parameter in a uniform liquidand the hexagonal phase is treated as a crystal. The model isthen used to consider a number of problems including meltingand freezing [3,9,10] and grain boundary effects [1,2,4].In this “standard” PFCmodel [Eqs. (1)–(3)] the hexagonallyarranged bumps are considered to be particles or colloids in acrystalline structure. In this context, the physical interpretationof the striped and hexagonally ordered hole structures isunclear and, as such, these phases are commonly ignored.The conjecture that the ordered bumps represent crystallineparticle structures can be extended by including a “vacancyterm” in the free energy [13,14], which strongly breaks thehole-bump (φ → −φ) symmetry of Eq. (3):F [φ] =∫dx[f (φ) + fvac(φ)]. (4)Using the free energy (4), it is possible to obtain structureswhich contain a mixture of bumps and vacant areas (the “va-cancies” are areas where the order parameter is approximatelyuniform around the value φ ≈ 0). As the system parametersare varied, one can observe states where the number of bumpschanges and the proportion of the system in which φ ≈ 0varies. In particular, one can go from states where the density061408-11539-3755/2012/85(6)/061408(22) ©2012 American Physical SocietyROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)of bumps is high, vacancies are absent, and the bumps arehexagonally ordered, much as in the standard PFC model.However, for other parameter values, the density of bumpsis significantly lower, with more vacant φ ≈ 0 regions andno long range order in the location of the bumps [13]. Thearrangement of the bumps in these states resembles snapshotsof the particles in a fluid. The interpretation that the bumpsrepresent particles and that this “vacancy phase-field crystal”(VPFC) may be used to model particle configurations on acoarse-grained scale of both crystalline solids and disorderedfluids is the central assumption of Ref. [13]. We return to theissue of the precise interpretation of the nature of the orderparameter field in the conclusion.In this paper we investigate the thermodynamics andthe structures formed in the augmented conserved Swift-Hohenbergmodel with the free energy functional (4), hereafterthe VPFC model, and extend the theory to a two-componentgeneralization of this model. The vacancy term takes thefollowing form [13,14]:fvac(φ) = Hφ2(|φ| − φ), (5)where H is a constant. We use the value H = 1500, as inRefs. [13,14], although the precise value ofH is not significant.We have confirmed that choosing H in the range 500–3000makes no qualitative difference to the structures formed. Thisvacancy term acts as a piecewise function which is zero forpositive values of φ and takes an increasingly large valuewhen φ < 0. Hence, this term penalizes negative values ofφ. This leads to the VPFC model forming periodic structureswhich are somewhat different from those of the regular PFC.In addition, the VPFC model has a large region of parameterspace at small ¯φ, where spatially localized structures form.The time evolution of the order parameter φ is governed bythe conserved dynamics used in the standard PFC model (1).We begin in Sec. II by considering the phase behaviorof the model, investigating the transition between periodicand localized states. We focus on understanding the bifurca-tion diagrams connecting the various uniform, periodic, andlocalized states exhibited by the model. We then go on toconsider how individual localized states or particles interactwith one another. In Sec. III we extend the model to considera two-component system, and we determine how the particlesin the binary mixture interact with one another. We finda transition between hexagonal and square ordering of theparticles as the concentration changes. Our conclusions followin Sec. IV and include a discussion of the proper interpretationof the order parameter field φ.II. ONE-COMPONENT SYSTEMA. Linear stability of a homogeneous profileWe begin by considering the phase behavior of the VPFCmodel [Eqs. (1) and (4)]. We calculate the limit of linearstability for a homogeneous flat state using a linear stabilityanalysis. In the context of colloidal suspensions exhibitingmicrophase separation and fluids of charged particles, this limitof linear stability is referred to as a “λ line” [15–18]. Sincefvac is nondifferentiable at φ = 0 we treat it in a piecewisemanner, by treating the two cases ¯φ > 0 and ¯φ < 0 separately[in fact, if φ(x) takes the form of Eq. (6) and ¯φ > |ξ |, thenfvac = 0 everywhere and the thermodynamics of the VPFCmodel reduces to that of the regular PFC model]. We assumethat the order parameter φ takes the form of a flat profile plusan additional small-amplitude harmonic modulation:φ = ¯φ + δφ = ¯φ + ξeik·xeβt , (6)where ¯φ is the average value of the order parameter andthe amplitude |ξ | 1. Substituting this expression into thefunctional derivative of the free energy (4) we obtainδFδφ= (r + q4) ¯φ + 3H ¯φ(| ¯φ| − ¯φ) + ¯φ3+ [(k2 − q2)2 + ]δφ + O(δφ2), (7)where = r + 6H (| ¯φ| − ¯φ) + 3 ¯φ2. (8)Inserting this expression for the functional derivative (7) intothe dynamical equation (1) and then linearizing we arrive atthe following dispersion relation:β = −k2α[(k2 − q2)2 + ]. (9)When the growth rate β(k) > 0, any small amplitude modula-tion with wave number k = |k| will grow over time. There isa local maximum in β (which becomes the global maximumwhen the uniform state is unstable) at the wave number:km = 13√6q2 + 3√q4 − 3 . (10)Thus, if one takes an initially almost flat profile φ(x,t =0) = ¯φ + X (x), where X (x) is composed of a sum of a largenumber of small-amplitude harmonic modulations [cf. Eq. (6)]with different wave numbers k [in practice X (x) is generatedby adding a small random number to the discretized initialprofile], then as the system evolves in time φ(x,t) will developspatial modulation on the length scale 2πkm, since this scalecorresponds to the maximum growth rate βm ≡ β(km). Thislength scale has an inverse dependence on the value of q; thatis, increasing the value of q reduces the length scale of thestructures which are formed.The limit of linear stability is defined as the locus of pointsin parameter space where the maximum in the dispersionrelation (9) is at zero, that is, βm = 0. The conditions β =∂β∂k= 0, subject to the requirement that km = 0 yield = 0,km = ±q. Thus, in Eq. (8) can be considered as a measureof stability: When < 0 the system is linearly unstable andwhen > 0 the system is linearly stable. The magnitude of indicates how “far” we are from the limit of stability. Figure 1shows the dispersion relations β(k) for various values of . Inaccordance with Eq. (10), the maximum at km ≈ q disappearswhen > q43 ; in this case only themaximumat k = 0 remains.It is important to note that these results are identical to thoseof the regular PFC model [Eqs. (1) and (2)] for positive valuesof the order parameter ¯φ > 0.B. One-dimensional modelIn order to develop a better understanding of the effectof the “vacancy term” (5) we initially consider the phasediagram for the system in one spatial dimension. The regularPFC model [Eqs. (1) and (2)] in 1D exhibits two distinct061408-2MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012) 0 0kq = -0.2 = 0.0 = 0.2 = 0.4FIG. 1. (Color online) Dispersion relation curves for the VPFCmodel [Eqs. (1) and (4)] when q = 1. Four cases are shown, with(i) β(km) > 0 (red solid line), (ii) β(km) = 0 (blue dashed line), (iii)β(km) < 0 but km > 0 (green dotted line), and (iv) β 0 and km = 0(magenta dash-dotted line).phases [2]: a nonuniform state in which the order parameterprofile resembles a sinusoid and a uniform state in whichthe order parameter is a constant. The phase diagram of theregular PFC model is symmetric around ¯φ = 0 owing to thesymmetry of the free energy (2) with respect to φ → −φ. Thisis no longer the case when the vacancy term (5) is added.The phase diagram for the 1D VPFC model is shownin Fig. 2(a) and is very different from that of the regularPFC [2]. As with the regular PFC model, modulated profilesare present below the limit of linear stability = 0 (bluedashed line) provided ¯φ < √3/2q2. However, with the addedvacancy term (5) the lower limit for the presence of themodulated phase is at ¯φ 0 (H 1). The tricritical pointwith ¯φ > 0 (red dot) familiar from the PFC model remains.Above this point the phase transition between the periodicand homogeneous phases is of second order. Below this pointa periodic phase with ¯φ = ¯φp coexists with a homogeneousphase with ¯φ = ¯φh and the phase transition between thesephases is of first order. Figure 2(a) shows the coexisting phasesusing fixed temperature (horizontal) tie lines connecting ¯φpand ¯φh (solid red lines).We can calculate the location of the tricritical point asfollows: Since the wave number near = 0 is km ≈ q weassume that the order parameter profile takes the formφ = ¯φp + A cos qx + B cos 2qx + · · · (11)and compute the free energy F . When ¯φp > 0 the vacancyterm drops out and we obtain the following expression forthe free energy per unit length fp = F/L of the periodicallymodulated phase:fp = 12 (r + q4) ¯φ2p + 14 ¯φ4p + 14 rA2 + 14 (r + 9q4)B2+ 34 ¯φ2p(A2 + B2) + 34 ¯φpA2B + 332 (A4 + 4A2B2 + B4).(12)We refer to the ansatz (11) as the two-mode approximation.This approximation is reliable around and above the tricriticalpoint since the amplitude of the modulation in φ is small when| | 1.Moreover the two-mode approximation appears to beexact at the tricritical point, since the location of the tricriticalpoint is unaffected by the inclusion of cos 3qx and other higherorder modes. In contrast, the mode B cos 2qx must be retainedin order to obtain the correct value of the amplitude A in thevicinity of the tricritical point (see below).To demonstrate this we minimize fp in Eq. (12) withrespect to the amplitudes A and B, obtaining the followingtwo conditions:r + 3 ¯φ2p + 3 ¯φpB + 34A2 + 32B2 = 0, (13)(r + 9q4)B + 3 ¯φ2pB + 32 ¯φpA2 + 32A2B + 34B3 = 0. (14)Solving these for the amplitudes A and B and substitutinginto Eq. (12), we obtain an approximation for the free energydensity of the periodic phase fp. Linearizing Eq. (14) in B,we find thatB = −¯φpA26q4+ O( pA2,A4), (15)where p ≡ r + 3 ¯φ2p and hence, from Eq. (13), thatA = 2√− p3(1 − 2¯φ2p3q4)−1/2+ O( p). (16)The free energy density of the homogeneous phase fh, havingφ(x) = ¯φh, is obtained simply setting A = B = 0 in Eq. (12)to obtain fh = 12 (r + q4) ¯φ2h + 14 ¯φ4h. The chemical potential inthe homogeneous phase is μh = ∂fh/∂ ¯φh, and in the periodicphase μp = ∂fp/∂ ¯φp.To calculate the location of the tricritical point we recall thatat coexistence between the periodic state and the homogeneousstate we must have μp = μh. We write the average value of φin the periodic state ¯φp = ¯φh + C, where C is the differencebetween the average value of the order parameter in the twocoexisting phases, implying that the coexistence condition isμp( ¯φh + C) − μh( ¯φh) = 0, (17)or equivalently,(r + q4)C + 3 ¯φ2hC + 32 ¯φhA2 + O(A4,CA2,C2) = 0. (18)Since the amplitude A of the modulated phase at coexistenceis small when | h| 1, where h ≡ r + 3 ¯φ2h, Eqs. (14)and (18) yield, for | h| 1, the expressionsB = −¯φhA26q4+ O( hA2,A4),(19)C = −3¯φhA22q4+ O( hA2,A4).Equation (13) then yieldsA = 2√− h3(1 − 38¯φ2h3q4)−1/2+ O( h). (20)Taking the limit C → 0 now takes us to the tricritical point. Atthe tricritical point h = p = 0 and the chemical potentialsμp( ¯φh) and μh( ¯φh) are identical. Thus, from Eq. (20) we seethat the tricritical point occurs at ¯φ = √3/38q2 ≈ 0.281q2,r = −(9/38)q4 ≈ −0.237q4. These coordinates agree pre-cisely with the result obtained from the common tangentconstruction between the free energy of the periodic phaseand the free energy of the homogeneous phase to determinephase coexistence and also with our numerical simulations of061408-3ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012) 0 0.2 0.4 0 50 100x(b) 0 0.2 0.4 0 50 100x(c) 0 0.2 0.4 0 50 100x(d) 0 0.2 0.4 0 50 100x(e) 0 0.3 0.6 0 50 100x(f) 0 0.5 1 0 50 100x(g)(a)(a)FIG. 2. (Color online) The phase diagram for the 1D VPFC model [Eqs. (1) and (4)] is displayed in (a) for the case q = 1. The red solidlines are the coexistence curves between the periodic and uniform phases; the red circle is the tricritical point. The blue dashed line is the locus = 0, which is the limit of linear stability for uniform profiles. The green dot-dashed lines are a guide showing the parameter space wherelocal and periodic structures are formed. (b)–(g) Examples of order parameter profiles from numerical simulations corresponding to (local)minima of the free energy for the values of ¯φ and r indicated in (a). The parameter values are q = 1, r = −0.9, α = 1 and (b) ¯φ = 0.01, (c)¯φ = 0.05, (d) ¯φ = 0.1, (e) ¯φ = 0.175, (f) ¯φ = 0.3, and (g) ¯φ = 0.5.the VPFC model performed with q = 1. Thus, when q = 1the phase transition is of second order for r > −9/38 and offirst order for r < −9/38. As already mentioned, this resultis exact in the sense that it is unchanged if more modes areincluded in the ansatz (11) and improves on the predictionr = −1/4 obtained for the PFC model using a one-modeapproximation [2]. Note that the vacancy term (5) does notaffect the transition because φ is positive everywhere.As discussed above, the transition between the periodic anduniform phases is of first order below the tricritical point. As rdecreases, this region of the phase diagram is increasinglyaffected by the vacancy term (5), as the amplitude of thestructures becomes large enough to reach negative φ values.We observe that including the vacancy term (5) decreasesthe distance between the coexistence curves. This is becausethe vacancy term increases the free energy of the profiles inthe periodic phase, which decreases the difference between thefree energy of the periodic structures and the homogeneousstate and hence a common tangent construction between thetwo yields values which are closer to the linear stabilityline = 0. We calculate the coexistence values below thetricritical point by numerically solving for the order param-eter profiles because the two-mode approximation becomesinaccurate when r < −0.3, where the order parameter profiledevelops regions where φ < 0 and so the vacancy term makesa contribution to the free energy. Note that for the regular PFCmodel (i.e., when H = 0) the two-mode approximation forthe free energy works very well, agreeing to two significantfigures or more with the exact free energy for r −0.9.To calculate the coexisting phases numerically we select avalue of r and determine the order parameter profile along thisline for different values of ¯φ. The free energy at each of thesepoints is minimized with respect to the domain size, whicheffectively gives us the minimum free energy for the infinite061408-4MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012)system. A polynomial is then fitted to these values to producea continuous curve which gives the free energy of the periodicphase for the chosen value of r . We then make the commontangent construction between the free energy of the periodicphase and the uniform phase to calculate the two coexisting ¯φvalues at the chosen value of r . This process is then repeatedfor different values of r . The resulting coexistence curves areplotted as the solid red lines in Fig. 2(a).The periodic structures which are formed by the VPFCmodel [Fig. 2(f)] are qualitatively very similar to the structureswhich can be found in the regular PFC model. However,the amplitude of the modulations is restricted by the largepenalty in the free energy accumulated when φ < 0. Inside thecoexistence region between the periodic and uniform states, weobserve interesting structures where the amplitude of the peaksdoes not remain constant and a second length scale is visible inthe structures [Fig. 2(g)]. This is also an effect which is presentin the regular PFC model and will be discussed in detailin future work. What is most intriguing, and is perhaps themost appealing aspect of the VPFC model, is the appearanceof localized states for small positive values of ¯φ when themagnitude of r is sufficiently large (r −0.6). We obtainorder parameter profiles by numerically integrating forwardin time Eqs. (1) and (4) until a stationary solution is reached,starting from the initial profile φ(x,t = 0) = ¯φ + X (x), whereX is a small amplitude random noise profile with zero mean.A rich variety of different patterns is observed, includingperiodic structures mixed with almost flat regions [Fig. 2(d)]and individual isolated peaks [Figs. 2(b) and 2(c)]. In Fig. 2(a)the green dot-dashed curves indicate the boundary of theregion where one observes regular periodic structures andwhere the localized structures are formed. Note that theseare guidelines only and are not thermodynamic coexistencecurves. The lower-left dot-dashed curve roughly denotes thelinear stability limit of the regular periodic structures, suchas that in Fig. 2(f). This is determined numerically. We beginwith a periodic profile and reduce the value of ¯φ gradually,minimizing the free energy at each step, while keeping rconstant. The limit point is then defined as the value of ¯φ wherethe periodic profile becomes linearly unstable and a vacancyis introduced. In a similar way, we determine the upper-rightdot-dashed line, which is the limit of linear stability of thestructures with defects. This is found by starting with a profilecontaining a single vacancy and increasing ¯φ until the vacancydisappears. These two points are calculated for different valuesof r and then a best fit to this data is shown in Fig. 2(a). Thereis some hysteresis in the region between these two curves, withthe type of profile produced depending heavily upon the initialconditions.Within the localized state region of the phase diagram itis possible to obtain order parameter profiles with a varyingnumber of peaks for a given system of length L. Keeping r −0.6 constant and varying ¯φ allows us to control the number(density) of bumps, as shown in Figs. 2(b)–2(e). Beginningwith ¯φ = 0 we find isolated peaks in large vacant areas (whereφ is approximately uniform with φ ≈ 0). As ¯φ is increasedthe number of peaks increases until we return to the familiarregular periodic structures. The assumption of Ref. [13] isthat, unlike in the regular PFC, where the uniform phase isassociated with the liquid and the modulated phase with thecrystal, in the VPFC model one may associate each bump inφ(x) as corresponding to a particle, so the model can describefluids [Figs. 2(c) and 2(d)], crystals with vacancies and defects[Fig. 2(e)], and regular crystals [Fig. 2(f)]. Note that the precisevalue of the parameter H in the free energy has no qualitativeeffect on the structures that are formed. The main effect ofvarying H is to change slightly the value of φ in the vacantareas between the bumps. For example, for ¯φ = 0.02 and r =−0.9, when H = 500 we find that φ ≈ −5.3 × 10−3 in thevacant areas, when H = 1500 (the value used for much ofthe work presented here) we find φ ≈ −3.0 × 10−3, and whenH = 3000 we find that φ ≈ −2.1 × 10−3 in the vacant areas.The findings presented in Fig. 2 indicate the existence of ahysteretic transition between periodic and localized states andare a consequence of homoclinic snaking [19–22] in the presentsystem. In the standard homoclinic scenario such localizedstates are present within a part of the coexistence region calledthe pinning region. The localized states in the lower left part ofthe parameter plane ( ¯φ,r) in Fig. 2(a) correspond to the globalenergy minimum or to other deep but local energy minima.Families of such steady state solutions can be obtained for theVPFC model that we study here [Eq. (1) with Eqs. (3), (4),and (5)] by employing the path continuation techniquesbundled in the package AUTO07P [23]. As an example, in Figs. 3and 4 we show the characteristics of localized solutions alongcuts through the plane ( ¯φ,r). In particular, Figs. 3 and 4 giveresults for changing r (at constant ¯φ = 0.1) and changing ¯φ(at constant r = −0.9), respectively. All solutions are charac-terized by their L2 norm ‖δφ‖ ≡√(1/L) ∫ L0 (φ(x) − ¯φ)2dx,chemical potential μ = δF/δφ, mean free energy densitydifference (F [φ(x)] − F0)/L, where F0 = F [ ¯φ], and meangrand potential density ω ≡ F [φ(x)]/L − ¯φμ, and satisfyperiodic boundary conditions on the domain 0 x L.It turns out that there exist three types of localized steadystates. (i) The heavy solid black line consists of x → −xsymmetric localized states that have a maximum at the center;that is, the overall number of bumps within the structure is odd.(ii) The dashed red line also represents x → −x symmetriclocalized states but this time with a hole (minimum) at thecenter. (iii) The localized solutions of the third type arenot symmetric under x → −x and are called “asymmetricstates.” These reside on branches that connect (via pitchforkbifurcations) the two branches of symmetric localized states.These branches are included in the bifurcation diagrams asdotted green lines.Examples of order parameter profiles of types (i)–(iii) arepresented in Fig. 5, corresponding to the various solutionbranches displayed in Fig. 4. This sequence of profiles expandsupon the few examples shown in Figs. 2(b)–2(g). Recall,however, that the results in Figs. 2(b)–2(g) are obtained startingfrom an order parameter profile with a small amplitude randomnoise and so they do not always exactly agree with the steadystates at the same ¯φ resulting from the path continuation. TheL2 norm, chemical potential μ, mean free energy density(F − F0)/L, and the mean grand potential F/L − ¯φμ havebeen calculated for the profiles obtained from time simulations[Figs. 2(b)–2(f)] and are plotted as blue dots in Fig. 4. A closeinspection reveals that the energy of the time simulation resultsis often slightly higher than that from the continuation results,061408-5ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)(a)-2 -1.5 -1 -0.5 r0.10.20.3||δφ||n=15n=16n=14-2 -1.5 -1 -0.5 r-0.5-0.4-0.3-0.2-0.10μn=15n=16n=14(b)(c)-2 -1.5 -1 -0.5 r-0.03-0.02-0.010(F-F0 )/Ln=15n=16n=14-2 -1.5 -1 -0.5 r-0.01-0.00500.0050.010.015ωn=15n=16 n=14(d)FIG. 3. (Color online) Bifurcation diagram showing localized solutions of the VPFC (an augmented conserved Swift-Hohenberg equation)[Eqs. (1) and (5)] with H = 1500, as a function of the parameter r , for the mean order parameter ¯φ = 0.1 and a fixed domain size of L = 100.The various solution profiles are characterized by their (a) L2 norm, (b) chemical potential μ, (c) mean free energy density (F − F0)/L, and(d) mean grand potential density ω ≡ F/L − ¯φμ. The heavy black dotted line corresponds to the homogeneous solution φ(x) = ¯φ. Periodicsolutions with n = 15 bumps are shown as a thin blue dashed line, whereas the nearby thin black dash-dotted lines represent the n = 14 andn = 16 solutions as indicated in the plot. The heavy solid black and dashed red lines that bifurcate from the n = 15 periodic solution representsymmetric localized states with a maximum (odd states) and a minimum (even states) at the center, respectively. The green dotted lines thatconnect the two branches of symmetric localized states correspond to asymmetric localized states. Together the branches of localized statesform a slanted snakes-and-ladders structure.indicating that in these cases the time simulation converges to alocal and not the global energyminimum.This is to be expectedas the solutions shown in the bifurcation diagrams are onlythe “tip of the iceberg.” For instance, there exist many moresolutions, where not all the inner distances between the bumpsare identical. This is related to the fact that individual bumpshave oscillatory tails and the “locking of these tails” allows fordifferent equilibrium distances [24]. The solutions presented inFigs. 3 and 4 represent the solution having the lowest energy inthe respective class. However, the energy differences betweenthese and the “less symmetric” solutions are often tiny. Thus,it is not surprising that time simulations starting from randominitial profiles often converge to solutions with greater disorderand energies than those shown in Fig. 4. For instance, thesolution in Fig. 2(d) at ¯φ = 0.1 is a 9-bump solution similarto the odd symmetric localized states shown in the first twopanels of the second row of Fig. 5. The amplitudes agree welland although the arrangements of the 9 bumps are different,the free energy and norm still agree to <1%. However, at largeaverage order parameter values ¯φ the time simulation resultscan converge to metastable states with energies quite differentfrom the minimum energy states for domains of this sizeL = 100. For example, the periodic solution obtained fromthe time simulation [shown in Fig. 2(f)] when ¯φ = 0.3 has 18bumps. However, from Fig. 4(c) we observe that the energeticminimum is obtained by a periodic profile with 15 bumps, asshown by the steady state solution in Fig. 5. The convergenceto a different number of bumps in the time simulation may becaused by discretization effects or by the initial noise profileused. As one would expect, the free energy associated withthe 18-bump periodic structure is significantly larger than the15-bump profile.In Fig. 3 ( ¯φ = 0.1) the localized states bifurcate subcrit-ically from the periodic solution branch (that itself emergesfrom the trivial homogeneous solution that is displayedas the heavy black dotted line). Therefore, one expectshysteretic behavior as encountered in the time simulations. Amagnification (not shown) allows us to determine the thresholdvalues for the hysteretic transition. When decreasing r in theregion where periodic solutions are always found, one firstpasses rsn = −0.685 where the last two branches of localizedsolutions annihilate in a saddle-node bifurcation [Fig. 3(a)].Slightly below rsn, both the periodic solution and the localizedstate with a single bump are local energetic minima. Althoughthe periodic solution represents the globalminimum, particulartime simulations sometimes converge to the localized state.061408-6MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012)(a)0 0.1 0.2 0.3φ00.050.10.150.2||δφ||n=15n=16n=140 0.1 0.2 0.3φ-0.05-0.04-0.03-0.02-0.010μn=15n=16n=14(b)(c)0 0.1 0.2 0.3φ-0.015-0.01-0.0050(F-F0 )/Ln=15n=16n=140 0.1 0.2 0.3φ-0.008-0.006-0.004-0.00200.002ωn=15n=16n=14(d)FIG. 4. (Color online) Bifurcation diagram showing localized solutions of the VPFC as a function of the mean order parameter ¯φ, forr = −0.9 and a fixed domain size of L = 100. The various solution profiles are characterized by their (a) L2 norm, (b) chemical potential μ,(c) mean free energy density (F − F0)/L, and (d) mean grand potential density ω ≡ F/L − ¯φμ. The line styles are as in Fig. 3. Here, however,the heavy solid black and dashed red lines bifurcate at large ¯φ from the n = 15 and n = 14 periodic solutions, respectively. Typical profiles forall the branches of localized states are given in Fig. 5. The vertical dotted lines in (a) correspond to values of ¯φ for results in Fig. 5. The bluedots correspond to the five time simulation profiles shown in Figs. 2(b)–2(f).The differences in energy between the two is <1% in the caseof Fig. 3. When r is further decreased below ren = −0.700the energy of the even symmetric states becomes smaller thanthe one of the n = 16 periodic solution, which is, however,still linearly stable. The situation changes at rc = −0.749,where both symmetric localized branches bifurcate from then = 16 branch; that is, below rc the latter is linearly unstable.Furthermore, below rc the energy of all localized states rapidlybecomes much smaller than the energy of all periodic states[Fig. 3(c)]. The hysteresis range displayed in Fig. 2 providesa good approximation for the region between rc and rsn. Thisregion becomes larger as ¯φ is increased.The situation is very similar when ¯φ is changed for fixed r(Fig. 4). The resulting hysteresis range is between ¯φ = 0.150and 0.239 for symmetric localized states with an odd numberof maxima and between ¯φ = 0.202 and 0.265 for symmetricstates with an even number of maxima. Overall, one shouldtherefore expect a wide hysteresis region roughly between¯φ = 0.15 and ¯φ = 0.25. The hysteresis range obtained fromthe time simulations (indicated in Fig. 2) is roughly 0.19 <¯φ < 0.22. This is narrower than the range deduced from thepath continuation analysis of the localized steady states, butlies right in the middle of it.Before we move on to discuss the 2D case, we should com-ment on howour results fit into thewider context of research onlocalized states. Much research on localized states focuses onthe nonconserved Swift-Hohenberg equation [19–21]. There,such states can exist only if the primary bifurcation of periodicstates from the homogeneous base state is subcritical. Thelocalized states exist in a subrange of the existence rangeof the periodic states bounded on either side by the saddle-node bifurcations of the branches of symmetric localizedstates. In the nonconserved Swift-Hohenberg equation theseaccumulate exponentially rapidly toward the parameter valuescorresponding to first and last tangencies between the unstablemanifold of the homogeneous state in space and the stablemanifold of the periodic state. These tangencies define thepinning region containing the different localized structures. Incontrast, in the presence of a conserved quantity localizedstates may exist outside the existence region of periodicstates and may occur even in the supercritical case, and thesaddle-node bifurcations of the localized states are no longeraligned; that is, one finds slanted snaking [25]. This is typicallya finite size effect [26].For the regular PFC (conserved Swift-Hohenberg equation)[Eq. (1) with Eqs. (2) and (3)], localized states are brieflymentioned in Ref. [8]. However, no systematic results alongthe lines of those presented in Refs. [19] for nonconservedor [25] for conserved order parameter fields are available.The model used here is a special case because it includes the061408-7ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)00.20.40.60.800.20.40.600.20.40.600.20.40.6φ(x)00.20.40.60 0.2 0.4 0.6 0.8 100.20.40.60.2 0.4 0.6 0.8 1 x/L0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 10.01 0.05 0.05 0.10.1 0.1 0.2 0.20.2 0.01 0.05 0.10.1 0.1 0.2 0.20.2 0.05 0.1 0.10.1 0.2 0.2 0.3FIG. 5. (Color online) A selection of steady state profiles φ(x) for r = −0.9 at ¯φ = 0.01, 0.05, 0.1, and 0.2. From top left to bottom rightwe show first nine type (i) solutions, that is, symmetric localized states with an odd number of maxima (in black), then eight type (ii) solutions,that is, symmetric localized states with an even number of maxima (in red), followed by six type (iii) solutions, that is, asymmetric localizedstates (in green). The final image is the n = 15 periodic solution at ¯φ = 0.3 (in blue). The number in each panel indicates the correspondingvalue of ¯φ. The solutions from the symmetric branches are shown in the sequence that follows the respective branch in Fig. 4(a), starting fromthe left. The asymmetric states for identical ¯φ are shown in the order of decreasing norm.nonanalytic vacancy term (5). However, a similar bifurcationstructure is found for the classical conserved Swift-Hohenbergequation, that is, the regular PFCmodel.Wewill present resultsfrom our study of localized states in the regular PFC modelelsewhere [27].C. Two-dimensional modelWe now move on to consider how the VPFC model behavesin 2D. As with the regular PFC model [2], when we expandinto 2D we observe stripes [see Fig. 6(b)] and hexagonallyordered bumps or holes [see Fig. 6(c)]. In Fig. 6(a) wedisplay the phase diagram of the VPFC model in 2D andtypical time simulation results from the striped [Fig. 6(b)]and hole [Fig. 6(c)] phases, calculated on a regular gridwith grid spacing dx = 0.5. Square ordering of bumps orholes does not appear in the phase diagram because thesestructures always have a higher free energy. However, thiscan be changed through appropriate alterations to the freeenergy [5]. Square ordering can also occur when extendingto a two-component mixture (cf. Sec. III C below). Using thesame method as outlined above, we calculate the regions of thephase diagramwhere there is coexistence between hexagonallyordered holes and the uniform distribution, between holes andstripes, and between stripes and hexagonally ordered bumps.The vacancy term (5) shifts the modulated phases into thepositive ¯φ > 0 plane. The section of the phase diagram whereholes are observed is much smaller when compared to theregular PFC model and now extends beyond the limit oflinear stability of the flat state (at = 0). This means that forcertain values of ¯φ (where 0 < 1), hexagonally arrangedholes are energetically favorable but are only observed in timesimulations for certain initial conditions; that is, when startingwith an order parameter profile φ(x,t = 0) which already hasmodulations which are sufficiently large in amplitude. Asr is decreased (i.e., for larger |r|) it becomes increasinglydifficult to obtain structures with holes up to and inside of thecoexistence region between the hole and the uniform phases.This is a direct consequence of the limit of linear stabilityoccurring in the middle of the hole phase. Therefore, theaccuracy of results for the coexistence region between the holeand uniform phases decreases as |r| becomes larger. Thestripe phase occurs between the two hexagonal phases.In the order parameter profiles displayed in Figs. 6(b) and 6(c)we observe various defects and in Fig. 6(c) “grain” boundariesbetween regions with different orientations, which depend onthe initial conditions (our initial profile was a flat state withadditional small amplitude white noise). The true minimumprofile for case (b) is a series of parallel stripes which areidentical to the periodic profiles in the 1D system [shown inFig. 2(f)].061408-8MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012)(c)(a)(b)(a)FIG. 6. (Color online) The phase diagram for the 2D VPFC model [Eqs. (1) and (4)] is displayed in (a) for the case q = 1. The red solid linesare the coexistence curves between the various phases. The blue dashed line is the limit of linear stability for uniform profiles = 0. The greendash-dotted lines indicate the region where localized and hexagonally ordered bump structures coexist. Simulations of (b) stripes and (c) hexag-onally ordered holes are also shown. The parameter values for these simulations are q = 1, r = −0.9, α = 1 and (b) ¯φ = 0.4 and (c) ¯φ = 0.53.The most important portion of the phase diagram from thematerials modeling point of view, is the bump phase becausethe basic assumption is that each bump represents a particle.When r −0.4 or when ¯φ has a value close to that in thecoexistence region between bumps and stripes, we observehexagonally arranged bumps, similar to those in the regularPFC model. However, in a manner similar to the 1D system,we observe localized structures at small values of ¯φ when r −0.4. In the phase diagram [Fig. 6(a)] the green dot-dashedlines are numerically obtained estimates for the location in thephase diagram of the limits of linear stability of the uniformperiodic states (lower curve) and the localized (vacancy) states(upper curve). They are determined in the same manner asdiscussed above for the 1D system for a square system ofside length L = 25. It is important to note that the parameterrange where localized bumps coexist with regular periodicallyordered bumps is much broader for the 2D system, implying alarge amount of hysteresis.We now focus our discussion on the portion of the phasediagram where isolated bumps form. As can be seen in Fig. 7,these profiles resemble particle configurations in gases, liquidsand crystalline solids and so the VPFC may be a valuablemodel for describing materials on the microscale [13]. Thisregion of the phase diagram is full of complexity and manyvaried structures may be observed. However, here we foregoa full systematic study of this large region in parameter spaceand limit ourselves to showing representative results obtainedfor a single value of the undercooling parameter r = −0.9 forwhich there is a fairly large range in ¯φ with isolated bumps. Weset the initial order parameter profile to be a uniform state witha small amplitude noise φ(x,t = 0) = ¯φ + λ(Y − 0.5), whereY is a random real number uniformly distributed between 0and 1 and λ = 10−6 is the amplitude of the noise. We considerthree cases: (i) ¯φ = 0.01 and (ii) ¯φ = 0.1, where a disorderedarrangement of localized bumps forms, and (iii) ¯φ = 0.24,which is in the region where bumps are hexagonally ordered.We average over many simulations to calculate the two-pointcorrelation function for each of these cases. This is done bylocating all the maxima in the equilibrium profile φ(x), for agiven initial realization of the noise; that is, we locate theposition of all the bumps. From these sets of coordinateswe calculate the radial distribution function g(x) in the usualway [28].We display a simulation result for case (i) in Fig. 7(a)and the corresponding radial distribution function g(x) inFig. 7(d). We find that there is almost no correlation betweenthe bumps in this circumstance except for the core repulsionand a very small peak at x ≈ 16, indicating that there is aweak attraction between the bumps. Therefore, simulationswith these parameter values appear to qualitatively describegaslike formations of particles or colloids. In Figs. 7(b)and 7(e) we plot a typical order parameter profile and thecorresponding g(x) for case (ii). We observe a large increasein the number of bumps as compared to the previous case. Theradial distribution function shows that we have strong shortrange ordering, but without any long range order. This is veryreminiscent of the ordering in liquids. There is a very sharppeak in g(x) at around x = 7.5 (which is approximately thediameter of the bumps) and a smaller peak around x = 15.A similar example is also given in Ref. [13]. If we furtherincrease the value of ¯φ we eventually find the more familiarhexagonally structured array of bumps which is reminiscentof the ordering in simple crystalline solids. In Fig. 7(c) wedisplay an example of the order parameter profile for case(iii) and in Fig. 7(f) we show the corresponding g(x). For061408-9ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012) 1 3 5 0 15 30 45g(x)x(d) 1 3 5 0 15 30 45g(x)x(e) 1 3 5 0 15 30 45g(x)x(f)(a) (b) (c)FIG. 7. (Color online) (a)–(c) Typical steady state order parameter profiles obtained in time simulations for increasing ¯φ. (d)–(f) Thecorresponding radial distribution function g(x) calculated from multiple simulations. The parameter values are α = 1, q = 1, r = −0.9 and in(a),(d) ¯φ = 0.01, in (b),(e) ¯φ = 0.1, and in (c),(f) ¯φ = 0.24.this case we observe that g(x) is highly structured, indicatingthe system has very strong short range correlations with asignificant degree of long range ordering. We observe the splitsecond and third peaks, which is a classic sign of crystallineorder. These results indicate that the VPFC model may beused to model crystalline structures, much like the regular PFCmodel. The major difference between the two models is theexistence of the fluidlike configuration of bumps observable inthe VPFC model. In contrast, the fluid phase in the PFC modelcorresponds to the homogeneous state.The variation in the size and shape of the bumps that areformed is fairly small. In Fig. 8(a) we display a selectionof results for the order parameter profile through the centerof the bumps for the case when r = −0.9 and ¯φ = 0.01. Wedetermine the shape of the bumps by plotting the value of theorder parameter φ against the distance from the peak of eachbump (as shown by the data points). We can then fit functionswhich take the following form:θ (x) = β0e−β1x2−β2x4−β3x6 cos(β4x) + β5. (21)We fit this form to the data using a least squaresmethod. Theexponential part of θ (x) describes the decay of the modulationas the distance from the peak increases and the cosine functioncaptures the oscillatory tail of the modulations, which is animportant factor in their interaction with other bumps [29,30].Figure 8(a) displays two cases: the + points and red solid lineshow the case q = 1 and the × points and blue dashed lineshow the case where q = 1.1. The size of the bump is reducedas we increase the value of q. This is because increasing thevalue of q increases the typical wave number, which results ina smaller typical length scale [cf. Fig. 1 and Eq. (10)].The curves obtained from fitting the bump profile can beused to obtain an approximation for the effective pair potentialV (x) between two isolated bumps, where x is the distancebetween the centers of the bumps. We take a uniform systemwith the value of φ equal to that in the uniform areas betweenbumps found in simulations for ¯φ = 0.01, corresponding tothe results in Fig. 7(a). We then impose upon this the profilesfor two bumps using the fitted curves shown in Fig. 8(a). Wevary the distance between the superposed bumps and calculatethe free energy of the system. We assume thereby that the twobumps retain their shape when they are close, despite the factthat, in reality, the bump shapes become distorted as bumpsare pushed close together.In Fig. 8(b) we display the results for q = 1 (red solidline) and q = 1.1 (blue dashed line). We observe that thereis a shallow minimum in the potential at the distance x ≈ 7.5when q = 1 and at x ≈ 7when q = 1.1 [see inset of Fig. 8(b)].The minimum is at a smaller distance when q is larger becauseof the decreased diameter of the bumps; recall that q deter-mines the size of the bumps. The resulting weak attractionbetween the bumps may also be inferred from the radialdistribution function g(x) calculated for the low density case¯φ = 0.01 when q = 1 displayed in Fig. 7(d). We observe asecond minimum in the potentials at x ≈ 3.15 when q = 1and at x ≈ 2.3 when q = 1.1, where the former rather appearslike a “shoulder.” The order parameter profiles at the secondminima resemble the elongated almost elliptical shapes whichare observed in and around the coexistence region betweenbumps and stripes. See also Fig. 6(c) where we also observeelliptical holes along some of the grain boundaries.III. TWO-COMPONENT SYSTEMWe now extend the model to consider a binary mixture inperhaps the most simple way possible, by adding together freeenergies like Eq. (4) for two order parameter fields φa(x,t) andφb(x,t). We introduce a simple coupling term which allows061408-10MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012) 0 0.2 0.4 0.6-4 -2 0 2 4 6x(a)q = 1.0q = 1.1 0 8 16 0 4 8 12V(x)x(b)q = 1.0q = 1.1-0.03 0 0.03 6 8 10FIG. 8. (Color online) (a) Several sets of numerical results for the order parameter profile through the center of a bump (+) for q = 1 and(×) for q = 1.1, together with fits to the data displayed as solid red and dashed blue lines. These fits are then used to calculate the effective pairpotential between two bumps. These pair potentials are displayed in (b). The inset displays a magnification of the tails of V (x). The parametervalues are r = −0.9 and ¯φ = 0.01.the two components to interact with each other. This gives usthe following expression for the free energy:F =∫dx[f (φa(x,t)) + fvac(φa(x,t)) + f (φb(x,t))+ fvac(φb(x,t)) + ηφaφb], (22)where η is the coupling coefficient and the functions f andfvac are defined as before in Eqs. (3) and (5). The value ofr is set equal for both components. However, we allow thevalue of q to be different for each species, so we now referto these values as qa and qb, where the subscript denotes thecorresponding component. Setting different values for q inthe two components (i.e., qa = qb) results in an asymmetricsystem in which the size of the bumps or modulations inφa(x) differs from that in φb(x), as discussed further belowin Secs. III B and III C. Note that a different coupling termis used in Ref. [14]; a somewhat different two-componentPFC model is presented in Ref. [31]. In our two-componentVPFC model one may interpret the bumps in the profile φa(x)as representing particles of species a and the bumps in theother profile φb(x) as representing particles of species b. Thereare also vacant “vacancy” areas, where both order parameterprofiles take values close to zero.Here, just as for the one-component model, we assume thedynamics of the system is governed by the following pair ofequations [cf. Eq. (1)]:∂φa∂t= αa∇2 δFδφa,∂φb∂t= αb∇2 δFδφb. (23)We also assume that the two mobility coefficients are equal:αa = αb = α. The two components are coupled purely by theterm ηφaφb in the free energy. When η > 0, this couplingterm leads to a repulsion between the two species and sopenalizes structures which overlap or form on top of eachother. The value of the parameter η determines the “strength”of the coupling, and so the two-component model reduces totwo disconnected one-component models in the limit η → 0.A. Phase behaviorWhen the coupling coefficient is fairly large, η 0.1, thecoupling term has a significant impact on the phase behaviorof the model. In particular, the limit of linear stability andthe phase coexistence curves extend to much larger valuesof ¯φ = ¯φa + ¯φb than for the one-component model. We nowdetermine the linear stability of a flat state in the model. Weassume that the order parameter profiles of both componentstake the formφa = ¯φa + δφ = ¯φa + ξeikxeβt , (24)φb = ¯φb + χδφ = ¯φb + χξeikxeβt ,where the amplitude |ξ | 1 and the parameter χ is theratio between the amplitude of the modulations in the twocomponents. The sign ofχ indicateswhether instabilities are inphase (χ > 0) or antiphase (χ < 0) between the two coupledorder parameter fields. From themagnitude ofχ we can deducewhether the instability is initiated from species a (|χ | 1),species b(|χ | 1), or a combination of both (|χ | = O(1)).Wemake a Taylor series expansion of the functional derivativesof the free energy with respect to the two order parameters φaand φb, to obtainδFδφa= (r + q4a) ¯φa + 3H ¯φa(| ¯φa| − ¯φa) + ¯φ3a + η ¯φb+ [(k2 − q2a)2 + a + χη]δφ + O(δφ2),δFδφb= (r + q4b) ¯φb + 3H ¯φb(| ¯φb| − ¯φb) + ¯φ3b + η ¯φa+ [χ(k2 − q2b)2 + χ b + η]δφ + O(δφ2), (25)where is defined as before Eq. (8) and the subscriptdenotes the corresponding component. We substitute theseexpressions into the dynamical equations (23), yielding thematrix problem [32,33]β(1χ)= M(1χ), (26)061408-11ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)whereM = −k2α((q2a − k2)2 + a ηη(q2b − k2)2 + b).We can now determine the dispersion relation β(k) bycalculating the eigenvalues of M:β(k) = Tr(M)2±√Tr(M)24− |M|. (27)The resulting dispersion relation β(k) is a double-valuedfunction. However, since the growth rate along the + branchis always larger than that along the − branch, the limit oflinear stability can be determined from the + branch alone. Ifwe assume that qa = qb = q, the dispersion relation simplifiessignificantly, yieldingβ(k) = −αk22[2(k2 − q2)2 + a + b−√( a − b)2 + 4η2]. (28)There is a local maximum of this expression which occurs atthe positive wave number:km = 16[24q2 + 6(4q4 − 6( a + b)+ 6√( a − b)2 + 4η2) 12] 12. (29)Substituting this wave number back into the dispersionrelation (28), allows us to calculate the parameter values suchthat β(km) = 0 (i.e., the limit of linear stability of a flat state).We arrive at the following relation: a b = η2. (30)When the system is linearly unstable it is possible for β(k =0) to be a minimum or maximum (this transition occurs at = −q2 in the one-component model). This is equivalent tothe coefficient of k2 changing from a positive value (minimum)to a negative value (maximum). The sign of the coefficient ofk2 is determined by the sign of the following quantity:C2 = ∂2g∂φ2a∂2g∂φ2b−(∂2g∂φb∂φa)2,= (q4a + a)(q4b + b)− η2, (31)where g(φa,φb) = f (φa) + fvac(φa) + f (φb) + fvac(φb) +ηφaφb. When C2 is negative (positive) β(k = 0) is a minimum(maximum); this relation also holds for asymmetric systemswhere qa = qb. In Fig. 9 we display typical dispersionrelations when qa = qb = 1, a = b = , and η = 4. Weshow the case when (i) the system is linearly unstable andC2 [Eq. (31)] is negative (red solid line), (ii) the system islinearly unstable and Eq. (31) is positive (blue dashed line),(iii) the system is at the limit of linear stability [i.e., Eq. (30)holds] (green dotted line), and (iv) the system is linearly stable(magenta dash-dotted line). We observe that when qa = qb,the typical wave number km → qa as we approach the limitof linear stability a b − η2 → 0. In the more general case 0 0kq = 2.8 = 3.5 = 4.0 = 5.0FIG. 9. (Color online) Dispersion curves for the two-componentVPFCmodel [Eqs. (22) and (23)], when qa = qb = 1, a = b = ,and η = 4. Four cases are shown: (i) β(km) > 0 with β(k = 0) aminimum (red solid line), (ii) β(km) > 0 and β(k = 0) a maximum(blue dashed line), (iii) β(km) = 0 (green dotted line), and (iv)β(km) < 0 (magenta dash-dotted line).with qa = qb the dispersion relation may have two maxima atpositive values of k, neither of which occurs at qa and qb. Inthis case the stability boundary is defined by the vanishing ofthe growth rate β(km) = 0 of the larger of the two possiblemaxima of β.From Eq. (30) it is clear that depending on the value of thecoupling coefficient η, the region of parameter space wherethe system is linearly unstable can be greatly larger thanthat for the one-component system. For example, picking thevalue η = 4 when r = −0.9 and setting the average valueof both order parameters to be equal ¯φa = ¯φb = ¯φ, we findthat the limit of linear stability increases from ¯φ = 0.548(for the one-component case) to ¯φ = 1.278. As one wouldexpect, this also increases the region of the phase diagramwhere modulated structures are formed. Our focus here ison the regions of parameter space where bumps are formedas this is the regime relevant to modeling crystalline solids.However, before proceeding to this, we make a brief surveyof some of the structures which may be observed for largervalues of φa and φb which lie outside of the bump phase.For the parameter values r = −0.9, η = 4, qa = qb = 1, and¯φa = ¯φb = ¯φ we show in Fig. 10 a sequence of order parameterprofiles with increasing ¯φ, for values of ¯φ that lie abovethe region where bumps are observed (see later sectionsfor a detailed analysis of the bump structures found in thetwo-component model). In Fig. 10 we display scaled plots oforder parameter profiles which are stationary states obtainedfrom time simulations. We plot an order parameter defined asthe normalized difference between the φi(x) values of the twocomponents ˆφ(x) ≡ φa(x)/ ˆφa − φb(x)/ ˆφb, whereˆφi = φmaxa φmaxb − φmina φminbφmaxi + φmini(32)and where φmaxi and φmini are, respectively, the maximumand minimum values of φi(x). ˆφ is defined so as to takea value in the range [−1,1]. When ˆφ ≈ +1, then the localvalue of φa is high while the value of φb is low. Conversely,when ˆφ ≈ −1, then the local φa is low and φb is high. The061408-12MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012)FIG. 10. (Color online) The scaled order parameter ˆφ for the two-component model, corresponding to minima of the free energy. Thepeaks in species a are shown in orange (light gray), peaks in species b are shown as blue (dark gray), and white areas show regions whereφa ≈ φb. The parameter values are η = 4, r = −0.9, qa = qb = 1, ¯φa = ¯φb = ¯φ, where (a) ¯φ = 0.25, (b) ¯φ = 0.3, (c) ¯φ = 1, (d) ¯φ = 1.15, (e)¯φ = 1.2, and (f) ¯φ = 1.27.average order parameter values in Fig. 10 are (a) ¯φ = 0.25, (b)¯φ = 0.3, (c) ¯φ = 1, (d) ¯φ = 1.15, (e) ¯φ = 1.2, (f) ¯φ = 1.27.The most palpable change from the one-component model isthat the phase diagram is largely dominated by the stripedprofiles, with stripes appearing in the range 0.22 ¯φ 1.28.Just outside the range of ¯φ where bump structures are formed,we observe order parameter profiles which contain a mixtureof bumps and stripes [see Fig. 10(a)]; this value of ¯φ mustlie inside the coexistence region between the bump and stripephases. As we increase the value of ¯φ we enter the largeregion of parameter space where stripe structures are formed[Figs. 10(b) and 10(c)]. The only significant change as weincrease ¯φ from 0.3 to 1 is the decrease in the width of thestripes; this is due to the fact that the typical length scale in thesystem is 2π/km, where the wave number km given by Eq. (29)is inversely proportional to the average order parameter values¯φa and ¯φb. Increasing the value of ¯φ further, we continue toobserve striped profiles, but now there are points where thestripes of one species “connect” to stripes of the other species[see Fig. 10(d), where these “connections” appear as whitelines]. Increasing ¯φ further, we observe a mixture of holesand stripes [Fig. 10(e)]. Close to the instability curve (30) wefind interesting profiles where we observe a mixture of stripes,holes, and regions where the profile is approximately uniformφa ≈ φb ≈ ¯φ [Fig. 10(e)]. Various modulated structures areobserved over a large range of parameter values. It wouldbe possible to consider the structures formed for differentvalues of the coupling coefficient η and different values ofthe average order parameters, where ¯φa = ¯φb. However, herewe do not make a systematic study of the entire parameterspace and instead focus on the various bump configurations.These structures closely resemble the configurations of parti-cles or colloids in condensed matter systems and we believethat in this regime the model may be useful in understandingthe fluid and solid phases of such systems.B. Intermolecular interactionsFor the remainder of this paper, we pursue the idea that thebumps in this two-component model represent two differenttypes of molecules or colloidal particles suspended in a fluidmedium. We perform time simulations of the two-componentmodel choosing parameter values which result in the formationof bump structures. We run these simulations until the orderparameter profiles reach an (almost) stationary state, whichcorresponds to being at (near) an energetic minimum. Wethen determine the coordinates of the particles by locatingthe position of the maximum of each of the peaks. Theradial distribution functions are calculated by analyzing thesecoordinates. We also calculate the effective pair potentialsbetween the bumps. Later in Sec. III C we consider thenearest neighbor bond angles and the ordering in crystallineconfigurations.The bump phase in the two-component model appears tobehave in a manner similar to that of the one-componentmodel (cf. Fig. 6). We observe bump structures when theaverage value of the order parameters ¯φa and ¯φb are small.In particular, when ¯φa = ¯φb = ¯φ and r = −0.9 we observebumps within the range 0 ¯φ 0.15. We study and comparetwo different systems: the symmetric case where qa = qb = 1and the asymmetric case where qa = 1 and qb = 1.1. In thesymmetric case, interactions between bumps of the same type(aa and bb) are identical in both components, but the natureof the interaction between a bump in φa and a bump in φb (ab)061408-13ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)is determined by the coupling term in the free energy. In theasymmetric case, the different values of q mean that the sizeof the bumps are different in φa and φb; hence, all possibleinteractions aa, bb, and ab are different.We begin by considering how the two order parameterprofiles change as we alter their average values ¯φa = ¯φb = ¯φ.Thus, we keep the concentration of the mixture fixed atc = 0.5, wherec =¯φa¯φa + ¯φb . (33)We set the other parameter values to α = 1, r = −0.9, andη = 4. In Fig. 11 we display typical results. We plot thenormalized difference between the two order parameters ˆφ[as defined in Eq. (32)]. In the left column we show profilesfrom the symmetric case and in the middle column wedisplay the profiles from the asymmetric system. In the rightcolumn we present the radial distribution functions, which areobtained by averaging over at least 50 runs, each with differentrealizations of the initial noise. The solid lines show the radialdistribution functions for the symmetric case and the dashedlines show the asymmetric case. It is very apparent that thisregion of the parameter space shares many similarities withthe one-component model in both 1D and 2D. If we selecta small value of ¯φ we find localized peaks surrounded byvacant areas, as shown in Figs. 11(a) and 11(b). We observea tendency for bumps in φa and in φb to sit pairwise next toeach other resembling configurations occurring in mixtures ofoppositely charged colloidal particles [34–36].When ¯φa ≈ ¯φb,the arrangement of the bumps also resembles snapshots ofmonovalent salts. It is very difficult to differentiate between thestructures formed by the symmetric and asymmetric modelsfor small values of ¯φ. This is because, structurally, there isvery little difference between the two cases. If we examinethe radial distribution functions gij (where i,j = a,b) for thesymmetric and asymmetric systems [Fig. 11(c)] we observethat the average distance between the different bumps seemsto be independent of the q values (any difference betweenthe curves is of the same order of magnitude as the statisticalerror). This is due to the large vacant areas, which means thatthere are not many bumps which are close to one another,especially between bumps of the same species (aa and bb).As we increase the values of ¯φ we find that the number ofbumps of both species increases. In Fig. 11(d)–11(f) we showthe case where ¯φ = 0.04 and in Fig. 11(g)–11(i) we show thecase where ¯φ = 0.06. There is now a clear difference betweenthe symmetric [Figs. 11(d) and 11(g)] and the asymmetric[Figs. 11(e) and 11(h)] cases. We observe a larger number ofbumps in φb when qb = 1.1. This is because the larger valueof q reduces the length scale of the modulations, meaning thatmore bumps can be created before the value of φb becomessmall (and negative) in vacant areas. There is an optimumvalueof φa and φb in the vacant (uniform) areas which depends onthe parameter values. This explains why increasing the valueof ¯φ increases the number of bumps (i.e., moremodulations areneeded in order to reach the optimum value of φ in the vacantregions). These intermediate values of ¯φ produce profileswith bump configurations that resemble real fluid structures.However, in stark contrast to the one-component system[shown in Fig. 7(b)], we now find the formation of chainsof alternating bumps reminiscent of structures observed incharged fluids. The radial distribution functions in Figs. 11(f)and 11(i) show that the asymmetry induced by the differentvalues of q begins to take effect at these intermediate values of¯φ. We observe that statistically the bumps sit closer togetherin the asymmetric case, especially when two bumps in φb arenext to each other (bb, shown by green dashed line). This isdue to the decreased size of the bumps in φb, allowing them tosit slightly closer to their neighbors.Increasing the average order parameter values ¯φ further,we begin to observe the formation of crystalline structuresas shown by Figs. 11(j)–11(l). The interesting thing is thatnow we observe square ordering of the particles instead ofthe hexagonal ordering which is present in the regular PFCmodel and the one-component VPFC model. This implies thatas we increase the concentration of one of the species fromc = 0 (almost a pure one-component system) to c = 0.5 theremust be a transition from hexagonal to square ordering ofbumps; this is something we return to below in Sec. III C.Just as for the one-component system, we find that there aremore modulations in φb when qb = 1.1. The profiles obtainedwith these parameter values resemble a compound crystalstructure with vacancies and grain boundaries. The radialdistribution functions in Fig. 11(l) show that the smaller sizedifference of the φb bumps in the asymmetric mixture hasa large impact on the average position of the bumps in thestructure compared to the symmetric mixture. This is becausethe higher concentration of particles forces them all closertogether, resulting in all pairs of bumps aa, bb, and ab beingcloser together.In Figs. 12(a) and 12(b) we show the shape of the individualbumps inφa andφb obtained in the low density limit ¯φ → 0. Todetermine these radially symmetric profiles we fit functions ofthe form θ (x) as defined above in Eq. (21). The bumps in φa arevirtually identical for both the symmetric and the asymmetricsystems. The φa bump in the symmetric system and the φbbump in the asymmetric system decay to different values dueto the different values of φa and φb in the vacant areas of theasymmetric system. We observe that in this two-componentmodel, a bump in one order parameter profile coincides with asmall depression in the other order parameter profile. This iscaused by the coupling term,whichmeans that the combinationof a bump in one order parameter and a hole in the other orderparameter reduces the free energy of the system. In Fig. 12(c)and 12(d) we show the shape of the “holes” which form inone profile under the bumps in the other order parameter field.These are determined the same way as the bump profiles:by fitting a function of the form shown in Eq. (21) to datapoints obtained from simulations. The depth of the holes ismuch smaller than the height of the bumps. This is becausethe vacancy term prevents large negative values of φa or φb.Using the fitted functions shown in Figs. 12(a)–12(d) wecalculate effective pair potentials Vij (x) between the differentparticles in the system (i,j = a,b). We do this by determiningthe free energy for a system containing two bumps and theircorresponding holes at various distances apart. In Fig. 12(e)we display the effective pair potentials for both the symmetric(solid lines) and the asymmetric (dashed lines) systems. Theresults show that there is an attraction between all of the bumps,just as we found for the one-component system [Fig. 8(b)]. In061408-14MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012) 1 5 9 0 15 30g ij(x)x(c)ab, qb = 1.0aa, qb = 1.0ab, qb = 1.1aa, qb = 1.1bb, qb = 1.1 1 5 9 0 15 30g ij(x)x(f)ab, qb = 1.0aa, qb = 1.0ab, qb = 1.1aa, qb = 1.1bb, qb = 1.1 1 5 9 0 15 30g ij(x)x(i)ab, qb = 1.0aa, qb = 1.0ab, qb = 1.1aa, qb = 1.1bb, qb = 1.1 1 5 9 0 15 30g ij(x)x(l)ab, qb = 1.0aa, qb = 1.0ab, qb = 1.1aa, qb = 1.1bb, qb = 1.1(a) (b)(d) (e)(h)(j) (k)(g)FIG. 11. (Color online) In the left-hand column we display typical simulation results for ˆφ for the symmetric case where qa = qb = 1. Inthe middle column we display results from the asymmetric case where qa = 1 and qb = 1.1. The orange (light gray) regions show where thereis a φa bump, while the blue (dark gray) show the φb bumps which are slightly smaller in the asymmetric mixture. In the right-hand columnthe radial distribution functions gij (x) are shown for the symmetric case (solid lines) and the asymmetric case (dashed lines). The parametervalues are α = 1, η = 4, r = −0.9, ¯φa = ¯φb = ¯φ, where (a)–(c) ¯φ = 0, (d)–(f) ¯φ = 0.04, (g)–(i) ¯φ = 0.06, and (j)–(l) ¯φ = 0.15.both the symmetric and asymmetric cases we find that theattraction between two bumps from different species (ab)is stronger and occurs at a smaller value of x than that oftwo bumps of the same species (aa and bb). This explainsthe tendency for the bumps to form chains at intermediatevalues of φa and φb [Figs. 11(d), 11(e), 11(g), and 11(h)]and square ordered crystalline structures at larger values ofφa and φb [Figs. 11(j) and 11(k)]. This is also consistent withthe appearance of the large peak in gab(x), which occurs at asmaller x value than the main peaks in gaa(x) and gbb(x) [seeFigs. 11(c), 11(f), 11(i), and 11(l)]. The effective pair potentialVaa(x) is almost identical in the symmetric and the asymmetric061408-15ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012) 0 0.2 0.4 0.6-4 -2 0 2 4x(a)ab 0 0.2 0.4 0.6-4 -2 0 2 4x(b)ab-0.02-0.015-0.01-0.005-4 -2 0 2 4x(c)ab-0.02-0.015-0.01-0.005-4 -2 0 2 4x(d)ab 0 30 60 90 0 4 8 12V ij(x)x(e)ab, qb=1.0aa, qb=1.0ab, qb=1.1bb, qb=1.1-0.05 0 5 7 9FIG. 12. (Color online) Fits to the shape of individual bumps (cf. Fig. 8) in the (a) symmetric, when qb = 1, and (b) asymmetric, whenqb = 1.1, systems and the corresponding dips in the other order parameter profile which occur under the bumps in (c) the symmetric and (d)the asymmetric systems. In the symmetric cases (a) and (c) the φa and φb curves are identical everywhere. The bump profile in φa is virtuallyidentical in the cases where qb = 1 and qb = 1.1. These fits are then used to calculate the effective pair potential between bumps, which aredisplayed in (e). The inset displays a magnification of the tails of Vij (x). The resulting pair potential V (x) between two bumps in φa whenqb = 1.1 lies on top of the aa, qb = 1 curve. The parameter values are α = 1, r = −0.9, η = 4, and ¯φa = ¯φb = 0.systems. This suggests that the small hole which appears inφb has little effect on the interaction between the bumps.The major difference between the symmetric and asymmetricsystems is that in the asymmetric mixture the minima of thepair potentials Vab(x) and Vbb(x) are at smaller values of xthan in the symmetric mixture. This is due to the reduced sizeof the φb bumps in the latter. The minimum in Vbb(x) is at aslightly larger value of x than the minimum in Vaa(x) and theattraction is also much weaker (in fact, it is so much weakerthat the minimum is barely visible in this plot). This to someextent explains why the effect of the asymmetry is not visiblefor smaller values of φa and φb, but becomes apparent forlarger values of φa and φb, where the vacant areas becomesmaller and we observe a close packing of the particles.C. Bond angles and the transition between hexagonaland square orderingIn the 2D one-component model [Eqs. (1) and (4)] weobserve hexagonally ordered structures for certain parametervalues [Figs. 6(a) and 7(c)]. However, in the two-componentmodel when ¯φa = ¯φb, we instead observe a square orderedcrystalline structure which alternates between species a andspecies b [Figs. 11(j) and 11(k)]. Thus, as the composition ofthe mixture is varied we should see a transition or crossoverfrom hexagonal to square ordering. The number of bumpsobserved in each field φi depends on the respective averagevalue ¯φi . When the concentration c ≈ 0 or c ≈ 1, where c isdefined in Eq. (33), (i.e., when either ¯φb ¯φa or ¯φa ¯φb),then the resulting order parameter profile ˆφ(x) has manymore bumps of one type than of the other, and in thesetwo limits we again observe hexagonal ordering. Note thatc in Eq. (33) is not a bump concentration, but instead is aratio between the two average order parameter values. Asthe φi may take a negative value, for c = 0 there are stilla few bumps of a and similarly there are still some speciesb bumps when c = 1. When c = 0.5 the number of bumpsis roughly the same in both species for the symmetric case(qa = qb), but this is not necessarily true for the asymmetricsystem (qa = qb). When ¯φa = ¯φb and qa < qb there are moreb bumps than a bumps.In Figs. 13(a)–13(c) we show the order parameter ˆφ forvarying values of c for the symmetric mixture (qa = qb). Wefix the total “density” ¯φa + ¯φb = 0.24 and investigate how thecrystalline structures change as the concentration c is varied. InFig. 13(a),when c = 0weobserve a profilewhich is dominatedby species b bumps. The crystal is hexagonally ordered withsome defects (these tend to occur in the vicinity of the φabumps). There are only a few φa bumps, which means thebumps in b are usually sitting next to each other, resulting inthem ordering themselves in a similar manner to that observedin the one-component model [Fig. 7(c)]. Increasing the valueof c from 0 to 0.25, we observe a loss of crystalline structure,as shown in Fig. 13(b). The loss of long range order isclearly visible in the associated radial distribution functions(not shown). The profile in Fig. 13(b) shows a somewhatamorphous structure which appears to include both squareand hexagonal ordering in equal measure. Increasing theconcentration further to c = 0.5, we observe a similar squareordering of bumps as in Figs. 11(j) and 11(k) [in Fig. 11(j),qb = 1, whereas in Fig. 11(k), qb = 1.1; all other parametervalues are the same]. The crystalline structure in Fig. 13(c)at ¯φa = ¯φb = 0.12 contains more vacancies and defects thanthe one in Fig. 11(j) at ¯φa = ¯φb = 0.15, as both average orderparameter values are smaller. Owing to the symmetry inducedby choosing qa = qb (i.e., φa → φb as c → 1 − c), a case withconcentration c is equivalent to the case with concentration1 − c. Thus, Fig. 13(b) also shows the case c = 0.75 if oneinterchanges the orange (light gray) and blue (dark gray)bumps. For this reason values c > 0.5 are not shown.For the asymmetric system the c → 1 − c symmetry doesnot exist and we therefore show five cases for c varying from 0to 1 in Figs. 13(d) c = 0, 13(e) c = 0.25, 13(f) c = 0.5, 13(g)c = 0.75, and 13(h) c = 1. We again observe a transition061408-16MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012)FIG. 13. (Color online) Plots of the order parameter ˆφ, in which bumps in species a appear in orange (light gray) and bumps in species bappear in blue (dark gray), for a constant total order parameter value ¯φa + ¯φb = 0.24. In (a)–(c) we show results from the symmetric mixture(qa = qb = 1), where the concentration of species a is (a) c = 0, (b) c = 0.25, and (c) c = 0.5. In (d)–(h) we show the results from theasymmetric mixture (qa = 1 and qb = 1.1), where (d) c = 0, (e) c = 0.25, (f) c = 0.5, (g) c = 0.75, and (h) c = 1. The parameter values areη = 4 and r = −0.9.from hexagonal ordering in Fig. 13(d) to square orderingin Fig. 13(f) and back to hexagonal ordering in Fig. 13(h)as the value of c is increased from 0 to 1. In betweenthe highly structured states we observe the mixed orderedstates [Figs. 13(e) and 13(g)] that were also present in thesymmetric system. By eye, it is very difficult to pick outthe differences between the symmetric and the asymmetriccases. As previously discussed, the different value of qb in theasymmetric system changes the shape, size, and quantity ofb bumps. In order to characterize and better understand theorganization of the crystalline structures that are formed, werequire a measure which may be used to quantify the structuresand distinguish between hexagonal and square ordering in boththe symmetric and the asymmetric systems. To do this, we useDelaunay triangulation [37,38] to calculate the distribution ofthe bond angles p() between nearest neighbors. We couldhave used other measures from stochastic geometry [39],which were used to characterize the hexagon-square transitionin Be´nard convection [40].The Delaunay triangulation is a triangulation of points (inour case the coordinates of the peaks of the bumps in both orderparameter fields), whichmaximizes theminimumangles of ev-ery triangle (i.e., avoids “skinny” triangles). This triangulationcan be calculated from the Voronoi diagram [37,38] of anyset of points on a 2D plane. The Voronoi diagram is a set ofpolygons, where each polygon represents an area in 2D spacewhich is closer to a particular point than to any of the otherpoints (i.e., the locus of points contained in each polygon iscloser to the bump inside the polygon than any other bump). InFig. 14 we show an example of how we calculate the Delaunaytriangulation for a given order parameter profile. The exampleshows the triangulation for a one-component profile (as thepairing between bumps in the two-component model makesthe triangulation harder to see) but the process is applied inthe same manner to the two-component model. We take thecoordinates of all the bumps to be our points on a 2D plane. Wethen calculate the Voronoi diagram [shown as the light blue(light gray) lines in Fig. 14(b)], which can be used to calculatethe Delaunay triangulation [shown as the red (dark gray) linesin Fig. 14(b)]. This can be done using any of the algorithmsoutlined in Refs. [37,38]. For an efficientmethod of calculatingVoronoi diagrams and Delaunay triangulations, see Ref. [41].061408-17ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)FIG. 14. (Color online) An example Delaunay triangulation is shown for a simple one-component case. In (a) we display a typical orderparameter profile for the one-component model where we observe isolated peaks. The coordinates of the maxima are calculated; these areshown as black points in (b). In (b) we show the Voronoi diagram as a light blue (light gray) polygon network and the Delaunay triangulationas red (dark gray) triangles for this particular set of coordinates.(Note that Delaunay triangulation becomes degenerate whenpoints appear in certain lines of symmetry. However, the initialnoise added to the order parameter fields prevents bumpsfrom forming in perfect symmetry.) We use the statistics ofthe triangles in the Delaunay triangulation to characterize thestructures produced by the bumps.We extract three quantities from the triangulation: the areaof the triangles, the length of the sides, and the angles ineach of the triangles. This information is gathered for fivedifferent realizations of the initial noise profile for systemsof size 200 × 200 and the information is sorted into bins.From these bins we obtain the probability distribution functionfor each quantity. Comparing the different distributions forvarious values of c allows us to observe how the triangles inthe triangulation change as we go from hexagonal to squareordering. Here we concentrate on the probability distributionsof the angles in the triangulation to characterize the crystallinestructures. For results from the other measures, see Ref. [42].In Fig. 15 we display the probability distribution p() forthe triangle corner angles, as the concentration c is varied from0 to 1. We show results for the symmetric (solid red line) andthe asymmetric (dashed blue line) systems. The distribution ofthe angles of the triangles clearly shows the transition betweenhexagonal and square ordering. When c = 0 we observehexagonal ordering in both the symmetric and the asymmetricsystems, which results in the formation of roughly equilateraltriangles in the Delaunay triangulation. This produces angledistributions which have a single peak at 60◦, as shownin Fig. 15(a). As the value of c increases the structurechanges to square ordering, transforming the triangles intoright-angled triangles. Hence, the angle distribution changesand we observe a peak slightly above the value 45◦ andanother peak (half the size) slightly below 90◦, as shownin Fig. 15(c). Increasing the concentration further to c = 1restores the hexagonal ordering; hence, the angle distributionreturns to the single peak at 60◦ [Fig. 15(e)]. In between thepurely hexagonal and the purely square ordered structures weobserve states where the distribution of bond angles is moreevenly spread, with small peaks occurring just above 45◦, ataround 60◦ and just below 90◦ [Figs. 15(b) and 15(d)]. Theserepresent the somewhat amorphous structures which lack thelong range ordering which is present in the hexagonally andsquare ordered structures. The position of these peaks in thebond angle distributions p() does not depend on the quantityor size of the bumps and so the peaks occur in (almost) thesame position for the symmetric and the asymmetric systemsfor all values of c (this is not the case for the area or lengthdistributions). Thismakes the bond angle distributions ideal forcomparing the structure of bump formations in different sys-tems. On comparing the symmetric and the asymmetric caseswe observe that p() appears smoother in the symmetric case.The distribution function p() has a more jagged appearancefor the asymmetric mixture, which we believe is due to the factthat there are different sized bumps in this mixture, making itmore difficult for the bumps to organize themselves into regularstructures. In Figs. 15(a) and 15(e) the distributions appear verysimilar for the symmetric and the asymmetric cases; however,in the other distributions [in particular, Figs. 15(b) and 15(d)]we observe a distinct difference in the height of the three peaks.This suggests that the transition between the different orderedstates occurs differently in the symmetric and the asymmetricsystems.To examine more closely the transition from the hexagonalto the square ordered states we introduce an order parameter which is calculated from the distribution of the angles from theDelaunay triangulation. We integrate the angle distributionsover three regions which cover the three different peaks (theseregions are determined arbitrarily from close examination ofthe angle distributions in Fig. 15) and define the quantitiesR0 =∫ 5325p()d, R1 =∫ 7253p()d,R2 =∫ 11572p()d. (34)We then define the order parameter in the following way: = R0 + R2R1. (35)061408-18MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012) 0 0.05 0.1 30 45 60 75 90p()(a)qb=1.0qb=1.1 0 0.05 0.1 30 45 60 75 90p()(b)qb=1.0qb=1.1 0 0.05 0.1 30 45 60 75 90p()(c)qb=1.0qb=1.1 0 0.05 0.1 30 45 60 75 90p()(d)qb=1.0qb=1.1 0 0.05 0.1 30 45 60 75 90p()(e)qb=1.0qb=1.1FIG. 15. (Color online) The bond angle distribution p() for a constant total order parameter value ¯φa + ¯φb = 0.24. The concentration ofspecies a is (a) c = 0, (b) c = 0.25, (c) c = 0.5, (d) c = 0.75, and (e) c = 1 (corresponding to the simulation snapshots shown in Fig. 13).Results for the symmetric system are shown as red solid lines and the asymmetric system results are shown as blue dashed lines. The parametervalues are qa = 1, η = 4, and r = −0.9.When a structure consists of mainly hexagonal configu-rations of bumps the value of is small (since → 0 asR0 → 0 and R2 → 0) and when a profile is dominated bysquare ordering the value of is large (since → ∞ asR1 → 0). Calculating this quantity for the angle distributionsfor different values of c gives us a measure for the hexagonalvs square ordering of the bumps.In Fig. 16 we show how the order parameter changeswith the concentration c for the symmetric (solid red line) andthe asymmetric (dashed blue line) mixtures. Both curves showa smooth continuous transition from hexagonal ordering tosquare ordering and back again. The different sized bumps inthe asymmetric system break the symmetry around c = 0.5and we observe that the maximum (which corresponds tothe strongest square ordering) occurs at around c ≈ 0.6 and 0 1 2 3 4 0 0.25 0.5 0.75 1cqb=1.0qb=1.1FIG. 16. (Color online) Plot showing the order parameter [defined in Eq. (35)] as a function of the concentration c of speciesa, where ¯φa + ¯φb = 0.24. The red solid line and points show thesymmetric case and the blue dashed line and points show theasymmetric case. The parameter values are qa = 1, η = 4, andr = −0.9.is actually higher than the peak in the symmetric case.The transition to and from square ordering appears to beslightly sharper in the asymmetric case. Even though thereis a difference in the transition between the different orderedstates in the symmetric and the asymmetric mixtures, theyappear to be qualitatively similar. It may be the case that fora larger difference in the values of qa and qb a different typeof transition from hexagonal to square ordering might occur,for example, a discontinuous transition. However, the effect ofvarying the ratio qa/qb is not studied in detail here.IV. CONCLUSIONSIn this paper we have investigated the VPFC model andits application to materials modeling. We first considered theone-component model proposed by Chan et al. in Ref. [13].We determined the linear stability of the homogeneous stateand discussed the dispersion relation. We examined the phasebehavior in 1D and calculated the phase diagram, computingexactly the location of the tricritical point between thehomogeneous and periodic states, and identified the regionof phase space where localized structures occur. Focusingon the latter region of the phase diagram, we investigatedthe localized steady state profiles and discussed the slantedhomoclinic snaking which occurs in the bifurcation diagrams.The one-component model was also studied in 2D and wedetermined the phase diagram, radial distribution functions,and effective pair potentials from our simulation data. Someof the behavior we have identified—the presence of transitionsresembling transitions from a solid phase to a liquid phaseand then to a gaslike phase—replicates behavior observedin nonconserved systems [43]. In Sec. III of the paper, weextended the model to include two coupled order parameterfields. We have considered how the coupling affects thelinear stability of flat profiles and then briefly touched on061408-19ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)the phase behavior of this two-component model. We havefocused on the bump structures which form, considering botha symmetric mixture where the bumps are of equal size andan asymmetric system where one of the bump species isslightly smaller than the other species. The radial distributionfunctions and effective pair potentials for these systems aresomewhat similar to those in binary mixtures of oppositelycharged colloidal particles. We have investigated how varyingthe concentration c of the mixture produces a crossover fromhexagonal to square ordered crystalline structures and how thetransition differs between the symmetric and the asymmetricsystems.A key issue on which we should comment concerns thequestion of what precisely does the order parameter profileφ(x,t) in the VPFC model represent? In the regular PFCmodel, the phase with the uniform flat profile is taken torepresent the liquid phase, while the bump phase correspondsto the crystalline solid. This interpretation is underpinned bythe fact that the regular PFC can be derived from densityfunctional theory (DFT) [2] and dynamical density functionaltheory (DDFT) [10], which is a theory for the dynamics of asystem of interacting Brownian (colloidal) particles [44–47].DFT [48–50] is a statistical mechanics theory for the one-bodynumber density ρ(x) of a system of particles, where ρ(x) =〈ρˆ(x)〉 and where ρˆ(x) = ∑i δ(x − xi) is the density operatorand 〈·〉 denotes a statistical ensemble average [48]. The centralquantity in DFT is the Helmholtz free energy functional F [ρ]and the equilibrium fluid density profile ρ(x) is that whichminimizes the grand free energy [ρ] = F [ρ] − μ ∫ dxρ(x).The DDFT for Brownian particles [44–47] takes as input thisfunctional and so yields the correct equilibrium fluid densityprofile. Making a truncated gradient expansion approximationfor F [ρ], expanding the free energy around that of a referenceliquid state with uniform density ρ0, one can argue that thefree energy is approximately given by Eqs. (2) and (3), wherethe order parameter φ(x) ∝ ρ(x) − ρ0. Thus, it is clear that ina bulk liquid, where ρ(x) is a constant, so too is φ(x) a constantand in the solid phase, where ρ(x) consists of a periodicarray of density peaks, then φ(x) also contains periodicmodulations. However, there are some problems extending thisinterpretation to the VPFC. Consider, for example, Fig. 7(a),where we see a few isolated localized peaks surroundedby a uniform background where φ(x) ≈ 0. Maintaining theabove PFC interpretation, this would correspond to a fewindividual “frozen” particles, surrounded by a fluid of mobileparticles. One might be tempted to think of this as somesort of glass transition [50–52], but the glass transition is acollective phenomenon: In a glass one does see “dynamicalheterogeneity,” that is, regions where the particles are totallyjammed and other regions which are more mobile, but to ourknowledge one never sees a single particle that is jammed onits own surrounded by more mobile particles. Thus, it maybe possible to assume this interpretation may be maintainedfor the VPFC, that is, by considering the localized peakssurrounded by a uniform background to be a dynamicallyheterogeneous glassy system, but there are problems with thispoint of view.An alternative interpretation for the order parameter profilein the VPFC model is that φ(x) is related to a coarse-graineddensity profile (rather than an ensemble average densityprofile) for the system ρ˜(x,t), that is, φ(x,t) ∝ ρ˜(x,t). Follow-ing Ref. [47], we may define the temporally coarse-graineddensity profile for a system of Brownian colloidal particlesas ρ˜(x,t) = ∫ K(t − t ′)ρˆ(x,t)dt ′, where K(t) is a normalizedfunction of finite support which defines a time window overwhich the density is coarse grained. One can then argue [47]that the time evolution equations for ρ˜(x,t) must be verysimilar or even the same as the DDFT equations for the timeevolution of the ensemble average density ρ(x,t), as long as thewidth in time τ of K(t) is large enough. By choosing the timeτ so that it is large compared to the time between the colloidalparticles receiving Brownian “kicks” from the solvent, but isshort compared to the diffusive time scale, corresponding tothe typical time for a particle to diffuse a distance equal toits own diameter, then the coarse-grained density ρ˜(x,t) andthe order parameter φ(x,t) will be quantities which containpeaks, each of which correspond to an individual particle inthe system [53]. Thus, in a low density colloidal suspensionone should see isolated peaks in the coarse-grained density,surrounded by regions where φ(x,t) ≈ 0, corresponding to noparticles being present in that region of the system. This isthe justification for the interpretation made by Chan et al. inRef. [13], that the peaks in the order parameter correspond toparticles and the uniform background corresponds to a portionof solvent free of particles. In order to observe the long timeBrownianmotion of the particles in this description, one shouldadd a stochastic noise term to the dynamical equations for thesystem (23), which continuously drives the system (as opposedto the small amount of noise that is present in our initial orderparameter profiles). However, in numerical simulations therecan be problems with such an approach, because the particlescan become pinned in place by the discrete grid on which theyare defined, and so do not move. We did not make a detailedinvestigation of theVPFCmodel with additional noise. Furtherissues arise as the noise renormalizes the parameters of thecontinuum model.There are state points in the PFC and VPFC phase diagramwhere all possible interpretations of φ break down: Theseare the state points where the equilibrium state is the stripeor the hole phase, such as those displayed in Fig. 10.Systems of spherical particles do not have an ensembleaverage density profile ρ nor a coarse-grained density profileρ˜ with stripes or holes, unless the particles in the systeminteract via pair potentials containing competing attractiveand repulsive parts [18,54,55]. We must conclude that forthe parameter values corresponding to these state points, thegradient expansion that is implicit in the PFC and VPFCfree energy functionals has broken down and that these orderparameter profiles are unphysical.The radial distribution functions for the one-componentmodel displayed in Fig. 7 (see also Fig. 5 of Ref. [13])are very similar to those in real fluids. We observe staticcorrelations which are very similar to what one observesin fluids. Increasing the value of ¯φ increases the numberof bumps and close packing causes long range (crystalline)ordering of the bumps. Calculating the effective pair potentialbetween isolated pairs of bumps, we find a pair potentialhaving an attractive minimum at a pair separation distancewhich is slightly larger than the diameter of the bumps. Thus,the interactions and correlations between bumps share certain061408-20MODELING THE STRUCTURE OF LIQUIDS AND . . . PHYSICAL REVIEW E 85, 061408 (2012)features with some colloidal fluids [56]. We also extend themodel to consider a two-component mixture, with a simplerepulsive coupling between the two order parameter profiles.At low values of ¯φ the bumps commonly appear in pairs and atintermediate values they tend to form chains. At higher valuesof ¯φ the system exhibits crystalline ordering. The appearanceof these structures is somewhat reminiscent of the arrangementof the particles in a binary mixture of oppositely chargedcolloidal particles (see, e.g., Refs. [34,57] and referencestherein). The radial distribution functions and the effectivepair potentials show that there is a fairly strong attractionbetween bumps of the opposite species a and b. The minimumin the ab effective pair potential is at a shorter pair separationdistance than the minimum in the aa and bb pair potentialsand so we observe square ordering when the concentrationc ≈ 1/2 and ¯φ is high enough for the bumps to pack intoa crystalline structure. However, when c ≈ 0 or c ≈ 1, weobserve hexagonal ordering and so we observe a transitionfrom hexagonal to square ordering as the concentration c isvaried. We find that this transition occurs smoothly but canbecome skewed by changing the size of one of the species ofbumps (qa = qb).It would be interesting to further investigate the effect thatvarying the ratio qa/qb has on the structures which form. Inparticular, determining the range of values of qa and qb forwhich bump profiles form in the 2D system would allow oneto determine the range of size ratios of particles (bumps) thatcan be modeled. The transition between hexagonal and squarestructures could then be studied for systems with very differentsized bumps and if theVPFC in this regime continues to be ableto model mixtures of charged colloidal particles, then a widerange of different crystal structures should be observed [34].Note also that the localized structures that we observe arenot a unique property of the VPFC model but are, in fact,also present in the regular PFC model for a small range ofparameter values outside the limit of linear stability. This issomething that we will focus on in future work.ACKNOWLEDGMENTSThis work was supported by the EU via the ITN MUL-TIFLOW (PITN-GA-2008-214919). M.J.R. also gratefullyacknowledges support from EPSRC and A.J.A. thanks RCUKfor support.[1] K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, Phys.Rev. Lett. 88, 245701 (2002).[2] K. R. Elder and M. Grant, Phys. Rev. E 70, 051605 (2004).[3] G. Tegze et al., Soft Matter 7, 1789 (2011).[4] J. Mellenthin, A. Karma, and M. Plapp, Phys. Rev. B 78, 184110(2008).[5] K.-A. Wu, M. Plapp, and P. W. Voorhees, J. Phys.: Condens.Matter 22, 364102 (2010).[6] J. Swift and P. C. Hohenberg, Phys. Rev. A 15, 319 (1977).[7] E. Knobloch, Phys. Rev. A 40, 1549 (1989).[8] P. C. Matthews and S. M. Cox, Nonlinearity 13, 1293 (2000).[9] K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and M. Grant,Phys. Rev. B 75, 064107 (2007).[10] S. van Teeffelen, R. Backofen, A. Voigt, and H. Lo¨wen, Phys.Rev. E 79, 051404 (2009).[11] Z.-F. Huang, K. R. Elder, and N. Provatas, Phys. Rev. E 82,021605 (2010).[12] H. Ohnogi and Y. Shiwa, Physica D 237, 3046 (2008).[13] P. Y. Chan, N. Goldenfeld, and J. Dantzig, Phys. Rev. E 79,035701 (2009).[14] J. Berry and M. Grant, Phys. Rev. Lett. 106, 175702 (2011).[15] G. Stell, J. Stat. Phys. 78, 197 (1995).[16] A. Ciach, W. T. Go´zdz, and R. Evans, J. Chem. Phys. 118, 3702(2003).[17] A. J. Archer, C. N. Likos, and R. Evans, J. Phys.: Condens.Matter 16, L297 (2004).[18] A. J. Archer, D. Pini, R. Evans, and L. Reatto, J. Chem. Phys.126, 014104 (2007).[19] J. Burke and E. Knobloch, Phys. Rev. E 73, 056211 (2006).[20] J. Burke and E. Knobloch, Chaos 17, 037102 (2007).[21] J. Burke and E. Knobloch, Phys. Lett. A 360, 681 (2007).[22] E. Knobloch, Nonlinearity 21, T45 (2008).[23] E. Doedel et al., Caltech, technical report, 2011 (unpublished).[24] J. Burke and E. Knobloch, in Dynamical Systems, DifferentialEquations and Applications, Discrete and Continuous Dynami-cal Systems-Suppl. September, edited by X.-J. Hou et al. (AIMS,Springfield, 2009), pp. 109–117.[25] J. H. P. Dawes, SIAM J. Appl. Dyn. Syst. 7, 186 (2008).[26] D. Lo Jacono, A. Bergeon, and E. Knobloch, J. Fluid Mech. 687,595 (2011).[27] U. Thiele, A. J. Archer, M. J. Robbins, H. Gomez, andE. Knobloch (unpublished).[28] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids(Clarendon Press, Oxford, 1999).[29] G. I. To´th et al., J. Phys.: Condens. Matter 22, 364101 (2010).[30] L. Gra´na´sy, G. Tegze, G. I. To´th, and T. Pusztai, Philos. Mag.91, 123 (2011).[31] S. Muralidharan and M. Haataja, Phys. Rev. Lett. 105, 126101(2010).[32] A. Pototsky, M. Bestehorn, D. Merkt, and U. Thiele, J. Chem.Phys. 122, 224711 (2005).[33] M. J. Robbins, A. J. Archer, and U. Thiele, J. Phys.: Condens.Matter 23, 415102 (2011).[34] M. E. Leunissen et al., Nature (London) 437, 235 (2005).[35] A. P. Hynninen, C. G. Christova, R. van Roij, A. van Blaaderen,and M. Dijkstra, Phys. Rev. Lett. 96, 138308 (2006).[36] A. P. Hynninen, M. E. Leunissen, A. van Blaaderen, andM. Dijkstra, Phys. Rev. Lett. 96, 018303 (2006).[37] M. de Berg, M. van Kreveld, M. Overmars, andO. Schwarzkopf, Computational Geometry: Algorithms andApplications (Springer, Berlin, 1997).[38] Ø. Hjelle and M. Dæhlen, Triangulations and Applications(Springer, Berlin, 2006).[39] D. Stoyan, W. S. Kendall, and J. Mecke, Stochastic Geometryand Its Applications, 2nd ed. (Wiley, West Sussex, 1995).[40] U. Thiele and K. Eckert, Phys. Rev. E 58, 3458 (1998).061408-21ROBBINS, ARCHER, THIELE, AND KNOBLOCH PHYSICAL REVIEW E 85, 061408 (2012)[41] A. Bowyer, Comput. J. 24, 162 (1981).[42] M. J. Robbins, Ph.D. thesis, Loughborough University, 2012.[43] U. Bortolozzo, M. G. Clerc, and S. Residori, New J. Phys. 11,093037 (2009).[44] U. M. B. Marconi and P. Tarazona, J. Chem. Phys. 110, 8032(1999).[45] U. M. B. Marconi and P. Tarazona, J. Phys.: Condens. Matter12, A413 (2000).[46] A. J. Archer and R. Evans, J. Chem. Phys. 121, 4246 (2004).[47] A. J. Archer and M. Rauscher, J. Phys. A 37, 9325 (2004).[48] R. Evans, Adv. Phys. 28, 143 (1979).[49] R. Evans, in Fundamentals of Inhomogeneous Fluids, edited byD. Henderson, Chap. 3 (Dekker, New York, 1992).[50] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids,3rd ed. (Academic Press, London, 2006).[51] L. Berthier and G. Biroli, Rev. Mod. Phys. 83, 587 (2011).[52] A. Cavagna, Phys. Rep. 476, 51 (2009).[53] Taking this point of view, one might then conclude that theparameter q in the VPFC, which determines the size of thebumps, must be related to τ , the width of the coarse-grainingtimewindow.However, since theVPFC is not derived by startingfrom the equations of motion for particles and then followingconsistently this coarse-graining approach, this interpretationcannot be maintained, because varying q affects the wholethermodynamics of the model.[54] A. J. Archer, Phys. Rev. E 78, 031402 (2008).[55] A. J. Archer and N. B. Wilding, Phys. Rev. E 76, 031501 (2007).[56] J.-L. Barrat and J. P. Hansen, Basic Concepts for Simple andComplex Liquids (Cambridge University Press, Cambridge,2003).[57] P. Bartlett and A. I. Campbell, Phys. Rev. Lett. 95, 128302(2005).061408-22

crystals

Abstract

Extracted data

crystals

Abstract

Extracted data

Related items

Related items