Low-dimensional elemental materials are a large family of two-dimensional (2D) materials.^[1–3] Graphene^[4–8] was the first mono-layer ever isolated for carbon and 2D materials while graphdiyne^[9] is an allotrope of it. Layers of group IV elements, known as silicene,^[10–12] germanene,^[13] stanene,^[14,15] and layers of a group III element, i.e., borophenes,^[16–18] as well as group V few-layers, i.e., 2D P,^[19–23] As,^[24] Sb,^[25] and Bi,^[26] were subsequently predicted and synthesized or isolated. These mono- and few-layers of groups III, IV, and V elements were experimentally shown or were theoretically predicted to have tunable bandgap,^[7,12,23,24] high carrier mobility,^[19–21,25] strong light absorption and response in infrared and visible lights ranges,^[27,28] exceptional mechanical and frictional properties,^[5,17,22] catalysis activities,^[29] topological electronic states,^{[8,10,14,26,30–35]} superconductivity,^[36,37] and among the others.^[38] However, the few-layer forms of group VI elements are still ambiguous and are yet to be unveiled.

Tellurium few-layers are a category of emerging group VI 2D layers.^[39–49] The few-layer α phase, cleavable from its bulk counterpart, shows amazing electronic, optical, vibrational, and topological properties^[50–52] and can be synthesized using wet-chemistry methods.^[41] Therefore, the synthesis of it does not need a substrate while other 2D few-layers need either exfoliation from their bulk counterparts^{[4–6,19–21,53–55]} or substrates to stabilize^[17,56–60] the synthesis of layered samples. A striking feature of it lies in that it has, at least, four few-layer allotropes predicted by density functional theory (DFT) calculations, the number of which is comparable with that of carbon. A previous theory shows that meta-stable few-layer phases could be stabilized with charge doping.^[52] Another theory predicted topological states in Te nanostructures,^[61] however, those structures are highly unstable. It would be thus of interest to know if there are any new phases, preferably with topological states, yet to be unveiled and the reason why Te could offer so many allotropes. Answers to these questions would boost both fundamental research and device applications.

Here, we predicted two novel forms, i.e., ε and ζ phases, of Te few-layers, among which the ζ phase shows extraordinary stability that its monolayer is 29 meV/Te more stable than the most-stable γ monolayer and its bilayer is over 30 meV/Te more stable than the most-stable α bilayer. An energetic crossover between the ζ and α phases occurs at the four-layer (12-sublayers) thickness that the ζ phase is prone to transform into the α phase beyond that thickness, while either hole or electron doping stabilizes the ζ phase and pushes the crossover to thicker layers. The ε phase is less stable than the ζ phase, but has a comparable stability with the monolayer γ or few-layer α phase. Phonon dispersion calculations suggest that the free-standing forms both phases are stable and could be exfoliated from thicker layers or substrates. These two novel phases strongly promote subsequent studies on polytypism of Te few-layers and add two new members to the family of Te allotropes.

Density functional theory calculations were performed using the generalized gradient approximation for the exchange–correlation potential, the projector augmented wave method,^[62,63] and a plane-wave basis set as implemented in the Vienna ab initio simulation package (VASP)^[64] and Quantum Espresso (QE).^[65] Density functional perturbation theory was employed to calculate phonon-related properties, including vibrational frequencies at the gamma point (VASP) and phonon spectrum (QE). The kinetic energy cut-off for the plane-wave basis set was set to 700 eV for geometric and vibrational properties and 300 eV for electronic structures calculations. A k-mesh of 15 × 11 × 1 was adopted to sample the first Brillouin zone of the conventional unit cell of few-layer Te (α phase) in all calculations. The mesh density of k points was kept fixed when calculating the properties for few-layer Te. A q-mesh of 6 × 6 × 1 was used for phonon spectrum calculations. In optimization of geometry, van der Waals interactions were considered at the vdW-DF^[66,67] level with the optB88 exchange functional (optB88-vdW),^[68–70] which was found to show appropriate electronic structures revealed with the Heyd–Scuseria–Ernzerhof (HSE06) hybrid functional.^[71] For all calculation of the phonon spectrum, we used the optB86b exchange functional, which was proved to be accurate in describing the structural properties of layered materials.^[21,72–74] The shape and volume of each supercell were fully optimized and all atoms in the supercell were allowed to relax until the residual force per atom was less than 1 × 10⁻⁴ eV/Å. Electronic bandstructures were calculated using the PBE and HSE06 functionals with and without spin–orbit coupling (SOC). The SOC effect was solely considered on the atomic basis which does not account for orbital contributions from formed electronic bands.

Charge doping on Te atoms was realized with the ionic potential method,^[75] which was used to model the chare transfer from graphite substrates. For electron doping, electrons are removed from a 4d core level of Te and placed into the lowest unoccupied band. For hole doping, electrons were removed from the valence band by adding a negative potential into the 4d core level of those three Te atoms. This method ensures the doped charges being located around Te atoms. It also keeps the neutrality of the whole supercell without introducing background charge, which eliminates effects of compensating charges. This method was widely used in the literature.^[52,76]

The α-phase, comprised of helical chains bonded with covalent-like-quasi-bonds (CLQB) along inter-chain directions, is the most stable form in Te few-layers and bulk^[46,77] (Fig. 1(a)). In monolayers, however, the γ-phase (Fig. 1(c)) was believed the most stable phase and the α-phase becomes unstable and spontaneously transforms into the β-phase (Fig. 1(d)). The γ-phase contains rhomboid chains along two directions forming a network in the C3_v symmetry while the β monolayer is comprised of parallel rhomboid chains with an inclination angle of 29.3° to the xy-plane. Strong charge doping increases the inclination angle of the rhomboid chains in the β-monolayer and the angle eventually reaches 90° at a doping level of 0.50 e/Te, giving rise to a new phase ε (Fig. 1(e)). This phase has parallel aligned diamond chains, in which Te–Te distance is 3.04 Å and angles are 87° and 93°, respectively, while the interlayer CLQB length is 3.59 Å. The ε-monolayer is 6 meV/Te less stable than the γ-monolayer (Fig. 1(g)), when comparing surface energies of these two phases, the relative stability reversed (23 meV/Ų). This is due to the much smaller surface area of the ε-monolayer resulted from its perpendicularly tilted rhomboid chains.

An even more stable ζ phase (Fig. 1(f)) was found by relaxing atomic coordinates from laterally shifted ε layers. The ζ monolayer consists of three sublayers where the Te atoms form a square lattice, with Te–Te bond length of 3.15 Å, in each sublayer, leading to a structure with an ultra-low specific area and a high symmetry (P4/MMM). A similar, but strongly tilted, structure was previously found in bulk Te under a high pressure of over 8 GPa,^[78] which transforms into the ζ bulk phase after the pressure is fast released (see Supporting information: Fig. S1). Although the ζ bulk form is less stable than the α bulk phase, the ζ monolayer and bilayer are much stable with at least 29 meV/Te (93 meV/Ų) and 35 meV/Te (200 meV/Ų) energy gains from the previously believed most-stable γ monolayer and α bilayer, respectively. Phonon dispersion calculations confirm their thermal stability although a likely charge density wave transition was evidenced for the ζ phase at lower temperatures. It would be thus interesting to find the transition boundary of energetic stability between the α and γ multilayers.

The ζ few-layers prefer an AA stacking by at least 9 meV/Te, in which a Te atom of an upper sublayer sits right over another Te atom underneath (see Supporting information: Fig. S2 and Table S1). We thus adopted the AA stacking in following calculations. Figure 1(g) plots the total energies of the six known phases as a function of the number of sublayers, which shows that the ζ few-layer (blue square) is energetically more stable than other five phases before the thickness reaches 12 sublayers (four layers). Beyond this thickness, the structure of the ζ phase still holds but the α phase becomes the most stable phase; this is, most likely, ascribed to weakened surface effects as the bulk properties dominate the behavior of the ζ phase in thicker layers. We also plotted surface energies in Fig. 1(h). It shows the α phase is the easiest one to cleave and the β phase has a comparable surface energy. Other phases, except the ε phase, show slightly higher but reasonable surface energies.

The intra-layer bond lengths (lattice constant a/b) are 3.02 Å and 3.08 Å in a mono- and bi-atomically thick ζ sublayers, respectively, which are much smaller that the bulk value of 3.21 Å. Figure 2(a) shows the evolution of intra-sublayer and inter-sublayer bond lengths as a function of layer thickness. The increased thickness significantly varies both the inter-sublayer and the intra-sublayer distances, indicating a strong inter-sublayer interaction. The more the sublayers stacked together, the stronger the charge transfers from pz orbitals of Te atoms to in-plane p_x / p_y orbitals and intra-sublayer regions, leading to undercut intra-sublayer and reinforced inter-sublayer Te–Te bonds. The intra-sublayer lattice constant (blue rectangular), as a result, expands 3.02 Å (1-sublayer) to the bulk value of 3.21 Å at 12-sublayer while the inter-sublayer distance shrinks from 3.38 Å (2-sublayer) to 3.21 Å also at 12-sublayer. Both intra- and inter-sublayer bond lengths converge to the bulk values (3.21 Å) at 12 sublayers, consistent with energetic crossover and the order of stability of the α and ζ bulk forms.

It is exceptional that structural relaxations were found in ζ multilayers that they are prone to form dimers or trimers with adjacent sublayers along the interlayer z direction. We used the bulk bond length of 3.21 Å as a criterion. The sublayers with bond lengths smaller than this value were regarded as dimerized and trimerized sublayers. Figure 2(b) presented the detailed distributions of dimers or trimers from mono- to 12-sublayers. Sublayers dimerizing or trimerizing together are marked by red dotted rectangles and the directions of atomic relaxations are indicated with black arrows. A trimer first appears in the tri-sublayer (a ζ monolayer), which is, most likely, due to a Fermi surface nesting induced electronic structure and geometry instability. Dimers, trimers and their mixtures emerge in thicker ζ sublayers with the thickness up to 12 sublayers. We tested different combinations of the dimers and trimers confirming the configurations shown in Fig. 2(b) are the most stable ones (see Supporting information: Fig. S3). All of them show mirror and central inversion symmetries along the inter-sublayer direction. Besides that, dimers would not show up at the surface region, which is consistent with the non-dimeric 2-sublayer ζ. The reason why these relaxations occurs is another research topic than we will discuss it elsewhere.

Figures 3(a) and 3(b) show the bandstructures of the ζ tri-sublayer calculated using the PBE functional without and with SOC, respectively. We found several band inversions as confirmed by the orbital decomposed band structures shown in Fig. 3(a). The inversion occurred around the G point forming a nodal ring is of particular interest that the inversion point sits roughly at the Fermi level. Inclusion of SOC opens a bandgap of 0.38 eV around the Fermi level. We thus calculated the Z₂ topological invariant using Quantum Espresso (QE) to verify the topological characteristic of the ζ tri-sublayer. The Te atom in the ζ tri-sublayer is in a square network structure, which has both time reversal and space central inversion symmetries. Therefore, the Z₂ topological invariant can be obtained by multiplying parities of filled states at all time-reversal invariant points, as shown in Fig. 3(c). Our calculation revealed the Z₂ value of (–1), which indicates the ζ tri-sublayer to be with nontrivial characteristic. However as present in Fig. 3(d), the surface states were overcovered by the bulk states. Details of the SOC induced bandgap opening and inversions were available in Supporting information: Figs. S4 and S5. In addition, we also found that the monolayer ε and 1-sublayer ζ phases are topologically trivial, as summarized in Supporting information: Tables S2 and S3. Quantum well states were explicitly observed for the states along the z direction, the direction normal to the layer planes, as shown in Fig. S6 where shows the evolution of the band structures of ζ few-layers with different thicknesses.

In a recent work,^[52] we found direct charge doping could change the relative stability of Te phases and transform a certain phase to another. We thus considered the stability of both the ε and ζ phases under electron or hole doping. Figure 4(a) shows the total energies of these two phases, γ and β phases in Te monolayers, in which the most stable ζ phase was chosen as the reference zero. The diagram shows that the ζ phase keeps its exceptional stability, with at least 30 meV/Te to the β phase, under either electron or hole doping. This statement was double confirmed by the SCAN + rVV10 functional calculations (see Fig. S7). Further calculations show the ζ monolayer should be the most stable phase even when the doping level reaches 0.8 e/Te or 0.8 h/Te. Such values are beyond the capability of modern ionic liquid gating techniques, which, together with vibrational frequency calculations, guarantee the superior stability of the ζ monolayer. These results indicate that the ζ monolayer has little chance transforming into other forms of low-dimensional Te, if it is fabricated. The ε monolayer is slightly less stable than the γ phase at the neutral state, while it becomes more stable when adding or removing electrons, which is substantially different from the previous found δ and mixed phases.^[52] This structural phase transition manipulates the semiconducting β or γ phase transforming into the metallic ε phase. The detailed discussion for the electronic properties of ε phase and the origin of β–ε phase transition in the Te monolayer are available in Fig. S8. Figure 4(b) shows the relative energies of different phases in 2L. It, again, shows the superior stability of the ζ phase and the more stable ε phase than the γ phase under doping. These results suggest doping could further stabilize the ζ phase, which pushes the ζ–α crossover thickness beyond 4L under charge doping as shown in Fig. 4(c).

In summary, we predicted two new low-dimensional Te allotropes, i.e., ε and ζ, which, especially the ζ phase, yield extraordinary stability. It has strong vertical inter-sub-layer interaction that shows quantum well states along the direction normal to the layers. As the most stable few-layer phase found so far, it was surprising that this phase has not been synthesized yet; this is, most likely, due to its substantially different geometry from the helical bulk-like form or the lack of a square substrate lattice. We expected that the ζ phase might be prepared by molecular beam epitaxy, physical vapor deposition, laser or electron beam deposition or even chemical vapor deposition with precisely controlled dosing rates, temperatures, substrates or from a fast released high pressure phase. Unlike semiconducting α, β, and γ layers, the δ and ε chains and the ζ phase are metallic with high and tunable density of states and strong band dispersions. We identified a weaker electronic interaction, with a typical interacting distance of roughly 3.2 Å–3.3 Å (marked by red lines in Supporting information: Fig. S10) in the metallic δ, ε, and ζ phases. Such distance is ∼ 0.3 Å longer than the lengths (2.8 Å–3.1 Å) of typical covalent bonds found in Te allotropes. The longer distances weaken the attraction to electrons from Te cores and thus lead to more delocalized Te p electrons between the Te chains in the δ and ε phases or among those individual Te atoms in the ζ phase, which is believed to result in those highly dispersive states in those metallic phases. In addition, those metallic phases are ideal for applications of layered electrodes. Our results added two more allotropes to few-layer Te and open a new avenue for studying topological properties in group VI 2D layers.

Acknowledgment. Calculations were performed at the Physics Laboratory of High-Performance Computing of Renmin University of China and the Shanghai Supercomputer Center.