Grain Boundary Structures and Collective Dynamics of Inversion Domains in Binary Two-Dimensional Materials

TL;DR

This paper develops a phase field crystal model to study grain boundary structures and dynamics in binary two-dimensional materials, revealing new dislocation core structures and defect-mediated growth mechanisms related to inversion symmetry breaking.

Contribution

It introduces a novel phase field crystal model tailored for binary 2D materials, uncovering new dislocation structures and defect-driven growth processes not previously documented.

Findings

Discovery of new dislocation core structures at grain boundaries.

Identification of defect-mediated growth dynamics of inversion domains.

Linking inversion symmetry breaking to atomic migration at grain boundaries.

Abstract

Understanding and controlling the properties and dynamics of topological defects is a lasting challenge in the study of two-dimensional materials, and is crucial to achieve high-quality films required for technological applications. Here grain boundary structures, energies, and dynamics of binary two-dimensional materials are investigated through the development of a phase field crystal model that is parameterized to match the ordering, symmetry, energy and length scales of hexagonal boron nitride. Our studies reveal some new dislocation core structures for various symmetrically and asymmetrically tilted grain boundaries, in addition to those obtained in previous experiments and first-principles calculations. We also identify a defect-mediated growth dynamics for inversion domains governed by the collective atomic migration and defect core transformation at grain boundaries and junctions, a process that is related to inversion symmetry breaking in binary lattice.

Figures

Equations

$$\partial n_A/\partial t={\nabla }^2\frac{\delta \mathcal{F}}{\delta n_A},\partial n_B/\partial t=m_B{\nabla }^2\frac{\delta \mathcal{F}}{\delta n_B},$$

$$\mathcal{F}$$

$$\gamma =\frac{b{Y}_2}{8\pi }\theta \left[1+\ln (b/r_0)-\ln (2\pi \theta )\right],$$

Full text

Topological defects, such as dislocations and grain boundaries (GBs), are known to be pivotal in controlling material properties. It is challenging to effectively capture the complexity of defects, given the nonequilibrium nature of material growth and evolution processes. Recent progress in the study of two-dimensional (2D) hexagonal materials, such as graphene, hexagonal boron nitride (h-BN), and transition metal dichalcogenides (TMDs), provides an excellent platform for the investigation of defect properties and dynamics. This is driven by the demand for controllable fabrication and synthesis of large-scale, high-quality samples of these atomically thin systems, which mostly rely on vapor-phase heteroepitaxy techniques particularly chemical vapor deposition. Such large-area 2D epitaxial films are usually polycrystalline Yazyev and Chen (2014), with various types of defects found in both theoretical Yazyev and Louie (2010); Liu and Yakobson (2010); Liu et al. (2012); Zou et al. (2013) and experimental Huang et al. (2011); Gibb et al. (2013); Cretu et al. (2014); Li et al. (2015); van der Zande et al. (2013); Björkman et al. (2013) studies of 2D materials. Typical examples include penta-hepta ($5|7$) defects in graphene Huang et al. (2011) and either penta-hepta or square-octagon ($4|8$) defects in h-BN Gibb et al. (2013); Cretu et al. (2014); Li et al. (2015) and TMD van der Zande et al. (2013) sheets.

Compared to 2D single-component materials such as graphene, in binary hexagonal materials (e.g., h-BN and TMDs) the inversion symmetry is broken in the corresponding binary honeycomb lattice. A much richer variety of GB configurations can be identified, some of which can significantly alter system electronic properties, as predicted by first-principles calculations Liu et al. (2012); Zou et al. (2013) and found in experiments of h-BN Li et al. (2015), MoS2 van der Zande et al. (2013); Zhou et al. (2013) and MoSe2 Liu et al. (2014) epitaxial monolayers. Of particular interest are the $60^{\circ}$ GBs [i.e., inversion domain boundaries (IDBs)], a characteristic of inversion symmetry breaking. Depending on the detailed structures of dislocation cores, these $60^{\circ}$ boundaries can cause a reduction of band gap (as in h-BN Liu et al. (2012)), the appearance of mid-gap states (for GBs consisting of $4|8$ cores in MoS2 Zou et al. (2013); van der Zande et al. (2013)), or a transition from semiconducting to localized metallic modes (for $4|4$ or $8|8$ cores in MoS2 Zou et al. (2013); Zhou et al. (2013) and MoSe2 Liu et al. (2014)).

It is of great difficulty to effectively track or control the dynamics of defect formation over the relevant spatial and temporal scales, via either in situ experimental techniques or simulations. Experimentally the studies of defect dynamics mostly rely on the activation process of electron irradiation that generates migrating vacancies in the sample Cretu et al. (2014); Lin et al. (2015); Lehtinen et al. (2015). Most theoretical studies are based on atomistic methods particularly first-principles density functional theory (DFT) and molecular dynamics. While large progress has been made for identifying lowest-energy defect structures and their electronic properties Yazyev and Chen (2014); Yazyev and Louie (2010); Liu and Yakobson (2010); Liu et al. (2012); Zou et al. (2013), the atomistic techniques are usually limited by the restrictions of small length and time scales and the pre-constructed defect core configurations. It is thus important to develop and apply modeling methods that are able to access large system sizes and realistic time scales while still maintaining microscopic spatial resolution.

Here we first construct such a model for binary 2D materials, based on the phase field crystal (PFC) method Elder et al. (2002); *Elder04; Elder et al. (2007); Emmerich et al. (2012); Mkhonta et al. (2013) which can simultaneously address mesoscopic structural profiles at diffusive time scales and resolve crystalline microstructures Huang and Elder (2008); Mellenthin et al. (2008); Tóth et al. (2011); Schwalbach et al. (2013); Adland et al. (2013); Ofori-Opoku et al. (2013); Backofen et al. (2014); Hirvonen et al. (2016); Mkhonta et al. (2016); Huang (2016). For simplicity, here a 2D planar model is developed that does not include out-of-plane deformations, given the constraint of 2D monolayers during epitaxial growth. As described in Supplemental Material SM , the model parameters are chosen to match the symmetry, sublattice ordering, Young’s modulus, and atomic spacing of h-BN. This model is used to systematically study GB structures, energies and the spontaneous formation of defects in systems up to micron size, without any pre-determined setup of defect cores. This allows us to identify complex defect structures that are absent in previous research, and also to predict a growth mechanism of collective dynamics and boundary defect shape transformation for inversion domains.

Our model is motivated from the classical dynamic DFT, similar to the previous procedure for PFC models Elder et al. (2007); Huang et al. (2010). For a binary AB system, the dimensionless fields $n_{A}$ and $n_{B}$ of atomic number density variation are governed by the relaxational dynamics, i.e.,

[TABLE]

where $m_{B}=M_{B}/M_{A}$ with $M_{A(B)}$ the mobility of A(B) component, and the free energy functional $\mathcal{F}$ is given by

[TABLE]

Here the rescaled parameters are chosen as $\epsilon_{A}=\epsilon_{B}=0.3$, $q_{A}=q_{B}=q_{AB}=1$, $\alpha_{AB}=0.5$, $\beta_{AB}=0.02$, $g_{A}=g_{B}=0.5$, $w=u=0.3$, and $\beta_{B}=v=m_{B}=1$. The corresponding stability and phase diagrams are given in Fig. S2 of Supplemental Material SM . We parameterize the model within the binary honeycomb regime of the phase diagram and match to the energy and length scales of h-BN. Given the Young’s modulus $Y=810$ GPa Kudin et al. (2001), lattice constant $a_{0}=2.51$ Å, and vertical layer spacing 3.33 Å for h-BN, a length unit ($l=1$) of our model corresponds to 0.342 Å and a PFC energy unit corresponds to 2.74 eV. Two pivotal factors have been incorporated in our model construction, including the lattice symmetry and a factor specific for binary compounds, that is, the heteroelemental A-B neighboring should be energetically favored as compared to homoelemental A-A or B-B ones.

In our PFC study four types of GBs are examined (see Fig. 1), for two adjoining grains of orientations $\theta_{1}$ and $\theta_{2}$ with GB misorientation angle $\theta=\theta_{1}+\theta_{2}$. The first two types correspond to symmetrically tilted GBs with $\theta_{1}=\theta_{2}=\theta/2$, starting either from the armchair edge (Fig. 1a, sym-ac, keeping mirror symmetry between the two grains), or from the zigzag direction (Fig. 1b, sym-zz). Given that polycrystalline samples are characterized by a prevalence of asymmetrically tilted GBs (i.e., $\theta_{1}\neq\theta_{2}$), we also simulate the configurations in the limit of completely asymmetric GBs with $\theta_{1}=0$ and $\theta_{2}=\theta$, for grain rotation from both armchair (Fig. 1c, asym-ac) and zigzag (Fig. 1d, asym-zz) directions.

The system simulated contains two parallel GBs to satisfy periodic boundary conditions. The total system size ranges from 14.1 nm $\times$ 24.4 nm for large misorientations to 0.9 $\mu$m $\times$ 0.3 $\mu$m for very small $\theta$ (corresponding to $1.3\times 10^{4}$ to $9.1\times 10^{6}$ equivalent atomic sites). For each angle the two misoriented grains are connected initially via a narrow band (around 2-6 grid points) of supersaturated homogeneous phase which spontaneously solidifies, causing the grains to merge and form a GB that evolves to a steady state. We tested various connection conditions by varying the relative lattice translation of adjoining grains to obtain the lowest energy states presented in Fig. 2.

Our simulations produce well-stitched GBs consisting of various types of defect rings. At very small $\theta$ the GB dislocation cores are mostly composed of well-separated $5|7$, $4|6$, or $4|10$ pairs, with some examples ($5|7$ and $4|6$) shown in Fig. 2b. GB structures of connected $5|7$ arrays are found at larger misorientations (Fig. 2c), containing either a single type of homoelemental neighboring, or the alternation of both types of opposite polarity (i.e., alternating A-A and B-B) as observed in experiments of h-BN Gibb et al. (2013); Li et al. (2015). Some new lowest-energy defect structures are also obtained. These include $4|10$ pairs comprising heteroelemental A-B neighboring (at $\theta=10.08^{\circ}$ of sym-zz), $6|8$ or $4|6|8$ arrays (sym-zz $\theta=38.82^{\circ}$), and an interesting case of GB between armchair and zigzag edges (asym $\theta=30^{\circ}$) showing as a connected $5|7$ array mixed by $6|8$ “fly-heads”. This “fly-head” configuration is similar to that obtained in graphene Liu and Yakobson (2010) as a result of the release of armchair-zigzag mismatch stress, although here a new $7|6|8$ fly-head structure appears instead of $7|5|7$ in graphene, given the energetically favorable heteroelemental-only neighboring in the $6|8$ rings compared to the unfavorable homoelemental one in $5|7$.

We have computed GB energies per unit length, $\gamma$, for all four types of GBs setup. For each setup the lowest-energy results at $\theta$ from 0 to $60^{\circ}$ are given in Fig. 2a (noting the full range of $[0,120^{\circ}]$ is symmetric around $\theta=60^{\circ}$). All the $\gamma$ values calculated are within the range of previous first-principles DFT results Liu et al. (2012). At small $\theta$ the $\gamma$ data are well fitted to the Read-Shockley equation Read and Shockley (1950)

[TABLE]

where $Y_{2}$ is the 2D Young’s modulus (set as 271 N/m for h-BN Kudin et al. (2001)), $b$ is the magnitude of the Burgers vector (assumed as $b=a_{0}=2.51$ Å for h-BN, for the shortest Burgers vector), and $r_{0}$ is the dislocation core radius. The small-angle $\gamma$ values for sym-zz and asym-zz cases are very close due to similar GB structures, giving $r_{0}=0.54\pm 0.01$ Å from the fitting. For the sym-ac GBs, the fitting yields $r_{0}=0.90\pm 0.01$ Å (in comparison, previous first-principles calculations Yazyev and Louie (2010) gave $r_{0}=1.2$ Å for graphene armchair GBs, while no prior results for h-BN are available). In contrast, the asymmetric tilting from armchair direction (asym-ac) leads to more sinuous configurations with non-regular spacing of dislocation cores when compared to the sym-ac GBs (see Fig. 2b). This explains the higher values and larger variation of small-angle $\gamma$ for asym-ac GBs. A different value of $r_{0}=0.60\pm 0.06$ Å is also obtained from fitting.

For $\theta$ approaching $60^{\circ}$, $4|8$ dislocation cores appear more frequently in low-energy structures. This square-octagon configuration, consisting only of heteroelemental neighboring, is found to be the lowest energy state of $60^{\circ}$ IDB, agreeing with the first-principles result for h-BN Liu et al. (2012). We also identify 6 other IDB structures that possess higher energy (Fig. 3), including arrays of $8|8$, $4|4$, and tilted 6-membered rings (Z6-I, equivalent to the absence of one A-atom column in $4|4$ structure) for zigzag GBs, arrays of tilted 8-membered (AT8) or tilted 12-membered (AT12) rings for armchair GBs, and a high-energy state composed of compressed 6-6 pairs (Z6-II). The stable AT12 structure is obtained at lower densities $n_{A0}=n_{B0}=-0.3$, but it is unstable (transforming to $4|8$) at $n_{A0}=n_{B0}=-0.28$ used in all other GBs calculations. Among them the $4|4$, $8|8$, and Z6-I structures, although not being found in previous h-BN studies, have been obtained in some theoretical Zou et al. (2013) and experimental Enyashin et al. (2013); Zhou et al. (2013); Liu et al. (2014) work of MoS2 and MoSe2 monolayers.

To identify the reliability and parameter dependence of our results, we have conducted a large number of additional simulations particularly for IDBs, with varying model parameters. All of them gave similar results for the defect core structures and GB energy sequence. Details are described in Supplemental Material.

To further examine the emergence and dynamics of $60^{\circ}$ boundaries, we conduct a series of simulations of grain nucleation and growth from a supersaturated homogeneous phase. The dynamics is modeled via $\partial n_{A}/\partial t=-{\delta\mathcal{F}}/{\delta n_{A}}+\mu_{A}$ and $\partial n_{B}/\partial t=-m_{B}({\delta\mathcal{F}}/{\delta n_{B}}-\mu_{B})$, where $\mu_{A(B)}$ is the chemical potential of A(B) component. We set $\mu_{A(B)}$ slightly above the equilibrium value to enable and control the grain growth via a constant flux. Also $\mu_{A}=\mu_{B}$ is assumed, so that a single grain is of hexagonal shape with both A- and B-terminated zigzag edges, consistent with some previous experimental Tay et al. (2014) and theoretical Zhao et al. (2015); Zhang et al. (2016) results for h-BN.

Figure 4a gives a typical simulation setup, starting from four nuclei that are of two different sizes and $60^{\circ}$ misorientation. The nuclei evolve to grains of hexagon shape and grow individually until they merge and form $60^{\circ}$ IDBs. This grain coalescence leads to the embedding of polygon-shaped inversion domains within the large grain or matrix and the subsequent domain shrinking (Fig. 4b). The resulting IDBs can be composed not only of $4|8$ cores (the lowest-energy type), but also of all other types of defect structure given in Fig. 3. Many of them are transients and transform to lower-energy configurations (mostly $4|8$, with some $8|8$ or $4|4$) during the evolution, as seen in the supplemental movies SM .

Figures 4d–f show the shrinking process of a typical triangle-shaped inversion domain. The domain boundary lines consist of connected $4|8$ pairs, joined at three junctions via decagon heart-shaped defect, with all defect cores containing only heteroelemental neighboring. The 3 heart-shaped junctions will finally merge and annihilate. The atomic evolution around a domain corner is detailed in Fig. 4g–i, over a cycle associated with one atomic-step shrinkage. Each cycle starts with the simultaneous displacement of atoms along two opposite directions (indicated by the arrows in Fig. 4g) inside each defect ring labeled from 1 to 8. This leads to the shape transformation of defect rings at a later time (Fig. 4h), e.g., $4\rightarrow 8$ transformation for ring 1, 2, 5, 7 and $8\rightarrow 4$ for ring 3, 6, 8, mediated by the rearrangement of corner decagon (white-dashed heart defect) to two 6-membered rings. Further atomic displacements inside rings 1-8 lead to further ring shape transformation (Fig. 4i), e.g., $8\rightarrow 4$ for ring 2, 5, 7 and $4\rightarrow 8$ for ring 4, 6, which are shifted downward in comparison to their starting locations in Fig. 4g. Also the heart-shaped junction is reconstructed at the location of ring 1 which, compared to the junction at the beginning of this cycle (4g, dashed heart), is moved one-step inward of the triangle (of both 1 atomic-row downward and 1 atomic-column leftward). The whole triangle domain then shrinks one-step rigidly.

Although similar heart-shaped defect (for h-BN Cretu et al. (2014)) or triangle inversion domains (for MoSe2 Lin et al. (2015); Lehtinen et al. (2015)) have been observed in experiments, they were generated by electron irradiation and thus of qualitatively different structures and dynamics. The $4|8$ and heart defects inside a single domain of irradiated h-BN flakes contained energetically unfavorable homoelemental bonds Cretu et al. (2014), instead of heteroelemental-only ones obtained here. In exfoliated or as-grown MoSe2 sheets, irradiation-induced Se vacancies gave rise to small triangle domains with $4|4$ defects at both boundaries and junctions, and the expansion, but not shrinking, of the domain is driven by the new vacancy creation and the triggered atomic motion Lin et al. (2015); Lehtinen et al. (2015).

In contrast, here polygon-shape domains are formed through coalescing of $60^{\circ}$ misoriented grains during growth. The subsequent grain dynamics is characterized by the collective atomic displacements and simultaneous shape transformation of defect core rings along the GBs and junctions, and involves the glide of boundary dislocations but not climb. Also no net translation, rotation or shear-induced deformation has been observed, a scenario different from the Cahn-Taylor shear driven mechanism with the coupling between normal and tangential motions Cahn and Taylor (2004); *CahnActaMater06; McReynolds et al. (2016). Such a mechanism originates from the breaking and maintaining of lattice plane continuity across a GB, and would not apply to the case of inversion domains for which the lattice planes always remain continuous and a $60^{\circ}$ GB forms purely due to inversion symmetry breaking in the binary lattice but not lattice sites mismatch.

During grain growth the characteristic grain size was historically expected to follow a power law $L\propto t^{\alpha}$ (corresponding to $t^{2\alpha}$ for 2D grain area or atomic site number $N$), where the exponent $\alpha=1/2$ for the classical curvature-driven growth Mullins (1956); *re:mullins89 but is found to be $<1/2$ and dependent on system temperature, noise strength, and bulk dissipation Adland et al. (2013); Backofen et al. (2014) and also GBs roughness that affects grain growth stagnation Holm and Foiles (2010). To quantify the collective dynamics identified here, we simulate large systems of size 56.2 nm $\times$ 97.4 nm ($2\times 10^{5}$ atomic sites), with a single triangle inversion domain being initialized in the center and consisting of $4|8$ IDBs and heart-defect junctions. As shown in Fig. 4c, two linear regimes of grain area shrinkage can be identified, for which the crossover to the faster shrinking rate occurs at a late time stage due to the defect interaction between GBs. Although the exponent $\alpha=1/2$ is obtained here, the corresponding dynamics of inversion domain is expected to be different from the conventional curvature-driven mechanism, given the straight or weakly curved domain boundaries and also the rigid and diffusionless boundary motion described above. The domain shrinking can be significantly slowed down for some other types of defects (e.g., a mixture of $4|8$, $4|4$, and $8|8$ cores), leading to the existence of metastable polygon domains in the sample. However, this still yields two linear regimes of decreasing $N$ (i.e., $\alpha=1/2$) before the stagnation of inversion domain, as shown in Fig. S6 of Supplemental Material. More work is needed for further understanding this $\alpha=1/2$ behavior, which is important in identifying the growth mechanisms of binary 2D materials that are different from those of traditional single-component or alloying systems.

In summary, we have systematically studied the structural and dynamical properties of 2D material systems at large spatial and temporal scales that are not accessible to traditional atomistic methods. Our results provide predictions and understanding of plasticity in binary 2D materials, including some new defect structures of symmetrically and asymmetrically tilted GBs and GB energies across the full range of misorientations. We also predict a dynamic behavior of grain motion for inversion domains, where the domain dynamics is governed by the collective atomic displacements of the connected dislocation core pairs along the boundary lines mediated via the heart-shaped defects at the junctions. These findings as well as the modeling approach presented here provide further insights for understanding the complex mechanisms of grain growth in binary 2D materials.