Interlayer decoupling in twisted bilayers of $\beta$-phosphorus and arsenic: a computational study

TL;DR

This computational study explores how twisting bilayers of blue phosphorus and grey arsenene creates interlayer decoupling, significantly altering their electronic and magnetic properties, including increased bandgap and induced ferromagnetism.

Contribution

It demonstrates that twisting these bilayers reduces interlayer coupling and enhances magnetic and electronic tunability, revealing new potential for 2D material engineering.

Findings

Interlayer decoupling increases with twist angle.

Bandgap enlarges by up to 50% in twisted bilayers.

Ferromagnetism is induced by hole doping in Moiré superlattices.

Abstract

We investigate magnetism and band structure engineering in Moir\'e superlattice of blue phosphorus ($\beta$-P) and grey arsenene ($\beta$-As) bilayers, using \textit{ab initio} calculations. The electronic states near the valence and conduction band edges have significant $p_z$ character in both the bilayers. Thus, twisting the layers significantly reduce the interlayer orbital overlap, leading to a decrease in the binding energy (up to $\sim33\%$) and an increase in interlayer distance (up to $\sim10\%$), compared to the most stable AA-stacking. This interlayer decoupling also results in a notable increase (up to $\sim$25-50\%) of the bandgap of twisted bilayers, with the valance band edge becoming relatively flat with van-Hove singularities in the density of states. Thus, hole doping induces a Stoner instability, leading to ferromagnetic ground state, which is more robust in Moir\'e superlattices, than that of AA-stacked $\beta$-P and $\beta$-As.

Tables

Table 1. Table 1: Structural parameters and bandgap of AA-stacked bilayer β 𝛽 \beta -P and β 𝛽 \beta -As. Our calculations are in good agreement with earlier works on β 𝛽 \beta -P Ghosh  et al. ( 2015 ) and β 𝛽 \beta -As Kecik  et al. ( 2016a ) (see main text for details).

Crystal	r (Å)	$\theta$	h (Å)	d (Å)	Eg (eV)
$\beta$-P	2.26	93.06∘	1.24	3.35	1.17
$\beta$-As	2.51	92.14∘	1.39	3.2	0.86

Equations

$$\Delta E_{DFT}(d)=\frac{E_{bilayer}(d)-2\times E_{monolayer}}{A},$$

Full text

I Introduction

The era of 2D materials started with the discovery of graphene, hosting massless Dirac quasiparticles, with gate tunable electronic and optical properties.Geim and Novoselov (2007); Castro Neto et al. (2009); Zhong et al. (2017) This motivated the discovery of several other 2D semiconductors. Among them, transition metal dichalcogenides (such as MoS2, WS2 etc.)Chhowalla et al. (2013); Xu et al. (2014); Xia et al. (2017) and black phosphoreneLiu et al. (2014); Lin et al. (2017) are some of the predominantly explored candidates. Several useful applications including transistors based on graphene,Schwierz (2010); Paek et al. (2017) MoS2,Radisavljevic et al. (2011) and black phosphoreneLi et al. (2014) have already been demonstrated.Dai and Zeng (2014) There has also been significant interest in twisted bilayers forming Moiré superlattices, with the twist angle offering another handle on the tunability of electronic and optical properties.Tong et al. (2016); Kim et al. (2017); Cao et al. (2018) In particular, several recent studies on twisted bilayer graphene explore the presence of flat bands and the van-Hove singularities (vHs) resulting in exciting many body instabilities, on account of the relative rotation between the layers.Tan et al. (2016); Lopes dos Santos et al. (2007); Li et al. (2009); Brihuega et al. (2012); Yan et al. (2012); Padhi et al. (2018)

The properties of ultrathin layered materials are governed by the stacking sequence of the monolayers, which dictates the interlayer orbital overlap among different layers. Thus, understanding the influence of the stacking sequences on the nature of the interlayer interaction among monolayers has been of great fundamental interest.Ping and Fuhrer (2012); Berashevich and Chakraborty (2011); Liu et al. (2014); Yan et al. (2015); Kecik et al. (2016a); Zhang et al. (2015); Shulenburger et al. (2015) Motivated by this, in this paper we explore the nature of the interlayer coupling and the modulation of the electronic and magnetic properties in twisted bilayers of blue phosphorus ($\beta$-P)Zhu and Tománek (2014); Guan et al. (2014) and grey arsenene ($\beta$-As),Mardanya et al. (2016); Kamal and Ezawa (2015) both having similar crystal structure and physical properties. Among several possible 2D allotropes predicted for group-V elements,Nahas et al. (2017) $\beta$-P has already been reported experimentally,Zhang et al. (2016); Gu et al. (2017) along with a detailed atomistic study of its growth mechanism.Han et al. (2017) $\beta$-allotrope of group-V elements have a graphene like buckled honeycomb structure, as opposed to the puckered honeycomb crystal of black phosphorene or $\alpha$-P. In terms of electronic transport properties, both the allotropes are found to be comparable.Priydarshi et al. (2018) Interestingly, a hole doping induced ferromagnetic ground state has been reported in case of $\beta$-P.Fu et al. (2017); Zhou et al. (2017)

Our study reveals that the $p_{z}$ orbitals contribute significantly to the valance as well as conduction band edge of the bilayers of AA-stacked $\beta$-P and $\beta$-As. Thus, the interactions among the $p_{z}$ orbitals on different layers are impacted significantly by the relative rotation of the layers relative to each other. We find that, increasing the rotation angle leads to larger interlayer distances and reduced binding energy in both $\beta$-As and $\beta$-P bilayers. As a result, the bandgap increases in Moiré superlattices, and also shows more modulation with applied electric field and strain. Additionally, the highest valance band of the Moiré superlattices also becomes relatively flat with van Hove singularities in the density of states. This leads to Stoner instability induced ferromagnetic ground state with hole doping, which is found to be more stable in twisted structures, than that of AA-stacked bilayers.

The manuscript is organized as follows: In Sec. II, we describe the details of the density functional theory (DFT) calculations. The electronic properties and the interlayer decoupling of the Moiré superlattice are discussed in detail in the Sec. III; in the context of binding energy, electronic band structure and its modulation with vertical electric field and strain, and magnetic instability in Moiré superlattices. The results are summarized in Sec. IV.

II Methodology

Structural relaxations and electronic band structure calculations are performed using density functional theory (DFT), as implemented in the Quantum ESPRESSO package.Giannozzi et al. (2009) A plane wave basis set with kinetic energy cutoff of 30 Ry and projector augmented wave pseudopotentials are used. Exchange-correlation effects are included within the framework of generalized gradient approximations (GGA), as proposed by Perdrew-Burke-Ernzerhof (PBE).Perdew et al. (1996) Dispersive forces are taken into account by using the vdW-DF-obk8 van der Waals correction.Thonhauser et al. (2015, 2007); Langreth et al. (2009) A $k$-point mesh of $24\times 24\times 1$ is used for Brillouin zone integrations in case of smaller unit cells of AA-stacked bilayers and a coarser mesh is chosen for the Moiré superlattices in twisted bilayers, according to their size. Structural optimizations are carried out until the energy difference between two successive steps of ionic relaxation are less than $10^{-4}$ Ry and all three components of the force on each atom belonging to the unit cell are less than $10^{-3}$ Ry/Bohr. A vacuum of 20 Å  is applied perpendicular to the plane to avoid spurious interactions between the periodic replicas in the out-of-plane direction. The crystal structures are prepared by XCrysden software.Kokalj (2003)

III Results and Discussions

III.1 AA-stacked bilayers: crystal and electronic band structure

The $\beta$-allotrope of monolayer P and As are known to have a buckled honeycomb structure.Ghosh et al. (2015); Kecik et al. (2016a); Mardanya et al. (2016); Nahas et al. (2017) Among several possible stacking sequences, AA-stacked bi-layers are known to have the lowest energy configuration.Ghosh et al. (2015); Kecik et al. (2016a) The top and side view of the crystal are shown in Fig. 1(a). The values of the structural parameters marked in Fig. 1(a), like the bond length (r), bond angle ($\theta$), buckling height (h) and interlayer separation (d), are reported in Table 1. These are in good agreement with the earlier studies on $\beta$-P Ghosh et al. (2015) (reported values of r = 2.26 Å, $\theta=93.06^{\circ}$, h = 1.24 Å, d = 3.23 Å) and $\beta$-AsKecik et al. (2016a) (reported values of r = 2.51 Å, $\theta=91.87^{\circ}$, h = 1.4 Å, d = 3.27 Å).

Surveying the structural data available in the literature for these bilayers, we find that the interlayer distance is sensitive to the choice of van der Waals interaction parameters, and values ranging from 3.23 to 3.40 Å  [3.14 to 3.27 Å] are reported for AA-stacked $\beta$-P [$\beta$-As].Ghosh et al. (2015); Pontes et al. (2018); Kecik et al. (2016a, b) Unlike the inter-layer distance, other structural parameters like unit cell dimension, bond length and bond angle are independent of the choice of the van der Waals interaction, and the reported values lie within a very narrow range.Ghosh et al. (2015); Pontes et al. (2018); Kecik et al. (2016a, b) In our calculations, the interlayer distance lies within the range of values reported in the literature, and the other structural parameters match very closely (within 1-2%). This validates our choice of the pseudo-potential and other calculation parameters like kinetic energy cut-off and k-point mesh used in this work.

The calculated electronic band structures of bilayer AA-stacked $\beta$-P and $\beta$-As are shown in Fig. 1(b) and Fig. 1(c), respectively. Evidently, both of them have similar qualitatively features: an indirect bandgap, with conduction band minimum (CBM) and valance band maximum (VBM) located at some point along the $\Gamma$M line. In the valence band, there is another point comparable to that of the VBM energy, lying in between the $\Gamma$K line. Relative weight of the $p_{z}$ orbitals on different energy bands clearly shows their dominance near the valence band edge. On the other hand, while $p_{z}$ orbitals have notable contributions near the conduction band edge, $p_{x}$ and $p_{y}$ orbitals are found to have equally significant contributions. There is a large anisotropy in the band dispersion, in vicinity of both VBM and CBM. This will lead to highly anisotropic charge carrier effective mass, as reported in case of monolayer allotropes.Priydarshi et al. (2018); Nahas et al. (2017) The indirect bandgap values are reported in Table 1, and they are in good agreement with the values reported in the literature.Ghosh et al. (2015); Kecik et al. (2016a) However, the bandgap in GGA based calculations are known to be underestimated, and nearly two fold increase is reported with use of more accurate hybrid functionals.Ghosh et al. (2015); Kecik et al. (2016a)

III.2 Binding energy

In case of layered materials, monolayers are believed to be stacked on top of each other and held together by Van der Waals (VDW) interactions. Such forces are isotropic in nature and bonding between two monolayers is expected to be independent of their relative orientation. This hypothesis is tested in case of twisted bilayers of $\beta$-P and $\beta$-As (forming Moiré patterns), by comparing their binding energies with that of AA-stacked bilayers. We use the optimized structures of AA-stacked bilayers to generate the Moiré superlattices of $\beta$-P and $\beta$-As, by rotating one layer with respect to the other. We select five rotation angles, namely, 10.9*∘, 13.2∘, 16.1∘, 21.8∘* and 27.8*∘*. The Moiré superlattices for different relative rotation angles, each of them having a hexagonal unit cell, are shown in Fig. 2. Binding between two monolayers is characterized by calculating the difference between the energy of two “infinitely separated” monolayers and a bilayer with inter-layer separation $d$:

[TABLE]

where $A$ is the area of the unit cell and $d$ is varied over a range from 2.5 to 9 Å, beyond which $\triangle E$ is found to be independent of $d$. In case of twisted bilayers with Moiré patterns, the interlayer distance $d$ is calculated by first measuring the distance between an atom at layer 1 and its closest neighbor at layer 2 and then taking an average of the distances measured between all such pairs. The calculated values of $\Delta E_{DFT}(d)$ for AA stacked and twisted bilayers are shown by the symbols in Fig. 3 (a) and (b) for $\beta$-P and $\beta$-As, respectively.

Further analysis is carried out by fitting the data points obtained from DFT calculations with a Morse potential, which is known to reasonably describe the non-bonded interlayer interactions, particularly at large $d$. The Morse potential is given by $\Delta E_{Morse}(d)=D\left[1-e^{-a(d-d_{0})}\right]^{2}$, where $D$ is the potential depth, $d_{0}$ is the equilibrium separation and $a$ (taken to be 1 in our case) controls the width of the potential. The potential depth or $D$, is the difference between the energies at equilibrium separation $d_{0}$ (equal to 3.35 and 3.2 Å for AA-stacked P and As) and at a large value of $d$, where there is no interaction between the monolayers. As shown in Fig. 3, this energy is almost constant for interlayer distance $\sim$9 Å, and thus we will use 9 Å as the large $d$ value to calculate $D$ for the Morse fit. As shown in Fig. 3 (a)-(b), AA-stacked $\beta$-As has larger (0.30 J/m2) binding energy than that of $\beta$-P (0.16 J/m2). Same is true for the twisted structure with Moiré patterns, with $\beta$-As (0.21 J/m2) having nearly twice the binding energy than that of $\beta$-P (0.11 J/m2). Comparing with the AA-stacked bilayers, the binding energy decreases and interlayer separation increases in case of twisted structures.

The variation of interlayer distance and binding energy as a function of stacking sequence can be attributed to the steric effect: the repulsion between the overlapping electron clouds. Steric effect depends on the size and in-plane interatomic distance of the constituent atoms. For example, larger size and interatomic distance of S atom in MoS2 compared to that of the C atom in graphene, results in a 3 times stronger steric effect in MoS2.Liu et al. (2014) Since P and As atoms are comparable in size with S atom, large steric effect is expected in case of $\beta$-P and $\beta$-As as well.

Indeed, bilayers of $\beta$-P and $\beta$-As show a wide range of interlayer distances and binding energies for different stacking sequences. The most stable ground-state stacking in the bilayers is the AA-configuration.Zhang et al. (2015); Kecik et al. (2016a) The twisted bilayers with Moiré pattern can be considered to be a “mixture” of all possible stacking sequences. Thus the binding energy of the twisted structures is expected to be lower than that of AA-stacking. Our calculations reveal nearly 33% reduction of binding energy and 10% increase of interlayer separation in 21.8∘ rotated twisted structures, compared to AA-stacked bilayers. The large dependence of the binding energy on the orientation of the constituent monolayers, clearly highlights the role of steric effects in these bilayers. Thus the interlayer bonding in the bilayers of $\beta$-P and $\beta$-As is not purely of VDW type. Similar observation has also been reported for the other layered allotrope such as black phosphorus.Shulenburger et al. (2015)

III.3 Electronic band structure of twisted bilayers

In the previous subsection, we showed that the formation of the Moiré patterns weakens the interlayer interaction. Increased interlayer distance in Moiré pattern of $\beta$-P and $\beta$-As leads to the reduction in the overlap among out of the plane $p_{z}$ orbitals. Since both valence band maximum (VBM) and conduction band minimum (CBM) have significant contributions from the $p_{z}$ orbitals [see Fig. 1(b) and (c)], decreasing overlap among such orbitals is expected to have significant effects.

Two such band structures are shown for $\theta=21.8^{\circ}$ rotated Moiré superlattice of $\beta$-P and $\beta$-As, in Fig. 4(a) and (b), respectively. Interestingly, for this particular rotation, both the bilayers are transformed to a direct bandgap semiconductor, with the valence and conduction band edges located in the middle of the $\Gamma$-M line. Note that, since the unit cells of the Moiré superlattices are also hexagonal, the high symmetry points in the reciprocal space are the same as in Fig. 1(b). The magnitude of the bandgap is plotted as a function of rotation angle $\theta$ in Fig. 4 (c) and (d), for $\beta$-P and $\beta$-As, respectively. As shown in the figures, compared to AA-stacked bilayers, bandgap in Moiré superlattices can increase up to $\sim$25% and 50% in case of P and As, respectively. While the bandgap values fall within a relatively small window, ranging from 1.4 to 1.5 eV in case of $\beta$-P Moiré superlattices, it changes prominently as a function of $\theta$ in case of $\beta$-As, which opens the possibility of bandgap tuning in the latter, by rotating one layer relative to the other. In order to analyze the connection between bandgap magnitude and interlayer separation, we also plot the latter as a function of the rotation angle in Fig. 4 (c) and (d). Evidently, the bandgap enhancement is highly correlated with the increasing interlayer separation and the decreasing overlap among the $p_{z}$ orbitals of adjacent layers.

Other than the $\theta=21.8^{\circ}$, rest of the Moiré superlattices are indirect bandgap semiconductors as shown in Fig. 4(c) and (d). However, due to the flatness of the valence band, the difference between the direct and indirect bandgap is very small in most of the cases [see Fig. S1 and S2 in the Supporting Information]. The flatness of the bands in vicinity of the band edges also gives rise to a large density of states, which is conducive to magnetic and other many-body instabilities. A recent example is that of the ‘magic angle’ twisted bilayer graphene, which shows superconducting instability.Lopes dos Santos et al. (2007); Li et al. (2009); Brihuega et al. (2012); Yan et al. (2012)

III.4 Band structure modulation with electric field

We now explore the possibility of bandgap engineering in these twisted bilayers by means of a vertical electric field. In case of $\beta$-P and $\beta$-As, very large electric field (about 0.5 V/Å) is known to reduce the bandgap, ultimately leading to an insulator to metal transition around 0.7 to 1 V/Å.Ghosh et al. (2015); Mardanya et al. (2016) Moreover, an interesting topological transformation is also reported for monolayer $\beta$-As, at nearly 1 V/Å  electric field.Mardanya et al. (2016)

A similar trend is observed in case of Moiré superlattices, where bandgap magnitudes decrease sharply beyond 0.3 [0.4 V/Å] for bilayer $\beta$-P [$\beta$-As]. A comparison of bandgap modification (with increasing electric field strength) in AA-stacked bilayers and Moiré super-lattices is shown in Fig. 5 (a). Interestingly, the latter (originally a direct bandgap semiconductor) is found to be converted to an indirect bandgap semiconductor at high electric electric field, both in case of $\beta$-P and $\beta$-As. This transition happens because the conduction band edge is shifted to the $\Gamma$ point (see Fig. S3 and S4 in SI). The computationally predicted electric field values of 0.3-0.4 V/Å  can further be reduced if we account for the presence of ripples and other structural imperfections.Agnihotri et al. (2018)

III.5 Band structure modulation with strain

Monolayers of $\beta$-P and $\beta$-As also show significant bandgap modulation under the influence of external strain.Zhu and Tománek (2014); Mardanya et al. (2016); Kecik et al. (2016a) A similar trend is observed in case of AA-stacked bilayers and Moiré superlattices of $\beta$-P and $\beta$-As. As shown in Fig. 5(b), the bandgap decreases monotonically under tensile strain. On the other hand, under compression, bandgap initially increases, followed by a decrease at higher value of strain. Evolution of electronic band structures as a function of strain for AA-stacked $\beta$-P and $\beta$-As are shown in SI (see Fig. S5 and S6). We find that, not only the magnitude of bandgap decreases under strain, band edges also shift from their original location. For example, in case of AA-stacked $\beta$-P, CBM and VBM shifts to the K-point and $\Gamma$-point at 2% and 4% compressive strain, respectively (see Fig. S5 in SI). On the other hand, in case of AA-stacked $\beta$-As, VBM moves to the $\Gamma$-point under 3% compressive strain, while CBM shifts to the $\Gamma$-point at 4% tensile strain (see Fig. S6 in SI). However, we do not find any indirect to direct bandgap conversion within $\pm$5% strain in case of AA-stacked $\beta$-P and $\beta$-As.

Interestingly, for Moiré superlattices, the bilayers can also convert from a direct to an indirect bandgap semiconductor under strain. This is more prominent in case of $\beta$-As Moiré superlattice ($\theta=21.8^{\circ}$), which changes to an indirect bandgap semiconductor under relatively small tensile (2%), as well as compressive (1%) strain. On the other hand, $\beta$-P Moiré superlattice ($\theta=21.8^{\circ}$) undergoes such a transition only under compression (3%). A detailed analysis reveals that the band edges shift to different valleys as a function of strain, leading to a direct to indirect bandgap transformation. For example, in case of $\beta$-P Moiré superlattice, the valence band edge moves to the $\Gamma$ point and the conduction band edge moves to the K-point at 3% compression (see Fig. S7 in SI). On the other hand, in case of $\beta$-As Moiré superlattice, the valence band edge moves to the $\Gamma$ point at 1% compressive strain, while a 2% tensile strain triggers a similar shift to the conduction band edge. (see Fig. S8 in SI).

III.6 Magnetic instability

As shown Fig. 1(b)-(c) and Fig. 4(a)-(b), the valence band of bilayer $\beta$-P and $\beta$-As have a relatively flat shape near the VBM. This results in a van Hove singularity in the electronic density of states near the valence band edge, which can be explored via hole doping. The van Hove singularity is conducive to quantum many body instabilities such as magnetism. In fact, magnetism arising from the van Hove singularity induced Stoner instability in monolayers and bilayers of P and As has already been predicted over a over a broad range of hole doping.Fu et al. (2017); Zhou et al. (2017)

Here we explore the existence of a ferromagnetic instability in hole doped bilayers of twisted $\beta$-P and $\beta$-As. The calculated magnetic moment (per hole) in the Ferromagnetic ground state is shown in Fig. 6. Our calculations reveal that the magnetic moment is equal to the number of holes in the system up to certain value of doping, beyond which both the bilayers fall back to the non-magnetic state rapidly. The critical value of hole concentration at which the crossover takes place is nearly double in AA-stacked $\beta$-P, than that of $\beta$-As. This is in good agreement with the reported data in the literature.Fu et al. (2017); Zhou et al. (2017) Note that, prediction of a ferromagnetic ground state up to $10^{14}$ cm*-2* hole concentration is very encouraging, because it is within the range of current experimental capabilities, as demonstrated for gate-tunable MoS2 device.Ye et al. (2012)

As shown in the respective insets for Fig. 6, the spin density is predominantly localized in the $p_{z}$ orbitals of the P and As atoms. Thus, the increase in the interaction of the $p_{z}$ orbitals with decreasing interlayer separation in bilayers is likely to result in the weakening of the magnetic ground state. This is confirmed by calculating the magnetic moment by changing the interlayer separation $d$ (slightly above and below the equilibrium value) for AA-stacked bilayers in Fig. 6. As expected we find that the critical value of hole doping, up to which the AA-stacked bilayers remain in ferromagnetic ground state, gradually decreases as the monolayers are brought closer to each other. This is also reflected in the fact that the ferromagnetic ground state survives up to higher hole concentration in case of Moiré superlattices (see Fig. 6 for results for $\theta=21.8^{\circ}$), which has a larger interlayer distance and smaller interactions among the $p_{z}$ orbitals. In fact, the ferromagnetic ground state can exist up to nearly 50% higher hole doping in Moiré superlattices as compared to the AA-stacked bilayers, making the former more suitable for magnetic and spintronic applications.

IV Conclusion

In conclusion, we have investigated the evolution of interlayer coupling in twisted $\beta$-P and $\beta$-As bilayers. Our calculations reveal nearly 33% reduction in binding energy, leading to 10% increase of interlayer separation in Moiré superlattices of $\beta$-P and $\beta$-As, compared to AA-stacked bilayers. We show that the band edges (particularly the VBM) are dominated by the out of plane $p_{z}$ orbitals in both $\beta$-P and $\beta$-As. Thus, increasing the interlayer separation reduces the overlap among $p_{z}$ orbitals from adjacent planes. This interlayer decoupling in Moiré superlattices leads to significant changes in the electronic, as well as magnetic (with hole doping) properties of $\beta$-P and $\beta$-As.

We have shown that compared to AA-stacked bilayers, bandgap increases by 25-50% in Moiré superlattices. Furthermore, on account of a relatively flat valence band, the difference between the direct and indirect bandgap is very small. In some cases, twisted bilayers are found to be direct bandgap semiconductors, while the AA-stacked bilayers are indirect bandgap semiconductors. Similar to other 2D materials, bandgap can be tuned by applying strain, as well as an electric field normal to the plane of the twisted bilayers.

Due to the Stoner instability originating from the van Hove singularity of density of states near the valence band edge, a ferromagnetic ground state is observed in bilayer $\beta$-P and $\beta$-As. The critical value of hole concentration, up to which the ferromagnetic ground state is sustained, increases by $\sim$50% in Moiré superlattices, compared to AA-stacked bilayers. Thus, twisted bilayers are better candidates than AA-stacked $\beta$-P and $\beta$-As for the purpose of magnetic and spintronic applications.

V Acknowledgments

The authors acknowledge funding from the Ramanujan fellowship research grant, the DST Nanomission project and SERB (EMR/2017/004970). The authors also thank the computer center of IIT Kanpur for providing HPC facility.