Simulations of twisted bilayer orthorhombic black phosphorus

We identified, by means of coincidence site lattice theory, an evaluative stacking phase with a wavelike Moiré pattern, denoted as $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$, from all potentially twisted bilayer orthorhombic black phosphorus. Such a twisted stacking comes with a low formation energy of $-162.8\mathrm{meV}$, very close to existing AB stacking, according to first-principles calculations. Particularly, classic molecular dynamic simulations verified that the stacking can be directly obtained in an in situ cleavage. The stability of $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$ stacking can be directly attributed to the corrugated configuration of black phosphorus leading to the van der Waals constraining forces, where the top layer can get stuck to the bottom when one layer rotates in plane relative to the other by $\sim 70.5^{\circ }$. Tribological analysis further revealed that the interlayer friction of $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$ stacking reaches up to 1.3 nN, playing a key role in the origin of $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$.

Figure 1

Lattice structures of monolayer black phosphorus (or phosphorene) in top and side views. Orthorhombic basis vectors ${\mathit{a}}_1$ and ${\mathit{a}}_2$ are along the zigzag and armchair directions in the unit cell, respectively, which can generate two supercell lattice vectors, ${\mathit{b}}_1$ and ${\mathit{b}}_2$.

Figure 2

Forming process of a commensurate twisted bilayer phosphorenes (BP). $R$ indicates a 2D rotation operator, and $T$ a rigid translation operator. ${\mathit{b}}_1$ and ${\mathit{b}}_2$ are two superlattice basis vectors in the twisted BP produced. $I$ is an inversion operator commuting with the coordinates $\mathit{X}$ mapped by the superlattice basis vectors.

Figure 3

Superlattice structure of $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$ stacking. Optimized superlattice structure in (a) top, (b) front, and (c) side views, together with unit-cell structural parameters. The interlayer constraining forces of (d) $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$ stacking and (e) AB stacking are denoted by red and black wavy arrows, whose in-plane components are presented in (a) and (f), respectively.

Figure 4

Computational model with top and front views for a piece of cleavage phosphorene drifting on a black phosphorus substrate surface. The average radii of the pie and plate are 1.73 and 6.32 nm, respectively. The small pie is initially set to be 0.35 nm above the big plate, and the intersection angle between their armchair directions is defined as a twist angle $\phi$.

Figure 5

Thermodynamic features of drifting-induced $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$/AB settlings at one tenth of room temperature (RT). Drifting trajectories of the mass center for the small pie including six typical moments marked for (a) $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$ stacking with an initial direction of ${60}^{\circ }$ and (c) AB stacking with an initial direction of ${45}^{\circ }$. Dynamic twist angle as a function of time and six typical Moiré pattern snapshots for (b) $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$ stacking and (d) AB stacking, the times of which are marked in the angle-time curve. The interlayer friction acting on the pie changes with the odometer for (e) AB stacking and (f) $2\mathrm{O}\text{−}\mathrm{t}\alpha \mathrm{P}$ stacking.