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.

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.
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 π-stack ) and between the O= and H-terminations (d OH ) 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.
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.
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 π-stack and d OH parameters are reported in Fig. S2 in SI.
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 OH = 1.76 and d π-stack = 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 tinybelow 0.1 eV -the values approaching 0.2 eV from the DFT and 0.5 eV from the hybrid-DFT were reported 18,19 . Scientific REPORts | 7: 9888 | DOI:10.1038/s41598-017-09954-z 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 σ xy plane, containing the stacking axis and the non-periodic direction, and (2) the σ xy 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. 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 π-stack and d OH , and are reported in Fig. 6. If the molecular stacking distance is fixed at 2.3 Å, and d OH 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 OH in the range [1.56, 1.76] Å, the edge states cross the Fermi level. A change of d π-stack 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.
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 π-stack and d OH 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 π-stack = 2.3 Å and d OH = 1.76 Å, shows the localization of the edge states at the both surfaces of the 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.
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 π-stack = 2.7 Å and d OH = 1.6 Å) are presented in Fig. S5 in SI.
Interestingly, the gapless edge-states with the surface-localized orbitals were found for the intermolecular distance d π-stack = 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 CH 2 CN 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 CH 2 CN group possessed special transport properties along the stacking axis 50 . The coincidence  . 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 π-stack = 2.58 Å and d OH = 1.76 Å. 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. 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. of the results for the stacking distances suggests, that the CH 2 CN 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 π-stack 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 A ) and slightly away  (k B and k C ). 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 A point at the band gap. However, this swapping is more pronounced between k A and k C (towards the Z-point) than between k A and k B (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 W 2 HfC 2 O 2 19 . 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 val = 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 ξ i 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.

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 π-stack and d OH . 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 G 0 W 0 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 G 0 W 0 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.