Phase field crystal model for heterostructures

TL;DR

This paper introduces an efficient phase-field crystal model for simulating atomistic configurations in 2D heterostructures, enabling exploration of material properties across large parameter spaces with realistic interface modeling.

Contribution

The paper develops a flexible phase-field crystal model for multi-species heterostructures, validated through benchmarking and applied to graphene-boron nitride systems.

Findings

Consistent elastic properties and lattice constants in simulations.

Realistic interface and crystal shape modeling.

Zigzag interfaces have lowest formation energy.

Abstract

Atomically thin 2-dimensional heterostructures are a promising, novel class of materials with groundbreaking properties. The possiblity of choosing the many constituent components and their proportions allows optimizing these materials to specific requirements. The wide adaptability comes with 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 and length scales related to physics of polycrystalline materials poses a challenge for computational studies. To this end, we present an efficient and flexible phase-field crystal model to describe the atomic configurations of multiple atomic species and phases coexisting in the same physical domain. We extensively benchmark the model for two-dimensional binary systems in terms of their elastic properties and phase boundary configurations and their energetics. As a concrete example, we demonstrate modelling lateral heterostructures of graphene and hexagonal boron nitride. We consider both idealized bicrystals and large-scale systems with random phase distributions. We find consistent relative elastic moduli and lattice constants, as well as realistic continuous interfaces and faceted crystal shapes. Zigzag-oriented interfaces are observed to display the lowest formation energy.

Figures

Tables

Table 1. Table 1: Summary of model parameters. They are listed below and their significance is explained. While not explicitly written in Eqs. ( 1 ) or ( 2 ), the average densities n ¯ i subscript ¯ 𝑛 𝑖 \bar{n}_{i} are important in controlling the relative stability of different phases.

$N$	number of density fields in the model
${\overline{n}}_i$	average density; controls the relative stability of
	different phases
${\alpha }_i$	temperature-related parameter; controls the
	diffuseness and facetedness of structures
${\beta }_i$	controls the elastic moduli
${\gamma }_i$	acts similarly to a chemical potential
${\delta }_i$	usually set to unity in PFC models
${\nu }_i$	wavenumber that sets the length scale
${\alpha }_{ij}$	controls the alignment of two lattices at their mutual
	interface but also introduces so-called "weak
	oscillations" (see Sec. III.2)
${\beta }_{ij}$	needed for h-BN Taha et al. (2017)
${\gamma }_{ij}$	needed for h-BN Taha et al. (2017)
${ϵ}_{ij}$	couples smoothed densities; controls phase separation
${\nu }_{ij}$	needed for h-BN Taha et al. (2017)
$\sigma$	spectral spread of the Gaussian smoothing kernel $G$

Table 2. Table 2: Average densities for lattice-mismatched bicrystalline heterostructures. The mismatch is indicated by λ 2 subscript 𝜆 2 \lambda_{2} .

${\lambda }_2$	${\overline{n}}_1^{\text{(c)}}$	${\overline{n}}_1^{\text{(d)}}$	${\overline{n}}_2^{\text{(c)}}$	${\overline{n}}_2^{\text{(d)}}$
1.0	0.19	0.51	0.19	0.51
1.05	0.16	0.55	0.12	0.55
1.1	0.24	0.47	0.21	0.48
1.2	0.28	0.43	0.23	0.45

Table 3. Table 3: Set of parameters for lateral heterostructures of graphene and h-BN. Note that α i = α i ​ j subscript 𝛼 𝑖 subscript 𝛼 𝑖 𝑗 \alpha_{i}=\alpha_{ij} when i = j 𝑖 𝑗 i=j and similarly for the other parameters.

$N$
3
${\alpha }_{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
${\beta }_{ij}$	$i=1$	$i=2$	$i=3$
$j=1$	2.25	–	–
$j=2$	0.0	1.0	–
$j=3$	0.0	0.02	1.0
${\gamma }_{ij}$	$i=1$	$i=2$	$i=3$
$j=1$	0.0	–	–
$j=2$	0.0	0.0	–
$j=3$	0.0	0.3	0.0
${\delta }_{ij}$	$i=1$	$i=2$	$i=3$
$j=1$	2.25	–	–
$j=2$	–	1.0	–
$j=3$	–	–	1.0
${ϵ}_{ij}$	$i=1$	$i=2$	$i=3$
$j=1$	–	–	–
$j=2$	-0.8	–	–
$j=3$	-0.8	0.0	–
${\lambda }_{ij}=1/{\nu }_{ij}$	$i=1$	$i=2$	$i=3$
$j=1$	1.0	–	–
$j=2$	0.0	1.018	–
$j=3$	0.0	1.018	1.018
${\overline{n}}_i^{\left(j\right)}$	$i=1$	$i=2$	$i=3$
$j=$ c	0.31	-0.32	-0.32
$j=$ d	0.66	-0.65	-0.65

Equations

$$\begin{array}{r}F=\int \mathrm{d}\mathbf{r}(\sum _{i=1}^N(\frac{{\alpha }_i}{2}{n}_i^2+\frac{{\beta }_i}{2}{n}_i{\left({\nu }_i^2+{\nabla }^2\right)}^2{n}_i\\ +\frac{{\gamma }_i}{3}{n}_i^3+\frac{{\delta }_i}{4}{n}_i^4)\\ +\sum _{i=1}^{N-1}\sum _{j=i+1}^N({\alpha }_{ij}{n}_i{n}_j+{\beta }_{ij}{n}_i{\left({\nu }_{ij}^2+{\nabla }^2\right)}^2{n}_j\\ +\frac{{\gamma }_{ij}}{2}\left(n_i^2{n}_j+n_i{n}_j^2\right)+{ϵ}_{ij}{\eta }_i{\eta }_j)).\end{array}$$

$$\begin{array}{rl}\frac{\partial n_i}{\partial t}=& {\nabla }^2\frac{\delta F}{\delta n_i}\\ +{\nabla }^2& ({\alpha }_i{n}_i+{\beta }_i{\left({\nu }_i^2+{\nabla }^2\right)}^2{n}_i+{\gamma }_i{n}_i^2+{\delta }_i{n}_i^3\\ & +\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N({\alpha }_{ij}{n}_j+{\beta }_{ij}{\left({\nu }_{ij}^2+{\nabla }^2\right)}^2{n}_j\\ & +\frac{{\gamma }_{ij}}{2}\left(2n_i{n}_j+n_j^2\right)+{ϵ}_{ij}G\ast {\eta }_j))\end{array}$$

$$f_{\text{e}}=\frac{C_{11}}{2}\left({\epsilon }_x^2+{\epsilon }_y^2\right)+C_{12}{\epsilon }_x{\epsilon }_y,$$

$$B=\frac{C_{11}+C_{12}}{2},$$

$$\mu =\frac{C_{11}-C_{12}}{2},$$

$$Y_{\text{2D}}=\frac{4B\mu }{B+\mu },$$

$$\nu =\frac{B-\mu }{B+\mu },$$

$$C_{11}=9\sum _i{\beta }_i{\varphi }_i^2$$

$$\begin{array}{rl}{\varphi }_i=& \frac{1}{15{\delta }_i}({\gamma }_i-3{\delta }_i{\bar{n}}_i\\ & -\sqrt{{\gamma }_i^2-15{\alpha }_i{\delta }_i+12{\delta }_i{\bar{n}}_i\left(2{\gamma }_i-3{\delta }_i{\bar{n}}_i\right)})\end{array}$$

$$C_{12}=C_{11}/3.$$

$$F=f{L}_{\perp }{L}_{\parallel }=f^{\ast }{L}_{\perp }{L}_{\parallel }+2\gamma L_{\parallel }+4\delta ,$$

$$f{L}_{\perp }=f^{\ast }{L}_{\perp }+2{\gamma }^{\ast },$$

$$\gamma =2{\gamma }_{\text{ZZ}}\sin \left(\theta \right)+2{\gamma }_{\text{AC}}\sin \left(30^{\circ }-\theta \right),$$

$$F_{\text{ct}}={\int }_{\Omega }\mathrm{d}\mathbf{r}\left[{ϵ}_{\mathrm{AB}}{\eta }_A(\mathbf{r}){\eta }_B(\mathbf{r})\right]$$

$$X_i=\{\mathbf{r}\in \Omega :{\eta }_X(\mathbf{r})={\eta }_X^{(i)}\},$$

$$\begin{array}{rl}{F}_{\text{ct}}& ={ϵ}_{\mathrm{AB}}\left[{\int }_{A_{\text{c}}\cap B_{\text{c}}}\mathrm{d}\mathbf{r}{\eta }_A^{\text{(c)}}{\eta }_B^{\text{(c)}}+{\int }_{A_{\text{c}}\cap B_{\text{d}}}\mathrm{d}\mathbf{r}{\eta }_A^{\text{(c)}}{\eta }_B^{\text{(d)}}+{\int }_{A_{\text{d}}\cap B_{\text{c}}}\mathrm{d}\mathbf{r}{\eta }_A^{\text{(d)}}{\eta }_B^{\text{(c)}}+{\int }_{A_{\text{d}}\cap B_{\text{d}}}\mathrm{d}\mathbf{r}{\eta }_A^{\text{(d)}}{\eta }_B^{\text{(d)}}\right]\\ & ={ϵ}_{\mathrm{AB}}\left[V(A_{\text{c}}\cap B_{\text{c}}){\eta }_A^{\text{(c)}}{\eta }_B^{\text{(c)}}+V(A_{\text{c}}\cap B_{\text{d}}){\eta }_A^{\text{(c)}}{\eta }_B^{\text{(d)}}+V(A_{\text{d}}\cap B_{\text{c}}){\eta }_A^{\text{(d)}}{\eta }_B^{\text{(c)}}+V(A_{\text{d}}\cap B_{\text{d}}){\eta }_A^{\text{(d)}}{\eta }_B^{\text{(d)}}\right].\end{array}$$

$$V(Y)=V(Y\cap X_l)+V(Y\cap X_s)$$

$$V(A_{\text{c}}\cap B_{\text{d}})=V(A_{\text{c}})-V(A_{\text{c}}\cap B_{\text{c}})$$

$$V(A_{\text{d}}\cap B_{\text{c}})=V(B_{\text{c}})-V(A_{\text{c}}\cap B_{\text{c}}).$$

$$\begin{array}{rl}V(A_{\text{c}}\cap B_{\text{d}})+V(A_{\text{c}}\cap B_{\text{c}})& =V(A_{\text{c}})=V(B_{\text{d}})\\ & =V(A_{\text{d}}\cap B_{\text{d}})+V(A_{\text{c}}\cap B_{\text{d}}),\end{array}$$

$$V(A_{\text{d}}\cap B_{\text{d}})=V(A_{\text{c}}\cap B_{\text{c}}).$$

$$\begin{array}{rl}{F}_{\text{ct}}& ={ϵ}_{\mathrm{AB}}\left[V(A_{\text{c}}\cap B_{\text{c}})({\eta }_A^{\text{(c)}}{\eta }_B^{\text{(c)}}+{\eta }_A^{\text{(d)}}{\eta }_B^{\text{(d)}}-{\eta }_A^{\text{(c)}}{\eta }_B^{\text{(d)}}-{\eta }_A^{\text{(d)}}{\eta }_B^{\text{(c)}})+V(A_{\text{c}}){\eta }_A^{\text{(c)}}{\eta }_B^{\text{(d)}}+V(B_{\text{c}}){\eta }_A^{\text{(d)}}{\eta }_B^{\text{(c)}}\right]\\ & ={ϵ}_{\mathrm{AB}}\left[V(A_{\text{c}}\cap B_{\text{c}})({\eta }_A^{\text{(c)}}-{\eta }_A^{\text{(d)}})({\eta }_B^{\text{(c)}}-{\eta }_B^{\text{(d)}})+V(A_{\text{c}}){\eta }_A^{\text{(c)}}{\eta }_B^{\text{(d)}}+V(B_{\text{c}}){\eta }_A^{\text{(d)}}{\eta }_B^{\text{(c)}}\right].\end{array}$$

$${\eta }_A^{\text{(c)}}>{\eta }_A^{\text{(d)}},{\eta }_B^{\text{(c)}}>{\eta }_B^{\text{(d)}}$$

$${\eta }_A^{\text{(c)}}<{\eta }_A^{\text{(d)}},{\eta }_B^{\text{(c)}}<{\eta }_B^{\text{(d)}}.$$

$$F_{ct}\ge {ϵ}_{\mathrm{AB}}\left[V(A_{\text{c}}){\eta }_A^{\text{(c)}}{\eta }_B^{\text{(d)}}+V(B_{\text{c}}){\eta }_A^{\text{(d)}}{\eta }_B^{\text{(c)}}\right].$$

$$F_{\text{ct}}=\int \mathrm{d}\mathbf{x}\left[{ϵ}_{\mathrm{AB}}{\eta }_{\mathrm{A}}(\mathbf{x}){\eta }_{\mathrm{B}}(\mathbf{x})\right]$$

$${\eta }_{\mathrm{X}}(\mathbf{x})=\int \mathrm{d}\mathbf{y}\left[G(\mathbf{x}-\mathbf{y})n_{\mathrm{X}}(\mathbf{y})\right].$$

$$\hat{f}(\mathbf{k})=\int \mathrm{d}\mathbf{x}\left[e^{-i\mathbf{k}\cdot \mathbf{x}}f(\mathbf{x})\right]$$

$$f(\mathbf{x})=\frac{1}{4{\pi }^2}\int \mathrm{d}\mathbf{x}\left[e^{i\mathbf{k}\cdot \mathbf{x}}\hat{f}(\mathbf{x})\right].$$

$$F_{\mathrm{ct}}=\frac{1}{4{\pi }^2}\int \mathrm{d}\mathbf{k}\left[{ϵ}_{\mathrm{AB}}{\hat{\eta }}_{\mathrm{A}}^{\ast }(\mathbf{k}){\hat{\eta }}_{\mathrm{B}}(\mathbf{k})\right]$$

Full text

††preprint: APS/123-QED

I Introduction

Most scientifically and technologically important materials are composed of multiple atomic species and phases with different chemical compositions, lattice structures and elastic properties. Some everyday examples include wood, rock, metallic alloys, and concrete. Such three-dimensional (3D) materials have been used for thousands of years and efforts towards their development continue in the age of nanophysics. Some modern examples include e.g. fiber-reinforced polymers and semiconductor heterostructures. The past decade has seen the emergence of a completely new type of materials, the atomically thin two-dimensional (2D) materials. The extraordinary properties of single component 2D materials Mas-Ballesté et al. (2011); Xia et al. (2014); Wang et al. (2015); Liu et al. (2015); Mak and Shan (2016); Shahzad et al. (2016); Akinwande et al. (2017); Manzeli et al. (2017); Fan et al. (2017); Dong et al. (2018) can be widely enhanced and adjusted by considering their heterostructures that can either be stacked to form vertical multilayer heterostructures Novoselov et al. (2016); Li et al. (2016); Kang et al. (2017), or they can be grown within a single material sheet into a lateral heterostructure Levendorf et al. (2012); Liu et al. (2013); Han et al. (2013); Ling et al. (2016).

The properties of pure or single phase crystalline materials are determined by the complex networks of microscopic defects and grains. In contrast, for many multiphase composite materials, macroscopic continuum models may provide sufficiently accurate predictions of many of their properties. This suggests that the role of their microscopic structure is less important. Nevertheless, for semiconductor heterostructures, as well as for vertical and lateral 2D heterostructures, the atomic-level structure of their phase interfaces plays a major role as said structures are miniaturized to the nanoscale where interfacial effects are important.

Predicting the atomic-level structure between two or more orientationally, structurally and elastically mismatched phases is particularly difficult. The number of possible atomic configurations is essentially endless and conventional atomistic modelling techniques cannot simultaneously reach all the length and time scales relevant to their formation. The more recently developed phase field crystal (PFC) method allows examination of long, diffusive time scales corresponding to the slow evolution of microstructures, and offers atomic-level spatial resolution up to mesoscopic length scales Elder et al. (2002). PFC models 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\left[n_{i}\right]$ that is minimized by a periodic $n_{i}$. PFC models have been extensively applied to study various complex systems and processes such as grain boundaries, vacancy diffusion, coarsening of polycrystals, heteroepitaxial growth, yield strength and fracture Elder et al. (2002); Elder and Grant (2004); Emmerich et al. (2012); Gránásy et al. (2019). In particular, multicomponent PFC models that explicitly incorporate multiple density fields $n_{i}$ coupled together, have been developed and applied to study crystal structures composed of multiple atomic species Elder et al. (2007, 2010); Greenwood et al. (2011); Ofori-Opoku et al. (2013); Taha et al. (2017); Balakrishna and W. Carter (2018).

In this work we introduce spatially smoothed atomic density fields coupled to the atomic density fields $n_{i}$ that enables well-controlled phase separation and, therefore, facilitates modelling heterostructures and composite materials. Smoothed densities have been employed in PFC modelling recently for introducing a vapor phase Kocher and Provatas (2015) and for controlling liquid/solid interface energies Guo et al. (2016). Here we apply this modelling approach to 2D heterostructures composed of multiple elements. We carry out a systematic investigation by varying model parameters one by one to determine their influence on the general behavior of the model. More specifically, we introduce mismatch in both the elastic moduli and the lattice constants between the two materials, as well as experiment with different couplings between the two density fields. Finally, we assess the model’s suitability to study graphene–hexagonal boron nitride 2D heterostructures.

The rest of the paper is organized as follows: Section II lays out and discusses the heterostructure model and gives some practical details of our calculations. Section III presents our investigation of the general properties of the model using binary heterostructures. Section IV assesses 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 continuum methods for multiscale modelling of polycrystalline materials and their complex microstructures. PFC models allow simultaneous access to both atomic and mesoscopic length scales, as well as to long, diffusive time scales. Formation and evolution of microstructures take place in this time regime. Conventional PFC models use a smooth, periodic density field $n$ to describe crystalline systems. The length scale and lattice symmetries, as well as the elastic properties of the model can be matched with the target material. These properties are determined by a free energy functional $F\left[n\left(\mathbf{r}\right)\right]$ governing the energetics of the system. In the solid phase $F$ is minimized by a periodic $n$ whose symmetries depend on the formulation of $F$ and average density $\bar{n}$ Elder et al. (2002); Provatas and Elder (2011). A PFC model can incorporate multiple density fields coupled together to allow the study of more complex structures. Such models have been applied to study multicomponent materials such as 2D hexagonal boron nitride (h-BN) Taha et al. (2017).

In conventional multicomponent PFC models the periodically oscillating densities representing the solid phase of each component overlap and form interlocking, mixed lattices. To study heterostructures with well-controlled phase separation, we propose the following dimensionless free energy functional:

[TABLE]

Here, the first sum contains the ideal contributions of the $N$ density fields and the second, nested sum the contributions of the interactions between them. In the first sum, the quadratic and quartic terms comprise a double-well potential, the cubic term acts similarly to a chemical potential and the gradient term gives rise to periodic solutions and elastic behavior. We refer the reader to Refs. Elder et al. (2002); Provatas and Elder (2011); Jaatinen (2010); Heinonen (2016) for a more in-depth discussion of PFC formulation. In the second sum, the quadratic and cubic terms are local couplings between the different density fields, whereas the rest are nonlocal terms. The model parameters and their roles are summarized in Table 1.

The term $\int\epsilon_{ij}\eta_{i}\eta_{j}\mathop{}\!\mathrm{d}\mathbf{r}$ in Eq. (1) is the essential coupling responsible for controlled phase separation. This term can effectively drive $\bar{n}_{i}$ and $\bar{n}_{j}$ apart in the same physical domain such that one corresponds to the disordered and the other to a crystalline phase in the phase diagram. The fields $\eta_{i}$ are spatially smoothed $n_{i}$ where the atomic-level structures have been filtered out defined as $\eta_{i}=G\ast n_{i}$. Here the asterisk denotes a convolution and $G$ is a Gaussian smoothing kernel with the Fourier transform $\hat{G}\left(\mathbf{k}\right)=e^{-\left|\mathbf{k}\right|^{2}/\left(2\sigma^{2}\right)}$. In the present work we found that $\sigma=0.2$ corresponding to a length scale of approximately five lattice constants with $\nu_{i}=1$ sufficiently smooths out the atomic-level structure. To enable atomistically sharp interfaces and, moreover, to keep the number of parameters to be tuned to a minimum, we did not consider values $\sigma<0.2$. The influence of this coupling term is demonstrated in Sec. III.1 and its ability to drive phase separation is shown analytically in Appendix A.

In order to find the $n_{i}$ that minimize $F$, we used density conserving gradient descent given by

[TABLE]

where $\alpha_{ij}=\alpha_{ji}$ and similarly for the other parameters. The density conservation constraint is essential for stabilizing the heterostructures. Note that while $t$ is called sometimes time, it is a relaxation parameter that is not related to real dynamics in this work. We solved Eq. (2) numerically using a semi-implicit spectral method described in Ref. 30. This method allows computing the gradients and convolutions present in Eqs. (1) and (2) accurately and efficiently by using fast Fourier transform routines. Note that according to the convolution theorem, convolutions can be expressed as $f\ast g=\mathcal{F}^{-1}\left\{\hat{f}\hat{g}\right\}$, where $\mathcal{F}^{-1}$ and the carets indicate inverse and forward Fourier transforms, respectively. The following upper bounds for spatial and temporal discretizations were used for all the calculations: $\Delta x=0.75$, $\Delta y=0.75$, and $\Delta t=0.25$.

Finally, we also used a model system size optimization algorithm Hirvonen et al. (2016) to eliminate strain in our bicrystalline model systems of heterostructures. We did not apply the method to polycrystalline systems, since we did not attempt to extract equilibrium densities from them or to analyze them quantitatively here.

III Binary heterostructures

We begin by demonstrating some general properties of the present model for simple binary heterostructures with $N=2$ and denote the two density fields by $n_{1}$ and $n_{2}$. We vary certain model parameters to investigate their influence and will refer to the periodic or “crystalline regions” in $n_{i}$ by $n_{i}^{\textrm{(c)}}$ and, similarly, to the disordered regions by $n_{i}^{\textrm{(d)}}$. The “crystalline phase $i$” encompasses regions where $n_{i}^{\textrm{(c)}}$ and $n_{j}^{\textrm{(d)}}$ coincide. Similarly, the “mixed phase” (disordered phase) spans the regions where $n_{i}^{\textrm{(c)}}$ and $n_{j}^{\textrm{(c)}}$ ($n_{i}^{\textrm{(d)}}$ and $n_{j}^{\textrm{(d)}}$) coincide.

III.1 Influence of smoothed coupling

The crucial parameter here is $\epsilon_{12}$ in the coupling term for the smoothed density fields. With $\epsilon_{12}>0$, $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ repel each other, whereas with $\epsilon_{12}<0$, $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ attract each other; see Appendix A for an analytical treatise. For the other parameters, we chose $\left(\alpha_{i},\beta_{i},\nu_{i},\gamma_{i},\delta_{i},\alpha_{12},\beta_{12},\gamma_{12}\right)=\left(-0.3,1.0,1.0,0.0,1.0,0.0,0.0,0.0\right)$, for simplicity. This choice of model parameters is used throughout Sec. III unless stated otherwise. Note that this choice of model parameters is symmetric, i.e., $F\left(n_{1},n_{2}\right)=F\left(n_{2},n_{1}\right)$.

We used simple model systems where we initialized both $n_{1}$ and $n_{2}$ roughly 50% crystalline and 50% disordered with average densities $\bar{n}_{i}^{\textrm{(c)}}\approx 0.32$ and $\bar{n}_{i}^{\textrm{(d)}}\approx 0.38$. The initial structure for $n_{i}^{\textrm{(c)}}$ was obtained using an inverted hexagonal one-mode approximation Jaatinen and Ala-Nissila (2010). We also arranged $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ in partial overlap to force some changes in them during relaxation. Figure 1 illustrates the relaxed heterostructures obtained with $\epsilon_{12}=\pm 0.2$. Panels (a) - (c) visualize the systems and panels (d) - (f) plot corresponding density profiles. Panel (a) shows the initial state with $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ in partial overlap. Panel (b) gives the repulsive case with $\epsilon_{12}=0.2$ where $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ have pushed themselves apart to eliminate the mixed phase. Panel (c) depicts the attractive case with $\epsilon_{12}=-0.2$ where $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ have come to a full overlap forming a coexistence between a mixed and a disordered phase. Recall that $n_{1}$ and $n_{2}$ are coupled here only via $\eta_{1}$ and $\eta_{2}$ whereby the two atomic lattices do not interact. Consequently, the lattices can end up arbitrarily aligned such as here; see panel (c).

Next we considered polycrystalline heterostructures. The density fields $n_{1}$ and $n_{2}$ were initialized with white noise with $\bar{n}_{i}=0.35$. With $\epsilon_{12}=0.2$, a mixed phase emerges first, followed by delayed decomposition into the two separate crystalline phases. We, therefore, used $\epsilon_{12}=\pm 1.0$ to drive $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ apart ($+$) or together ($-$). Note, however, that if the coupling strength is increased significantly more, stripe phases Elder et al. (2002) may replace the crystalline ones as the most stable phase. Figure 2 demonstrates both the repulsive and attractive cases after 7 500 time units of relaxation. In the repulsive case shown in panel (a), $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ are well separated and there is no mixed or disordered phase. The interfaces between the two phases appear fuzzy and disordered as expected due to no interaction between the two underlying lattices. The attractive case is shown in panel (b) where $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ are in full overlap, yielding a patched coexistence between a mixed and a disordered phase. The arbitrary misorientations and translations between the lattices in $n_{1}^{\textrm{(c)}}$ and $n_{2}^{\textrm{(c)}}$ result in a multitude of Moiré patterns.

III.2 Influence of $\alpha_{ij}$

Next, we varied the quadratic coupling parameter $\alpha_{12}$ to study its influence on the interfaces between the two crystalline phases and on the heterostructures as a whole. Here, we fix $\epsilon_{12}=1.0$. We chose $\alpha_{12}<0$ to achieve commensurate alignment of the two crystalline lattices at their interface to ensure the continuity of the underlying honeycomb lattice there. A side effect of this coupling is that it causes $n_{i}^{\textrm{(c)}}$ to induce oscillations in $n_{j}^{\textrm{(d)}}$. The amplitude of such weak oscillations in $n_{j}^{\textrm{(d)}}$ should be constrained to keep the two crystalline phases from mixing together. Indeed, the weak oscillations can be viewed as slight intermixing of the different atomic species. Intermixing is common in metallic alloys and in doped semiconductors and has been observed in lateral heterostructures of graphene and hexagonal boron nitride as well Sutter et al. (2012).

While constrained by the amplitude of the weak oscillations induced, the magnitude of $\alpha_{12}$ should be maximized to ensure continuity even for lattice-mismatched or misoriented interfaces. We optimized $\alpha_{12}$ using bicrystalline heterostructures. We observed that the heterostructures are rendered unstable when $\alpha_{12}=-0.1$, but with $\alpha_{12}=-0.03$ they retain their stability while the amplitude of the weak oscillations remains negligible. Figure 3 illustrates the interface in a relaxed bicrystalline heterostructure with $\alpha_{12}=-0.03$, $\bar{n}_{i}^{\textrm{(c)}}=0.12$ and $\bar{n}_{i}^{\textrm{(d)}}=0.58$. It is clear both from the visualization of the heterostructure as well as from the density profiles below it that the honeycomb lattice is highly continuous from one phase to the other. Furthermore, the interface is atomically sharp with an approximate width of two lattice constants.

We further demonstrated the influence of $\alpha_{12}$ for polycrystalline heterostructures. We initialized $n_{1}$ and $n_{2}$ with white noise where $\bar{n}_{i}=0.35$. Figure 4 demonstrates a relaxed heterostructure. On a larger scale, the system resembles that shown in Fig. 2 (a), but here the interfaces between the two crystalline phases are better ordered and more continuous. The width of the interfaces appears small for all misorientations. Note also the weak oscillations in $n_{2}^{\textrm{(d)}}$ visible in panel (a).

III.3 Influence of $\beta_{i}$

The crystalline phase $i$ is present where $n_{i}^{\textrm{(c)}}$ and $n_{j}^{\textrm{(d)}}$ coincide. The elastic properties of said phase should be dictated by $n_{i}$, but $n_{j}$ can have a minor contribution as well. We demonstrate here to what extent the elastic properties of the two crystalline phases can be controlled via $\beta_{1}$ and $\beta_{2}$, the coefficients of the gradient terms in $F$ responsible for elastic contribution. The ability to control the elastic stiffness of both phases separately is essential when modelling realistic heterostructures. Note that we show in Appendix B that the smoothed densities $\eta_{i}$ have a negligible elastic contribution.

The contribution from a uniform elastic deformation to the free energy density is given by

[TABLE]

where $\varepsilon_{x}$ and $\varepsilon_{y}$ are the $x$ and the $y$ components of strain, and $C_{11}=C_{22}$ and $C_{12}=C_{21}$ are the stiffness coefficients. Furthermore, the bulk, shear and 2D Young’s moduli, as well as Poisson’s ratio are given by

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

respectively Chaikin et al. (1995).

We determined the elastic coefficients of the two crystalline phases separately by straining single-crystals of either phase in the small deformation limit. More precisely, we varied $-0.002\leq\varepsilon_{x}\leq 0.002$ and $-0.002\leq\varepsilon_{y}\leq 0.002$ independently. We fixed $\beta_{1}=1.0$, and varied $0.25\leq\beta_{2}\leq 4$. For $0.9\leq\beta_{2}\leq 1.3$, we used $\epsilon_{12}=1$, but, for $\beta_{2}\leq 0.7$ ($\beta_{2}\geq 2$), we had to adjust $0.5\leq\epsilon_{12}\leq 0.75$ ($1.5\leq\epsilon_{12}\leq 2$) to retain the stability of the heterostructures. We again used $\alpha_{12}=-0.03$ to include the weak oscillations. We determined the average densities $\bar{n}_{i}^{\textrm{(c)}}$ and $\bar{n}_{i}^{\textrm{(d)}}$ in equilibrium by relaxing bicrystalline heterostructures and by extracting the average densities from the middle of the crystalline phases.

Figure 5 shows the 2D Young’s modulus as a function of $\beta_{2}$ for both crystalline phases. For crystalline phase 1, the modulus is essentially unaffected by $\beta_{2}$, i.e., the corresponding linear fit has a negligible slope. In contrast, the modulus for the crystalline phase 2 is linearly proportional to $\beta_{2}$ with a slope of 0.17. Independent control of the elastic stiffness of either of the crystalline phases appears straightforward. In addition, for each value of $\beta_{2}$, we observed $C_{12}\approx C_{11}/3$, whereby $\nu\approx 1/3$. This is a feature common to many simple PFC models Hirvonen et al. (2016).

Finally, we compared the numerical results against an analytical prediction Provatas and Elder (2011) where for one crystalline phase

[TABLE]

with an amplitude

[TABLE]

and

[TABLE]

Here, $C_{11}$ is given simply as a sum over the individual density fields. The values obtained for $Y_{\textrm{2D}}$ using the analytical expressions above are also plotted in Fig. 5 for the crystalline phase 2. The numerical and analytical results are in very good agreement. Note that the amplitude of the weak oscillations in the crystalline phase 1 is roughly an order of magnitude lower than that of the oscillations in the crystalline phase 2. Since $C_{11}\propto\phi_{i}^{2}$, the influence of these oscillations is negligible.

III.4 Influence of $\nu_{i}$

As a final demonstration of our model for binary heterostructures, we introduced lattice mismatch between the two crystalline phases via $\nu_{i}=1/\lambda_{i}$. We fixed $\lambda_{1}=1.0$ and varied $\lambda_{2}$ through values 1.05, 1.1 and 1.2. We ensured the stability of the heterostructures by choosing $\epsilon_{12}=$ 0.9, 0.75 and 0.5, respectively, and by setting $\alpha_{12}=-$0.03. Here $\beta_{i}=1.0$ for simplicity.

We first considered bicrystalline model systems either with different permutations of armchair and zigzag edges along the interface or with symmetrically tilted crystals with a tilt angle $2\theta=\theta-\left(-\theta\right)$, where $5.5^{\circ}\leq 2\theta\leq 55.8^{\circ}$. We considered two different strains. For unstrained systems, the periodicities of both bicrystal halves were matched separately with the periodic domain whereby the lattice mismatch is accommodated by misfit dislocations. For strained systems, both bicrystal halves were initialized with an average lattice constant and were again matched with the domain whereby the lattice mismatch is accommodated via elastic deformation. The average densities for the different strain and mismatch cases are given in Table 2. Note that said densities were chosen to yield an approximate 1:1 coexistence between the two crystalline phases. For reference, we considered here also $\lambda_{2}=1.0$ with $\epsilon_{12}=1.0$.

Overall, the phase interfaces obtained displayed well-defined structures. While misorientation and lattice mismatch introduce defects, extensively fuzzy and ill-defined structures are rare. In addition, the vast majority of the highly strained systems remained stable during relaxation. Figure 6 offers a representative sample of the structures obtained and shows a comparison between different strain and mismatch cases.

The first row of panels in Fig. 6 demonstrates zero-misorientation armchair-armchair interfaces between lattices of varying mismatch. We observed perfect hexagonal order for the reference and strained cases. Indeed, the strained heterostructures retained their stability without experiencing any stress-relieving reconstructions such as via subsequent nucleation, creation and annihilation of dislocations. In the unstrained structures, we observed periodic arrays of point-like misfit dislocations along the interface. We obtained similar structures for zigzag-zigzag interfaces.

The second row of panels in Fig. 6 gives low-misorientation interfaces between symmetrically tilted lattices of varying mismatch. Here, the tilt angle $2\theta\approx 51.4^{\circ}$. All structures display fairly periodic arrays of dislocations. For corresponding graphene grain boundaries, alternatingly slanted dislocations are expected Yazyev and Louie (2010); Hirvonen et al. (2016). For the reference structure, the highly symmetric initial state used has lead to extended defect structures without such symmetry breaking. The low-strain structure with $\lambda_{2}=1.05$ indeed displays alternatingly slanted dislocations, whereas the high-strain structure with $\lambda_{2}=1.2$ again does not, due to having achieved some strain-relief via annihilations of dislocations. The unstrained structures appear very similar to the corresponding strained ones.

The third row of panels in Fig. 6 depicts high-misorientation armchair-zigzag interfaces between lattices of varying mismatch. Despite the extreme misorientation between the two crystalline phases, the present model performs well in stitching them together with rather well-defined atomic-level structures. In fact, the lattice-mismatched structures do not appear visibly fuzzier than the reference.

Finally, we simulated the growth of polycrystalline heterostructures with the aforementioned lattice mismatches. The density fields were initialized with white noise where $\bar{n}_{i}=[\bar{n}_{i}^{\textrm{(c)}}+\bar{n}_{i}^{\textrm{(d)}}]/2$.

Figure 7 gives examples of the lattice-mismatched polycrystalline heterostructures obtained. Despite the various misorientations between the two crystalline phases at their interfaces, the model performs well in localizing the mismatch into point-like dislocations. Some interfaces, especially those with larger mismatch, appear somewhat diffuse, but are well comparable to some single-phase PFC grain boundaries; cf. Ref. 38, for example. For the cases with $\lambda_{2}=1.05$ and $\lambda_{2}=1.1$, the mismatch is minor and the interfaces appear highly continuous; see panels (a) and (b), respectively. There are a number of dislocations along the interfaces in both heterostructures, but many are due to lattice misorientation. For the case with $\lambda_{2}=1.2$, the interfaces are still fairly continuous, but display several regions with somewhat fuzzy features; see panel (c). These regions seem to coincide with greater interfacial curvature. Last, it also appears that the mobility of the interfaces decreases with increasing mismatch. This is evident from the noticeably smaller domain sizes in the heterostructure with $\lambda_{2}=1.2$; note that all three have been relaxed for the same 250 000 time units.

IV Application to lateral graphene–hexagonal boron nitride heterostructures

In this section, we consider a three-component model for lateral 2D heterostructures between graphene and h-BN (G–h-BN). We focus here on demonstrating the suitability of the model to qualitative modelling of G–h-BN. Finding optimal model parameters for quantitatively accurate modelling of G–h-BN (involving, e.g., fitting to interfacial formation energies) will be presented in future work.

IV.1 Model requirements and parameters

To model G–h-BN, we set $N=3$ and chose to model the graphene phase with $n_{1}^{\textrm{(c)}},n_{2}^{\textrm{(d)}}$ and $n_{3}^{\textrm{(d)}}$ and the h-BN phase with $n_{1}^{\textrm{(d)}},n_{2}^{\textrm{(c)}}$ and $n_{3}^{\textrm{(c)}}$. The parameters for $n_{2},n_{3}$ and their mutual couplings were adopted from Ref. 26. The other parameters were chosen by trial and error by varying them one at a time. We use the following criteria, guiding principles and simplifying assumptions:

• •

Start by choosing parameters for $n_{1}$ that are similar to those for $n_{2}$ and $n_{3}$.

• •

Assume that boron and nitrogen atoms are interchangeable, i.e., $F\left(n_{1},n_{2},n_{3}\right)=F\left(n_{1},n_{3},n_{2}\right)$.

• •

Match the relative Young’s moduli and lattice constants to the experimentally and theoretically determined values $Y_{\textrm{h-BN}}/Y_{\textrm{G}}\approx 0.87$ Lee et al. (2008); Falin et al. (2017) and $a_{\textrm{h-BN}}/a_{\textrm{G}}\approx 1.018$ Giovannetti et al. (2007); Li et al. (2008); Wang (2010), respectively.

• •

Require sharp, continuous interfaces and faceted crystal shapes Sutter et al. (2012); Liu et al. (2013); Gao et al. (2013); Drost et al. (2014); Lu et al. (2014); Liu et al. (2014a); Drost et al. (2015).

• •

Retain the stability of the heterostructures and limit the amplitude of weak oscillations.

Table 3 gives a set of model parameters that was found to satisfy the criteria listed above. Most importantly, this choice of parameters yielded $Y_{\textrm{h-BN}}/Y_{\textrm{G}}=0.84$ and $a_{\textrm{h-BN}}/a_{\textrm{G}}=1.021$ in fair agreement with the target values. In the following, we demonstrate how the model behaves and how it fulfills the other criteria above.

IV.2 Atomic configurations

Figure 8 demonstrates a zigzag-oriented interface in a bicrystalline G–h-BN system. Panel (a) gives both the distribution of the two phases and a representation of the same structure by $m=n_{1}+n_{2}+n_{3}$ for a clearer illustration of the atomic-level structure. Note that in the latter the h-BN phase appears brighter facilitating identification of the two phases in such figures. The interface displays perfect hexagonal order and is again atomistically sharp. Panel (b) gives the profiles of $n_{i}$ and $\eta_{i}$ along the horizontal periodic edge of the system and shows that the amplitude of the weak oscillations is small. Note that in contrast to the binary heterostructures considered in Sec. III, here the average densities $\bar{n}_{2}$ and $\bar{n}_{3}$ are negative.

Figure 9 demonstrates a large polycrystalline G–h-BN system grown from white noise where $\bar{n}_{i}=[\bar{n}_{i}^{\textrm{(c)}}+\bar{n}_{i}^{\textrm{(d)}}]/2$. In panel (a), an overview of an approx. 50 $\times$ 50 nm2 system is given by a coarse-grained representation where the graphene (h-BN) phase appears cyan (red). After a relaxation of $2.5\times 10^{6}$ time units, the heterostructure assumes configurations typical to spinodal decomposition in binary systems. Coarsening is slow because we have strived here for stable sharp crystalline structures instead of diffuse high-temperature ones. Panels (b) - (e) show a blow-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, excluding few fuzzy patches. Panels (d) and (e) show $n_{1}$ and $n_{2}$, respectively. Both appear faceted and have sharp interfaces with a primary (secondary) preference for the zigzag (armchair) direction. Weak oscillations are also visible in both panels.

The coarsening of G–h-BN was found slow with the present model and set of parameters. Moreover, concurrent nucleation and growth is not how said heterostructures are produced in practice Levendorf et al. (2012); Liu et al. (2013); Han et al. (2013); Ling et al. (2016). We demonstrated preparing more realistic model systems of random polycrystalline G–h-BN with larger phase domain and grain sizes. For initialization, we used Voronoi grain structures with random seed points, crystal orientations and phases Dong et al. (2018) and relaxed for 25 000 time units for local relaxation of the interfaces and grain boundaries. Figure 10 gives an overview of one such system where the initial, large-scale Voronoi structure has remained essentially unchanged as shown by the coarse-grained depiction of the system in panel (a). Panels (b) - (e) show the total density $m$ illustrating the atomic level structure of selected interfaces and boundaries. Panel (b) displays two triple junctions, one within h-BN and the other between graphene and h-BN, connected by an inversion boundary within h-BN. In h-BN, an inversion boundary is formed between two crystals with a misorientation of $60^{\circ}$ as the ordering of boron and nitrogen becomes inverted in one crystal with respect to the other Taha et al. (2017); Dong et al. (2018). The interfaces between graphene and h-BN have small-to-intermediate misorientations, whereas the grain boundaries within h-BN are large-angle boundaries. The interfaces and grain boundaries appear disordered but have fairly well-defined atomic level structures. The inversion boundary is formed by a perfect 4|8 chain as expected Liu et al. (2012); Taha et al. (2017). Panels (c) - (e) show longer interfaces between graphene and h-BN: (c) a large-misorientation interface, (d) a small-misorientation zigzag interface and (e) a small-misorientation armchair interface. While the large-angle interface shown in panel (c) appears disordered, all interfaces display well-defined atomic level structures.

IV.3 Interface energies

Finally, we investigated the relative stability of G–h-BN interfaces in different lattice directions by studying their formation energies. We considered 12 different interface angles $0^{\circ}\leq\theta\leq 30^{\circ}$, where $\theta=0^{\circ}$ corresponds to armchair and $\theta=30^{\circ}$ to zigzag interfaces. For simplicity, we considered here only strained configurations with perfect honeycomb order and no misfit dislocations along the interfaces (see Fig. 11 for an example) to avoid very large system sizes. We assumed all interfaces to be composed of zigzag and armchair segments as shown in Fig. 11. The model appears to yield stepped interfaces, typically with minimal segment lengths $L_{\textrm{ZZ}}$ and $L_{\textrm{AC}}$. Interfaces initialized with longer segments are also at least metastable.

Assuming a sufficiently large bicrystalline heterostructure with dimensions $L_{\bot}$ perpendicular and $L_{\parallel}$ parallel to the two interfaces, the total formation energy of the system can be written as

[TABLE]

where $f$ is the free energy density per unit area obtained by evaluating Eq. (1) and by dividing by the total area, and $f^{\ast}=xf_{\textrm{G}}+\left(1-x\right)f_{\textrm{h-BN}}$ is the effective free energy density per unit area given by the equilibrium bulk free energy densities of the two phases weighted by their area fractions $x$ and $1-x$. Note that we fixed $x\approx 0.5$ by fixing $\bar{n}_{i}=\pm 0.485$ and by initializing the two phases with equal or very close to equal areas. Here $\gamma$ is the average formation energy of the interface per unit length (two interfaces appear here; one has C-B and the other C-N bonds along its zigzag segment) and $\delta$ is the average formation energy of a vertex (each system has a total of four vertices; two are convex and the other two concave with respect to one of the two constituents).

The energy terms $f^{\ast}$, $\gamma$ and $\delta$ can be obtained by fitting Eq. (11) to simulation data using the method of least squares. For each interface angle considered, we varied $L_{\bot}=nL_{\bot}^{0}$ and $L_{\parallel}=mL_{\parallel}^{0}$, where $L_{\bot}^{0}$ is roughly 10 nm, $L_{\parallel}^{0}$ is the minimal $L_{\parallel}$ that satisfies periodic boundary conditions, $n=1,2,\dotsc,5;$ and $m=1,2,\dotsc,8$. For each combination of $L_{\bot},L_{\parallel}$, we initialized a corresponding bicrystalline heterostructure with segmented interfaces, relaxed it, extracted the final $f,L_{\bot}$ and $L_{\parallel}$111The energy $F$ needs to be relaxed with respect to $L_{\bot}$ and $L_{\parallel}$ to relieve possible mechanical stresses, and fitted Eq. (11) to these data. Outliers in the data and visually divergent configurations were excluded from the analysis. We determined $\gamma$ and $\delta$ for long segment lengths with large $m$. In practice, we fitted to data points unaffected by nonlinear finite-size effects where typically $n,m>1$. In addition, we considered minimal segment lengths with $m=1$ (and $L_{\parallel}$ hence a constant), i.e., interfaces with the maximal packing density of vertices. In this case the vertex energy needs to be absorbed into the interface energy as they cannot be separated without varying $L_{\parallel}$. This gives the scaling relation

[TABLE]

where $\gamma^{\ast}=\gamma+2\delta/L_{\parallel}$. This scaling relation is one-dimensional. Note that for long segment lengths the energy contribution of the vertices is negligible $\gamma\gg 2\delta/L_{\parallel}$, and hence $\gamma^{\ast}\approx\gamma$.

In the limit of long segment lengths, one can derive an analytical expression for $\gamma$ as a sum of the two segments’ individual contributions. Since the angle between the zigzag and armchair segments is 150∘, and as we know $\theta$ and $L_{\parallel}$, simple geometrical considerations give Liu et al. (2010)

[TABLE]

where $\gamma_{\textrm{ZZ}}=\gamma\left(\theta=30^{\circ}\right)$ and $\gamma_{\textrm{AC}}=\gamma\left(\theta=0^{\circ}\right)$ are the interface energies of pure zigzag and armchair interfaces, respectively.

Figure 12 shows the dimensionless interface $\gamma$ and vertex energies $\delta$ as a function of the interface angle $\theta$. All error bars are given by two-sigma confidence intervals. In the limit of long segment lengths, the numerical data for $\gamma$ agrees well with the analytical expression given by Eq. (13) predicting a maximal interface energy at $\theta\approx 10^{\circ}$. One should note that both the zigzag and armchair interfaces have a locally minimal interface energy with respect to the interface angle. In this case the zigzag interface has the lower energy. This is consistent with a number of experimental findings Gao et al. (2013); Drost et al. (2015); Liu et al. (2014b); Sutter et al. (2012); Liu et al. (2014a); Lu et al. (2014), reporting a preference towards zigzag interfaces. The dominance of zigzag interfaces can be explained by the growth process: zigzag-faceted crystals of one phase are typically formed first and serve as seeds for the subsequent growth of the second phase. The thermodynamic stability of the zigzag interfaces has also been verified computationally using density functional (DFT) theory Gao et al. (2013); Drost et al. (2015); Liu et al. (2011). A contradictory preference for armchair interfaces has also been reported by some DFT studies Bhowmick et al. (2011); Zhang et al. (2016).

Figure 12 also gives the interface energy for interfaces with minimal segment lengths ($m=1$), i.e., interfaces with the maximal packing density of vertices. Again, the zigzag direction yields the lowest energy. However, here $\gamma$ is approximately constant for $\theta\leq 15^{\circ}$ and decreases with $\theta$ for $\theta\geq 30^{\circ}$. There are two data points with significantly larger error bars. For these cases, the scaling is not as linear as for the other interface angles despite visually ideal configurations. It appears that interfaces with minimal segment lengths have lower formation energies in general.

Figure 12 also shows that the vertex energy $\delta$ is zero for both zigzag and armchair directions where there are no vertices. For the intermediate directions, $\delta$ is found slightly negative and roughly constant. This explains why the formation energy is generally lower for interfaces with more vertices. A negative $\delta$ is possible, since the vertices cannot exist independently of the segments whose $\gamma>0$. Furthermore, interfaces with long segment lengths proved stable as the highly symmetric initial states provided insufficient driving force to overcome the energy barriers for nucleating more vertices. Related point defects, triple junctions between grain boundaries, have also been shown to display negative formation energies King (1999); Srinivasan et al. (1999); Hirvonen et al. (2017).

V Summary and outlook

We have introduced an efficient and flexible phase field crystal model intended for studying heterostructures or composite materials for the general case of $N$ atomic species. This model allows well-controlled phase separation via the use of smoothed density fields. The lattice symmetries and the length scale, as well as the elastic properties, of the individual phases can be controlled readily. This model offers a straightforward approach to modelling systems with multiple ordered phases.

We have carried out a comprehensive demonstration of the model’s properties using simple 2D binary systems. More specifically, we have varied several model parameters independently to investigate their influence on phase separation, on the atomic-level structure and order of the phase interfaces and on the elastic properties of the two phases. A lattice constant mismatch between the two phases results in disordered but generally well-connected interfacial configurations. The elastic properties of the two phases can be controlled independently and robustly. We have also demonstrated the applicability of this model by considering graphene–hexagonal boron nitride lateral heterostructures (G–h-BN). We have shown that the model can reproduce many of the features of G–h-BN, such as the relative lattice constants and Young’s moduli of the two phases, as well as continuous interfaces with a preference for zigzag and armchair directions. We have also demonstrated how to model large, complex G–h-BN microstructures.

One obvious extension to this study would be to further optimize the model parameters used for G–h-BN. While we have matched the relative lattice constants and Young’s moduli approximately to their experimental values, most model parameters were either adopted from previous works or were chosen on qualitative grounds. Especially the coupling coefficients in the model could be fitted by matching the structure and formation energy of phase interfaces to corresponding results from atomistic calculations. The different chemical affinities between carbon and boron, and carbon and nitrogen, in particular, could be incorporated to the model via these parameters. In addition, structural or other more sophisticated PFC models could be incorporated to the model by replacing the terms in the energy $F$ proportional to $\beta_{i}$ and $\beta_{ij}$ with convolution kernels for more accurate material description or to allow a broader range of lattice symmetries. Although here we have focused on the 2D case for conceptual simplicity, the model is also applicable to 3D problems, where various nanoscale heterostructures and mesoscopic multiphase microstructures or composite materials could be considered. Constraining the coupling between the different lattices to their mutual interfaces could facilitate eliminating the occasional fuzzy structures without amplifying the weak oscillations. Achieving this without making the equations of motion significantly more complicated is a topic for future work.

VI Acknowledgments

This work has been supported in part by the Academy of Finland through its QFT Center of Excellence Program grant (no. 312298). We acknowledge the computational resources provided by the Aalto Science-IT project and the CSC IT Center for Science, Finland. P.H. acknowledges financial support from the Vilho, Yrjö and Kalle Väisälä Foundation of the Finnish Academy of Science and Letters. K.R.E. acknowledges financial support from the National Science Foundation under Grant No. DMR-1506634 and from the Aalto Science Institute (ASCI).

Appendix A Phase separation

In this Appendix we examine the smoothed phase separating part of the energy

[TABLE]

appearing in Eq. (1). Here $\Omega$ is the domain that can be thought of as a box with a finite volume and periodic boundary conditions (flat torus). We neglect the contribution of the disorder–crystalline boundary and assume that $\eta_{A}$ and $\eta_{B}$ take constant values specific to the phase. Let

[TABLE]

where $X$ is the component label $A$ or $B$ and $i\in\{\text{d},\text{c}\}$ corresponding to disordered and crystalline phases. We define $V(X)$ as the volume (area) of set $X$. A diagram of the setup is shown in Fig. 13. The system is set up such that $V(A_{\text{d}})=V(B_{\text{c}})$. Since $X_{\text{d}}$ and $X_{\text{c}}$ split $\Omega$ perfectly, also $V(A_{\text{c}})=V(B_{\text{d}})$. This can be done by choosing the overall number of constituents $A$ and $B$ correctly.

Now the energy contribution from Eq. (14) is

[TABLE]

Any set $Y\subset\Omega$ can be divided such that

[TABLE]

because $X_{l}$ and $X_{s}$ split $\Omega$. From this it follows that

[TABLE]

and

[TABLE]

Moreover,

[TABLE]

from which it follows that

[TABLE]

Now

[TABLE]

The components $A$ and $B$ have a similar phase diagram in the sense that either

[TABLE]

or

[TABLE]

From this it follows that $(\eta_{A}^{\text{(c)}}-\eta_{A}^{\text{(d)}})(\eta_{B}^{\text{(c)}}-\eta_{B}^{\text{(d)}})>0$, which implies that

[TABLE]

Therefore, at the ground state $V(A_{\text{c}}\cap B_{\text{c}})=0$. This shows that $F$ is minimized when different phases of the different components appear together. Even if $V(A_{\text{c}})\neq V(B_{\text{d}})$, the areas $V(A_{\text{c}}\cap B_{\text{c}})$ and $V(A_{\text{d}}\cap B_{\text{d}})$ would be minimized. If $\epsilon_{\rm AB}<0$, $V(A_{\text{c}}\cap B_{\text{c}})=0$ is maximized and the crystalline phases of the constituents overlap.

Appendix B Elastic effects due to smoothed density fields

In this Appendix we will study the non-local effects due to the introduction of the smoothed number density fields $\eta_{\rm X}$ (${\rm X}\in\{{\rm A,B}\}$ ) in $F$ (Eq. (1) ) The smoothed number densities appear in the term

[TABLE]

that might contribute to excess elastic energy if the system is deformed. Throughout this section we assume $\nu_{i}=1$. This sets the length scale of the bulk oscillations of the density fields.

The smoothed fields are defined as

[TABLE]

This convolution gives rise to non-local self-interactions. Let

[TABLE]

be the Fourier transform of $f(\mathbf{x})$. Now the inverse 2D transform is

[TABLE]

The energy $F_{\rm ct}$ can be written in terms of the Fourier transforms as

[TABLE]

by using the Plancherel theorem. Here $\hat{\eta}_{\rm A}^{*}$ is the complex conjugate of $\hat{\eta}_{\rm A}$. The fields $\hat{\eta}_{\rm X}$ ($\rm X\in\{\rm A,\rm B\}$ ) can be easily expressed using the convolution theorem as

[TABLE]

where $\hat{G}$ is the Fourier transform of the Gaussian convolution kernel, also a Gaussian

[TABLE]

Here $k=|\mathbf{k}|$ and $\gamma$ gives the length scale of the smoothing. Notice that $\hat{G}$ is real. Now

[TABLE]

We consider deformations of the form $\mathbf{k}-\mathbf{p}(\mathbf{k})$, with $|\mathbf{p}|\ll|\mathbf{k}|$. As an example, for a uniform compression of 5% , $\mathbf{p}\approx 0.05\mathbf{k}$. Now $F_{\rm ct}$ becomes

[TABLE]

Making a change of variables $\mathbf{k}\to\mathbf{k}+\mathbf{p}$ gives

[TABLE]

where $\nu$ is the change in the volume element that is given by the determinant of the Jacobian $\mathbf{I}+\nabla\mathbf{p}$.

The fields $\hat{n}_{\rm X}$ have non-zero structure at the nearest neighbor length scale (PFC fluctuations) and close to $\mathbf{k}=0$ (order–disorder boundaries). The length scale given by $1/k=1$ corresponds to nearest neighbor distance of the PFC lattice and $\gamma$ is chosen such that $\hat{G}(k=1)\ll 1$ implying that deformations at this length scale do not contribute to $F_{\rm ct}$. We will investigate the other important regime, where $\mathbf{k}$ is small. Let $F_{\rm ct}=\int\mathop{}\!\mathrm{d}\mathbf{k}\,f_{\rm ct}.$ Expanding $\hat{G}$ around $\mathbf{k}$ gives

[TABLE]

Here we have used the fact that for small $\mathbf{p}$, $\nu\approx 1+\nabla\cdot\mathbf{p}$ and assume that the part proportional to $\nabla\cdot\mathbf{p}$ is much smaller than unity. Also, the system is initially in equilibrium implying that $\hat{G}(\mathbf{k}+\mathbf{p})^{2}$ has to be expanded up to second order. Now the excess part of the energy at $\mathbf{k}$ is

[TABLE]

In order to calculate the contribution to the excess elastic energy due to an interface, we assume that $\hat{n}_{\rm A}$ and $\hat{n}_{\rm B}$ vary only in one direction and are peaked around $\mathbf{k}=0$. Furthermore, we can estimate $\hat{n}_{\rm A}^{*}\hat{n}_{\rm B}<\phi^{2}$, where $\phi$ is the amplitude of the one-mode oscillations in the crystal. More precisely

[TABLE]

The Fourier amplitudes due to the interfaces should be significantly smaller than the Fourier amplitude of the bulk oscillations. The excess energy $\Delta f_{\rm ct}$ is maximized for parallel $\mathbf{k}$ and $\mathbf{p}$. Let us assume that $\mathbf{p}(\mathbf{k})<\delta k_{x}$, where $k_{x}$ is the component of $\mathbf{k}$ parallel to the interface and $\delta$ is small. Now we can estimate the contribution of the interface per interface length as

[TABLE]

We can compare $\Delta f_{\rm ct}^{\text{int}}$ to the bulk elastic excess energy. The bulk elastic energy density due to component $A$ is

[TABLE]

where $\Omega$ is a compact domain that is taken to infinity and $V_{\Omega}$ is its area. We consider the contribution of the bulk oscillations and set

[TABLE]

where $\mathbf{q}_{j}$ are the principal reciprocal lattice vectors with $q_{j}=1$. Now

[TABLE]

We evaluate $f_{\rm el}$ in Fourier space as

[TABLE]

where

[TABLE]

is the Fourier transform of the pattern forming operator $(1+\nabla^{2})$. We repeat the earlier calculation by replacing $\hat{G}$ with $\hat{\mathcal{L}}$. We use

[TABLE]

Now

[TABLE]

Up to second order in $\mathbf{p}$

[TABLE]

Let us consider small linear deformations with $\mathbf{p}(\mathbf{k})=\delta\mathbf{J}\mathbf{k}$, with some matrix $\mathbf{J}$ with squared eigenvalues $\lambda_{1}^{2}+\lambda_{2}^{2}=\operatorname{Tr}(\mathbf{J}^{2})=1$ and some small $\delta$. For a hexagonal lattice, the reciprocal lattice vectors $\mathbf{q}_{j}$ form a star with 6-fold symmetry. It can be shown Heinonen et al. (2014) that

[TABLE]

In order to compare with $\Delta f_{\rm ct}^{\text{int}}$ of Eq. (28), $\Delta f_{\rm el}$ needs to be multiplied by the thickness of the interface, which we assume to be two lattice constants $a=4\pi/\sqrt{3}$. This gives

[TABLE]

Inserting $\beta_{\rm A}=1$, $\epsilon_{\rm AB}=1$ we get an estimate

[TABLE]

which proves that the contribution of the smoothing term to the elastic excitation energies is insignificant.