Phase-field crystal model for heterostructures

Atomically thin two-dimensional heterostructures are a promising, novel class of materials with ground-breaking properties. The possibility of choosing many constituent components and their proportions allows optimization of 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 modeling 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 modeling 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.

Figure 1

Influence of the coupling parameter ${\epsilon }_{12}$ on the heterostructures. (a) The initial state with $n_1^{(}\mathrm{c})$ and $n_2^{(}\mathrm{c})$ in partial overlap. (b) The relaxed heterostructure for the repulsive case where ${\epsilon }_{12}=0.2$. Here, the crystalline phases 1 and 2 are shown in cyan and red, respectively. (c) The relaxed coexistence between a mixed and a disordered phase for the attractive case where ${\epsilon }_{12}=-0.2$. Here, the mixed phase appears white due to the coincidental alignment of the structures in $n_1^{(}\mathrm{c})$ and $n_2^{(}\mathrm{c})$, and the disordered phase appears black. Note also a slight excess of disordered phase due to not having chosen the average densities to correspond to exactly equal parts of either phase. (d) Profiles of the smoothed densities ${\eta }_1$ (cyan) and ${\eta }_2$ (red) along the periodic edge in the horizontal direction of the initial state from (a). (e) Profiles of the densities (solid lines) $n_1$ (cyan) and $n_2$ (red) and of the smoothed densities (dashed lines) ${\eta }_1$ (cyan) and ${\eta }_2$ (red) along the periodic edge in the horizontal direction of the relaxed heterostructure from (b). (f) Same profiles for (c).

Figure 2

Influence of ${\epsilon }_{12}$ on polycrystalline heterostructures. (a) A close-up of a larger system for the repulsive case where ${\epsilon }_{12}=1.0$ after 7500 time units. The left-hand side of the panel reveals the distribution of the two phases and the right-hand side represents the heterostructure by $m=n_1+n_2$ for a clearer illustration of the atomic-level structure. (b) A close-up of a larger system for the attractive case where ${\epsilon }_{12}=-1.0$ after 7500 time units. Moire overlap patterns between the two lattices are clearly visible.

Figure 3

Influence of ${\alpha }_{12}$ and the corresponding coupling on a bicrystalline heterostructure. (a) A close-up of a relaxed heterostructure with ${\alpha }_{12}=-0.03$. (b) Profiles of the densities (solid lines) $n_1$ (cyan) and $n_2$ (red) and of the smoothed densities (dashed lines) ${\eta }_1$ (cyan) and ${\eta }_2$ (red) along the periodic horizontal edge of the relaxed heterostructure from (a).

Figure 4

Influence of ${\alpha }_{12}$ on a polycrystalline heterostructure. (a) A close-up of a larger system after a relaxation of 25 000 time units. The left-hand side of the panel reveals the distribution of the two phases and the right-hand side demonstrates $n_2$ with weak oscillations in $n_2^{(}\mathrm{d})$. (b) The heterostructure from (a) represented by $m=n_1+n_2$ for a clearer illustration of the atomic-level structure.

Figure 5

Two-dimensional Young's modulus $Y_{\text{2D}}$ as a function of the gradient term coefficient ${\beta }_2$ for both crystalline phases. The markers give the actual data and the lines are optimal linear fits. The slope for the second fit is 0.17. The analytical prediction obtained for the crystalline phase 2 using the analytical expressions for the stiffness coefficients $C_{11}$ and $C_{12}$ [Eqs. (8)–(10)] is plotted using solid black markers.

Figure 6

Collage of interface structures for different lattice mismatches, strains and misorientations. The left-hand side of each panel reveals the distribution of the two crystalline phases, and the right hand side represents the heterostructure by $m=n_1+n_2$ for a clearer illustration of the atomic level structure. Note that in many cases only a small part of the total length of the interface modeled is shown. The first column of panels gives reference (R) structures with no lattice mismatch between the two phases. The next two columns give strained (S) structures where the mismatch is accommodated by elastic deformation. The last two columns give unstrained (U) structures where the mismatch is accommodated by misfit dislocations. The mismatch for each column is indicated via ${\lambda }_2=1/{\nu }_2$. The first row of panels gives structures with zero-misorientation armchair-armchair (AC-AC) interfaces. The second row depicts low-misorientation tilt interfaces with $2\theta \approx 51.4^{\circ }$. The third row demonstrates high-misorientation armchair-zigzag (AC-ZZ) interfaces. In panels where feasible, we indicate one dislocation by $\perp$ and adjacent lines in the armchair direction.

Figure 7

Lattice-mismatched polycrystalline heterostructures. The panels offer close-up view of larger systems after a relaxation of 250 000 time units. The top halves of the panels show the distribution of the two phases and the bottom halves represent the heterostructures by $m=n_1+n_2$ for a clearer illustration of the atomic-level structure. We have fixed ${\lambda }_1=1/{\nu }_1=1$ and have varied ${\lambda }_2=1.05$ (a), 1.1 (b), and 1.2 (c).

Figure 8

Zigzag interface from a bicrystalline G–h-BN lateral heterostructure. (a) A visualization of the heterostructure in the top half and in the bottom half the same structure represented by $m=n_1+n_2+n_3$ for a clearer illustration of the atomic level structure. In the top half, graphene appears cyan, whereas boron and nitrogen are in magenta and yellow. In the bottom half, graphene (h-BN) appears darker (brighter). (b) Profiles of the densities (solid lines) $n_1$ (cyan), $n_2$ (magenta), and $n_3$ (yellow) and of the smoothed densities (dashed lines) ${\eta }_1$ (cyan), ${\eta }_2$ (magenta), and ${\eta }_3$ (yellow) along the periodic edge in the horizontal direction of the relaxed heterostructure.

Figure 9

Large random polycrystalline graphene–h-BN lateral heterostructure relaxed from white noise for $2.5\times {10}^6$ time units. The sides of the system are approx. 50 nm in length. (a) A coarse-grained representation of the large-scale structure where graphene appears cyan and h-BN red. [(b)–(e)] Close-ups of the region indicated by the blue square in (a). The width (height) of the region shown in the close-ups is approximately 9 nm. (b) A visualization explained in Fig. 8, (c) the total density $m=n_1+n_2+n_3$, (d) $n_1$ and (e) $n_2$ in the close-up.

Figure 10

Large random polycrystalline graphene–h-BN lateral heterostructure from a random Voronoi grain structure. A side of the system is approximately 50 nm long. (a) A coarse-grained representation of the large-scale structure where graphene appears cyan and h-BN red. [(b)–(e)] Close-ups of the atomic level structure of the regions indicated by the blue squares in (a). The regions shown in the close-ups are 6 nm wide. (b) A collection of graphene–h-BN interfaces and h-BN grain and inversion boundaries. (c) A large-misorientation interface. Small-misorientation (d) zigzag and (e) armchair interfaces.

Figure 11

Examples of strained bicrystalline heterostructures used to study the formation energy of G–h-BN interfaces with perfect honeycomb order. The interface angle is $\theta \approx 16.1^{\circ }$ for both structures and is indicated by the white wedge in (a). In (a), $m=2$ and in (b) $m=4$. For both, $n=1$ and the total width is $L_{\perp }\approx 12$ nm (see text for the definition of $m$ and $n$). The zigzag and armchair segments of the right-hand side interfaces are traced in red and blue, respectively. Their lengths $L_{ZZ}$ and $L_{AC}$ are also indicated. In (a), $L_{ZZ}$ is 4 lattice constants and $L_{AC}$ is $2\sqrt{3}$ lattice constants. In (b), both are twice as long. Note that the system shown in panel (a) can be decomposed into two identical fields by cutting the domain in half perpendicular to the interface. Here we consider such subdomains with four vertices.

Figure 12

Dimensionless interface $\gamma$ (in blue on left axis) and vertex $\delta$ (in red on right axis) energies as a function of the interface angle $\theta$; see Eq. (11). Open markers correspond to $\gamma$ in the limit of long segment lengths (large $m$) and solid markers to $\gamma$ in the limit of minimal segment lengths ($m=1$); see Eq. (12). The solid curve gives the analytical expression for $\gamma$ from Eq. (13).

Figure 13

A pictorial showing an example of disordered and ordered regions of the two components A and B. As a consequence of $V\left(A_{\text{d}}\right)=V\left(B_{\text{c}}\right)$ [and $V\left(A_{\text{c}}\right)=V\left(B_{\text{d}}\right)$], the areas $V(A_{\text{c}}\cap B_{\text{c}})$ and $V(A_{\text{d}}\cap B_{\text{d}})$ are equal. See text for details.