Kolmogorov-Crespi Potential For Multilayer Transition Metal Dichalcogenides: Capturing Structural Transformations In Moir\'e Superlattices

TL;DR

This paper develops a Kolmogorov-Crespi potential for multilayer transition metal dichalcogenides, enabling efficient structural relaxation of moiré superlattices with results comparable to DFT, facilitating studies of their electronic properties.

Contribution

It introduces a new interlayer potential for TMDs and a mixing rule for heterostructures, allowing accurate and efficient modeling of moiré superlattices.

Findings

KC potential accurately reproduces DFT relaxations

Relaxed structures yield similar electronic band structures to DFT

Enables efficient study of structural and electronic features in TMD moiré systems

Abstract

We develop parameters for the interlayer Kolmogorov-Crespi (KC) potential to study structural features of four transition metal dichalcogenides (TMDs): MoS$_2$, WS$_2$, MoSe$_2$ and WSe$_2$. We also propose a mixing rule to extend the parameters to their heterostructures. Moir\'e superlattices of twisted bilayer TMDs have been recently shown to host shear solitons, topological point defects and ultraflatbands close to the valence band edge. Performing structural relaxations at the DFT level is a major bottleneck in the study of these systems. We show that the parametrized KC potential can be used to obtain atomic relaxations in good agreement with DFT relaxations. Furthermore, the moir\'e superlattices relaxed using DFT and the proposed forcefield yield very similar electronic band structures.

Figures

Tables

Table 1. Table 1: Parameters for interlayer Kolmogorov-Crespi interactions.

	$z_0$	C0	C2	C4	C	$\delta$	$\lambda$	A
	(Å)	(meV)	(meV)	(meV)	(meV)	(Å)	(Å-1)	(meV)
S2-S3	3.857	7.074	2.624	0.024	24.859	0.982	2.762	55.637
Se2-Se3	3.992	7.886	6.355	-0.011	26.762	0.935	2.753	60.829
Mo1-S3, Mo2-S2	4.049	-2.718	-2.719	-0.325	16.459	0.616	2.564	14.239
Mo1-Se3, Mo2-Se2	5.080	-1.120	-0.761	-0.117	4.624	0.777	1.996	6.209
W1-S3, W2-S2	4.331	-1.344	-1.370	-0.162	8.372	0.587	2.082	7.078
W1-Se3, W2-Se2	5.079	-1.109	-0.762	-0.120	4.627	0.783	1.999	5.112
Mo1-Mo2	9.236	1.471	1.126	0.131	1.928	1.183	0.924	3.439
W1-W2	10.545	0.254	0.139	-0.028	0.506	1.010	0.944	1.356

Table 2. Table 2: Comparison of the shear (SM) and layer breathing modes (LBM) of bilayer TMDs computed using SW+KC with experimental measurements.

	LBM (cm-1)
	SW+KC	Expt.
BLMoS2	38.6	40 63, 41.6 64
BLWS2	30.5	33.8 64
BLMoSe2	31.0	34.3 64, 29 65
BLWSe2	27.2	29.1 64, 27 65
	SM (cm-1)
	SW+KC${}_{M-M}{}^{}{}_{}{}^{X-X}$	Expt.
BLMoS2	20.9	22 63, 24.2 64
BLWS2	18.0	19.6 64
BLMoSe2	18.3	21 64, 18 65
BLWSe2	15.7	17.7 64, 17 65

Equations

$$V_{ij}=e^{-\lambda (r_{ij}-z_0)}{V}_{\rho }-A{\left(\frac{r_{ij}}{z_0}\right)}^{-6}$$

$$V_{\rho }=[C+f({\rho }_{ij})+f({\rho }_{ji})],$$

$$f(\rho )=e^{(-\rho /\delta {)}^2}\sum _{n=0}^2{C}_{2n}(\rho /\delta {)}^{2n}$$

Full text

Indian Institute of Science] Center for Condensed Matter Theory, Department of Physics, Indian Institute of Science, Bangalore 560012, India

1 Introduction

The growing family of layered materials and the ’Lego set’ 1 of possible heterostructures is an attractive field of research 2, 3. Layered materials are composed of two-dimensional (2D) atomic layers weakly held together by van der Waals (vdW) forces. The widely used approach to theoretically simulate 2D materials and their heterostructures is through vdW corrected density functional theory 4, 5, 6, 7 (DFT). This quantum mechanics based approach is accurate but computationally intractable for applications that require large scale simulations of 2D materials with more than 10,000 atoms. One such application is the study of moiré patterns in 2D materials.

On introducing a small-angle twist between the two layers of a bilayer system leads to the formation of a large scale moiré superlattice (MSL) 8, 9, 10. MSLs have interesting properties at the electronic as well as the structural level. The moiré pattern is composed of various local high-symmetry stackings. 9, 8, 10 Structural reconstructions lead to the formation of shear strain solitons at stacking boundaries and topological point defects in twisted bilayer graphene (tBLG) 11, 12, 8, 13, 14. Unconventional superconductivity was recently observed in tBLG at a ’magic’ 15 twist angle due to the formation of ultraflatbands close to the Fermi level 16, 17, 18, 19. Recent DFT calculations show that small-angle twisted bilayer TMDs can also host ultraflatbands, shear solitons and topological point defects 9. However, angles smaller than 3.5∘ could not be explored in this study 9 due to higher computational cost. Obtaining the relaxed coordinates of the atoms in the MSL is a major bottle-neck in these calculations.

Classical force-field based methods have been used to replace expensive DFT calculations to study the structural properties of 2D materials 20, 21, 14, 8, 22. An important ingredient to classical force-field based methods is the proper modelling of the vdW interaction between layers. 23, 24, 25 The Lennard-Jones (LJ) model has been widely used to account for the interlayer vdW interactions. Sliding one layer of a 2D material with respect to the other leads to different stacking configuration of the atoms, which have different binding energies. This stacking dependence of the binding energy is not captured by the LJ model. Tribology 26, 27 and the study of moiré patterns in 2D materials are particularly sensitive to the stacking dependence 28, 29.

To overcome this drawback, the Kolmogorov-Crespi29, 30 (KC) potential has been developed for bilayer graphene. This potential includes an additional stacking dependent term coupled with the short-range interaction. It has been successfully applied to graphene 8, 31, 32 and BN 33 to explore solitons in moiré patterns 8, 14, dislocations, shear modes, self-retraction 34, 35, 36 and lubricity of graphene flakes 37, 38, 39, 27 and carbon nanotubes 22, 40, 41, 29, 27, etc. However, it has been restricted to these materials and not been extended to the vast family of transition metal dichalcogenides, and their heterostructures 42, 43. We note that other methods like the continuum approach 44, 45 and a Gaussian-LJ 46 method have been developed to treat shear in bilayer MoS2.

In this work, we parametrize the KC potential for multilayers of four transition metal dichalcogenides: MoS2, MoSe2, WS2 and WSe2. The Stillinger-Weber (SW) 47, 48 force field is used to model the constituent single layers and the KC potential determines the interlayer interactions. The KC parametrization is performed to fit the binding energy of the various bilayer stackings computed using DFT. We also propose a mixing rule for these parameters to simulate heterostructures using these materials. The layer breathing mode and shear mode frequencies computed using these potentials are in good agreement with experiments. We use the SW and fitted KC parameters to perform the structural relaxation of the 3.5∘ and 56.5∘ twisted bilayer MoS2. We find that the atomic relaxations using these forcefield calculations are in good agreement with relaxations performed using DFT. Furthermore, the electronic band structure computed using the forcefield relaxed structure is in good agreement with that computed using the DFT relaxed structure. Indicating that this method can replace the computationally expensive DFT relaxation steps while studying the electronic properties of moiré superlattices.

2 Computational Details

The vdW corrected DFT calculations are carried out using the plane-wave pseudopotential package Quantum Espresso 49. We use ultrasoft 50 pseudopotentials and the local density approximation 51 (LDA) for the exchange-correlation functional. The van der Waals interactions are computed using the van der Waals density functional 52, 53 along with Cooper 54 exchange (vdW-DF-C09) 55. We use a plane-wave energy cut-off of 50 Ry for the wavefunctions and a cut-off of 500 Ry for the charge density. The Brillouin zone is sampled with a 12$\times$12$\times$1 k-point mesh. The DFT calculations are performed with the optimized geometry of the cell for each TMD. The forcefield calculations are performed using the LAMMPS 56, 57 package. The SW 48 potential is used for intralayer interactions and KC potential for interlayer interactions. The minimizations are performed using conjugate gradient method (force tolerance of 10*-6* eV/$\mathrm{\AA}$) as implemented in LAMMPS with the DFT lattice constant. Parametrization of the KC potential to fit the DFT binding energy is performed with help of the Dakota 58, 59 package. Commensurate moiré superlattices of twisted bilayer MoS2 are constructed using the Twister code 60. KC potential file and LAMMPS input files for the bilayer TMDs are available with the Twister package 60.

3 Kolmogorov-Crespi potential

Fig. 1 shows the various possible high-symmetry stackings of a bilayer TMD system obtained by sliding one layer with respect to the other 9. The AA and AB stackings are obtained when atoms in the top layer are exactly above atoms in the bottom layer. In the case of AA stacking (Fig. 1 (a)), similar atoms are on top of one another (M on M and X on X), while in AB stacking (Fig. 1 (d)) dissimilar atoms are on top of one another (M on X and X on M). The other stackings are Bernal type (Fig. 1 (b), (c), (e) and (f)), ie. one atom in the top layer is in the hexagonal cavity of the bottom layer. In our notation, BM/X implies a Bernal (B) stacking and that atom M in the top layer is directly above atom X in the bottom layer (Fig. 1 (c)). On starting with the AA stacking, we obtain BX/M and BM/X on sliding one layer with respect to the other in the $x$ direction as shown in Fig. 1. We denote this shear by $\mathrm{SH}_{\mathrm{AA}}$. Note that BX/M and BM/X are equivalent stackings. On sliding, starting with AB stacked layers, two non-equivalent high-symmetry stackings are obtained: BM/M and BX/X (Fig. 1). We denote this shear by $\mathrm{SH}_{\mathrm{AB}}$. There are thus five unique high-symmetry stackings possible with TMDs in the 2H-phase.

We aim to model the binding energy (BE) landscape of the bilayer system as a function of shear between the two layers using classical forcefields. The BE landscape computed using vdW corrected DFT for bilayer MoS2 with fixed interlayer spacing of 6.35 $\mathrm{\AA}$ is shown in Fig. 2 (a) and (d). The Kolmogorov-Crespi is a reliable interlayer potential since it contains an explicit registry dependent term. The form of the potential is given by:

[TABLE]

The potential, $V_{ij}$, is defined between atom $i$ in one layer and atom $j$ in the adjacent layer. The potential is set to zero after a cut-off radius, $r_{cut}$. The KC potential includes a stacking dependent term, $V_{\rho}$, multiplying the short-range repulsive interaction.

[TABLE]

where, $\rho_{ij}^{2}=r_{ij}^{2}-(\mathbf{n_{i}}\cdot\mathbf{r_{ij}})^{2}$, $\rho_{ji}^{2}=r_{ij}^{2}-(\mathbf{n_{j}}\cdot\mathbf{r_{ij}})^{2}$ and

[TABLE]

$V_{\rho}$ determines the energy barrier to shear one layer with respect to the other in a bilayer system. The KC potential thus consists of 8 parameters for each type of interaction. A taper function 24, 23, 61, 21 is often used to ensure the interlayer potential goes to zero smoothly. Using a large $r_{cut}$ = 14 $\mathrm{\SIUnitSymbolAngstrom}$ does not affect the results significantly in the absence of the taper function. See Supplementary Information (SI) for more details. $\mathbf{n_{i}}$ and $\mathbf{n_{j}}$ are surface normals at site $i$ and $j$, respectively. To compute the surface normal at the position of an atom $i$, we find the neighbours of the same type in a radius $r_{n}$ around atom $i$. This is illustrated in Fig. 4. $r_{n}$ is chosen to accommodate the six nearest neighbours of the same type. Six normals are then constructed from consecutive pairs of neighbouring atoms and averaged to obtain $\mathbf{n_{i}}$. See SI for more details. We refer to interactions that take into account the normals as KC-n. A simplification to the present form of the potential can be introduced by setting $\mathbf{n_{i}}=\mathbf{n_{j}}=\hat{z}$. This approximation only works for multilayers whose bending leads to normals that do not deviate significantly from $\hat{z}$. For the case of twisted bilayers, we find this to be a good approximation since the normals are close to $\hat{z}$. We will refer to this approximation as KC-z. This is discussed further in subsequent sections.

Steric effects govern the relative energy of the stackings. The stackings with chalcogen atoms in the top layer directly above chalcogen atoms in the bottom layer, ie. AA and BX/X, are unfavourable and highest in energy 9, 62. The other three stackings are relatively lower in energy. It would thus seem reasonable to consider the KC interactions between the X2 and X3 (Fig. 3 (a)) atoms alone. On fitting the KC parameters of this interaction, KCX-X, to DFT BE leads to a good representation of shear starting with the AA stacking ($\mathrm{SH}_{\mathrm{AA}}$), as shown in Fig. 3 (b). But fails to completely represent the BE of shear starting with the AB stacking (Fig. 3 (c)). In particular, the difference in BE between AB and BM/M stackings is not captured. This is expected, because in the presence of only X2-X3 interactions, AB stacking is indistinguishable from the BM/M stacking.

To differentiate the BE of the AB and BM/M stackings in $\mathrm{SH}_{\mathrm{AB}}$ we further introduce the M1-X3, M2-X2 and M1-M2 interlayer interactions (see Fig. 3 (a)). We introduce these additional interactions and refit all the KC parameters. In the supplementary information we propose two other approaches to obtain the BE landscape with reasonable accuracy. The first includes only X2-X3 and M1-M2 interactions, KC${}^{\mathrm{X-X}}_{\mathrm{M-M}}$, and the second approach includes only X2-X3, M1-X3 and M2-X2 interactions, KC${}^{\mathrm{X-X}}_{\mathrm{M-X}}$. We have also fit KC${}^{\mathrm{X-X}}_{\mathrm{M-M}}$ parameters including the taper function discussed above. KC${}^{\mathrm{X-X}}_{\mathrm{M-M}}$ and KC${}^{\mathrm{X-X}}_{\mathrm{M-X}}$ parameters are provided and further discussed in the supplementary information. For the rest of the article, we will show the performance of the parametrized KC potential when all interlayer interactions, X2-X3, M1-X3, M2-X2 and M1-M2, are included. These parameters are provided in Table 1. The BE landscape for MoS2 with SW+KC is shown in Fig. 2 (b) and (e). Fig. 2 (c) and (f) show the deviation of the landscape computed using SW+KC from DFT.

The KC potential parameters are obtained by fitting the BE as a function of interlayer spacing (out-of-plane separation between M atoms of top and bottom layer) for the high-symmetry stackings and as a function of shear between the two layers, $\mathrm{SH}_{\mathrm{AA}}$ and $\mathrm{SH}_{\mathrm{AB}}$. The fitting to the shear is performed at the average interlayer spacing of the high-symmetry stackings. Fitting shear at the average interlayer spacing is sufficient to reproduce the shear at the minimum and maximum interlayer spacings as well. Fig. 5 (a), (c) and (e) compare the BE computed using SW+KC with the corresponding DFT BE for bilayer MoS2. Fig. 5 (b), (d) and (f) compares SW+KC BE with DFT for bilayer WSe2. Fig. 6 similarly shows the performance of the KC parameters for bilayer WS2 and MoSe2. Furthermore, the shear and layer breathing modes 66 of the bilayer TMDs computed using SW+KC are in good agreement with experimental measurements (Table 2). See supplementary information for the performance of the other set of parameters: KC${}^{\mathrm{X-X}}_{\mathrm{M-M}}$ and KC${}^{\mathrm{X-X}}_{\mathrm{M-X}}$.

4 Heterostructures: mixing rule

Mechanical exfoliation 67 and chemical vapour deposition 68, 69, 70 techniques can be used to construct vertical stacks 71 of different TMDs. These heterostructures generally form a Type II heterojunction 67, which makes them suitable for applications in nano- and opto-electronics71. The KC parameters developed above for the bilayer TMDs can be extended to heterostructures by the use of a simple mixing rule. The cross-interaction parameters for S-Se are obtained by taking the arithmetic mean (AM) of S-S and Se-Se parameters. Similarly the Mo-W interaction are computed as an AM of Mo-Mo and W-W parameters. We tried other mixing rules like the geometric mixing rule and find that the AM works best to reproduce the DFT BE.

We apply the mixing rule to simulate the MoS2/MoSe2 and MoS2/WS2 heterostructure. MoS2 and MoSe2 are not lattice matched. For the purpose of comparison, we strain the layers to their average lattice constants in the DFT and forcefield calculations. We compute the BE as a function of interlayer spacing and shear as shown in Fig. 7 for MoS2/MoSe2 and MoS2/WS2 heterostructure. AA stacking in the context of heterostructures implies transition metal (chalcogen) atom in the top layer is above a transition metal (chalcogen) atom in the bottom layer. AB stacking implies transition metal (chalcogen) atom in the top layer is above a chalcogen (transition metal) atom in the bottom layer. The BE computed using SW+KC is in good agreement with DFT (Fig. 7).

5 Moiré superlattices

Twisted bilayer TMDs form distinct moiré patterns for twist angles close to 0∘ and 60∘ 9, 72, 73. Fig. 8 shows the MSL for 5.1∘ (M5.1) and 54.9∘ (M54.9) twists. The MSL is composed of the various high-symmetry stackings (Fig. 1)9. The regions O, A and B are AA, BX/M and BM/X respectively. The regions O’, A’ and B’ are AB, BM/M and BX/X respectively. Relaxing these moiré patterns in DFT starting from the rigidly twisted structure leads to significant in-plane and out-of-plane displacements of the atoms 9. Because of large bending rigidity 66, the out-of-plane displacements smoothly vary across the MSL surface. The ILS is largest for the AA and BS/S stacking regions as shown in Fig. 9 (e) for the case of 3.5∘ twisted bilayer MoS2. We also use the SW+KC to perform the relaxations of the same MSL, starting with the rigidly twisted structure. The ILS distribution obtained using the forcefield approach (Fig. 9 (f)) is in good agreement with the DFT results. We have performed these relaxation using KC-n, ie. including normals in the interlayer interaction.

The undulations in each of the layers is smooth. The out-of-plane displacements of each layer, while significant, vary over a large area. We thus find that the normals do not deviate significantly from $\hat{z}$ for small or large twist-angles. The distribution of $\phi$ in KC-n relaxed structures is shown in Fig. 11 for a relatively large (7.3∘) and small (1∘) twist angle MSL. For angles greater than 7.3∘ and smaller than 52.7∘, the out-of-plane displacements are small. $\phi$ takes a maximum value of 1.2∘ in Fig. 11, and cos(1.2∘) $\approx$ 1. Furthermore, we find the energy computed using KC-n and KC-z for these angles differ by less than 0.01 meV/atom. KC-z is thus sufficient to capture relaxations in moiré superlattices. We also compare the forces in the MSL computed using DFT, KC-n and KC-z for unrelaxed 7.3∘ twisted bilayer MoS2 in Fig. 12. The out-of-plane forces are in excellent agreement between DFT and the proposed forcefields.

The in-plane displacements of the MoS2 units in the moiré pattern lead to an increase in the area of low-energy stacking with respect to the rigidly twisted MSL 9. The stackings in the relaxed structure can be identified by means of an order-parameter (OP). The OP is defined 9 as the shortest displacement vector that takes any given stacking in the moiré pattern to the highest energy stacking. We define $\vec{u}$ for the moiré pattern formed by twist angles close to 0∘ and $\vec{v}$ for those formed by twist angles close to 60∘. Hence $\vec{u}$ is the shortest displacement vector that takes any given stacking in M3.5 to AA stacking and $\vec{v}$ is the shortest displacement vector that takes any given stacking in M56.5 to BS/S 9. Fig. 9 (a) shows the OP distribution for the rigidly twisted MSL, M3.5. The distribution of the OP in the DFT relaxed MSL (Fig. 9 (b)) is in good agreement with that in the SW+KC${}^{\mathrm{S-S}}_{\mathrm{Mo-S}}$ relaxed MSL (Fig. 9 (c)). This indicates that the registry of atoms are similar in the MSL relaxed within DFT and the forcefield. OP and ILS using the other set of parameters, KC${}^{\mathrm{S-S}}_{\mathrm{Mo-Mo}}$ and KC${}^{\mathrm{S-S}}_{\mathrm{Mo-S}}$, are also in good agreement with DFT (see supporting information).

The electronic structure of MSLs in bilayers 9 as well as heterostructures 74 has been shown to host flat bands close to the valence band edge. We demonstrate that the electronic structure calculations can replace the time-consuming relaxation steps in DFT with the forcefield relaxations. To this effect, we compute the band structure of 7.3∘, 9.4∘, 50.6∘, 52.7∘ twisted MSL using the DFT relaxed structure and the structure from SW+KC relaxations. We find that these band structures are in good agreement as shown in Fig. 13.

In conclusion, we parametrize the KC potential for multilayer TMDs and propose a mixing rule to extend the parameter to TMD heterostructures. The parametrization is performed to fit the BE landscape computed using van der Waals corrected DFT. We show that the forcefield can be used to simulate structural transformations in MSLs of twisted bilayer MoS2. These structural reconstructions are in good agreement with those computed using DFT. Furthermore, the electronic band structure computed using the forcefield relaxed structure is in good agreement with band structure computed using the DFT relaxed structure. This indicates that the computationally expensive DFT relaxation steps can be replaced by the proposed forcefield method. Apart from predicting the structure of MSLs, several tribological properties can be studied using these KC potentials. The lubricity of the various bilayers and bilayer heterostructures can be computed. Moreover, the dependence of lubricity on the twist angle between the bilayers can be explored. The effect of finite temperature on the tribological properties and on solitons in the MSLs can also be studied with molecular dynamics simulations using these potentials.

{acknowledgement}

We thank the Supercomputer Education and Research Centre (SERC) at IISc for providing the computational facilities.

{suppinfo}

Supporting Information. Verification of computation of normals with peridically modulated MoS2 sheet; alternate parameter sets KC${}^{\mathrm{X-X}}_{\mathrm{M-M}}$ and KC${}^{\mathrm{X-X}}_{\mathrm{M-X}}$; performance of the alternate parameter sets; table containing number of atoms in the MSL and MSL dimension; ILS, OP and electronic band structure of the MSL using SW+KC${}^{\mathrm{S-S}}_{\mathrm{Mo-Mo}}$ for twisted bilayer MoS2; refitted parameters with the taper function for KC${}^{\mathrm{X-X}}_{\mathrm{M-M}}$.