Tunable Dirac cones in single-layer selenium

Dirac cone, one of the main characters of topological materials, provides us an approach to explore topological phase transitions and topological states. Single-element 2D-Xenes are prominent candidates for hosting Dirac cones. Till now, the multiple Dirac cones, Dirac-like cones, and semi-metal Dirac point have been discovered in them. However, it is still difficult to realize the tunable Dirac cones due to the lack of appropriate materials. Using first-principles calculations, this paper proposes that monolayer selenium with square lattice could achieve tunable Dirac cones and a topological phase transition. Double structural phases of the monolayer selenium can be distinguished according to strain applied, i.e., buckled square and buckled rectangular phases, which have rich Dirac physics. There exist four anisotropic Dirac cones in the buckled square phase, owing to fourfold symmetry. The buckled rectangular phase hosts a topological phase transition from a 2D topological insulator with double Dirac cones to a simple insulator, with a Dirac semi-metal having single Dirac point as the phase transition point. Moreover, the topological insulator has a global band gap of 0.16 eV, suggesting its potential utilizations in room-temperature devices. These studies will greatly promote the development of the Dirac physics and widen the application ranges of 2D-Xenes.


Introduction
In the past decade, the 2D-Xenes are being extensively explored for potential applications in nanoelectronic and spintronic devices [1][2][3][4][5]. Many extraordinary properties originate from the Dirac-cone-type band structure which is immune to non-magnetic impurities and disorder [6][7][8][9][10]. Dirac cone, protected by crystal symmetry or accidental degeneracy, is featured as linear band dispersion around the central degenerate point in momentum space. After its discovery in graphene, the investigations on single-element 2D-Xenes (X = Si, Ge, Sn, B, P, Bi, etc.) are inspired. Various types of Dirac cones have been reported: 6 Dirac cones [11][12][13][14], 4 Dirac cones [15,16], 2 Dirac cones [17][18][19], Dirac-like cones [13] and semi-metal Dirac point without spin-orbit coupling (SOC) interaction [11,12]. In order to enrich the physics of topological phase transition and to enlarge the application range of 2D-Xenes in devices, it is in urgent need to design a new system which can manipulate the features of Dirac carriers, such as density, mass, velocity, and so on.
In this work, using first-principles calculations, we report that the monolayer selenium with square lattice hosts tunable Dirac cones and a topological phase transition both via the strain engineering. Under different strain levels, two structural phases are identified, i.e., buckled square (BS) phase and buckled rectangular (BR) phase. When SOC interaction is taken into consideration, four anisotropic Dirac cones having the direction-dependent electronic conductivities emerge in the BS structure. In addition, the BR structure experiences a topological phase transition from a 2D topological insulator (TI) with two Dirac cones to a simple insulator, where the phase transition point is actually a Dirac semi-metal (DSM) containing single Dirac point. Furthermore, tuning the strain makes it possible to achieve the various mass of Dirac carriers by controlling its direct band gap. Therefore, the monolayer selenium provides us an encouraging candidate for Dirac cone engineering to achieve ideal Dirac carrier features.

Methods
Electronic properties of the 2D selenium system is explored by employing the Vienna ab initio simulation package (VASP) [20,21] to perform the first-principles calculations based on density functional theory. The projected augmented wave pseudopotential is used in the calculation [22,23] and exchange-correlation interaction is treated within the generalized gradient approximation in the form of Perdew-Burke-Ernzerhof (PBE) [24]. To obtain the optimized structure, the residual force is carefully checked to ensure it is less than 0.001 eV Å −1 , while the vacuum layer thickness is set to be 20 Å. In the process of calculating self-consistent electronic structure, k-point mesh of 24 × 24 × 1 is used. The nonlocal Heyd-Scuseria-Ernzerh of screened hybrid functional (HSE06) [25,26] is executed to prove the band structure results. To verify the topological properties, maximally-localised Wannier functions [27,28] are generated for constructing tight-binding models by using WannierTools [29].

Crystal structures
Up to now, ultra-thin two-dimensional (2D) materials can be synthetized in a variety of ways, including mechanical cleavage [30], molecular beam epitaxy [31], molten-salt-assisted chemical vapor deposition [32], and so on. Here, a monolayer selenium system can be obtained through cutting from an experimental bulk counterpart existing in high pressure phase [figure 1(a)] [33]. The optimized structure is a monolayer BR structure as shown in figure 1(b). Via the strain engineering, the system is found to be an affluent ground of topology. The 2D BR selenium hosts two inequivalent atoms in a unit cell since the middle atom deviates from the center. The optimized lattice constant is a 0 = 3.67 Å and the corresponding band gap is 0.12 eV. Additionally, with the strain change, the energy of the monolayer selenium varies accordingly, as shown in figure 1(c). The equiaxial strain is specific as (a − a 0 )/a 0 × 100% (positive strain means expansion, while negative one means compression). From figure 1(c), we can see that when the strain shifts from −6.0% to 9.0%, the change of energy is about 0.20 eV. Meanwhile, an interesting topological phase transition happens and it will be discussed later. The 2D selenium system possesses two structural phases depending on the applied strain. Firstly, the system with strain in the range from −13% to −8.0% is investigated and shows BS structure, which belongs to space group P 4 /nmm (No. 129). As shown in figure 2(a), one of the two inequivalent Se is presented in the middle. Then, in the large strain range from −8.0% to 13%, the system maintains the BR structure, belonging to space group C 2 /m (No. 12) [figure 2(b)]. Unlike the structure of BS structure, the middle atom of BR structure shifts away from the center. Figure 2(c) displays the 2D Brillouin zone (BZ) correspondingly.
In order to explain the mechanism of structural phase transition, we first calculate the strain dependence of angle θ between Se-Se bond and z axis (appendix A figure A1). We can see that when the strain is −8.0%, there appears a sudden decrease of angle, indicating a structural phase transition. In addition, there exist three kinds of bonds in selenium system. The strain dependence of bond lengths is then studied (appendix A figure A2). With the increase of strain, the differences among the three bond lengths in BR structure turn into larger and larger, while those in BS structure become zero. These above results imply that the phase transition is of the first order. It is also worth noting that when the strain is greater than 8.6%, there appears a small reduce both in angle and short bond length and a small increase both in medium and long bond lengths, suggesting a change of electronic property, which will be discussed later in detail.

Band structures and a topological phase transition
To understand the impact of various structure on topology, the energy band structures are further calculated. It is worth noticing that the topology related to band touching or band inversion is impacted through pressure control.
We firstly explore the topological features of BS structure. When SOC is not considered, the band structure (red lines) is shown in figure 3(a). Along M 1 to Γ and Γ to M 2 directions, the band intersections with linear dispersion display near the Fermi level. These band crossings actually construct a nodal line in xy-plane, protected by PT (product of P and T) symmetry. However, when SOC is included, the system transits into 2D TI phase due to the gap opening along the nodal line. For example, when the strain = −8.0%, the system is at the phase transition point, the SOC interaction gaps the Dirac points, driving the system turn into a 2D TI with a small energy gap of 0.27 meV, as shown in figure 3(a). To be specific, the Dirac carriers acquire a small finite mass proportional to half the energy gap [34]. With the further decrease of strain ( < −8.0%), the effective Dirac mass of the BS structure becomes larger and larger, keeping staying in the 2D TI phase. The results are also confirmed by WIEN2k calculations (appendix A figure A4).
Moreover, the topology hosted in BR structure is studied. When SOC is neglected, the BR structure maintains two linear band intersections along the M 1 -Γ direction in the whole strain range (−8.0% < < 13%), as shown in figures 3(b)-(d). The wave-function analysis implies the band intersections are protected by the C 2 point group [35], indicating two Dirac points are formed in the M 1 -Γ direction. However, along the Γ-M 2 direction, the band gap opens when the strain is smaller than 8.6%, due to the hybridization between two nearest Se atom. When the strain is 8.6%, there exists a band gap of 7.8 meV at Γ point [right upper panel in figure 3(c)]. In addition, when the strain is 13%, another intersection appears in the Γ-M 2 direction near the Fermi level [ figure 3(d)]. Overall, during the increase process of the strain from −8.0% to 13%, the direct band gap value gradually increases firstly and then decreases to 0 eV.
The introduction of SOC interaction opens the global gap of the BR structure, but only closes the gap at Γ point when the strain = 8.6% [figures 3(b)-(d)], implying a topological phase transition. Interestingly, when the strain = 5.0%, the global energy gap of the system achieves the maximum value of 0.16 eV, as shown in figure 3(b). In addition, when the strain further increases to 8.6%, there exists an energy band intersection at Γ point, around which the band disperses linearly, as shown in figure 3(c). Therefore, the structure is actually a DSM. The fourfold degenerate Dirac fermion is protected by C 2h point group [35], having different dispersion along M 1 -Γ and Γ-M 2 .
To further shed light on the mechanism of topological phase transition in BR structure, we explore orbital-resolved band structures with SOC interaction. The components of the band were projected to p x,y,z The red circle and green square represent the (p x + p y ) orbitals and p z orbitals, respectively. atomic orbitals, as shown in figures 3(e)-(h). From figure 3(e), we can see that lowest conduction band crosses highest valence band, forming an inverted band ordering in the vicinity of Γ point. Thus, this leads to a TI with a large Γ point inversed band gap (GPIBG) ∼ 6.5 eV. Here, (p x + p y ) states are much higher than p z states. With the increase of strain from −8.0% to 8.6% by stretching the lattice, we can see the GPIBG reduces gradually to 0 eV and a DSM appears, due to the reduction of p-p hybridization. However, when the strain further increases up to 13%, the band gap opens again, where p z states are higher than (p x + p y ) states, resulting in a positive energy gap and a simple insulator is formed.

Tunable Dirac cones and topological properties
In order to give a holistic picture of the process of topological phase transition, we plot in figure 4(a) the strain-driven phase diagram with SOC interaction in detail, including 2D TI, DSM and trivial insulator phases. Interestingly, the 2D TI phase can be further classified as four Dirac cones and two Dirac cones subphases. The blue dots in figure 4(a) show the stain dependency of direct energy gap when considering SOC. In BS structure, when the strain decreases ( < −8.0%), the direct gap linearly increases correspondingly.
When the strain = −8.0%, the energy gap approaches the minimum value ∼0.27 meV, as shown in figure 3(a). Hybrid functional HSE06 is performed and the calculations verify the results. This indicates that the topology protected four gapped Dirac cones with mass ∼0.135 meV are formed. It is obviously from figure 4(b) that the four Dirac cones located in diagonal relation are identical due to the fourfold rotation symmetry. To provide a visual of the band dispersion near the Dirac cones, one of them is enlarged and shown in figure 4(c). From the picture, we can see that the band dispersion near the Dirac cone is in the shape of wild goose and anisotropic. In addition, the Dirac cone has two different velocity branches along the k x and k y directions. The highest Fermi velocity is 1.08 × 10 6 m s −1 , which is comparable with that of graphene (∼1.0 × 10 6 m s −1 ) [36].
With the further increase of strain, the system becomes a 2D TI in BR structure, which hosts two equivalent Dirac cones due to C 2 symmetry. To prove the quantum spin Hall effect (QSHE) in BR structure, the Z 2 invariant is calculated by means of tracing the evolution of WCC in an adiabatic pumping process. The Z 2 invariant is defined as the odd number crossings of WCC (black lines) of the reference line (red line). Take the case of the strain 5.0% for example, the Z 2 invariant index (crossing number modulo 2) is calculated to be 1, as revealed in figure 4(d). This means the BR structure is indeed a TI, presenting nontrivial helical edge states. Especially, in the case of the strain 5.0%, the global energy gap of structure achieves its maximum value ∼0.16 eV among the whole TI phase, as illustrated in figure 3(b). Hybrid functional HSE06 calculations also verify the results, even with a larger band gap. Then, the edge states are performed, as presented in figure 4(e). Nontrivial edge states (red lines) emerge in the band gap, connecting the conduction bands and the valence bands. The big gap makes it feasible for the QSHE applied in room temperature. Moreover, as shown in the middle panel of figure 4(a), the direct gap of the BR structure gradually increases first and then drop to 0 eV after its peak value ( = 7.0%), when the strain changes from −8.0% to 8.6%. In other words, the strain-driven tunable Dirac mass varies from 0 eV to 0.084 eV, indicating the wide application range of the monolayer selenium in devices.
When the strain is 8.6%, the structure is at the phase transition point from 2D TI phase to simple insulator phase, and becomes a DSM. To demonstrate the topological property of the DSM, the fourfold degenerate Dirac point at Γ point is calculated and presented in figures 4(f) and (g). From the pictures, we can see away from the Γ point, the band dispersion satisfies twofold rotation symmetry, while it is anisotropic when near the Γ point. Moreover, the robustness of topology of monolayer selenium, which maintains in a wide range of the strain from −8.0% to 8.6%, implies its potential application under an abroad mechanical condition. Furthermore, when the strain is above 8.6%, the structure is a trivial insulator. The direct energy gap enlarges with the increase of strain, as shown in figure 4(a).
Finally, we notice that the studies of phonon structure and molecular dynamics simulations have verified the stability of freestanding monolayer selenium [37]. In addition, monolayer selenium with planar square structure deposited on a substrate Au (100) has been observed in experiment [38]. Also, we propose that the insulating substrate SrTiO 3 [39][40][41] can not only stabilize the strained monolayer selenium with BS lattice, but also maintains the Dirac cones (appendix A figure A5). These results indicate the potentially experimental realization of the monolayer selenium in this paper. Therefore, the tunability of number or mass of the Dirac cones in monolayer selenium provides a potential candidate for technological applications in spintronics and opens a pathway to the development of novel devices.

Conclusion
In this paper, using first-principles calculations, a monolayer 2D square selenium is proposed to achieve Dirac cones and a topological phase transition. The computational results show that the monolayer selenium displays two structural phases: the BS phase and BR phase. In the BS phase, the structure is a 2D TI having four Dirac cones with small masses. With the decrease of the strain, the structure stays in the TI phase and the mass of Dirac carriers become larger and larger. The BR phase maintains in the TI phase in the large strain range from −8.0% to 8.6%, but with two Dirac cones having tunable Dirac mass. These topological features are verified by hybrid Wannier function technique. When the strain is 5.0%, the maximum global band gap ∼0.16 eV is obtained in BR phase, enabling the QSHE to be observed at room temperature. When the strain becomes 8.6%, the structure turns into a DSM with single Dirac point located at Γ point. Then, with the further increase of strain, the structure transits into a simple insulator without Dirac cone. It is notably that in the whole strain range, the Dirac mass is continuously tuned in a controllable way. Therefore, our study reveals that the monolayer selenium has high tunability of Dirac carriers, which could be applied for potential design of novel tunable topological materials and devices.