Insight into the stacking effect on shifted patterns of bilayer phosphorene: a comprehensive first-principles study

It is crucial to deeply understand how the interlayer interaction acts on controlling the structural and electronic properties of shifted patterns of bilayer phosphorene. A comprehensive first-principles study on the bilayer phosphorene through relative translation along different directions has revealed that there is a direct correlation between the potential energy surface and the interlayer equilibrium distance. The shorter the interlayer distance, the lower the potential energy surface. The shifted patterns with the most stable state, the metastable state, and the transition state (with energy barrier of ∼1.3 meV/atom) were found associated with the AB, the Aδ, and the TS stacking configurations, respectively. The high energy barriers, on the other hand, are ∼9.3 meV/atom at the AA stacking configuration along the zigzag pathway, ∼5.3 meV/atom at the AB′ stacking configuration along the armchair pathway, and ∼11.2 meV/atom at the AA′ stacking configuration along the diagonal pathway, respectively. The character of electronic bandgap with respect to the shifting shows an anisotropic behavior (with the value of 0.69–1.22 eV). A transition from the indirect to the direct bandgap occurs under the shifting, implying a tunable bandgap by stacking engineering. Furthermore, the orbital hybridization at the interfacial region induces a redistribution of the net charge (∼0.002–0.011 e) associated with the relative shifting between layers, leading to a strong polarization with stripe-like electron depletion near the lone pairs and accumulation in the middle of the interfacial region. It is expected that such interesting findings will provide a fundamental reference to deeply understand and analyze the complex local structural and electronic properties of twisted bilayer phosphorene and will make the shifted patterns of bilayer phosphorene promising for nanoelectronics as versatile shiftronics materials.


Introduction
The discovery of the unconventional superconductivity and Mott insulation in twisted graphene bilayer at a magic angle [1,2] has opened a new door for exploring moiré patterns induced by relative twisting or shifting between two-dimensional (2D) layered nanomaterials [3], which are of great importance in discovering promising applications to modify the spectrum of 2D electron systems and, ultimately, to create tailored electronic, optical, and magnetic properties in a van der Waals (vdW) system.Various unique properties, such as the commensurate-incommensurate phase transition through the formation of a network of domain walls in bilayer graphene by stretching one layer relative to another layer or by rotating between graphene and hexagonal boron nitride layers [4][5][6][7][8], the self-retraction of graphene bilayers [9,10], the Hofstadter's butterfly and the fractal quantum Hall effect [11,12], as well as the intrinsic quantum anomalous Hall effect in twisted bilayer graphene aligned to hexagonal boron nitride [13], the tunable exciton funnel in twisted bilayer graphene/boron nitride with regard to the certain stacking environments [14], a nearly dispersion less bands with strong orbital anisotropy from destructive interference in twisted bilayer MoS 2 [15], and a localization feature of interlayer moiré excitons in twisted MoS 2 and WS 2 heterostructures [16], etc, have been reported from the experimental observations and theoretical calculations.Apparently, the vdW interaction between 2D layers is responsible for those unique phenomena [3].
Mainly due to the direct bandgap nature of black phosphorene [17] coupled with its unique anisotropic properties in the electron conductivity [18], the thermal conductance [19], and the thermoelectricity [20], it has attracted tremendous attention in material science community to explore the unique physical properties of moiré patterns of bilayer black phosphorene and to understand how their anisotropic behaviors could be modulated by shifting or twisting one layer relative to the other through commensurate or incommensurate angles.One of the interesting findings from constructing moiré patterns of bilayer black phosphorene is that any twisted pattern consists of several specific domains with certain stacking symmetries.For instance, a moiré pattern with either 3.8°or 5.4°twist angle was identified by the combination of four local stacking domains, namely, the AA, AB, AA′, and AB′ stacking configurations [21][22][23][24][25][26].Furthermore, a systematic theoretical study for small twist angles found that there are three different twist angle regimes where the conduction and valence states take different geometries [24].The first regime is called Hubbard regime (θ < 2°), where electrons (holes) localize at AB′ (AB) stacking sites to form mesoscale rectangular lattices (obeys Hubbard Model Physics).The second regime is called Tomonaga-Luttinger regime (2°< θ 10°), where electrons and holes delocalize along the long axis of the unit cell to form arrays of quasi 1D states.The third regime is called Ballistic regime (θ > 10°), where electrons and holes fully delocalize in 2D, exhibiting anisotropic dispersions efficiently modulated by the interlayer twist angle.Thus, various stacking configuration domains will coexist as the building blocks of the moiré patterns of twisted bilayer phosphorene and should play a crucial role in determining their local structural and electronic properties.
To deeply understand the whole aspect of the stacking effect on the moiré patterns of bilayer phosphorene, we conducted a comprehensive study of the stacking effect on the shifted bilayer phosphorene focusing on the correlations between the stacking configurations and the vdW interlayer interaction to the anisotropic structural and the electronic properties of phosphorene bilayers.Specifically, we concentrated on the feature of the potential energy surface (PES), which reflects the strength of the interlayer interaction and determines the commensurate-incommensurate phase transition, the self-retraction, and the operation of nanoelectromechanical systems under a relative motion between two layers [27].Understanding the PES of bilayer phosphorene is crucial for exploring its anisotropic structural and dynamical behavior.Recent molecular dynamic study reported [28] that the typical PESs are strongly correlated with the friction properties of the bilayer interface.The potential corrugation and hence the frictional force typically increase when the interlayer spacing is reduced, consistent with their moiré pattern [29][30][31].In this work, we have revealed that there is a direct correlation between the PES and the interlayer equilibrium distance under shifting.The shorter the interlayer distance, the lower the PES.Furthermore, a polarization pattern of net charge distribution was found in the interfacial region, and the bandgap value as well as the direct/indirect bandgap nature can be tuned by shifting one layer with respect to the other layer.Such findings will provide a fundamental reference to deeply understand and analyze the complex local structural and electronic properties of twisted bilayer phosphorene.

Computational details
We employed the density functional theory (DFT) [32,33], as implemented in the Vienna ab initio simulation package (VASP) [34], and performed the structure optimization, the dynamic stability analysis, and the PES and electronic properties calculations.The electron-ion interactions were described by the projector augmented wave (PAW) [35], while electron exchange-correlation interactions were treated by the generalized gradient approximation (GGA) [36] in the scheme of Perdew-Burke-Ernzerhof (PBE) [37].The vdW interactions were considered by employing the zero damping DFT-D3 method [38][39][40].The structural relaxation was carried out using Congregate-Gradient algorithm [41].The periodic boundary conditions were chosen in the layered plane with a vacuum space of at least 15 Å between adjacent bilayers to avoid any mirror interactions.An energy cutoff was set to be 400 eV for the plane wave basis in all calculations, and the criteria for the convergences of energy and force in relaxation processes were set to be 10 -5 eV and 10 −4 eV Å −1 , respectively.A 1 × 1 rectangular primitive cell with lattice constants a and b was chosen to study the moiré patterns of bilayer black phosphorene by shifting, and the Brillouin zones were sampled by 15 × 15 × 2 k-point meshes generated in accordance with the Monkhorst-Pack scheme [42] in the optimization and band structure calculations.
The dynamic stability of the systems was examined through the phonon dispersion calculations based on Parlinski-Li-Kawazoe method [43,44].The temperature dependence of free energy, entropy, and heat capacity are calculated by based on the phonon dispersion as follows: The band structures were calculated at the DFT_PBE level and corrected by a hybrid functional to the DFT_HSE06 level [45].The Bader analysis scheme [46][47][48][49][50] was employed to analyze charge redistribution and charge transfer between layers, and the electronic properties at the interface region were analyzed from the differential electron charge density (DCD), defined as Δρ = ρ total -ρ top -ρ bottom , where ρ total is the total electron charge density of the bilayer system, ρ top and ρ bottom are the electron charge densities associated with the top and bottom layers in the bilayer system, respectively.The plane-averaged differential charge density (Δρ) along zdirection normal to the phosphorene bilayer is calculated by the integral of x y z x y , , d d , and the amount of net charge electrons transferred far from the bottom of the bilayer up to z level is analyzed by The bilayer black phosphorene with AA stacking configuration (space group Pmna) was first carried out under a fully relaxed process with the degrees of freedom on the atomic positions, the lattice symmetry, and the volume of the unit cell until the pressure on the supercells is decreased to the value less than 0.17 kB.The optimized lattice constants for a and b are 3.298 Å and 4.515 Å, respectively, and the optimized equilibrium interlayer distance is 3.497 Å, consistent with other theoretical calculations [21,22,[51][52][53].In order to systematically study how the interlayer interaction acts on the structural and electronic properties of shifted patterns of bilayer phosphorene through relative translation, we built up a series bilayer structures, starting from the AA stacking configuration and shifting the top layer with respect to the bottom layer on a 2D grid of 15 × 15 points with the steps of 0.11 Å and 0.15 Å in the zigzag and armchair directions, respectively (see figure S1).Clearly, the high symmetric moiré patterns with AB, Aδ, AB′, TS, and AA′ configurations (as defined in the insets of figure 1) are automatically obtained when shifting the top layer along the zigzag direction by a/2 (i.e. the AB stacking configuration with space group Pbcm) or along the armchair direction by b/ 4 (i.e. the Aδ stacking configuration with space group p2/m) and by b/2 (i.e. the AB′ stacking configuration with space group Pccm), or along the diagonal direction by (a + b)/5 ( i.e. the TS stacking configuration with space P1) and by (a + b)/2 (i.e. the AA′ stacking configuration with space group Pmma), respectively.There are 225 stacked bilayer structures generated by shifting the top layer relative to the bottom layer on the 2D grid, and each bilayer structure was fully relaxed allowing the degrees of freedom in vertical direction to optimize interlayer distance between two layers.

Structural properties
The structural properties of optimized bilayer phosphorene under relative shifting are analyzed in terms of the bond lengths, the bond angles, and the equilibrium interlayer distances, as shown in figure 2 and listed in table 1.It was found that the in-plane bond lengths (b 2 , b 3 ) and bond angle (θ in ) are almost independent of the stacking configuration, and the outof-plane bond lengths (b 1 ) and bond angle (θ out ) show a slight variation (<0.05 Å and <0.47°).But the equilibrium interlayer distance (d) is sensitive to the stacking arrangement, which is the shortest at the AB stacking configuration (3.22 Å) and longest at the AA′ stacking configuration (3.77Å), reflecting the role played by the vdW interaction between two layers.The relative energy (E r ) with respect to the energy at AB stacking configuration for high symmetry shifted patterns (i.e.E r = E stacking − E AB , where E stacking and E AB are total energies of the combined systems for a given stacking configuration and the total energy of the combined system with the AB configuration, respectively) are listed in the last column of table 1.It was found that the most energetically favorable structure is the bilayer system with the AB stacking configuration, followed by the Aδ, TS, AB′, AA stacking configurations, and the AA′ stacking configuration has the highest relative energy.The small difference in relative energy (in the range of a few meV/atom) indicates that the structure transition from one stacking configuration to another can be realized experimentally by sliding or can even happen during synthesis of bilayer phosphorene [53,54].
The dynamic stability of optimized bilayer black phosphorene with high symmetric shifted patterns has been studied from the analysis of phonon dispersion.As shown in figure 3, the bilayer system with the AB, Aδ, AB′, and TS stacking configurations has no imaginary frequency for all optical branches, confirming its dynamical stability with these stacking configurations [53].There are few small negative values (<20 cm −1 ) near the Γ edge point in the AA and AA′ stacking configurations, which were confirmed that they are not from the imaginary frequencies of transition states but are numerical errors due to the symmetry broken during the phonon spectrum calculations [43,44,[55][56][57].Thus, bilayer phosphorene with AA and AA′ stacking configurations are also dynamically stable but energetically are less preferential as compared with AB, Aδ, TS, and AB′ stacking configurations (see E r in table 1).The phonon dispersion bands at Γ point (including three Raman active modes characterizing the out-of-plane (A g 1 ) and the in-plane (B , 2 ) vibrations) show a tiny change in different stacking ordering (within the range of ∼2 cm −1 ).The low frequency phonon dispersion branches, on the other hand, show similarity between AB and AB′, or Aδ and AA′, or AA and TS, respectively.These features indicate the correlation between the interlayer coupling and  the stacking configuration.Furthermore, the phonon dispersion (reflecting the lattice vibration) also plays a significant role in understanding sound velocity, thermal, elastic, and optical properties of materials.The temperature dependence of the free energy, entropy, and heat capacity (C v ) of the shifted patterns of bilayer phosphorene with various stacking configurations are calculated (equations ( 1)-( 3)) and presented in figure S2.They are almost independent of the stacking ordering, i.e. the free energy (entropy) monotonically decreases (increases) with the elevating temperature, and the heat capacity first increases and then saturates to 200 (J/K/ mol) after 800 K regardless the stacking effect.

Potential energy surface (PES) and energy barrier
An important characteristic of interlayer interaction in 2D bilayer systems is the PES (also called the interlayer interaction energy [29,30]) which represents the energy landscape of a system as a function of its atomic coordinates.In the case of bilayer phosphorene, the PES describes the energy variation as the positions of the phosphorus atoms in both layers and provides insights into the structural changes, phase transitions, and chemical reactions of the material.The Van der Waals interaction between the two layers of phosphorene plays a crucial role in determining the PES. Figure 4 shows the contour (left) and 2D profile (middle) of the relative PES values (with respect to that at AB stacking configuration) under various shifting between the top and the bottom phosphorene layers.An anisotropic behavior of PES was observed associated with minimum/maximum values at high symmetrical stacked structures (see the symbols on the PES).The global minimum or the most energetically preferential point (i.e. the dark blue area on the left/middle panels of figure 4) is found when sliding the top layer with  Stacking  .Therefore, the intrinsic resistance is stronger when sliding along the diagonal direction, but weaker when sliding along the MEP.Furthermore, the small transition barrier on PES along MEP indicates that the bilayer system at the metastable state with Aδ stacking can be self-retracted to the most stable state with AB stacking even below the room temperature.Thus, bilayer phosphorene can undergo structural transitions between different stacking configurations or phases.These transitions can be driven by temperature, strain, or other external factors.
Interestingly, a linear relationship between the PES and the equilibrium interlayer distance surface was identified by comparing the 2D profiles of PES (the middle panel on figure 4) and the interlayer distance (the right panel of figure 4).As seen in figure 4, the shorter the interlayer distance, the lower the PES (e.g. the dark blue areas around the AB stacking configuration), and the weaker the intrinsic resistance for sliding.Oppositely, the larger the interlayer distance, the higher the PES (e.g. the dark red areas around the AA′ stacking configuration), and the stronger the intrinsic resistance for sliding.Apparently, such straightforward correlation between the interlayer distance and the PES directly identify the role played by the attractive vdW interaction and the interfacial hybridization between two phosphorene layers, since the interlayer distance, primarily governed by van der Waals forces, affects the preferred stacking configurations of the bilayer phosphorene and influences the energy barriers for structural transitions between different stacking configurations.As the interlayer distance changes, the energy required to transition from one stacking arrangement to another can vary.

Electronic band structures
The electronic band structures and total density of states (DOS) of high symmetric shifted patterns of bilayer black phosphorene are calculated at the DFT-PBE level (figures S3 and S4) and corrected by a hybrid functional at the DFT-HSE06 level [45] (figure 6).They show a semiconduction behavior with the bandgap ranging from 0.09 (0.69) eV at AB′ stacking configuration to 0.61 (1.22) eV at TS stacking configuration at DFT-PBE (DFT-HSE06) level [24,52].The underestimated bandgaps at DFT-PBE level are corrected at DFT-HSE06 level with a rigid shift of ∼0.6 eV (see the 2nd and 3rd columns in table 2), while the nature of the bandgap keeps the same at both levels.It was found that the nature of the bandgap depends on the stacking configuration (the 4th column in table 2).In most cases, the bandgap shows indirect behavior having the conduction band minimum (CBM) located at the Γ point and the valence band maximum (VBM) located along the Γ and X edge points (e.g.AA and AA′ stacking configurations) or along the Γ and Y edge points (e.g.Aδ and TS stacking configurations).But in the cases of AB and AB′ stacking configurations, the direct bandgap feature with both VBM and CBM located at the Γ point was found.The direct to indirect band gap transition is found mainly controlled by the curvature of VBM near the Γ point (figure 6).A small deviation of the VBM from the Γ point (less than 10%) along the Γ-X/Y path leads to a transition from the direct (e.g.AB or AB′ stacking configuration) to the indirect bang gap (e.g.AA or AA′ stacking configuration).From the partial band composition (figure S3), the partial density of states (DOS) (figure S4), and the partial charge density distribution of VBM and CBM (figure S5), we found that the feature of VBM band is mainly dominated by p z orbitals near the lone pairs in the interfacial region.Such orbital hybridization between layers strongly depends on the stacking arrangement: the p z orbitals at the top layer are staggered to those at the bottom layer in the AB (AB′) stacking arrangement, but they face each other in the AA (AA′) stacking arrangement (figure S5).Since the phosphorene monolayer possesses a direct bandgap, the stacking symmetry dependence of the bandgap nature in bilayer phosphorene directly reflects the effect of the interlayer coupling between layers.
Furthermore, it was found that the degrees of slopes of the top valence band and the bottom of conduction band, characterizing the charge carriers' dispersion, are sensitive to the dispersion direction in the k-space.Calculated the effective mass for electron (m e * ) and hole (m h * ) carriers (5-8th columns in table 2) are about an order in magnitude heaver along the Γ-X direction than along the Γ-Y direction, retaining the anisotropy of bilayer phosphorene in band structures and charge carriers [51,54,59].In addition, due to the flat nature of the VBM band along the Γ-X direction near the Γ edge point, the m h * is heavier than m e * and disperses slowly.
The landscape of bandgap character associated with the sliding was analyzed from the band structure calculations for all shifting on the 2D grid (shown on figure S1 with 225 shifts) and plotted on the 2D profile (figure 7).The color at the right-side column in figure 7  The partial charge density distribution of VBM and CBM (with the isosurface of 4.2 × 10 −3 e Å −3 ) for AA, Aδ, AB′, AB, TS, and AA′ stacking configurations (figure S5) show localized behavior and evenly distributed pattern on both layers, resulting in a strong interlayer interaction [54].The bonding nature and lone pairs, on the other hand, are clearly seen in VBM (as example, indicated in AA stacking case), and the antibonding nature are demonstrated in CBM.As the top layer shifts from AA along either the zigzag, or the armchair, or the diagonal direction, the lone pair orbital in CBM states in the interfacial region is enhanced.

Interfacial hybridization
Each phosphorus atom has five valence electrons.When they form black phosphorene monolayer, three of the five valence electrons form the sp 3 type of bonding and the remaining two are lone pair protruding away from the layer, which can interact with other phosphorus atoms at adjacent layers and play an important role in the electronic properties of phosphorene.When two phosphorene layers come into proximity forming bilayer phosphorene, the electronic states of each layer will couple, mostly induced by lone pairs, leading to a formation of interlayer hybridized states.Even though the Table 2. Calculated bandgaps of high symmetric shifted patterns of bilayer phosphorene (with AA, AB, Aδ, AB′, TS, and AA′ stacking configurations) at DEF-PBE (2nd column) and DFT-HSE06 (3rd column) levels.The nature of the bandgap is indicated in the 4th column together with the k points at VBM/CMB position (in the parentheses).The effective mass for electron (m e * ) and hole (m h * ) carriers along Γ-X and Γ-Y directions (DFT-HSE06 level) are listed in the 5th−8th columns in the unit of free electron effective mass (m 0 ).interlayer interaction is primarily governed by van der Waals forces, which are relatively weak compared to covalent bonds, these weak interactions can still significantly affect the electronic structure of the bilayer system, and associated interlayer hybridization can result in a creation of new electronic bands or a modification of existing bands and significantly influence the electronic and optical properties of the system.Consequently, the interlayer hybridization can lead to a renormalization of the electronic band structures (e.g.altering their effective masses and energy dispersions), introduce additional optical transitions (e.g.broadening the absorption spectrum and influencing the excitonic properties of the material), and influence the charge carrier transport behavior (e.g. carrier scattering processes, carrier mobility, and the overall conductivity of the material), etc.Our Bader analysis [46][47][48][49][50] found that a net charge transfer between two layers at the interfacial region occurs, mainly from the hybridization of lone pairs.Such charge transfer is sensitive to the stacking configuration with the smallest transfer (∼0.002 e) at the AB′ stacking configuration and largest transfer (∼0.011 e) at the TS stacking configuration (see table S1).This stacking dependence of charge transfer at the interfacial region directly demonstrates the character of the interlayer binding and reveals that the relative alignment of the two layers of phosphorene plays a crucial role in the interfacial hybridization.
To track how the net charge is redistributed in the interfacial region when the two phosphorene layers forming a shifted pattern and deeply understand the feature of the interfacial hybridization, we analyzed the differential electron charge density (DCD) of bilayer black phosphorene under various stacking configurations.Figure 8 presents the projected (top), the side view from front (middle), and the side view from left (bottom) of the DCD contours of the system (with the isosurface of 1.9 × 10 −4 e Å −3 ) with high symmetric shifted patterns (i.e. the AA, AB, Aδ, AB′, TS, and AA′ stacking configurations).The yellow color denotes the electron accumulation, and the blue color denotes the electron depletion, respectively.It was found that electrons deplete along the valley of the top layer and the crest of the bottom layer (except the AB′ and AA′ stacking configurations with antisymmetric symmetry).Additional electron depletion was found near the lowermost P atoms in the top layer and the uppermost P atoms in the bottom layer (i.e.around loan pairs).They align as depletion strips.On the other hand, the strong interlayer lone-pairs hybridization leads to a covalentlike charge localization at the midpoint of the interlayer region, forming accumulation strips.Furthermore, it was found that the shape and the location of such electron depletion/accumulation strips vary with the relative shifting between two layers in different translating directions.Apparently, such stacking symmetry dependence of the electron accumulation/depletion is originated from the interlayer hybridization, which is sensitive to the relative positions of the atoms in the interfacial region, hence strongly depends on the average interlayer distance (d) due to the different stacking configuration of the atoms.In addition, some isolated electron depletion and accumulation clouds are found located on the uppermost P atoms in the top layer and the lowermost P atoms in the bottom layer in the cases of AA, Aδ, and TS stacking configurations, indicating that the interlayer hybridization also spreads into the layers.But such phenomenon disappears in the cases of AB, AB′, and AA′ stacking configurations [61,62].
The plane-averaged differential charge density (Δρ) along z-direction normal to the phosphorene bilayer (black curves) and the amount of net charge electrons (ΔQ) transferred from bottom bilayer (at z = 0) up to z level (red curves) of high symmetric shifted patterns of bilayer black phosphorene with AA, AB, Aδ, AB′, TS, and AA′ stacking configurations are calculated and plotted in figure 9.The yellow color denotes the electron accumulation, and the blue color denotes the electron depletion, respectively.The space between two black-dashed lines in each panel denotes the interfacial region.A fluctuation of the net charge transfer between two layers was found.Except the cases of AA and AA′ stacking configurations where the fluctuation is small (∼|0.008-0.01|e), a larger fluctuation (∼|0.02|e) was found in other cases of AB, AB′, Aδ, and TS stacking configurations, indicating the strong orbital hybridization between layers in these stacking configurations.Consistent with the DCD analysis at the interfacial region, it was found that the net electrons accumulate at the bottom (top) surface of the top (bottom) layer (i.e.yellow peaks at the two black-dashed lines), followed by the net electron depletion near the lone pairs (i.e.blue peaks between the two black-dashed lines), and then a strong net electron accumulation located at the center of the region (i.e.yellow peaks at the middle point between the two black-dashed lines).They form alternative net positive (yellow peaks) and negative (blue peaks) charge strips parallel to the phosphorene bilayer, leading to a polarization in the interfacial region, totally different from the graphene bilayer system where no net charge was transferred between layers.Such polarization is mainly induced by the hybridization of lone pairs in the interlayer space and triggered by the vdW interlayer interaction.In addition, it was also found that, at the interfacial region, the distribution of the total amount of net charge electrons (ΔQ) along z-direction show a negative peak near the top surface of the bottom layer and a positive peak near the bottom surface of the top layer, forming a thin 'dipole-like' polarization layer.Exception is found in AA stacking configuration, where an oscillation of the positive/negative peaks appear inside the interfacial region.
To gain more insight into the stacking symmetry dependence of interlayer interaction, we calculate the electrostatic potential (i.e. the Coulomb potential between electrons and ions and between electrons) and depict the profiles of electrostatic potential along the z-direction for high symmetric shifted patterns of bilayer phosphorene in figure S7.Interestingly, different from the twisted moiré patterns of bilayer phosphorene [51] where the potential energy drops differently at two layers, the profiles of electrostatic potential along the z-direction are identical in shifted patterns of bilayer phosphorene.The potential drop from the uppermost P atoms in the bottom layer to the lowermost P atom in the top layer is almost zero, indicating the effective potentials are identical in adjacent layers and the interlayer coupling is strong which can affect the electronic structure, charge transport, and optical properties of phosphorene bilayers.Similarly, the work function varies slightly between 4.437 and 4.687 eV under a relative translation between layers.

Conclusion
A comprehensive first-principles study on the shifted patterns of bilayer phosphorene has been carried out to deeply understand the role played by the interlayer interaction in controlling the structural and electronic properties through relative translation.Structurally, the attractive vdW interaction stabilizes the bilayer with the shorter equilibrium interlayer distance by a shifting of a/2 from the AA stacking along the zigzag direction (i.e. the AB stacking), followed by a shifting of b/2 from the AA stacking along the armchair direction (i.e. the Aδ and AB′ stackings).The largest equilibrium interlayer distance is found by a shifting of (a + b)/2 from the AA stacking along the diagonal direction (i.e. the antisymmetric AA′ stacking).Energetically, the most stable state was found in the case of the AB stacking configuration, the metastable state was found in the case of Aδ stacking configuration, and a transition state between them was found at the saddle point with the TS stacking configuration.Other high energy states are AA, AB′, and AA′ respectively, indicating a direct correlation between the equilibrium interlayer distance and the relative energy connected by the vdW interaction.Calculated PES shows a corrugation pattern with the lowest energy barrier of ∼1.3 meV/atom along MEP, an energy barrier of ∼9.3 meV/atom along the zigzag pathway, an energy barrier of ∼5.3 meV/atom along the armchair pathway, and the highest energy barrier of ∼11.2 meV/atom along the diagonal pathway, respectively.Such results provide insights into the stability, reactivity, and mechanical properties of bilayer phosphorene, facilitating its potential applications in electronics, optoelectronics, and beyond.The electronic bandgap is in the range of 0.69-1.22eV and shows a transition from the indirect bandgap to the direct bandgap under the translation, implying a tunable bandgap can be applied by stacking engineering.Furthermore, the orbital hybridization due to the lone pairs at the interfacial region induces a redistribution of the net charge associated with the sliding top layer with respect to the bottom layer, leading to a strong polarization along the zigzag direction with alternative arrangement of stripe-like electron depletion and accumulation and a net charge transfer of ∼|0.002-0.011|e between layers.The interplay between the interlayer coupling and the intrinsic property of each layer gives rise to unique electronic, optical, and transport phenomena in bilayer phosphorene.These novel findings will provide a fundamental reference to deeply understand and analyze the complex local structural and electronic properties of twisted bilayer phosphorene.Furthermore, our calculations also provide worthful references for friction application and bandgap engineering on phosphorene bilayers, indicating that bilayer phosphorene will be promising as versatile shiftronics material.

Figure 1 .
Figure 1.Schematic illustration of shifting the top layer with respect to the bottom layer of black phosphorene bilayer (starting from AA stacking configuration) on the 2D grid of 15 × 15 points.The high symmetric shifted patterns are automatically obtained by shifting the top layer along the zigzag, armchair, and diagonal directions, respectively.The insets are the top (left) and side (right) views of those high symmetric shifted patterns along with corresponding stacking and shifting notations (x, y).

Figure 2 .
Figure 2. The top and side views of optimized high symmetric shifted patterns of bilayer phosphorene with (a) AA, (b) AB, (c) Aδ, (d) AB′, (e) TS, and (f) AA′ stacking configurations, respectively.The equilibrium interlayer distances are indicated by the values.(g) The schematic illustration of the bond lengths and angles in optimized black phosphorene bilayer.The in-plane bond lengths and bond angle are denoted by b 2 , b 3 , and θ in , and the out-of-plane bond length and angle are denoted by b 1 and θ out , respectively.The indices P# denote atom number, and the directions along the lattice vectors are represented by colored arrows.

Figure 3 . 1 )
Figure 3.The phonon dispersion spectra of bilayer phosphorene with AA, AB, Aδ, AB′, TS, and AA′ stacking configurations.The out-ofplane (A g 1 ) and in-plane (B g 2 and A g 2 ) vibration modes at Γ edge point are indicated.

Table 1 .
Optimized bond lengths (b 1 , b 2 , b 3 ), bond angles (θ in , θ out ), equilibrium interlayer distances (d), and the relative energy (E r ) with respect to the total energy of AB stacking configuration of the shifted patterns of bilayer black phosphorene.The space group for each pattern is indicated in the corresponding parenthesis.
respect to the bottom layer by a/2 along the zigzag direction (i.e. the AB stacking configuration)[21,53,54,58].A local minimum or energetically metastable point (see the light blue area on left/middle panels of figure4), on the other hand, was found when sliding the top layer with respect to the bottom layer by b/4 along the armchair direction (i.e. the Aδ stacking configuration).A minimum energy path (MEP)[59] was also identified between the energetically metastable (Aδ) configuration and the energetically most stable (AB) configuration with a transition barrier at the saddle point (indicated by the black-dashed line on the middle panel of figure4).Such saddle point corresponds to the stacking configuration by shifting the top layer with respect to the bottom layer along the diagonal direction by (a + b)/5 and is denoted as TS symmetry (left/middle panels of figure4).These results are consistent with other DFT calculations[60].The corrugation of the PES is found not only produced by the vdW interlayer interaction but also by the electrostatic interaction induced by the net charge redistribution due to the relative atomic positions between layers (see the discussion in 3.4).Apparently, such corrugation results in the intrinsic resistance during sliding[16] and the maximum energy that can be dissipated by frictional mechanisms.[31,60].The corrugation of the PES indicates that a certain energy is required for the bilayer system to transition between different stacking structures.To analyze the intrinsic resistance generated by the PES corrugation and determine the energy barriers, we investigated the potential energy profile along different relative translational pathways, as shown in figure5.There is no potential barrier along the zigzag direction (i.e.X-shift) if the top layer slides from AA to AB stacking configurations.Oppositely, if the top layer slides from the global minimum at AB to AA, it needs ∼9.3 meV/atom to climb over the potential peak at AA (see figure5(a)).Along the armchair direction (Y-shift) between AA and AB′ stacking configurations, on the other hand, the potential energy decreases from AA and reaches the local minimum at Aδ, then increases to AB′.The corresponding energy barrier along this pathway from Aδ climbing over AA is 5.3 meV/ and that from Aδ climbing over AB′ is 4.6 meV/atom, respectively (see figure 5(b)).Similarly, there is a potential valley at the saddle point along the diagonal direction between AA and AA′ stacking configurations, and the corresponding energy barriers are 4.1 meV/atom from TS climbing over AA and 11.2 meV/atom from TS climbing over the highest potential energy point at AA′, respectively (see figure 5(c)).The corresponding transition barrier along the MEP path from Aδ to AB is 1.3 meV/atom (see figure 5(d))

Figure 4 .
Figure 4.The contour (left) and the 2D profile (middle) of PES, as well as the 2D profile of the equilibrium interlayer distance (right) of shifted patterns of bilayer phosphorene.The colors at the right-side column of each panel present the scales of the relative energy E r (meV/ atom) and the equilibrium interlayer distance d (Å), respectively, with the lowest values in blue and highest values in red colors.The shifting paths along the zigzag and the armchair directions, as well as the global, local minima and the saddle point are indicated by black arrows.The stacking configurations are denoted by symbols.

Figure 5 .
Figure 5.The potential energy profile of shifted patterns of bilayer phosphorene along different pathways: (a) in the zigzag direction (X-shift) between AA and AB stacking configurations, (b) in the armchair direction (Y-shift) between AA and AB′ stacking configurations, (c) in the diagonal direction between AA and AA′ stacking configurations, and (d) along the MEP path between Aδ and AB stacking configurations.The corresponding energy barriers along various paths are presented.
presents the scale of bandgap with the lowest value in blue and highest value in red colors.The direct (indirect) bandgap nature is indicated by D (I), and the stacking configurations are denoted by symbols.Clearly, the bandgap is sensitive to the relative translation of one layer relative to the other and shows in an anisotropic pattern.It fluctuates in a small range of ∼0.05 eV with relative translation from AA to AB along the zigzag direction, while increases from AA to Aδ by about 0.25 eV and then decreases by ∼0.52 eV from Aδ to AB′ with the relative translation along the armchair direction (see figure S6).Similarly, an anisotropic bandgap nature was observed on 'bandgap' surface with most regions having indirect bandgap nature except specific domains around AB and AB′ stacking configurations which have the direct bandgap nature.Thus, a transition from indirect to the direct bandgap happens when sliding the top layer with respect to the bottom layer from AA (TS and AA′) to AB stacking configurations (see figures 7 and S6).

Figure 6 .
Figure 6.Calculated band structure of high symmetric shifted patterns of bilayer phosphorene (with AA, AB, Aδ, AB′, TS, and AA′ stacking configurations) at the DFT-HSE06 level.The values of bandgap (E g ) are presented for each stacking ordering and the direct (indirect) bandgap nature was indicated by the red arrows.

Figure 7 .
Figure 7.The 2D profile of the bandgap character (DFT-PBE level) of the shifted patterns of bilayer phosphorene.The color at the rightside column presents the scale of bandgap energy with the lowest value in blue and highest value in red colors.The direct (indirect) bandgap nature is indicated by D (I), and the stacking configurations are denoted by symbols.

Figure 8 .
Figure 8.The projected (top), side view from front (middle), and side view from left (bottom) of the DCD contours (with the isosurface of 1.9 × 10 −4 e Å −3 ) of high symmetric shifted patterns of bilayer phosphorene with AA, AB, Aδ, AB′, TS, and AA′ stacking configurations, respectively.The yellow color denotes the electron accumulation, and the blue color denotes the electron depletion.The pink color denotes the phosphorus atoms.

Figure 9 .
Figure 9.The plane-averaged differential charge density (Δρ) along z-direction normal to the phosphorene bilayer (black curves) and the amount of net charge electrons (ΔQ) transferred from bottom of the bilayer (at z = 0) up to z (red curves) of high symmetric shifted patterns of black phosphorene bilayer with AA, AB, Aδ, AB′, YS, and AA′ stacking configurations, respectively.The yellow color denotes the electron accumulation, and the blue color denotes the electron depletion.The space between two black-dashed lines in each panel denotes the interfacial region.