Stacked bilayer phosphorene: strain-induced quantum spin Hall state and optical measurement

Bilayer phosphorene attracted considerable interest, giving a potential application in nanoelectronics owing to its natural bandgap and high carrier mobility. However, very little is known regarding the possible usefulness in spintronics as a quantum spin Hall (QSH) state of material characterized by a bulk energy gap and gapless spin-filtered edge states. Here, we report a strain-induced topological phase transition from normal to QSH state in bilayer phosphorene, accompanied by band-inversion that changes number from 0 to 1, which is highly dependent on interlayer stacking. When the bottom layer is shifted by 1/2 unit-cell along zigzag/armchair direction with respect to the top layer, the maximum topological bandgap 92.5 meV is sufficiently large to realize QSH effect even at room-temperature. An optical measurement of QSH effect is therefore suggested in view of the wide optical absorption spectrum extending to far infra-red, making bilayer phosphorene a promising candidate for opto-spintronic devices.

QSH material should have a large bandgap and be easily fabricated. Graphene has a superior high carrier mobility but its bandgap is extremely small [31][32][33] , while bilayer bismuth (111) is potentially large bandgap material but its fabrication is usually difficult due to its strong interlayer bonding 34,35 . Bilayer phosphorene has not only a high carrier mobility and a large bandgap but also easy fabrication. Hence, bilayer phosphorene can be a potential excellent candidate for future applications in QSH devices, such as topological field-effect transistors and spin valve devices.
In this work, we will investigate a quantum spin Hall state of bilayer phosphorene by tunning strain and stacked order, then find a topological phase transition from normal to QSH state accompanied by a band inversion that changes  2 number from 0 to 1. Furthermore, the direct-indirect bandgap and semiconductor-metal transitions can be observed by adjusting the weight of vdW interaction. Our results reveal that the tuning of topological behavior in bilayer phosphorene is dependent on interlayer stacking order and the direction of the applied in-plane strain that should be along either armchair or zigzag direction. When the bottom layer is shifted by 1/2 unit cell along the zigzag or the armchair direction with respect to the top layer, the maximum topological bandgap can reach up to 92.5 meV, which is sufficiently large to realize the quantum spin Hall effect even at room temperature. The optical property of bilayer phosphorene is examined based on a real-space and real-time time-dependent density functional theory, which shows that the optical absorption spectrum of bilayer phosphorene becomes wide and is even extended to far-infra-red region in QSH state. Such improvement in optical response is indispensable for the broadband photodetection. An optical experimental setup is therefore improved to measure the QSH effect in bilayer phosphorene.

Results
Crystal structures and electronic structures of bilayer phosphorene for four stacking orders. The crystal structure of phosphorene is given in Fig. 1. Bilayer phosphorene can be viewed as cleaved from the (0001) surface of black phosphorus (space group Cmca), which allows various stacking orders due to its weak interlayer van der Waals (vdW) interaction, as shown in Fig. 1(a,b). We consider four possible high-symmetry stacking orders: (Ta) the top layer is stacked vertically on the bottom layer (space group Pmna), (Tb) the bottom layer is shifted by 1/2 unit cell along x (zigzag) or y (armchair) direction with respect to the top layer (space group Pbcm), (Tc) the bottom layer is shifted by one unit cell along x or y direction with respect to the top layer, and thus the top and bottom layers are mirror images of each other (space group Pmma), and (Td) the bottom layer is shifted by 3/2 unit cells along x or y direction with respect to the top layer (space group Pccm), as shown in Fig. 1(c-f). Table 1 gives the optimized lattice constants, the bond lengths and the other structural parameters of the bilayer phosphorene for the four stacking orders, which is well consistent with the previously theoretical data 26 with errors lower than 0.5%. For the different stacking orders, the bond length R 1 and ′ R 1 are almost same, being shorter than the bond length R 2 . The bond angles α is smaller than β. The smallest interlayer distance d int is 3.503 Å (Ta), 3.085 Å (Tb), 3.739 Å (Tc) and 3.291 Å (Td). Among the four stacking orders, the Tb stacking order has the lowest cohesive energy. Therefore, the phonon spectrum of the Tb stacking order is further calculated, which shows the excellent structural stability (see Supplementary information). The cohesive energies of the Ta, Tc, Td stacking orders are slight larger than the Tb stacking order with the difference less than 10 meV, which indicates that the Ta, Tc and Td stacking orders also have the excellent structural stability.
The electronic structure of the bilayer phosphorene is given in Fig. 2. Shown in Fig. 2(a-d) are the band structures of the four different stacking orders obtained by using Perdew-Burk-Ernzerhof (PBE) functional 38 , which gives the bandgap E g = 0.434 eV (Ta), 0.442 eV (Tb), 0.264 eV (Tc) and 0.002 eV (Td) at Γ point. Since our hybrid Heyd-Scuseria-Emzerhof (HSE06) 39,40 calculations indicate that the PBE E g values are underestimated about 0.56 eV, the E g of the bilayer phosphorene should be larger than 0.56 eV for the four different stacking orders, which is well consistent with the previous HSE06 data 26 . The iso-surfaces of the band-decomposed charge density of the valence band maximum (VBM) and the conduction band minimum (CBM) at Γ point show that the four stacking orders have the similar charge density distribution and bonding character, as shown in Fig. 2(e-h). Basically, the charge density distributions of the VBM and CBM in one nonplanar monolayer have no effect on the stacking order: an anti-bonding-like feature in one sublayer and a bonding-like feature between two sublayers are visible for the VBM, while the bonding-like features in one sublayer and the anti-bonding-like features between two sublayers are found for the CBM. On the contrary, the charge density distributions of the VBM and CBM between two layers are dependent on the stacking order: the sign of the chemical bond between two layers is not visible for the VBM in the four stacking orders, but a chemical bonding-like character between two layers for the CBM is found in the Tc and Td stacking orders, which is absent in Ta and Tb stacking orders. The bonding features for the CBM in the Tc and Td stacking orders push their CBM down to the low energies, which leads to the smaller bandgaps than the Ta and Tb stacking orders. Therefore, it is confirmed that the different chemical bonding in the interfacial area between two layers is related to the bandgap E g for the four stacking orders.
Quantum spin Hall state induced by in-plane strain in bilayer phosphorene. Figure 3 illustrates the variation of the bandgap of the bilayer phosphorene versus the in-plane strain for the four stacking orders, where σ x , σ y and σ xy are the in-plane strains along x-direction, y-direction, and both xand y-directions, respectively. When the in-plane strain is applied, the lattice constants along the other directions are relaxed fully until the residual force on each atom is less than 0.01 eV/Å in order to ensure the complete relaxation of the crystal structure for the strained bilayer phosphorene. In our work, the in-plane uni-axial strain is defined as σ = ( − )/ a a a x 0 are lattice constants along the x,y directions under strain, respectively, and a 0 ,b 0 are the corresponding equilibrium lattice constants without strain. The in-plane strain σ xy is loaded synchronously in the x-and y-directions. In experiment, the in-plane strain on the bilayer phosphorene can be realized by bending its flexure substrate similar to graphene, where the amount of the strain is proportional to 2D mode position of the bilayer phosphorene 41 . Bandgaps of the four stacking orders are more sensitive to σ xy , as compared with σ x , σ y . When σ xy > 0 (tensile), the bandgap increases with σ xy , then turns to decrease at a critical value, being 4.0% (Ta), 3.0% (Tb), 5.0% (Tc and Td). When σ xy < 0 (compression), the bandgap decreases with σ xy and becomes zero at σ xy = − 3.0% (Ta and Tc), − 2.77% (Tb) and 0.02% (Td). However, to our surprised, after the bandgaps are closed, E g turns to increase again for the Ta, Tb, and Td stacking orders [see the inset where filled and opened circles are two different sublattices, R 1 and ′ R 1 are P-P bond lengths in up and bottom sublattices, and α is the angle between two R 1 bonds. One unit cell has four atoms, as denoted by the blue shadow. (c-f) The top and side views of the four different stacking orders: (Ta) the top layer is stacked vertically on the bottom layer, (Tb) the bottom layer is shifted by half of one unit cell along x or y with respect to the top layer, (Tc) the bottom layer is shifted by one unit cell along x or y direction with respect to the top layer, and thus the top and bottom layers are mirror images of each other, and (Td) the bottom layer is shifted by one and a half of one unit cell along x or y direction. Here, R 2 is the P-P bond length between up and bottom sublattices of a nonplanar monolayer, β is the angle between the R 2 and ′ R 1 bonds, and d int is the smallest interlayer distance.
Next, the electronic structure of the Tb-stacked bilayer phosphorene is given in Fig. 4. At Γ point under σ xy = 0, the conduction band (CB) and the valence band (VB) have the even ('+ ') and odd ('− ') parities, respectively. Shown in Fig. 4(a-c) are the electronic band structures for the Tb stacking order under different σ xy values. The CB and VB at Γ point tend to approach together as σ xy increases, then become overlapped under σ xy = − 2.77%. For larger compression strain (σ xy = − 3.0%), the CB and VB at Γ point are separated because of the repulsion between them. Most remarkably, the parities of the CB and VB at Γ are exchanged when the compression strain σ xy < − 2.77%, as shown in Fig. 4(d), which shows the band inversion of the VB and CB at Γ . Such band-inversion character is also observed in the density of states (DOS) and the orbital-projected band structures when σ xy = − 3.0% [see Fig. 4(e-g)]. The p orbital makes a significant contribution to the total DOS. The conduction band near E F mainly comes from the p y,z obitals, but the bottom of the conductor band originates mostly from the p z . The valence band near E F mainly comes from the p z obital, but the top of the valence band originates mostly from the p y orbital, which indicates an band inversion process when the compression strain σ xy is increased.
Further, we calculate the  2 number of the Tb-stacked bilayer phosphorene when σ xy = 0 and − 3.0%, by using the method of Fu and Kane 42 . Such method is valid since the Tb stacking order has both spatial invention and time reversal symmetries (four time reversal invariant points in the 2D Brillouin zone). Inversion center in the crystal ensures ε ε is the electron energy for the n-th band with spin index α at k wave vector in the Brillouin zone. The time reversal symmetry makes , where α is the spin opposite to α. The calculated parities of all occupied bands at four time-reversal invariant momenta are summarized in Table 2. We can find that the product of parities of occupied bands contributes to a + 1 parity at the four time-reversal invariant momenta when σ xy = 0, yielding a trivial topological invariant =  0 2 . As the strain is increased up to σ xy = − 3.0%, the product of parities of occupied bands is − 1 at Γ but + 1 at the three other time-reversal invariant momenta like X, Y, and S. The Tb-stacked bilayer phosphorene under σ xy = − 3.0% is identified as topological insulators The results shown in Fig. 4(d) suggest that the Tb-stacked bilayer phosphorene maintain to be the QSH state under the compression strain σ xy = − 2.77 ~ − 7.0%, where the maximum topological bandgap E g = 92.5 meV is obtained when σ xy = − 5.0%.
The PBE method may underestimate the bandgap compared with reality, which leads to a low evaluation of the critical strain needed to produce a quantum spin Hall state. In order to avoid such problem, we have carried out the band structure calculations of the possible topological insulators (Ta and Tb stacking orders) under the in-plane strains by the HSE06 method which has been proven reliable for few layer phosphorene systems 26 . Fig. 5(a,b) present the P-p z orbital-projected band structures of the Tb-stacked bilayer phosphorene by the HSE06 method when σ xy = − 4.0% and − 6.0%, respectively. We can find that the top of the valence band and the bottom of the conduction band for the Tb stacking order exchange the weight of P-p z orbital when the strain σ xy = − 6.0%, which shows a band inversion. Actually, the Tb-stacked bilayer phosphorene has already been in quantum spin Hall state when the strain is beyond − 4.8%, as shown in Fig. 5(c). The corresponding strain range where the Tb-stacked bilayer phosphorene maintain to be the QSH state is − 4.8 ~ − 10.0% by the HSE06 method. It is worth  60 , where a and b are the lattice constants along x and y directions, respectively, d int is the smallest interlayer distance, R 1 and ′ R 1 are the P-P bond lengths in up and bottom sublayers of a nonplanar monolayer, respectively, R 2 is the P-P bond length between up and bottom sublayers of a nonplanar monolayer, α is the angle between two R 1 bonds , β is the angle between the R 2 and ′ R 1 bonds, E coh is the cohesive energy, and ΔE coh is the relative difference of the cohesive energy for Ta, Tc, and Td stacking orders and the Tb stacking order.
Eg=0.264 noting that the lattice structure of the Tb-stacked bilayer phosphorene is still stable under the compression strain up to − 5.0% by calculating its phonon spectrum (see Supplementary information), which is critical for the practical applications. The topological phase transition for the Ta stacking order induced by in-plane strain is proven by the same method. Like the Tb stacking order, the Ta stacking order also undergoes the bandgap closed and reopened process with the increasing σ xy , as shown in Fig. 5(d), which shows a band inversion. Shown in Fig. 5(e) is the evolution of the bandgap with σ xy , which indicates that the PBE bandgap is underestimated compared with the HSE06 bandgap. Hence the critical strain σ xy (about − 5.0%) obtained by the HSE06 method, where the Ta-stacked bilayer phosphorene is in quantum spin Hall state, is higher than that (about − 3.0%) by the PBE method. For the Td stacking order, the VBM and CBM have the same parity, which rules out the possibility of topological phase transition.

Semiconductor-metal transition induced by interlayer interaction in bilayer phosphorene.
The electronic band structure as function of the weight of van der Waals interaction between the top and bottom layers of the bilayer phosphorene is shown in Fig. 6. The weight of van der Waals interaction is incorporated into our DFT calculations by using a semi-empirical van der Waals approach, known as DFT-D2 method 43 . The total energy is defined as the Kohn-Sham DFT total energy, E vdW is the vdw energy, and W is a scaling factor in order to consider the weight of van der Waals (WvdW) interaction. The variation of the bandgap E g of the bilayer phosphorene versus the WvdW interaction is illustrated in Fig. 6(a). The E g tends to decrease with the increasing WvdW interaction for all four stacking orders. When the increasing WvdW interaction is beyond 2.70 t (t is the real value for the WvdW interaction), E g for each stacking order will decrease gradually to zero, which indicates that the VBM and CBM of the bilayer phosphorene have been overlapped, yielding a semiconductor to metal transition. For the Td stacking order, the decrease of E g versus WwdW is much stronger than the Ta, Tb, and Tc stacking orders. Very specifically, we find a bandgap closed and reopened process for the Td stacking order, as shown in the inset of the Fig. 6(a), whereas this phenomenon is absent for the Ta, Tb and Tc stacking orders. The P-p z orbital-projected band structures of the Td-stacked bilayer phosphorene are presented in   Fig. 6(c)], which can be due to the reason that the top and the bottom layers of the Td-stacked bilayer phosphorene tend to be bonded as companied with decrease of d int . Furthermore, as accompanied by the bandgap closed and reopened process, the VBM and CBM for the Td stacking order exchange the weight of P-p z orbital, indicating the band inversion of the VBM and CBM. However, as mentioned above, the VBM and CBM for the Td stacking order have the same parity, which rules out the possibility of topological phase transition. Hence, despite that the WvdW interaction causes the decrease of E g and the break of the band degeneracy for the four stacking orders, but does not cause the topological phase transition, as shown in our computations. Our calculations establish the basic conditions for materials to exhibit strain-induced topological phase transition. First, the direction of the applied strain on materials should be along the bonding direction, i.e., the direction with the maximum atomic wave function overlap. For the bilayer phosphorene,  Optical response of bilayer phosphorene and optical measurement of quantum spin Hall effect. As the most primary structure inheriting from the bulk black phosphorus, the Tb stacking order has the highest structural stability and is most facile to be fabricated in experiment among the four stacking orders. Therefore, the Tb-stacked bilayer phosphorus would be the best phosphorus candidate for the future optical electronic devices. The optical response property of the Tb-stacked bilayer phosphorene is investigated, as shown in Fig. 7. The optical response calculation is conducted for 4 × 4 supercell of Tb-stacked bilayer phosphorene [see Fig. 7(a)]. In Fig. 7(b), we compare the optical response of the Tb-stacked bilayer phosphorene in the normal state (σ xy = 0) and the QSH state (σ xy = − 3.0%), where the structural stability in QSH state is proved by calculating the phonon spectrum (see Supplementary information). Our results show that the photonic band-gap is decreased in the QSH state greatly, for example, the obtained photonic band-gap (PBG) is 0.5 eV in the normal state but only 0.01 eV in the QSH state for the case of the incident light polarized along the armchair direction. There exists two adsorption peaks in the optical adsorption spectrum. The first adsorption peak appears in the region that the optical strength of the incident light is less than 2 eV and the corresponding adsorption strength is also very low. When the incident light is polarized along the armchair direction, the first adsorption peak shows a red-shift in the QSH state [see red arrow in inset Fig. 7(b)], which leads to wider absorption zone that is highly desirable in broadband photodetection. The second adsorption peak occurs at the region of the optical strength of the incident light around 10 eV, which indicates that the bilayer phosphorene occurs a σ plasmon resonance. When the incident light is polarized along the zigzag direction, the second absorption peak becomes smaller and shows a blue shift in the QSH state, as compared with the normal state [see the black and blue arrows in Fig. 7(b)]. The blue shift is because that the increased interatomic spacing in the zigzag direction causes the larger resonance level spacing. The mechanism of the optical absorption in low-energy resonance behind the resonance phenomena is further analyzed in terms of the induced charge response in real-time propagation. The induced electron and hole charge densities tends to be separated and finally located at two edges. The induced charge density in the QSH state is larger than that in the normal state, emerging in the center areas in addition to the two edges. This is attributed to that the shielding effect becomes stronger in the QSH state when the interatomic spacing becomes smaller due to the in-plane strain (see Supplementary information).
We propose an experimental setup [see Fig. 7(c,d)] to measure QSH effect of strained Tb-stacked bilayer phosphorene by using Scanning Kerr rotation microscopy. Our setup has the similar working principle with the previous optical measurement of local spin accumulation at the sample edges in n-GaAs and n-InGaAs [44][45][46] , but combines with the electronic measurement. Experimentally, because the transverse spin currents do not lead to net charge imbalance across the sample in nonmagnetic systems, the spin Hall effect is difficult to be observed by the simple electrical measurement. Hence, we employ the Scanning Kerr rotation microscopy to image the spin polarization at the edges of the samples that has 320-μm-long and 100-μm-wide, as shown in Fig. 7(c). The probe beam is bent 90° by a prism mirror placed in the He-flow cryostat, then is incident upon the samples through an objective lens with 0.73 μm lateral spatial resolution, which ensures that the beam is focusing on a circular spot with a full width at half-maximum of 1.1 mm 44 . An image from the Scanning Kerr rotation microscopy is taken by moving the cryostat integrated on a piezo nanoscanning stage, just like in the case of GaAs 46 . The Tb-stacked bilayer phosphorene is grown epitaxially on a substrate of comparable thickness, which makes the layers strained due to a mismatch in lattice parameter between the film and substrate materials. We set the center of the sample ("O" point) as the origin of the coordinate, and x and y axes are taken to be across and along the sample channel, respectively. This channel is defined by the standard photolithography and the wet etching. Kerr rotation is measured as a function of magnetic field (B ext ) and position (x, y), which is at a maximum at the two the channel edges when B ext = 0 mT and falls off rapidly with distance from the edge 44 . This indicates an accumulation of electron spins polarized at the edges, which is also a strong signature of the spin Hall effect because the spin polarization is out-of plane and changes sign for opposing edges [47][48][49] . Secondly, the two-terminal conductance G SD can be measured by using the electronic measurement, as shown in Fig. 7(d). When sample is a normal insulator, G SD should vanish at low temperatures, whereas G SD should equal a value ~ne 2 /h (n is the integer, e is the electronic charge and h is the Planck constant). It is worth noting that the electronic voltage V 1 and V 2 vanish between M 1 and M 2 and between M 3 and M 4 , which is one of the remarkable manifestations of the QSH state 50 .

Discussion
We have proven that the Ta-and Tb-stacked bilayer phosphorene are excellent topological materials when the compression strain loaded synchronously in the x-and y-directions is larger than about 3%. The quantum spin Hall (QSH) state of the bilayer phosphorene is revealed due to the band inversion between the valence band maximum of p z orbital and the conduction band maximum of p y and p z orbitals. Our results show that the strain manipulation is very effective for the bilayer phosphorene.
In addition to the semiconductor-to-metal transition, strain can also induce the nontrivial topological band structure of the bilayer phosphorene. The obtained maximum topological bandgap is high up to 92.5 meV in the Tb-stacked bilayer phosphorene, implying that the QSH effect of bilayer phosphorene is observable even at room temperature. On the other hand, as compared with the chemical decoration and the electric field, the strain manipulation is more feasible in experiment, considering the current microfabrication technique, for example, heterepitaxy of the bilayer phosphorene grown on substrates of smaller lattices. Our results confirm that the bilayer phosphorene is an ideal material to study novel quantum states, showing new potential for future applications in spintronics devices such as topological field-effect transistors and spin valve devices. Furthermore, the strain is also an excellent way to tuning optical property besides those aforementioned. The optical absorption spectrum of the strained bilayer phosphorene becomes broadened, and is even extended to the far-infra-red region and leads to a wider range of brightness. An improvement optical experimental equipment is designed therefore to measure the QSH effect in strained bilayer phosphorene, which shows that the strain-induced topological bilayer phosphorene can also be used as an outstanding opto-spintronic device. Electronic structure calculation. We calculate the lattice configurations as well as electronic band structures of bilayer phosphorene with four different stacking orders based on the density functional theory (DFT) implemented in the Vienna Ab-initio Simulation Package (VASP) 51,52 . The projector augmented wave (PAW) 53 method and the Perdew-Burk-Ernzerhof (PBE) 38 exchange-correlation functional are adopted. Long-range dispersion corrections have been taken into account within a semi-empirical van der Waals approach proposed by Grimme known as the DFT-D2 method 43 (where D2 stands for the second generation of this method). The kinetic energy cutoff for the plane wave basis set is chosen to be 600 eV, and the reciprocal space is meshed at 14 × 10 × 1 using Monkhorst-Pack method 54  vacuum space of at least 25 Å along the z direction is used to separate the bilayer systems in order to avoid spurious interactions due to the nonlocal nature of the correlation energy 55 . In order to correct the PBE bulk bandgaps, we apply a hybrid Heyd-Scuseria-Emzerhof (HSE) 39,40 functional in which the exchange potential is separated into a long-range and a short-range part, where 1/4 of the PBE exchange is replaced by the Hartree-Fock exact exchange and the full PBE correlation energy is added. Hence the HSE functional corrects the GGA bandgaps 40 significant by partially correcting the self-interaction. Spin-orbit coupling (SOC) is included in the calculations after the structural relaxations 56 .

Methods
Optical response. We calculate the optical responses of the Tb-stacked bilayer phosphorene under the strains σ xy = 0.0 and − 3.0%, based on a real-space and real-time time-dependent density functional theory (TDDFT) as implemented in the OCTOPUS code 57 . The Hartwigsen-Goedecker-Hutter pseudopotentials 58 and Generalized Gradient Approximation (GGA) with PBE functional for the exchange-correlation are used to calculate both the ground state and excited state. We consider a 4 × 4 supercell of Tb-stacked phosphorene [see Fig. 7(a)] in order to ignore the interaction between the two boundaries. The simulation zone is defined by assigning a sphere around each atom with a radius of 6 Å and a uniform mesh grid of 0.3 Å. The induced density plane paralleled to the bilayer phosphorene plane, and located at the middle of the top pucker-layer in the vertical direction. In the real time propagation, excitation spectrum is extracted by Fourier transform of the dipole strength induced by an impulse excitation 59 . In the real-time propagation, the electronic wave packets are evolved for typically 6000 steps with a time step of Δ t = 0.005 Å/eV. The obtained optical bandgap is close to the electronic bandgap calculated by PBE implemented in VASP, which means that OCTOPUS is reliable on calculation for the optical property.