Phase field crystal model for heterostructures
Petri Hirvonen, Vili Heinonen, Haikuan Dong, Zheyong Fan, Ken R. Elder, Tapio Ala-Nissila
PPhase field crystal model for heterostructures
Petri Hirvonen, ∗ Vili Heinonen, Haikuan Dong, Zheyong Fan, Ken R. Elder, and Tapio Ala-Nissila
1, 5 QTF Centre of Excellence, Department of Applied Physics,Aalto University School of Science, P.O. Box 11000, FIN-00076, Aalto, Espoo, Finland Department of Mathematics, Massachusetts Institute of Technology,77 Massachusetts Avenue, Cambridge, MA 02139-4307, USA School of Mathematics and Physics, Bohai University, Jinzhou 121000, China Department of Physics, Oakland University, Rochester, MI 48309, USA Interdisciplinary Centre for Mathematical Modelling and Department of Mathematical Sciences,Loughborough University, Loughborough, Leicestershire LE11 3TU, UK (Dated: August 16, 2019)Atomically thin 2-dimensional heterostructures are a promising, novel class of materials withgroundbreaking properties. The possiblity of choosing the many constituent components and theirproportions allows optimizing these materials to specific requirements. The wide adaptability comeswith a cost of large parameter space making it hard to experimentally test all the possibilities.Instead, efficient computational modelling is needed. However, large range of relevant time andlength scales related to physics of polycrystalline materials poses a challenge for computationalstudies. To this end, we present an efficient and flexible phase-field crystal model to describethe atomic configurations of multiple atomic species and phases coexisting in the same physicaldomain. We extensively benchmark the model for two-dimensional binary systems in terms of theirelastic properties and phase boundary configurations and their energetics. As a concrete example,we demonstrate modelling lateral heterostructures of graphene and hexagonal boron nitride. Weconsider both idealized bicrystals and large-scale systems with random phase distributions. Wefind consistent relative elastic moduli and lattice constants, as well as realistic continuous interfacesand faceted crystal shapes. Zigzag-oriented interfaces are observed to display the lowest formationenergy.
I. INTRODUCTION
Most scientifically and technologically important mate-rials are composed of multiple atomic species and phaseswith different chemical compositions, lattice structuresand elastic properties. Some everyday examples in-clude wood, rock, metallic alloys, and concrete. Suchthree-dimensional (3D) materials have been used forthousands of years and efforts towards their develop-ment continue in the age of nanophysics. Some mod-ern examples include e.g. fiber-reinforced polymers andsemiconductor heterostructures. The past decade hasseen the emergence of a completely new type of mate-rials, the atomically thin two-dimensional (2D) materi-als. The extraordinary properties of single component2D materials can be widely enhanced and adjustedby considering their heterostructures that can either bestacked to form vertical multilayer heterostructures ,or they can be grown within a single material sheet intoa lateral heterostructure .The properties of pure or single phase crystalline ma-terials are determined by the complex networks of micro-scopic defects and grains. In contrast, for many multi-phase composite materials, macroscopic continuum mod-els may provide sufficiently accurate predictions of manyof their properties. This suggests that the role of theirmicroscopic structure is less important. Nevertheless, forsemiconductor heterostructures, as well as for verticaland lateral 2D heterostructures, the atomic-level struc-ture of their phase interfaces plays a major role as said structures are miniaturized to the nanoscale where inter-facial effects are important.Predicting the atomic-level structure between two ormore orientationally, structurally and elastically mis-matched phases is particularly difficult. The number ofpossible atomic configurations is essentially endless andconventional atomistic modelling techniques cannot si-multaneously reach all the length and time scales rele-vant to their formation. The more recently developedphase field crystal (PFC) method allows examinationof long, diffusive time scales corresponding to the slowevolution of microstructures, and offers atomic-level spa-tial resolution up to mesoscopic length scales . PFCmodels describe crystalline matter in terms of smooth,classical density fields n i of the different atomic species.The essential thermodynamic quantity is the free energy F [ n i ] that is minimized by a periodic n i . PFC mod-els have been extensively applied to study various com-plex systems and processes such as grain boundaries, va-cancy diffusion, coarsening of polycrystals, heteroepitax-ial growth, yield strength and fracture . In partic-ular, multicomponent PFC models that explicitly incor-porate multiple density fields n i coupled together, havebeen developed and applied to study crystal structurescomposed of multiple atomic species .In this work we introduce spatially smoothed atomicdensity fields coupled to the atomic density fields n i thatenables well-controlled phase separation and, therefore,facilitates modelling heterostructures and composite ma-terials. Smoothed densities have been employed in PFCmodelling recently for introducing a vapor phase and a r X i v : . [ c ond - m a t . m e s - h a ll ] A ug for controlling liquid/solid interface energies . Here weapply this modelling approach to 2D heterostructurescomposed of multiple elements. We carry out a system-atic investigation by varying model parameters one byone to determine their influence on the general behaviorof the model. More specifically, we introduce mismatch inboth the elastic moduli and the lattice constants betweenthe two materials, as well as experiment with differentcouplings between the two density fields. Finally, we as-sess the model’s suitability to study graphene–hexagonalboron nitride 2D heterostructures.The rest of the paper is organized as follows: SectionII lays out and discusses the heterostructure model andgives some practical details of our calculations. SectionIII presents our investigation of the general propertiesof the model using binary heterostructures. Section IVassesses the model’s suitability to modelling graphene–hexagonal boron nitride lateral heterostructures. Finally,Sec. V summarizes the work. II. HETEROSTRUCTURE MODEL
Phase field crystal (PFC) models are a family of contin-uum methods for multiscale modelling of polycrystallinematerials and their complex microstructures. PFC mod-els allow simultaneous access to both atomic and meso-scopic length scales, as well as to long, diffusive timescales. Formation and evolution of microstructures takeplace in this time regime. Conventional PFC models usea smooth, periodic density field n to describe crystallinesystems. The length scale and lattice symmetries, as wellas the elastic properties of the model can be matched withthe target material. These properties are determined bya free energy functional F [ n ( r )] governing the energet-ics of the system. In the solid phase F is minimized by aperiodic n whose symmetries depend on the formulationof F and average density ¯ n . A PFC model can incor-porate multiple density fields coupled together to allowthe study of more complex structures. Such models havebeen applied to study multicomponent materials such as2D hexagonal boron nitride (h-BN) .In conventional multicomponent PFC models the peri-odically oscillating densities representing the solid phaseof each component overlap and form interlocking, mixedlattices. To study heterostructures with well-controlledphase separation, we propose the following dimensionlessfree energy functional: F = (cid:90) d r (cid:32) N (cid:88) i =1 (cid:18) α i n i + β i n i (cid:0) ν i + ∇ (cid:1) n i + γ i n i + δ i n i (cid:19) + N − (cid:88) i =1 N (cid:88) j = i +1 (cid:16) α ij n i n j + β ij n i (cid:0) ν ij + ∇ (cid:1) n j + γ ij (cid:0) n i n j + n i n j (cid:1) + (cid:15) ij η i η j (cid:17) (cid:33) . (1)Here, the first sum contains the ideal contributions of the N density fields and the second, nested sum the contri-butions of the interactions between them. In the firstsum, the quadratic and quartic terms comprise a double-well potential, the cubic term acts similarly to a chemi-cal potential and the gradient term gives rise to periodicsolutions and elastic behavior. We refer the reader toRefs. for a more in-depth discussion of PFC for-mulation. In the second sum, the quadratic and cubicterms are local couplings between the different densityfields, whereas the rest are nonlocal terms. The modelparameters and their roles are summarized in Table I.The term (cid:82) (cid:15) ij η i η j d r in Eq. (1) is the essential cou-pling responsible for controlled phase separation. Thisterm can effectively drive ¯ n i and ¯ n j apart in the samephysical domain such that one corresponds to the disor-dered and the other to a crystalline phase in the phasediagram. The fields η i are spatially smoothed n i wherethe atomic-level structures have been filtered out definedas η i = G ∗ n i . Here the asterisk denotes a convolutionand G is a Gaussian smoothing kernel with the Fouriertransform ˆ G ( k ) = e −| k | / ( σ ) . In the present work wefound that σ = 0 . corresponding to a length scale of ap-proximately five lattice constants with ν i = 1 sufficientlysmooths out the atomic-level structure. To enable atom-istically sharp interfaces and, moreover, to keep the num-ber of parameters to be tuned to a minimum, we did notconsider values σ < . . The influence of this couplingterm is demonstrated in Sec. III A and its ability to drivephase separation is shown analytically in Appendix A.In order to find the n i that minimize F , we used densityconserving gradient descent given by ∂n i ∂t = ∇ δFδn i + ∇ α i n i + β i (cid:0) ν i + ∇ (cid:1) n i + γ i n i + δ i n i + N (cid:88) j =1 j (cid:54) = i (cid:16) α ij n j + β ij (cid:0) ν ij + ∇ (cid:1) n j + γ ij (cid:0) n i n j + n j (cid:1) + (cid:15) ij G ∗ η j (cid:17) (2)where α ij = α ji and similarly for the other parameters.The density conservation constraint is essential for sta-bilizing the heterostructures. Note that while t is calledsometimes time , it is a relaxation parameter that is notrelated to real dynamics in this work. We solved Eq. (2)numerically using a semi-implicit spectral method de-scribed in Ref. 30. This method allows computing thegradients and convolutions present in Eqs. (1) and (2)accurately and efficiently by using fast Fourier transformroutines. Note that according to the convolution theo-rem, convolutions can be expressed as f ∗ g = F − (cid:110) ˆ f ˆ g (cid:111) ,where F − and the carets indicate inverse and forwardFourier transforms, respectively. The following upperbounds for spatial and temporal discretizations were usedfor all the calculations: ∆ x = 0 . , ∆ y = 0 . , and ∆ t = 0 . .Finally, we also used a model system size optimiza-tion algorithm to eliminate strain in our bicrystallinemodel systems of heterostructures. We did not apply themethod to polycrystalline systems, since we did not at-tempt to extract equilibrium densities from them or toanalyze them quantitatively here. III. BINARY HETEROSTRUCTURES
We begin by demonstrating some general properties ofthe present model for simple binary heterostructures with N = 2 and denote the two density fields by n and n .We vary certain model parameters to investigate theirinfluence and will refer to the periodic or “crystalline re-gions” in n i by n (c) i and, similarly, to the disordered re-gions by n (d) i . The “crystalline phase i ” encompasses re-gions where n (c) i and n (d) j coincide. Similarly, the “mixedphase” (disordered phase) spans the regions where n (c) i and n (c) j ( n (d) i and n (d) j ) coincide. TABLE I. Summary of model parameters. They are listed be-low and their significance is explained. While not explicitlywritten in Eqs. (1) or (2), the average densities ¯ n i are impor-tant in controlling the relative stability of different phases. N number of density fields in the model ¯ n i average density; controls the relative stability ofdifferent phases α i temperature-related parameter; controls thediffuseness and facetedness of structures β i controls the elastic moduli γ i acts similarly to a chemical potential δ i usually set to unity in PFC models ν i wavenumber that sets the length scale α ij controls the alignment of two lattices at their mutualinterface but also introduces so-called "weakoscillations" (see Sec. III B) β ij needed for h-BN γ ij needed for h-BN (cid:15) ij couples smoothed densities; controls phase separation ν ij needed for h-BN σ spectral spread of the Gaussian smoothing kernel G A. Influence of smoothed coupling
The crucial parameter here is (cid:15) in the couplingterm for the smoothed density fields. With (cid:15) > , n (c) and n (c) repel each other, whereas with (cid:15) < , n (c) and n (c) attract each other; see Ap-pendix A for an analytical treatise. For the otherparameters, we chose ( α i , β i , ν i , γ i , δ i , α , β , γ ) =( − . , . , . , . , . , . , . , . , for simplicity. Thischoice of model parameters is used throughout Sec. IIIunless stated otherwise. Note that this choice of modelparameters is symmetric, i.e. , F ( n , n ) = F ( n , n ) .We used simple model systems where we initializedboth n and n roughly 50% crystalline and 50% disor-dered with average densities ¯ n (c) i ≈ . and ¯ n (d) i ≈ . .The initial structure for n (c) i was obtained using an in-verted hexagonal one-mode approximation . We alsoarranged n (c) and n (c) in partial overlap to force somechanges in them during relaxation. Figure 1 illustratesthe relaxed heterostructures obtained with (cid:15) = ± . .Panels (a) - (c) visualize the systems and panels (d) - (f)plot corresponding density profiles. Panel (a) shows theinitial state with n (c) and n (c) in partial overlap. Panel(b) gives the repulsive case with (cid:15) = 0 . where n (c) and n (c) have pushed themselves apart to eliminate themixed phase. Panel (c) depicts the attractive case with (cid:15) = − . where n (c) and n (c) have come to a full overlapforming a coexistence between a mixed and a disorderedphase. Recall that n and n are coupled here only via η and η whereby the two atomic lattices do not interact.Consequently, the lattices can end up arbitrarily alignedsuch as here; see panel (c).Next we considered polycrystalline heterostructures.The density fields n and n were initialized with white Distance (lattice constants) D e n s i t y ( d i m e n s i o n l e ss )
40 45 50 55 60
Distance (lattice constants) -0.6-0.4-0.200.20.40.60.8 D e n s i t y ( d i m e n s i o n l e ss )
15 20 25 30 35
Distance (lattice constants) -0.6-0.4-0.200.20.40.60.8 D e n s i t y ( d i m e n s i o n l e ss ) (a)(b)(c) (d) (f)(e) phase 1 phase 2mixed phase disord. phase FIG. 1. Influence of the coupling parameter (cid:15) on the heterostructures. (a) The initial state with n (c) and n (c) in partialoverlap. (b) The relaxed heterostructure for the repulsive case where (cid:15) = 0 . . Here, the crystalline phases 1 and 2 are shownin cyan and red, respectively. (c) The relaxed coexistence between a mixed and a disordered phase for the attractive case where (cid:15) = − . . Here, the mixed phase appears white due to the coindicental alignment of the structures in n (c) and n (c) , and thedisordered phase appears black. (d) Profiles of the smoothed densities η (cyan) and η (red) along the periodic edge in thehorizontal direction of the initial state from panel (a). (e) Profiles of the densities (solid lines) n (cyan) and n (red) and ofthe smoothed densities (dashed lines) η (cyan) and η (red) along the periodic edge in the horizontal direction of the relaxedheterostructure from panel (b). (f) Same profiles for panel (c). noise with ¯ n i = 0 . . With (cid:15) = 0 . , a mixed phaseemerges first, followed by delayed decomposition intothe two separate crystalline phases. We, therefore, used (cid:15) = ± . to drive n (c) and n (c) apart ( + ) or together( − ). Note, however, that if the coupling strength is in-creased significantly more, stripe phases may replacethe crystalline ones as the most stable phase. Figure 2demonstrates both the repulsive and attractive cases af-ter 7 500 time units of relaxation. In the repulsive caseshown in panel (a), n (c) and n (c) are well separated andthere is no mixed or disordered phase. The interfaces be-tween the two phases appear fuzzy and disordered as ex-pected due to no interaction between the two underlyinglattices. The attractive case is shown in panel (b) where n (c) and n (c) are in full overlap, yielding a patched co-existence between a mixed and a disordered phase. Thearbitrary misorientations and translations between thelattices in n (c) and n (c) result in a multitude of Moirépatterns. B. Influence of α ij Next, we varied the quadratic coupling parameter α to study its influence on the interfaces between the twocrystalline phases and on the heterostructures as a whole.Here, we fix (cid:15) = 1 . . We chose α < to achievecommensurate alignment of the two crystalline latticesat their interface to ensure the continuity of the underly-ing honeycomb lattice there. A side effect of this couplingis that it causes n (c) i to induce oscillations in n (d) j . Theamplitude of such weak oscillations in n (d) j should be con-strained to keep the two crystalline phases from mixingtogether. Indeed, the weak oscillations can be viewed asslight intermixing of the different atomic species. Inter-mixing is common in metallic alloys and in doped semi-conductors and has been observed in lateral heterostruc-tures of graphene and hexagonal boron nitride as well .While constrained by the amplitude of the weak oscilla-tions induced, the magnitude of α should be maximizedto ensure continuity even for lattice-mismatched or mis-oriented interfaces. We optimized α using bicrystallineheterostructures. We observed that the heterostructuresare rendered unstable when α = − . , but with α = − . they retain their stability while the amplitude of (a) (b) FIG. 2. Influence of (cid:15) on polycrystalline heterostructures. (a) A blow-up of a larger system for the repulsive case where (cid:15) = 1 . after 7 500 time units. The left-hand side of the panel reveals the distribution of the two phases and the right-handside represents the heterostructure by m = n + n for a clearer illustration of the atomic-level structure. (b) A blow-up of alarger system for the attractive case where (cid:15) = − . after 7 500 time units. Moiré overlap patterns between the two latticesare clearly visible.
15 20 25 30 35
Distance (lattice constants) -0.8-0.6-0.4-0.200.20.40.60.8 D e n s i t y ( d i m e n s i o n l e ss ) (a)(b) FIG. 3. Influence of α and the corresponding coupling on abicrystalline heterostructure. (a) A blow-up of a relaxed het-erostructure with α = − . . (b) Profiles of the densities(solid lines) n (cyan) and n (red) and of the smoothed den-sities (dashed lines) η (cyan) and η (red) along the periodichorizontal edge of the relaxed heterostructure from panel (a). the weak oscillations remains negligible. Figure 3 illus-trates the interface in a relaxed bicrystalline heterostruc-ture with α = − . , ¯ n (c) i = 0 . and ¯ n (d) i = 0 . . Itis clear both from the visualization of the heterostruc-ture as well as from the density profiles below it that thehoneycomb lattice is highly continuous from one phase tothe other. Furthermore, the interface is atomically sharpwith an approximate width of two lattice constants. We further demonstrated the influence of α for poly-crystalline heterostructures. We initialized n and n with white noise where ¯ n i = 0 . . Figure 4 demonstratesa relaxed heterostructure. On a larger scale, the sys-tem resembles that shown in Fig. 2 (a), but here theinterfaces between the two crystalline phases are betterordered and more continuous. The width of the inter-faces appears small for all misorientations. Note also theweak oscillations in n (d) visible in panel (a). C. Influence of β i The crystalline phase i is present where n (c) i and n (d) j coincide. The elastic properties of said phase should bedictated by n i , but n j can have a minor contribution aswell. We demonstrate here to what extent the elasticproperties of the two crystalline phases can be controlledvia β and β , the coefficients of the gradient terms in F responsible for elastic contribution. The ability to controlthe elastic stiffness of both phases separately is essentialwhen modelling realistic heterostructures. Note that weshow in Appendix B that the smoothed densities η i havea negligible elastic contribution.The contribution from a uniform elastic deformationto the free energy density is given by f e = C (cid:0) ε x + ε y (cid:1) + C ε x ε y , (3)where ε x and ε y are the x and the y components of strain,and C = C and C = C are the stiffness coef-ficients. Furthermore, the bulk, shear and 2D Young’smoduli, as well as Poisson’s ratio are given by (a) (b) FIG. 4. Influence of α on a polycrystalline heterostructure. (a) A blow-up of a larger system after a relaxation of 25 000time units. The left-hand side of the panel reveals the distribution of the two phases and the right-hand side demonstrates n with weak oscillations in n (d) . (b) The heterostructure from panel (a) represented by m = n + n for a clearer illustration ofthe atomic-level structure. B = C + C , (4) µ = C − C , (5) Y = 4 BµB + µ , (6)and ν = B − µB + µ , (7)respectively .We determined the elastic coefficients of the two crys-talline phases separately by straining single-crystals ofeither phase in the small deformation limit. More pre-cisely, we varied − . ≤ ε x ≤ . and − . ≤ ε y ≤ . independently. We fixed β = 1 . , and varied . ≤ β ≤ . For . ≤ β ≤ . , we used (cid:15) = 1 , but,for β ≤ . ( β ≥ ), we had to adjust . ≤ (cid:15) ≤ . ( . ≤ (cid:15) ≤ ) to retain the stability of the heterostruc-tures. We again used α = − . to include the weak os-cillations. We determined the average densities ¯ n (c) i and ¯ n (d) i in equilibrium by relaxing bicrystalline heterostruc-tures and by extracting the average densities from themiddle of the crystalline phases.Figure 5 shows the 2D Young’s modulus as a functionof β for both crystalline phases. For crystalline phase 1,the modulus is essentially unaffected by β , i.e., the cor-responding linear fit has a negligible slope. In contrast, the modulus for the crystalline phase 2 is linearly pro-portional to β with a slope of 0.17. Independent controlof the elastic stiffness of either of the crystalline phasesappears straightforward. In addition, for each value of β , we observed C ≈ C / , whereby ν ≈ / . This isa feature common to many simple PFC models .Finally, we compared the numerical results against ananalytical prediction where for one crystalline phase C = 9 (cid:88) i β i φ i (8)with an amplitude φ i = 115 δ i (cid:16) γ i − δ i ¯ n i − (cid:113) γ i − α i δ i + 12 δ i ¯ n i (2 γ i − δ i ¯ n i ) (cid:19) (9)and C = C / . (10)Here, C is given simply as a sum over the individualdensity fields. The values obtained for Y using the ana-lytical expressions above are also plotted in Fig. 5 for thecrystalline phase 2. The numerical and analytical resultsare in very good agreement. Note that the amplitude ofthe weak oscillations in the crystalline phase 1 is roughlyan order of magnitude lower than that of the oscillationsin the crystalline phase 2. Since C ∝ φ i , the influenceof these oscillations is negligible. FIG. 5. Two-dimensional Young’s modulus Y as a functionof the gradient term coefficient β for both crystalline phases.The markers give the actual data and the lines are optimallinear fits. The slope for the second fit is 0.17. The analyticalprediction obtained for the crystalline phase 2 using the an-alytical expressions for the stiffness coefficients C and C [Eqs. (8) - (10)] is plotted using solid black markers.TABLE II. Average densities for lattice-mismatched bicrys-talline heterostructures. The mismatch is indicated by λ . λ ¯ n (c) ¯ n (d) ¯ n (c) ¯ n (d) D. Influence of ν i As a final demonstration of our model for binary het-erostructures, we introduced lattice mismatch betweenthe two crystalline phases via ν i = 1 /λ i . We fixed λ = 1 . and varied λ through values 1.05, 1.1 and 1.2.We ensured the stability of the heterostructures by choos-ing (cid:15) = α = − β i = 1 . for simplicity.We first considered bicrystalline model systems eitherwith different permutations of armchair and zigzag edgesalong the interface or with symmetrically tilted crystalswith a tilt angle θ = θ − ( − θ ) , where . ◦ ≤ θ ≤ . ◦ .We considered two different strains. For unstrainedsystems, the periodicities of both bicrystal halves werematched separately with the periodic domain wherebythe lattice mismatch is accommodated by misfit dis-locations. For strained systems, both bicrystal halveswere initialized with an average lattice constant and wereagain matched with the domain whereby the lattice mis-match is accommodated via elastic deformation. Theaverage densities for the different strain and mismatchcases are given in Table II. Note that said densities werechosen to yield an approximate 1:1 coexistence betweenthe two crystalline phases. For reference, we consideredhere also λ = 1 . with (cid:15) = 1 . .Overall, the phase interfaces obtained displayed well-defined structures. While misorientation and lattice mis- match introduce defects, extensively fuzzy and ill-definedstructures are rare. In addition, the vast majority ofthe highly strained systems remained stable during re-laxation. Figure 6 offers a representative sample of thestructures obtained and shows a comparison between dif-ferent strain and mismatch cases.The first row of panels in Fig. 6 demonstrates zero-misorientation armchair-armchair interfaces between lat-tices of varying mismatch. We observed perfect hexago-nal order for the reference and strained cases. Indeed, thestrained heterostructures retained their stability withoutexperiencing any stress-relieving reconstructions such asvia subsequent nucleation, creation and annihilation ofdislocations. In the unstrained structures, we observedperiodic arrays of point-like misfit dislocations along theinterface. We obtained similar structures for zigzag-zigzag interfaces.The second row of panels in Fig. 6 gives low-misorientation interfaces between symmetrically tiltedlattices of varying mismatch. Here, the tilt angle θ ≈ . ◦ . All structures display fairly periodic arrays ofdislocations. For corresponding graphene grain bound-aries, alternatingly slanted dislocations are expected .For the reference structure, the highly symmetric initialstate used has lead to extended defect structures withoutsuch symmetry breaking. The low-strain structure with λ = 1 . indeed displays alternatingly slanted disloca-tions, whereas the high-strain structure with λ = 1 . again does not, due to having achieved some strain-reliefvia annihilations of dislocations. The unstrained struc-tures appear very similar to the corresponding strainedones.The third row of panels in Fig. 6 depicts high-misorientation armchair-zigzag interfaces between lat-tices of varying mismatch. Despite the extreme misori-entation between the two crystalline phases, the presentmodel performs well in stitching them together withrather well-defined atomic-level structures. In fact,the lattice-mismatched structures do not appear visiblyfuzzier than the reference.Finally, we simulated the growth of polycrystallineheterostructures with the aforementioned lattice mis-matches. The density fields were initialized with whitenoise where ¯ n i = [¯ n (c) i + ¯ n (d) i ] / .Figure 7 gives examples of the lattice-mismatchedpolycrystalline heterostructures obtained. Despite thevarious misorientations between the two crystallinephases at their interfaces, the model performs well in lo-calizing the mismatch into point-like dislocations. Someinterfaces, especially those with larger mismatch, appearsomewhat diffuse, but are well comparable to some single-phase PFC grain boundaries; cf. Ref. 38, for example.For the cases with λ = 1 . and λ = 1 . , the mismatchis minor and the interfaces appear highly continuous; seepanels (a) and (b), respectively. There are a numberof dislocations along the interfaces in both heterostruc-tures, but many are due to lattice misorientation. Forthe case with λ = 1 . , the interfaces are still fairly con- A C - ZZ λ = 1.0 (R) λ = 1.05 (S) λ = 1.2 (S) λ = 1.05 (U) λ = 1.2 (U) A C - A C θ = . ° FIG. 6. Collage of interface structures for different lattice mismatches, strains and misorientations. The left hand side ofeach panel reveals the distribution of the two crystalline phases, and the right hand side represents the heterostructure by m = n + n for a clearer illustration of the atomic level structure. Note that in many cases only a small part of the totallength of the interface modeled is shown. The first column of panels gives reference (R) structures with no lattice mismatchbetween the two phases. The next two columns give strained (S) structures where the mismatch is accommodated by elasticdeformation. The last two columns give unstrained (U) structures where the mismatch is accommodated by misfit dislocations.The mismatch for each column is indicated via λ = 1 /ν . The first row of panels gives structures with zero-misorientationarmchair-armchair (AC-AC) interfaces. The second row depicts low-misorientation tilt interfaces with θ ≈ . ◦ . The thirdrow demonstrates high-misorientation armchair-zigzag (AC-ZZ) interfaces. tinuous, but display several regions with somewhat fuzzyfeatures; see panel (c). These regions seem to coincidewith greater interfacial curvature. Last, it also appearsthat the mobility of the interfaces decreases with increas-ing mismatch. This is evident from the noticeably smallerdomain sizes in the heterostructure with λ = 1 . ; notethat all three have been relaxed for the same 250 000time units. IV. APPLICATION TO LATERALGRAPHENE–HEXAGONAL BORON NITRIDEHETEROSTRUCTURES
In this section, we consider a three-component modelfor lateral 2D heterostructures between graphene and h-BN (G–h-BN). We focus here on demonstrating the suit-ability of the model to qualitative modelling of G–h-BN.Finding optimal model parameters for quantitatively ac- (a) (b) (c)
FIG. 7. Lattice-mismatched polycrystalline heterostructures. The panels offer blow-ups of larger systems after a relaxation of250 000 time units. The top halves of the panels show the distribution of the two phases and the bottom halves represent theheterostructures by m = n + n for a clearer illustration of the atomic-level structure. We have fixed λ = 1 /ν = 1 and havevaried λ = 1 . in (a), λ = 1 . in (b) and λ = 1 . in (c). curate modelling of G–h-BN (involving, e.g., fitting tointerfacial formation energies) will be presented in futurework. A. Model requirements and parameters
To model G–h-BN, we set N = 3 and chose to modelthe graphene phase with n (c) , n (d) and n (d) and the h-BNphase with n (d) , n (c) and n (c) . The parameters for n , n and their mutual couplings were adopted from Ref. 26.The other parameters were chosen by trial and error byvarying them one at a time. We use the following criteria,guiding principles and simplifying assumptions: • Start by choosing parameters for n that are similarto those for n and n . • Assume that boron and nitrogen atoms are inter-changeable, i.e., F ( n , n , n ) = F ( n , n , n ) . • Match the relative Young’s moduli and latticeconstants to the experimentally and theoreticallydetermined values Y h-BN /Y G ≈ . and a h-BN /a G ≈ . , respectively. • Require sharp, continuous interfaces and facetedcrystal shapes . • Retain the stability of the heterostructures andlimit the amplitude of weak oscillations.Table III gives a set of model parameters that was foundto satisfy the criteria listed above. Most importantly,this choice of parameters yielded Y h-BN /Y G = 0 . and a h-BN /a G = 1 . in fair agreement with the target val-ues. In the following, we demonstrate how the modelbehaves and how it fulfills the other criteria above. TABLE III. Set of parameters for lateral heterostructures ofgraphene and h-BN. Note that α i = α ij when i = j andsimilarly for the other parameters. N α ij i = 1 i = 2 i = 3 j = 1 -1.4 – – j = 2 -0.04 -0.3 – j = 3 -0.04 0.5 -0.3 β ij i = 1 i = 2 i = 3 j = 1 j = 2 j = 3 γ ij i = 1 i = 2 i = 3 j = 1 j = 2 j = 3 δ ij i = 1 i = 2 i = 3 j = 1 j = 2 – 1.0 – j = 3 – – 1.0 (cid:15) ij i = 1 i = 2 i = 3 j = 1 – – – j = 2 -0.8 – – j = 3 -0.8 0.0 – λ ij = 1 /ν ij i = 1 i = 2 i = 3 j = 1 j = 2 j = 3 ¯ n ( j ) i i = 1 i = 2 i = 3 j = c 0.31 -0.32 -0.32 j = d 0.66 -0.65 -0.65 B. Atomic configurations
Figure 8 demonstrates a zigzag-oriented interface in abicrystalline G–h-BN system. Panel (a) gives both the0
15 20 25 30 35
Distance (lattice constants) -1.5-1-0.500.511.5 D e n s i t y ( d i m e n s i o n l e ss ) (a)(b) FIG. 8. Zigzag interface from a bicrystalline G–h-BN lateralheterostructure. (a) A visualization of the heterostructurein the top half and in the bottom half the same structurerepresentedy by m = n + n + n for a clearer illustrationof the atomic level structure. In the top half, graphene ap-pears cyan, whereas boron and nitrogen are in magenta andyellow. In the bottom half, graphene (h-BN) appears darker(brighter). (b) Profiles of the densities (solid lines) n (cyan), n (magenta) and n (yellow) and of the smoothed densities(dashed lines) η (cyan), η (magenta) and η (yellow) alongthe periodic edge in the horizontal direction of the relaxedheterostructure. distribution of the two phases and a representation of thesame structure by m = n + n + n for a clearer illustra-tion of the atomic-level structure. Note that in the latterthe h-BN phase appears brighter facilitating identifica-tion of the two phases in such figures. The interface dis-plays perfect hexagonal order and is again atomisticallysharp. Panel (b) gives the profiles of n i and η i along thehorizontal periodic edge of the system and shows that theamplitude of the weak oscillations is small. Note that incontrast to the binary heterostructures considered in Sec.III, here the average densities ¯ n and ¯ n are negative.Figure 9 demonstrates a large polycrystalline G–h-BNsystem grown from white noise where ¯ n i = [¯ n (c) i + ¯ n (d) i ] / .In panel (a), an overview of an approx. 50 ×
50 nm system is given by a coarse-grained representation wherethe graphene (h-BN) phase appears cyan (red). After arelaxation of . × time units, the heterostructure as-sumes configurations typical to spinodal decompositionin binary systems. Coarsening is slow because we havestrived here for stable sharp crystalline structures insteadof diffuse high-temperature ones. Panels (b) - (e) show ablow-up of the region indicated by a blue square in panel(a). Panel (b) visualizes the region and panel (c) presents m for a clearer illustration of the atomic-level structure.Despite the various misorientations present in the system,the structure of the interfaces is overall well-defined, ex- cluding few fuzzy patches. Panels (d) and (e) show n and n , respectively. Both appear faceted and have sharpinterfaces with a primary (secondary) preference for thezigzag (armchair) direction. Weak oscillations are alsovisible in both panels.The coarsening of G–h-BN was found slow with thepresent model and set of parameters. Moreover, concur-rent nucleation and growth is not how said heterostruc-tures are produced in practice . We demonstratedpreparing more realistic model systems of random poly-crystalline G–h-BN with larger phase domain and grainsizes. For initialization, we used Voronoi grain struc-tures with random seed points, crystal orientations andphases and relaxed for 25 000 time units for local re-laxation of the interfaces and grain boundaries. Figure10 gives an overview of one such system where the ini-tial, large-scale Voronoi structure has remained essen-tially unchanged as shown by the coarse-grained depic-tion of the system in panel (a). Panels (b) - (e) showthe total density m illustrating the atomic level struc-ture of selected interfaces and boundaries. Panel (b)displays two triple junctions, one within h-BN and theother between graphene and h-BN, connected by an in-version boundary within h-BN. In h-BN, an inversionboundary is formed between two crystals with a mis-orientation of ◦ as the ordering of boron and nitro-gen becomes inverted in one crystal with respect to theother . The interfaces between graphene and h-BNhave small-to-intermediate misorientations, whereas thegrain boundaries within h-BN are large-angle boundaries.The interfaces and grain boundaries appear disorderedbut have fairly well-defined atomic level structures. Theinversion boundary is formed by a perfect 4|8 chain asexpected . Panels (c) - (e) show longer interfaces be-tween graphene and h-BN: (c) a large-misorientation in-terface, (d) a small-misorientation zigzag interface and(e) a small-misorientation armchair interface. Whilethe large-angle interface shown in panel (c) appears dis-ordered, all interfaces display well-defined atomic levelstructures. C. Interface energies
Finally, we investigated the relative stability of G–h-BN interfaces in different lattice directions by studyingtheir formation energies. We considered 12 different in-terface angles ◦ ≤ θ ≤ ◦ , where θ = 0 ◦ correspondsto armchair and θ = 30 ◦ to zigzag interfaces. For sim-plicity, we considered here only strained configurationswith perfect honeycomb order and no misfit dislocationsalong the interfaces (see Fig. 11 for an example) to avoidvery large system sizes. We assumed all interfaces tobe composed of zigzag and armchair segments as shownin Fig. 11. The model appears to yield stepped inter-faces, typically with minimal segment lengths L ZZ and L AC . Interfaces initialized with longer segments are alsoat least metastable.1 (a) (d) (e)(b) (c) FIG. 9. Large random polycrystalline graphene–h-BN lateral heterostructure relaxed from white noise for . × time units.The sides of the system are approx. 50 nm in length. (a) A coarse-grained representation of the large-scale structure wheregraphene appears cyan and h-BN red. (b) - (e) Blow-ups of the region indicated by the blue square in panel (a). The width(height) of the region shown in the blow-ups is approximately 9 nm. (b) A visualization explained in Fig. 8, (c) the total density m = n + n + n , (d) n and (e) n in the blow-up. (e)(d) (c)(b)(a) (b) (c)(d)(e) FIG. 10. Large random polycrystalline graphene–h-BN lateral heterostructure from a random Voronoi grain structure. A sideof the system is approximately 50 nm long. (a) A coarse-grained representation of the large-scale structure where grapheneappears cyan and h-BN red. (b) - (e) Blow-ups of the atomic level structure of the regions indicated by the blue squares inpanel (a). The regions shown in the blow-ups are 6 nm wide. (b) A collection of graphene–h-BN interfaces and h-BN grainand inversion boundaries. (c) A large-misorientation interface. Small-misorientation (d) zigzag and (e) armchair interfaces. (a)(b) L ZZ L AC FIG. 11. Examples of strained bicrystalline heterostructuresused to study the formation energy of G–h-BN interfaces withperfect honeycomb order. The interface angle is θ ≈ . ◦ forboth structures and is indicated by the white wedge in panel(a). In panel (a), m = 2 and in panel (b) m = 4 . For both, n = 1 and the total width is L ⊥ ≈ nm (see text for thedefinition of m and n ). The zigzag and armchair segmentsof the right-hand side interfaces are traced in red and blue,respectively. Their lengths L ZZ and L AC are also indicated.In panel (a), L ZZ is 4 lattice constants and L AC is √ latticeconstants. In panel (b), both are twice as long. Note thatthe system shown in panel (a) can be decomposed into twoidentical fields by cutting the domain in half perpendicularto the interface. Here we consider such subdomains with fourvertices. Assuming a sufficiently large bicrystalline heterostruc-ture with dimensions L ⊥ perpendicular and L (cid:107) parallelto the two interfaces, the total formation energy of thesystem can be written as F = f L ⊥ L (cid:107) = f ∗ L ⊥ L (cid:107) + 2 γL (cid:107) + 4 δ, (11)where f is the free energy density per unit area obtainedby evaluating Eq. (1) and by dividing by the total area,and f ∗ = xf G + (1 − x ) f h-BN is the effective free energydensity per unit area given by the equilibrium bulk freeenergy densities of the two phases weighted by their areafractions x and − x . Note that we fixed x ≈ . byfixing ¯ n i = ± . and by initializing the two phaseswith equal or very close to equal areas. Here γ is theaverage formation energy of the interface per unit length(two interfaces appear here; one has C-B and the otherC-N bonds along its zigzag segment) and δ is the averageformation energy of a vertex (each system has a total offour vertices; two are convex and the other two concavewith respect to one of the two constituents).The energy terms f ∗ , γ and δ can be obtained by fit-ting Eq. (11) to simulation data using the method of leastsquares. For each interface angle considered, we varied L ⊥ = nL ⊥ and L (cid:107) = mL (cid:107) , where L ⊥ is roughly 10 nm, L (cid:107) is the minimal L (cid:107) that satisfies periodic boundaryconditions, n = 1 , , . . . , and m = 1 , , . . . , . For eachcombination of L ⊥ , L (cid:107) , we initialized a correspondingbicrystalline heterostructure with segmented interfaces, relaxed it, extracted the final f, L ⊥ and L (cid:107) , and fittedEq. (11) to these data. Outliers in the data and visuallydivergent configurations were excluded from the analysis.We determined γ and δ for long segment lengths withlarge m . In practice, we fitted to data points unaffectedby nonlinear finite-size effects where typically n, m > .In addition, we considered minimal segment lengths with m = 1 (and L (cid:107) hence a constant), i.e. , interfaces with themaximal packing density of vertices. In this case the ver-tex energy needs to be absorbed into the interface energyas they cannot be separated without varying L (cid:107) . Thisgives the scaling relation f L ⊥ = f ∗ L ⊥ + 2 γ ∗ , (12)where γ ∗ = γ + 2 δ/L (cid:107) . This scaling relation is one-dimensional. Note that for long segment lengths the en-ergy contribution of the vertices is negligible γ (cid:29) δ/L (cid:107) ,and hence γ ∗ ≈ γ .In the limit of long segment lengths, one can derive ananalytical expression for γ as a sum of the two segments’individual contributions. Since the angle between thezigzag and armchair segments is 150 ◦ , and as we know θ and L (cid:107) , simple geometrical considerations give γ = 2 γ ZZ sin ( θ ) + 2 γ AC sin (30 ◦ − θ ) , (13)where γ ZZ = γ ( θ = 30 ◦ ) and γ AC = γ ( θ = 0 ◦ ) are theinterface energies of pure zigzag and armchair interfaces,respectively.Figure 12 shows the dimensionless interface γ and ver-tex energies δ as a function of the interface angle θ . Allerror bars are given by two-sigma confidence intervals.In the limit of long segment lengths, the numerical datafor γ agrees well with the analytical expression givenby Eq. (13) predicting a maximal interface energy at θ ≈ ◦ . One should note that both the zigzag and arm-chair interfaces have a locally minimal interface energywith respect to the interface angle. In this case the zigzaginterface has the lower energy. This is consistent with anumber of experimental findings , reporting apreference towards zigzag interfaces. The dominance ofzigzag interfaces can be explained by the growth process:zigzag-faceted crystals of one phase are typically formedfirst and serve as seeds for the subsequent growth of thesecond phase. The thermodynamic stability of the zigzaginterfaces has also been verified computationally usingdensity functional (DFT) theory . A contradictorypreference for armchair interfaces has also been reportedby some DFT studies .Figure 12 also gives the interface energy for interfaceswith minimal segment lengths ( m = 1 ), i.e., interfaceswith the maximal packing density of vertices. Again, thezigzag direction yields the lowest energy. However, here γ is approximately constant for θ ≤ ◦ and decreaseswith θ for θ ≥ ◦ . There are two data points with sig-nificantly larger error bars. For these cases, the scaling is3 FIG. 12. Dimensionless interface γ (in blue on left axis) andvertex δ (in red on right axis) energies as a function of theinterface angle θ ; see Eq. (11). Open markers correspondto γ in the limit of long segment lengths (large m ) and solidmarkers to γ in the limit of minimal segment lengths ( m = 1 );see Eq. (12). The solid curve gives the analytical expressionfor γ from Eq. (13). not as linear as for the other interface angles despite visu-ally ideal configurations. It appears that interfaces withminimal segment lengths have lower formation energiesin general.Figure 12 also shows that the vertex energy δ is zerofor both zigzag and armchair directions where there areno vertices. For the intermediate directions, δ is foundslightly negative and roughly constant. This explainswhy the formation energy is generally lower for interfaceswith more vertices. A negative δ is possible, since the ver-tices cannot exist independently of the segments whose γ > . Furthermore, interfaces with long segment lengthsproved stable as the highly symmetric initial states pro-vided insufficient driving force to overcome the energybarriers for nucleating more vertices. Related point de-fects, triple junctions between grain boundaries, have alsobeen shown to display negative formation energies . V. SUMMARY AND OUTLOOK
We have introduced an efficient and flexible phase fieldcrystal model intended for studying heterostructures orcomposite materials for the general case of N atomicspecies. This model allows well-controlled phase sepa-ration via the use of smoothed density fields. The lat-tice symmetries and the length scale, as well as the elas-tic properties, of the individual phases can be controlledreadily. This model offers a straightforward approach tomodelling systems with multiple ordered phases.We have carried out a comprehensive demonstrationof the model’s properties using simple 2D binary sys- tems. More specifically, we have varied several modelparameters independently to investigate their influenceon phase separation, on the atomic-level structure andorder of the phase interfaces and on the elastic prop-erties of the two phases. A lattice constant mismatchbetween the two phases results in disordered but gener-ally well-connected interfacial configurations. The elas-tic properties of the two phases can be controlled in-dependently and robustly. We have also demonstratedthe applicability of this model by considering graphene–hexagonal boron nitride lateral heterostructures (G–h-BN). We have shown that the model can reproduce manyof the features of G–h-BN, such as the relative latticeconstants and Young’s moduli of the two phases, as wellas continuous interfaces with a preference for zigzag andarmchair directions. We have also demonstrated how tomodel large, complex G–h-BN microstructures.One obvious extension to this study would be to fur-ther optimize the model parameters used for G–h-BN.While we have matched the relative lattice constants andYoung’s moduli approximately to their experimental val-ues, most model parameters were either adopted fromprevious works or were chosen on qualitative grounds.Especially the coupling coefficients in the model couldbe fitted by matching the structure and formation en-ergy of phase interfaces to corresponding results fromatomistic calculations. The different chemical affinitiesbetween carbon and boron, and carbon and nitrogen, inparticular, could be incorporated to the model via theseparameters. In addition, structural or other more sophis-ticated PFC models could be incorporated to the modelby replacing the terms in the energy F proportional to β i and β ij with convolution kernels for more accurate mate-rial description or to allow a broader range of lattice sym-metries. Although here we have focused on the 2D casefor conceptual simplicity, the model is also applicable to3D problems, where various nanoscale heterostructuresand mesoscopic multiphase microstructures or compositematerials could be considered. Constraining the couplingbetween the different lattices to their mutual interfacescould facilitate eliminating the occasional fuzzy struc-tures without amplifying the weak oscillations. Achievingthis without making the equations of motion significantlymore complicated is a topic for future work. VI. ACKNOWLEDGMENTS
This work has been supported in part by the Academyof Finland through its QFT Center of Excellence Pro-gram grant (no. 312298). We acknowledge the computa-tional resources provided by the Aalto Science-IT projectand the CSC IT Center for Science, Finland. P.H. ac-knowledges financial support from the Vilho, Yrjö andKalle Väisälä Foundation of the Finnish Academy of Sci-ence and Letters. K.R.E. acknowledges financial sup-port from the National Science Foundation under GrantNo. DMR-1506634 and from the Aalto Science Institute(ASCI).4 ∗ [email protected] R. Mas-Ballesté, C. Gómez-Navarro, J. Gómez-Herrero,and F. Zamora, Nanoscale , 20 (2011). F. Xia, H. Wang, D. Xiao, M. Dubey, and A. Ramasub-ramaniam, Nature Photonics , 899 (2014). F. Wang, Z. Wang, Q. Wang, F. Wang, L. Yin, K. Xu,Y. Huang, and J. He, Nanotechnology , 292001 (2015). H. Liu, Y. Du, Y. Deng, and P. D. Ye, Chemical SocietyReviews , 2732 (2015). K. F. Mak and J. Shan, Nature Photonics , 216 (2016). F. Shahzad, M. Alhabeb, C. B. Hatter, B. Anasori,S. Man Hong, C. M. Koo, and Y. Gogotsi, Science ,1137 (2016). D. Akinwande, C. J. Brennan, J. S. Bunch, P. Egberts,J. R. Felts, H. Gao, R. Huang, J.-S. Kim, T. Li, Y. Li,K. M. Liechti, N. Lu, H. S. Park, E. J. Reed, P. Wang,B. I. Yakobson, T. Zhang, Y.-W. Zhang, Y. Zhou, andY. Zhu, Extreme Mechanics Letters , 42 (2017). S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev,and A. Kis, Nature Reviews Materials , 17033 (2017). Z. Fan, P. Hirvonen, L. F. C. Pereira, M. M. Ervasti, K. R.Elder, D. Donadio, A. Harju, and T. Ala-Nissila, NanoLetters , 5919 (2017). H. Dong, P. Hirvonen, Z. Fan, and T. Ala-Nissila, PhysicalChemistry Chemical Physics , 24602 (2018). K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H.Castro Neto, Science (2016), 10.1126/science.aac9439. M.-Y. Li, C.-H. Chen, Y. Shi, and L.-J. Li, MaterialsToday , 322 (2016). K. Kang, K.-H. Lee, Y. Han, H. Gao, S. Xie, D. A. Muller,and J. Park, Nature , 229 (2017). M. P. Levendorf, C.-J. Kim, L. Brown, P. Y. Huang, R. W.Havener, D. A. Muller, and J. Park, Nature , 627(2012). Z. Liu, L. Ma, G. Shi, W. Zhou, Y. Gong, S. Lei, X. Yang,J. Zhang, J. Yu, K. P. Hackenberg, A. Badakhani, J.-C.Idrobo, R. Vajtai, J. Lou, and P. M. Ajayan, Nature Nan-otechnology , 119 (2013). G. H. Han, J. A. Rodríguez-Manzo, C.-W. Lee, N. J. Ky-bert, M. B. Lerner, Z. J. Qi, E. N. Dattoli, A. M. Rappe,M. Drndic, and A. T. C. Johnson, ACS Nano , 10129(2013). X. Ling, Y. Lin, Q. Ma, Z. Wang, Y. Song, L. Yu, S. Huang,W. Fang, X. Zhang, A. L. Hsu, Y. Bie, Y.-H. Lee, Y. Zhu,L. Wu, J. Li, P. Jarillo-Herrero, M. Dresselhaus, T. Pala-cios, and J. Kong, Advanced Materials , 2322 (2016). K. R. Elder, M. Katakowski, M. Haataja, and M. Grant,Physical Review Letters , 245701 (2002). K. R. Elder and M. Grant, Physical Review E , 051605(2004). H. Emmerich, H. Löwen, R. Wittkowski, T. Gruhn, G. I.Tóth, G. Tegze, and L. Gránásy, Advances in Physics ,665 (2012). L. Gránásy, G. I. Tóth, J. A. Warren, F. Podmaniczky,G. Tegze, L. Rátkai, and T. Pusztai, Progress in MaterialsScience , 100569 (2019). K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, andM. Grant, Physical Review B , 064107 (2007). K. R. Elder, Z.-F. Huang, and N. Provatas, Physical Re-view E , 011602 (2010). M. Greenwood, N. Ofori-Opoku, J. Rottler, andN. Provatas, Physical Review B , 064104 (2011). N. Ofori-Opoku, V. Fallah, M. Greenwood, S. Esmaeili,and N. Provatas, Physical Review B , 134105 (2013). D. Taha, S. K. Mkhonta, K. R. Elder, and Z.-F. Huang,Physical Review Letters , 255501 (2017). A. R. Balakrishna and W. C. W. Carter, Physical ReviewE , 043304 (2018). G. Kocher and N. Provatas, Physical Review Letters ,155501 (2015). C. Guo, J. Wang, Z. Wang, J. Li, Y. Guo, and Y. Huang,Soft Matter , 4666 (2016). N. Provatas and K. Elder,
Phase-field methods in materialsscience and engineering (John Wiley & Sons, 2011). A. Jaatinen,
Modeling materials with phase field crystalmodels , Ph.D. thesis, Aalto University (2010). V. Heinonen,
Phase field crystal models and fast dynamics ,Ph.D. thesis, Aalto University (2016). P. Hirvonen, M. M. Ervasti, Z. Fan, M. Jalalvand, M. Sey-mour, S. M. Vaez Allaei, N. Provatas, A. Harju, K. R.Elder, and T. Ala-Nissila, Physical Review B , 035414(2016). A. Jaatinen and T. Ala-Nissila, Journal of Physics: Con-densed Matter , 205402 (2010). P. Sutter, R. Cortes, J. Lahiri, and E. Sutter, Nano Letters , 4869 (2012). P. M. Chaikin, T. C. Lubensky, and T. A. Witten, PhysicsToday , 82 (1995). O. V. Yazyev and S. G. Louie, Physical Review B ,195420 (2010). G. M. L. Boissonière and R. Choksi, Modelling and Sim-ulation in Materials Science and Engineering , 035001(2018). C. Lee, X. Wei, J. W. Kysar, and J. Hone, Science ,385 (2008). A. Falin, Q. Cai, E. J. Santos, D. Scullion, D. Qian,R. Zhang, Z. Yang, S. Huang, K. Watanabe, T. Taniguchi, et al. , Nature Communications , 15815 (2017). G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly,and J. van den Brink, Physical Review B , 073103(2007). J. Li, G. Gui, and J. Zhong, Journal of Applied Physics , 094311 (2008). S. Wang, Journal of the Physical Society of Japan ,064602 (2010). Y. Gao, Y. Zhang, P. Chen, Y. Li, M. Liu, T. Gao, D. Ma,Y. Chen, Z. Cheng, X. Qiu, W. Duan, and Z. Liu, NanoLetters , 3439 (2013). R. Drost, A. Uppstu, F. Schulz, S. K. Hämäläinen, M. Er-vasti, A. Harju, and P. Liljeroth, Nano Letters , 5128(2014). J. Lu, L. C. Gomes, R. W. Nunes, A. H. Castro Neto, andK. P. Loh, Nano Letters , 5133 (2014). M. Liu, Y. Li, P. Chen, J. Sun, D. Ma, Q. Li, T. Gao,Y. Gao, Z. Cheng, X. Qiu, Y. Fang, Y. Zhang, and Z. Liu,Nano Letters , 6342 (2014). R. Drost, S. Kezilebieke, M. M. Ervasti, S. K. Hämäläinen,F. Schulz, A. Harju, and P. Liljeroth, Scientific reports ,16741 (2015). Y. Liu, X. Zou, and B. I. Yakobson, ACS Nano , 7053(2012). A d ∩ B c A c ∩ B c A d ∩ B d A c ∩ B d FIG. 13. A pictorial showing an example of disordered andordered regions of the two components A and B. As a conse-quence of V ( A d ) = V ( B c ) (and V ( A c ) = V ( B d ) ), the areas V ( A c ∩ B c ) and V ( A d ∩ B d ) are equal. See text for details. The energy F needs to be relaxed with respect to L ⊥ and L (cid:107) to relieve possible mechanical stresses. Y. Liu, A. Dobrinsky, and B. I. Yakobson, Physical ReviewLetters , 235502 (2010). L. Liu, J. Park, D. A. Siegel, K. F. McCarty, K. W. Clark,W. Deng, L. Basile, J. C. Idrobo, A.-P. Li, and G. Gu,Science , 163 (2014). Y. Liu, S. Bhowmick, and B. I. Yakobson, Nano Letters , 3113 (2011). S. Bhowmick, A. K. Singh, and B. I. Yakobson, The Jour-nal of Physical Chemistry C , 9889 (2011). J. Zhang, W. Xie, X. Xu, S. Zhang, and J. Zhao, Chem-istry of Materials , 5022 (2016). A. H. King, Interface Science , 251 (1999). S. G. Srinivasan, J. W. Cahn, H. Jónsson, and G. Kalonji,Acta Materialia , 2821 (1999). P. Hirvonen, Z. Fan, M. M. Ervasti, A. Harju, K. R. Elder,and T. Ala-Nissila, Scientific reports , 4754 (2017). V. Heinonen, C. V. Achim, K. R. Elder, S. Buyukdagli,and T. Ala-Nissila, Phys. Rev. E , 032411 (2014). Appendix A: Phase separation
In this Appendix we examine the smoothed phase sep-arating part of the energy F ct = (cid:90) Ω d r [ (cid:15) AB η A ( r ) η B ( r )] (A1)appearing in Eq. (1). Here Ω is the domain that can bethought of as a box with a finite volume and periodicboundary conditions (flat torus). We neglect the contri-bution of the disorder–crystalline boundary and assumethat η A and η B take constant values specific to the phase.Let X i = { r ∈ Ω : η X ( r ) = η ( i ) X } , (A2)where X is the component label A or B and i ∈ { d , c } corresponding to disordered and crystalline phases. Wedefine V ( X ) as the volume (area) of set X . A diagram ofthe setup is shown in Fig. A. The system is set up suchthat V ( A d ) = V ( B c ) . Since X d and X c split Ω perfectly,also V ( A c ) = V ( B d ) . This can be done by choosing theoverall number of constituents A and B correctly.Now the energy contribution from Eq. (A1) is F ct = (cid:15) AB (cid:20)(cid:90) A c ∩ B c d r η (c) A η (c) B + (cid:90) A c ∩ B d d r η (c) A η (d) B + (cid:90) A d ∩ B c d r η (d) A η (c) B + (cid:90) A d ∩ B d d r η (d) A η (d) B (cid:21) = (cid:15) AB (cid:104) V ( A c ∩ B c ) η (c) A η (c) B + V ( A c ∩ B d ) η (c) A η (d) B + V ( A d ∩ B c ) η (d) A η (c) B + V ( A d ∩ B d ) η (d) A η (d) B (cid:105) . (A3)Any set Y ⊂ Ω can be divided such that V ( Y ) = V ( Y ∩ X l ) + V ( Y ∩ X s ) because X l and X s split Ω . From this it follows that V ( A c ∩ B d ) = V ( A c ) − V ( A c ∩ B c ) and V ( A d ∩ B c ) = V ( B c ) − V ( A c ∩ B c ) . Moreover, V ( A c ∩ B d ) + V ( A c ∩ B c ) = V ( A c ) = V ( B d )= V ( A d ∩ B d ) + V ( A c ∩ B d ) , from which it follows that V ( A d ∩ B d ) = V ( A c ∩ B c ) . Now6 F ct = (cid:15) AB (cid:104) V ( A c ∩ B c )( η (c) A η (c) B + η (d) A η (d) B − η (c) A η (d) B − η (d) A η (c) B ) + V ( A c ) η (c) A η (d) B + V ( B c ) η (d) A η (c) B (cid:105) = (cid:15) AB (cid:104) V ( A c ∩ B c )( η (c) A − η (d) A )( η (c) B − η (d) B ) + V ( A c ) η (c) A η (d) B + V ( B c ) η (d) A η (c) B (cid:105) . (A4)The components A and B have a similar phase diagramin the sense that either η (c) A > η (d) A , η (c) B > η (d) B or η (c) A < η (d) A , η (c) B < η (d) B . From this it follows that ( η (c) A − η (d) A )( η (c) B − η (d) B ) > ,which implies that F ct ≥ (cid:15) AB (cid:104) V ( A c ) η (c) A η (d) B + V ( B c ) η (d) A η (c) B (cid:105) . (A5)Therefore, at the ground state V ( A c ∩ B c ) = 0 . Thisshows that F is minimized when different phases of thedifferent components appear together. Even if V ( A c ) (cid:54) = V ( B d ) , the areas V ( A c ∩ B c ) and V ( A d ∩ B d ) would beminimized. If (cid:15) AB < , V ( A c ∩ B c ) = 0 is maximized andthe crystalline phases of the constituents overlap. Appendix B: Elastic effects due to smoothed densityfields
In this Appendix we will study the non-local effectsdue to the introduction of the smoothed number densityfields η X ( X ∈ { A , B } ) in F (Eq. (1) ) The smoothednumber densities appear in the term F ct = (cid:90) d x [ (cid:15) AB η A ( x ) η B ( x )] (B1)that might contribute to excess elastic energy if the sys-tem is deformed. Throughout this section we assume ν i = 1 . This sets the length scale of the bulk oscillationsof the density fields.The smoothed fields are defined as η X ( x ) = (cid:90) d y [ G ( x − y ) n X ( y )] . (B2)This convolution gives rise to non-local self-interactions.Let ˆ f ( k ) = (cid:90) d x (cid:2) e − i k · x f ( x ) (cid:3) (B3)be the Fourier transform of f ( x ) . Now the inverse 2Dtransform is f ( x ) = 14 π (cid:90) d x (cid:104) e i k · x ˆ f ( x ) (cid:105) . The energy F ct can be written in terms of the Fouriertransforms as F ct = 14 π (cid:90) d k [ (cid:15) AB ˆ η ∗ A ( k )ˆ η B ( k )] (B4)by using the Plancherel theorem. Here ˆ η ∗ A is the complexconjugate of ˆ η A . The fields ˆ η X ( X ∈ { A , B } ) can beeasily expressed using the convolution theorem as ˆ η X = ˆ G ˆ n X , (B5)where ˆ G is the Fourier transform of the Gaussian convo-lution kernel, also a Gaussian ˆ G ( k ) = e − ( γk ) . (B6)Here k = | k | and γ gives the length scale of the smooth-ing. Notice that ˆ G is real. Now F ct = (cid:15) AB π (cid:90) d k (cid:104) ˆ G ( k ) ˆ n ∗ A ( k )ˆ n B ( k ) (cid:105) . (B7)We consider deformations of the form k − p ( k ) , with | p | (cid:28) | k | . As an example, for a uniform compression of5% , p ≈ . k . Now F ct becomes F ct = (cid:15) AB π (cid:90) d k (cid:104) ˆ G ( k ) ˆ n ∗ A ( k − p )ˆ n B ( k − p ) (cid:105) . Making a change of variables k → k + p gives F ct = (cid:15) AB π (cid:90) d k ν ( k ) (cid:104) ˆ G ( k + p ) ˆ n ∗ A ( k )ˆ n B ( k ) (cid:105) , (B8)where ν is the change in the volume element that is givenby the determinant of the Jacobian I + ∇ p .The fields ˆ n X have non-zero structure at the nearestneighbor length scale (PFC fluctuations) and close to k =0 (order–disorder boundaries). The length scale given by /k = 1 corresponds to nearest neighbor distance of thePFC lattice and γ is chosen such that ˆ G ( k = 1) (cid:28) im-plying that deformations at this length scale do not con-tribute to F ct . We will investigate the other importantregime, where k is small. Let F ct = (cid:82) d k f ct . Expanding ˆ G around k gives f ct ≈ (cid:15) AB π ν ˆ n ∗ A ˆ n B (cid:88) i,j δ i δ j ∂ ij ˆ G ( k ) ≈ (cid:15) AB π ˆ n ∗ A ˆ n B (cid:2) p · k ) γ − | p | γ (cid:3) ˆ G ( k ) . Here we have used the fact that for small p , ν ≈ ∇ · p and assume that the part proportional to ∇ · p is7much smaller than unity. Also, the system is initially inequilibrium implying that ˆ G ( k + p ) has to be expandedup to second order. Now the excess part of the energy at k is ∆ f ct := f ct − f ct | p =0 = (cid:15) AB γ π (cid:2) γ ( p · k ) − | p | (cid:3) ˆ G ( k ) ˆ n ∗ A ˆ n B . (B9)In order to calculate the contribution to the excesselastic energy due to an interface, we assume that ˆ n A and ˆ n B vary only in one direction and are peaked around k = 0 . Furthermore, we can estimate ˆ n ∗ A ˆ n B < φ , where φ is the amplitude of the one-mode oscillations in thecrystal. More precisely ˆ n ∗ A ( k )ˆ n B ( k ) = 2 πδ ( k y )ˆ n ∗ A ( k x )ˆ n B ( k x )( k ) < πδ ( k y ) φ . The Fourier amplitudes due to the interfaces should besignificantly smaller than the Fourier amplitude of thebulk oscillations. The excess energy ∆ f ct is maximizedfor parallel k and p . Let us assume that p ( k ) < δk x ,where k x is the component of k parallel to the interfaceand δ is small. Now we can estimate the contribution ofthe interface per interface length as ∆ f int ct < (cid:15) AB δ γ φ π (cid:90) ∞−∞ d k x (cid:2) γ k x − k x (cid:3) ˆ G ( k x ) = (cid:15) AB δ φ √ πγ . (B10)We can compare ∆ f int ct to the bulk elastic excess en-ergy. The bulk elastic energy density due to component A is f el = lim V Ω →∞ V Ω (cid:90) Ω d x (cid:20) β A n A (1 + ∇ ) n A (cid:21) , (B11)where Ω is a compact domain that is taken to infinityand V Ω is its area. We consider the contribution of thebulk oscillations and set n A = φ (cid:88) j e i q j · x , (B12)where q j are the principal reciprocal lattice vectors with q j = 1 . Now ˆ n A = 2 πφ (cid:88) j δ ( k − q j ) . (B13) We evaluate f el in Fourier space as f el = lim V Ω →∞ β A π V Ω (cid:90) d k | ˆ n A | ˆ L ( k ) , (B14)where L ( k ) = (1 − k ) (B15)is the Fourier transform of the pattern forming operator (1 + ∇ ) . We repeat the earlier calculation by replacing ˆ G with ˆ L . We use lim V Ω →∞ δ ( k − q j ) δ ( k − q i ) V Ω = δ ij δ ( k − q j ) . Now f el = β A φ (cid:90) d k L ( k + p ) (cid:88) j δ ( k − q j ) . (B16)Up to second order in p ∆ f el = f el − f el | p =0 = 2 β A φ (cid:88) j [ q j · p ( q j )] . (B17)Let us consider small linear deformations with p ( k ) = δ Jk , with some matrix J with squared eigenvalues λ + λ = Tr( J ) = 1 and some small δ . For a hexagonallattice, the reciprocal lattice vectors q j form a star with6-fold symmetry. It can be shown that ∆ f el = 32 δ β A φ (cid:2) Tr( J ) + (Tr J ) + Tr ( J T J ) (cid:3) > β A φ δ ( λ + λ ) = 3 β A φ δ . (B18)In order to compare with ∆ f int ct of Eq. (B10), ∆ f el needsto be multiplied by the thickness of the interface, whichwe assume to be two lattice constants a = 4 π/ √ . Thisgives a ∆ f el > √ πβ A φ δ . (B19)Inserting β A = 1 , (cid:15) AB = 1 we get an estimate a ∆ f el ∆ f int ct > ,,