Giant edge state splitting at atomically precise graphene zigzag edges

Zigzag edges of graphene nanostructures host localized electronic states that are predicted to be spin-polarized. However, these edge states are highly susceptible to edge roughness and interaction with a supporting substrate, complicating the study of their intrinsic electronic and magnetic structure. Here, we focus on atomically precise graphene nanoribbons whose two short zigzag edges host exactly one localized electron each. Using the tip of a scanning tunnelling microscope, the graphene nanoribbons are transferred from the metallic growth substrate onto insulating islands of NaCl in order to decouple their electronic structure from the metal. The absence of charge transfer and hybridization with the substrate is confirmed by scanning tunnelling spectroscopy, which reveals a pair of occupied/unoccupied edge states. Their large energy splitting of 1.9 eV is in accordance with ab initio many-body perturbation theory calculations and reflects the dominant role of electron–electron interactions in these localized states.

GNRs. The circle area is proportional to the electron density, while grey/black indicates the sign of the wave function. c, Highest occupied molecular orbital (HOMO, "bonding") and lowest unoccupied molecular orbital (LUMO, "antibonding") of (7,12) GNR.  obtained from the PBE Kohn-Sham gap and from the G 0 W 0 gap. GNR length is specified by the number of zigzag lines n, as defined in Figure 1a of the main text.

Supplementary Note 1: Transfer of GNRs onto NaCl monolayer islands
Individual (7,n) GNRs are physisorbed on Au(111), where they can be moved laterally by an STM tip.
However, direct pushing of ribbons onto NaCl monolayer islands was found impossible. The apparent height of NaCl (2.2 Å) is higher than that of GNRs (1.8 Å) and, moreover, the NaCl islands themselves and (4) lateral positioning of the GNR fully onto the NaCl monolayer.
Step 3 was introduced because releasing a GNR from the tip was found much easier when one end of the GNR is still adsorbed on Au(111). Transfer of GNRs onto bilayer NaCl has also been attempted, but was found to be much more challenging due to a lack of stable adsorption of the GNRs on NaCl bilayers in step 3.

Supplementary Note 2: Finite and infinite zigzag edges within tight binding
Single-orbital nearest-neighbor tight binding is a simple, yet instructive model that provides qualitative insights into graphene's π-electronic structure 3 . Within this model, the edge states of the (7,∞) GNR directly derive from those found at infinite zigzag edges, as is discussed in the following.
The corresponding Hamiltonian is given by For finite (7,n) GNRs, the edge states of the two termini overlap. This gives rise to a finite energy splitting ∆ !! between a bonding and an anti-bonding linear combination of edge states. Supplementary   Figure 5c depicts these states for the (7,12) GNR. In the length range ≥ 12 investigated experimentally, the weak overlap between states at the two termini gives rise to negligible splittings ∆ !! < 5 meV.

Supplementary Note 3: Computational details
Electronic structure calculations were performed within the framework of density functional theory (DFT), using the PBE generalized-gradient approximation to the exchange-correlation functional as implemented in the Quantum ESPRESSO package 5 . The Kohn-Sham orbitals were expanded on a plane-wave basis set with energy cutoff at 150 Ry, using norm-conserving pseudopotentials.
For the finite (7,n) GNRs, the simulation cell was chosen at least twice as large as the 10 -5 /a 0 3 isosurface of the charge density in order to enable the Coulomb-cutoff technique 6 in the GW calculations.
For the (7, ∞) GNR, 18 Å of vacuum were introduced in the directions perpendicular to the GNR axis, while 16 k-points were used to sample the first Brillouin zone. The lattice parameter of the (7,∞) GNR was determined to be 4.285 Å and atomic positions were relaxed until the forces acting on the atoms were below 3 meV/Å.
Quasi-particle corrections were computed within the framework of many-body perturbation theory, using the G 0 W 0 approximation to the self-energy as implemented in the BerkeleyGW package 7,8 . The electronic structure from DFT was recalculated using 60 Ry plane-wave cutoff (and 64 k-points in the first Brillouin zone for the (7, ∞) GNR). We have computed sufficient numbers of empty states to cover the energy range up to 2.1 Ry (infinite GNR) and 1.6 Ry (finite GNRs) above the highest occupied band. The static dielectric matrix ε was calculated in the random phase approximation with 8 Ry cutoff for the plane-wave basis. ε !! was extended to the real frequency axis using the generalized plasmonpole model by Hybertsen and Louie 7 . A rectangular Coulomb-cutoff was employed along the aperiodic dimensions as described in reference 6 . In the calculation of the self-energy, the static remainder approach was used to speed up the convergence with respect to the number of empty bands 9 .
With regard to the magnetic structure, the PBE functional predicts ultrashort GNRs of lengths n=2 (anthracene) and n=4 (bisanthene) to have a spin-unpolarized ground state. From length n=6 onwards, PBE finds a ground state with spin-polarized, antiferromagnetically coupled edge states. The energy difference between the antiferromagnetic (spin-singlet) and the ferromagnetic (spin-triplet) solution, however, decays exponentially with GNR length, reaching ≈1 meV already for n=12, the shortest GNR length investigated experimentally. Similar findings have been obtained previously by the complete active space self-consistent field (CASSCF) method