Gapless edge states in (C,O,H)-built molecular system with p-stacking and hydrogen bonds.

The gapless edge states have been found in a 2D molecular system built with light atoms: C,O,H. This prediction is done on the basis of combined density functional theory (DFT) and tight-binding calculations. The system does not exhibit any effect of the spin-orbit coupling (SOC), neither intrinsic nor Rashba type. The band structure and the edge states are tuned with a strength of the p-stacking and O...H interactions. The elementary cell of this noncovalent structure, does not have the 3D inversion or rotational symmetry. Instead, the system transforms via a superposition of two reflections: with respect to the xz and xy mirror planes, both containing the non-periodic direction. This superposition is equivalent to the inversion in the 2D subspace, in which the system is periodic. The energy gap obtained with the DFT method is 0.11 eV, and largely opens (above 1 eV) with the GW and hybrid-DFT approaches. The bands inversion is partial, i.e. the Bloch states are mixed, with the "swapping" and "non-swapping" atomic contributions.

Introduction

Topological insulators (TIs) attracted much attention, starting from the work by Kane and Mele[1] in 2005 to the Nobel Prize for David J. Thouless, F. Duncan M. Haldane and J. Micheal Kosterlitz in 2016. This is due to many plausible applications of their metallic surface states, with the linear dispersions, propagating electrons without the elastic scattering[2]. First experimentally found TIs were HgTe/CdTe[3] and InAs/GaSb[4] quantum wells, followed by the 3D crystals with the time-reversal symmetry, such as Bi2Sb3, Bi2Te3 and Sb2Te3 [5, 6]. Further, the Heusler compounds containing rare-earth elements[7], (Tl, Hg, Pb, Cu, Ag)-chalcogenides[8, 9], Zintl compounds[10], antiperovskites[11] and strongly correlated oxides[12] have been found. All contained heavy elements, and some of them were under pressure[9]. The spin-orbit coupling in these materials is a key-parameter, because it is responsible for the bands inversion and opening of the gap[13]. The 2D TIs, such as a single and multilayer Bi[14, 15], silicene on the silver surface[16], electrically controlled multilayer germanene[17], oxygen-functionalized MXenes[18, 19], arsenene monolayers[20], SiTe[21] and porous MoS2 [22] were discovered. Then, a few organic networks were theoretically predicted to possess the topological edge states. Their list contains: triphenyl-lead[23], s-triazines[24], nickel-bis-dithiolene[25] honeycomb lattices, and ligand complexes of Pd and Pt[26].

In 2011, a new class of TIs was theoretically postulated, namely topological crystalline insulators (TCIs)[27]. They were supposed to be built with light atoms and not show any effect of SOC. These materials have the doubly degenerate (due to spin) surface states with quadratic dispersions. They are protected by time-reversal and discrete rotational symmetry[27-29] (e.g. C 6 and C 4 in 3D). The electron’s orbital degrees of freedom should play there a role similar to spin. In 2014, the group of TCIs was theoretically extended to the 3D systems without the time-reversal symmetry, but possessing the rotational symmetries with the mirror planes (i.e. C for n = 3, 4, 6, and not 2)[30].

Early reports about TCIs, based on numerical simulations, predicted the gapless edge states in PbSe, PbTe, PbS[31] and pyrochlore oxides A2Ir2O7 (A = rare-earth element)[32]. The strong ionic character and SOC were responsible for the bands inversion in these materials. Later, SnSe and SnS were theoretically proposed as TCIs[33]. Although the SOC was still present, it did not play a main role in the topological mechanism of the bands inversion, but it just opened the gap. The experimental reports confirmed the TCI phase in SnTe[34, 35] and Pb1−SnSe/Te alloys[36-38]. The 2D TCIs have been predicted in TlS and TlSe[39], SnTe[40], PbSe[41], containing heavy atoms. The 3D heterostructure of 2D materials PbSe/h-BN was postulated, on the base of the mirror symmetry consideretions, to possess the TCI phase[42]. It was theoretically suggested, that the topological order and Rashba SOC can be tuned with the hydrostatic pressure, epitaxial strain, and additionally controlled with the electric field, in the antiferroelectric ABC-compounds[43]. In contrast, up to date, there is probably only one theoretical report on the 2D TCI system with a negligible SOC, namely the twisted graphene multilayers[44].

In this work, the gapless edge states were found in a 2D molecular crystal composed of phenalene derivative. The chemical structure of this molecule is presented in Fig. 1. It is not aromatic, because a complete termination of one C-ring with oxygens fixes a pattern of the double bonds. Similarity of the studied molecule to 1H-phenalene-1,3(2H)-dione[45] makes an expectation that it could be synthesized.

Figure 1

The chemical structure of the studied phenalene derivative.

In a plausible TCI reported here, the molecules are π-stacked with a given order. Their columns are periodically repeated in such a way, that the hydrogen bonds form between them along one of the planar directions. The edge states are tuned with two intermolecular distances: between the molecular planes (d ) and between the O= and H-terminations (d ) of the neighbouring molecules.

Results

Topological band effects show up when non-aromatic molecules are π-stacked in a given geometrical order. By applying a pressure along the molecular columns, one can pass from the insulating to metallic state. In some cases the gapless edge states, between two “normal” phases, appear.

Figure 2 shows five orders of the 1D stacking: “on-top”, “rotated”, “rotated-shifted”, “flipped” and “flipped-shifted”. The corresponding band structures of the wires, at three chosen intermolecular distances, are also displayed. It is interesting to note, that for two orders, namely “on-top” and “rotated”, the metallic phase forms under pressure. While for the third one, called the “rotated-shifted” structure, two different “normal” insulating phases are separated by the “topological” phase. These states are followed by a metal phase with the negative gap, at a very high pressure. The band structure of the “flipped” order has a semiconducting character. While for the “flipped-shifted” order, the semiconducting phase transforms into the metallic phase at quite low pressure. Thus, further in this work, only the “rotated-shifted” order is examined, as a potential topological insulator.

Figure 2

The geometries of five molecular wires, called accordingly to the stacking patterns: “on-top”, “rotated”, “rotated-shifted”, “flipped” and “flipped-shifted”, as well as their band structures, for three intermolecular stacking distances. The pressures are given above the corresponding panels. The Fermi level is at zero energy.

The band gap obtained with the DFT method for the molecular wire with the “rotated-shifted” order at  the stacking distance of 2.3 Å is 0.11 eV. This fundamental gap largely opens to 0.4 or 1.16 eV with the GW approach, depending on a polarization of the electric field - parallel or perpendicular to the wire, respectively. Strongly anisotropic GW results for the molecular wires are not surprising[46]. Figure 3 presents the band structures.

Figure 3

Band structures of a wire with the “rotated-shifted” molecular order, for a chosen intermolecular distance of 2.3 Å. The results were obtained with the DFT method and the GW approach for two electric-field polarizations: perpendicular and parallel to the axis of the π-stack. The calculated points are marked with circles and the band-lines are interpolated.

The 2D structure is constructed via a repetition of the π-stacks along the chosen planar direction - the one which forms the O…H bonds. This is displayed in Fig. S1 in the supporting information (SI). The band structure, plotted for the intermolecular distances of 2.58 and 1.76 Å in the stacking and planar O…H-bond directions, respectively, is presented in Fig. 4. While changes of the band structure with a variation of the d and d parameters are reported in Fig. S2 in SI.

Figure 4

The geometric parameters of the 2D structure formed by a π-stacked wires with the “rotated-shifted” molecular order. The planar intermolecular hydrogen-bonds connect the columns. The band structure was obtained with the DFT method for d  = 2.58 Å and d  = 1.76 Å.

In order to take into account an effect of the exact exchange on the band structure of the 2D case, the hybrid-DFT calculations were performed for the intermolecular distances d  = 1.76 and d  = 2.3 Å. The fundamental gap obtained with this method is 1.2 eV. A comparison of the hybrid-DFT energy gaps at the symmetry points with the DFT band structure is presented in Fig. S3 in SI.

Most of gaps in TIs open due to the SOC effect. As one would expect, for very light atomic systems, an effect of the spin-orbit coupling is vanishing. The lack of the SOC in the system studied here has been checked with the fully relativistic DFT method. It was necessary, since the low-dimensional systems might be surprising. For example, graphene weakly doped with hydrogen opens the gap to 10−2 eV, which is a few orders of magnitude larger than that of pure graphene, 10−6 eV[47, 48]. Although the SOC bandgaps in most of the TIs are rather tiny - below 0.1 eV - the values approaching 0.2 eV from the DFT and 0.5 eV from the hybrid-DFT were reported[18, 19].

The 2D system with “rotated-shifted” molecular order does not possess the rotational symmetry or the inversion in the 3D space. Instead, it transforms with a superposition of two reflections: with respect to (1) the σ plane, containing the stacking axis and the non-periodic direction, and (2) the σ plane, which separates the molecules across the π-staking axis. The mirror planes are depicted in Fig. 5. This superposition of reflections is equivalent to the inversion in the 2D subspace spanning the periodic directions of the system.

Figure 5

Two mirror planes in the 2D structure formed by the π-stacked wires with the “rotated-shifted” molecular order, which are repeated along the hydrogen intermolecular bonds.

For the finite system in the stacking direction, the topological crystalline insulating phase is supposed to appear, due to the gapless edge states. They are doubly degenerate with spin components, exhibit the parabolic dispersions, and propagate along the O…H intermolecular bonds. These edge states are tunable with the intermolecular distances, d and d , and are reported in Fig. 6. If the molecular stacking distance is fixed at 2.3 Å, and d varies between 1.36 and 1.96 Å, the edge states touch the Fermi line from the lower energies at the Γ-point for short O…H, or from the higher energies at the π-point for larger distances. While for d in the range [1.56, 1.76] Å, the edge states cross the Fermi level. A change of d to smaller or larger values does not move the edge states away from the Fermi level. In Fig. 6, the finite size along the π-stack was set to fifty elementary cells. An evolution of the edge states with a thickness of the stack is presented in Fig. S4 in SI. The pressures corresponding to the studied intermolecular distances are collected in Table 1.

Figure 6

The edge states in the 2D structure presented in Fig. 4. The thickness of the π-stacked layers was 50 elementary cells along the z-axis. The y-axis along O-H bonds was treated periodically. The intermolecular distances are given above the corresponding plots.

Table 1

Pressures (in GPa), obtained with the DFT method for the 2D structure presented in Fig. 4.

d π-stack	P π-stack	d OH	P OH
2.15	13.7	1.36	2.5
0.30	8.5	1.56	0.3
0.45	5.2	1.60	0.0
0.58	3.3	1.76	−0.9
0.70	1.8	1.96	−1.7

The d distance (in Å was varied when the d distance was fixed to 1.76 Å. The dependence on the d distance (in Å) was obtained with the fixed d  = 2.3 Å.

The localization of the wavefunctions of the edge states can be easily visualized using the Wannier-function[49] based tight-binding parametrization of the Hamiltonian. In this approach, the periodic problem is solved with constraints for a finite dimension - in this case along the π-stack axis. The largest square components of the chosen eigenfunctions are plotted in Fig. 7, for the same d and d distances which are reported in Fig. 6. There is no exact coincidence between the gapless edge states and the localization at the surface or in the interior. Only one case, d  = 2.3 Å and d  = 1.76 Å, shows the localization of the edge states at the both surfaces of the finite stack. For most of the cases, the edge orbitals localize at one side of the stack. Only for very high- or very low-pressure cases, the edge-state orbitals are delocalized in the interior. For a sake of completeness, the bands, edge states and edge orbitals at ambient pressure (d  = 2.7 Å and d  = 1.6 Å) are presented in Fig. S5 in SI.

Figure 7

Visualization of the highest-energy occupied eigenvectors of the tight-binding Hamiltonian of the 2D structure presented in Fig. 4. The π-stacking and O-H distances were varied. The matrix elements involving the Wannier functions were included when .

Interestingly, the gapless edge-states with the surface-localized orbitals were found for the intermolecular distance d  = 2.3 Å. And this stacking distance, taken twice and equal to 4.6 Å, was obtained during an optimization of the ferroelectric organic solar cell[50]. The mentioned system was built by benzene rings terminated with the COOH and CH2CN dipole groups. These dipoles - of about 0.07 and 0.6 Debye, respectively - tend to orient ferroelectrically in the stacking direction. The COOH group was used to connect molecules in the planar directions, while the CH2CN group possessed special transport properties along the stacking axis[50]. The coincidence of the results for the stacking distances suggests, that the CH2CN dipole groups could be utilized for a chemical stabilization of the TCI proposed here - if we use them for a connection of every second molecule.

Since the DFT is not a trustworthy method for the orbital order and localization, the pseudopotential self-interaction correction approach[51-53] (pSIC) has been applied. The pSIC scheme, by its construction, is similar to the LDA + U method. The difference is that all atomic shells (not only d- or f-type) are corrected, and there is no ‘ad hoc’ parameter such as Hubbard-U. This parameter-free method is an easier alternative to the hybrid-DFT approach, where an amount of the exact exchange can be varied with respect to the applied density-functional exchange. The band structure, edge states and orbitals obtained with the pSIC are presented in Fig. S6 in SI, for the “rotated-shifted” case. Generally, the results qualitatively agree with these from the DFT. However, the pSIC edge states obtained at chosen d parameters seem to be more correlated with the corresponding DFT results at smaller distances.

In order to check the bands inversion property, firstly, the elementary cell is divided into “upper” (M1) and “lower” (M2) molecular components. Figure S7 in SI shows the bands projected at the Wannier functions[54] localized at M1 and M2. The highest occupied Bloch states (HOBS), through whole BZ except the Z-M line, and the lowest unoccupied Bloch states (LUBS) are composed of mixed M1- and M2-centered states in 1:1 ratio. Additionally, in Fig. 8, a composition of HOBS and LUBS with respect to the atomic contributions is mapped at two slices, just below M1 and above M2. This is done for three k-points: at the bandgap (k ) and slightly away (k and k ). The scheleton of the M1 molecule is used as a reference (a guide for an eye) for projections at both surfaces. These plots demonstrate only partial bands inversion. The characteristic swap, between HOBS and LUBS for some atomic components, occurs on both sides of the k point at the band gap. However, this swapping is more pronounced between k and k (towards the Z-point) than between k and k (towards the G-point). Such incomplete, and not shaply marked in BZ, bands inversion is novel among known TIs and TCIs. In almost all topological insulators, the exchange of the bands character takes place between the topological k-points in BZ. However, there is an example in which bands inversion occurs exactly at the topological point and not away from it, namely W2HfC2O2 [19].

Figure 8

The case with gapless edge states: d  = 2.3 Å and d  = 1.76 Å. Squares of the highest occupied Bloch state (HOBS) and lowest unoccupied Bloch state (LUBS), at three chosen k-points: k = 0.3(3)π, k = 0.3 π, k = 0.36(6)π, are plotted for two chosen 2D slices - just below the upper molecule (M1) and above the lower molecule (M2) in the elementary cell. All plots are done with the same scale for the atomic coefficients. As a reference, a scheleton of the upper molecule is displayed in all panels.

Further, an analysis of the symmetry properties of the Bloch functions at the Γ-point and three TRIM (time-reversal invariant momentum) points is done. The manifold of the entangled valence bands is composed of N  = 50 states. However, at high pressures, the six lower bands merge with the upper fifty states. The energy gap which separates the bottom of the valence band manifold, for all studied cases, is presented in Fig. S8 in SI. The parities of all valence states and the lowest unoccupied state, for all studied cases, are collected in Tables T1 and T2, and T3 in SI. The product of Bloch’s parities, , is also given. It indicates that the 2D system with “rotated-shifted” order might be a weak topological crystalline insulator at some small pressures. The conclusive data are collected in Table 2.

Table 2

Characterization of the valence band parities at the Γ-point and three TRIM points: Z = (0, 0, 1/2), Y = (0, 1/2, 0), M = (0, 1/2, 1/2), for the 2D structures with various [d , d ] distances (in Å).

Case	N val	δ(Γ)	δ(Z)	δ(Y)	δ(M)	H/L(Γ)	H/L(Z)	H/L(Y)	H/L(M)	gapless edge st.	localization of edge st.
DFT
2.3, 1.36	50	+	−	+	−	−/−	−/+	−/−	−/+	No	one side
0.3, 1.56	50	+	−	+	−	−/−	−/−	−/−	−/−	Yes	one side
0.3, 1.76	50	+	−	+	−	−/−	−/−	−/−	−/+	Yes	both sides
0.3, 1.96	50	+	−	+	−	−/−	−/+	−/−	−/+	No	one side
2.15, 1.76	56	+	+	+	+	−/−	−/−	−/−	−/−	Yes	delocalized
0.45, 1.76	50	+	−	+	−	−/−	−/+	−/−	−/+	Yes	one side
0.58, 1.76	50	+	−	+	−	−/−	−/+	−/−	−/+	Yes	delocalized
2.70, 1.60	50	+	−	+	−	−/−	+/+	−/−	−/+	Yes	delocalized
pSIC
2.30, 1.76	56	+	+	+	+	−/−	−/+	−/−	+/+	Yes	both sides
0.58, 1.76	50	++	−	+	−	−/−	+/+	−/−	+/+	Yes	one side
0.70, 1.60	50	+	−	+	−	−/−	−/+	−/−	−/+	Yes	delocalized

N is a number of the entangled valence bands, δ is a product of the Bloch’s parities within the valence-band manifold, H/L means HOBS/LUBS. The rightmost two columns indicate whether the edge states cut the Fermi level or not, and report the localization of these states at the surface or in the interior.

Conclusion

In this work, a new plausible topological crystalline insulator has been found via the DFT-based tight-binding simulations. This is a 2D molecular system with π-stacks in one dimension and hydrogen bonds in the other. It is composed of light elements: C,O,H. The gapless edge states, with the quadratic dispersions, propagate through the hydrogen bonds within the molecular planes, for the finite stacking size. The character of these states - the parity and localization - are tunable with the intermolecular distances d and d . There is completely no effect of the spin-orbit interaction, neither intrinsic nor Rashba type. Only partial bands inversion occurs, and it is not sharply marked in a well defined BZ region. This makes the proposed system to be unique among known or predicted up-to-date TCIs.

Methods

The density functional theory[55, 56] calculations were performed employing the Quantum ESPRESSO (QE) suit of codes[57]. It is based on the plane waves and the pseudopotentials describing the atomic cores. The normconserving pseudopotentials were chosen in this work. The energy cutoff of 60 Ry for the plane waves was enough to converge the band structures. The uniform k-mesh of Monkhorst and Pack[58] with the 12-point grid along the stack and 8-point grid along the OH intermolecular bonds was sufficient for the self-consistent calculations. The exchange-correlation functional in the gradient corrected Perdew-Burke-Ernzerhof parametrization was chosen[59]. The vacuum space between the 2D slabs in the non-periodic direction was set to 30 Å.

In order to accurately interpolate the band structures, the wannier90 package[60] was used. It enables finding the maximally-localized Wannier functions for the composite bands[49, 54]. The Kohn-Sham Hamiltonian[55, 56] between the Wannier functions[54], which were obtained from the DFT Bloch functions, was used as a tight-binding model for the edge-states analysis. The cut-off for the nonvanishing matrix elements was set to 0.005 eV. For the tight-binding simulations, the pythTB[61] code was utilized. The parities were obtained with the bands.x utility from the QE package.

The GW calculations were performed with the Yambo code[62] in the non self-consistent G0W0 scheme[63, 64]. The frequency-dependent electronic screening was calculated within the plasmon-pole approximation[65]. A cutoff of 12 Ry was used for the exchange self-energy. The dielectric matrix was constructed summing up to 800 unoccupied bands and using a cutoff of 6 Ry. Up to 1200 unoccupied bands were taken for the G0W0 summation. For the hybrid-DFT calculations, within the PBE0 flavour[66], the uniform mesh of the q-vectors with the 3-point grid along the stack and 2-point grid along the OH-bonds was used. This is a commensurate subset of the k-points used in the DFT calculations.

Supplementary information