Disentangling the electronic structure of an adsorbed graphene nanoring by scanning tunneling microscopy

Graphene nanorings are promising model structures to realize persistent ring currents and Aharonov–Bohm effect at the single molecular level. To investigate such intriguing effects, precise molecular characterization is crucial. Here, we combine low-temperature scanning tunneling imaging and spectroscopy with CO functionalized tips and algorithmic data analysis to investigate the electronic structure of the molecular cycloarene C108 (graphene nanoring) adsorbed on a Au(111) surface. We demonstrate that CO functionalized tips enhance the visibility of molecular resonances, both in differential conductance spectra and in real-space topographic images. Comparing our experimental data with ab-initio density functional theory reveals a remarkably precise agreement of the molecular orbitals and enables us to disentangle close-lying molecular states only separated by 50 meV at an energy of 2 eV below the Fermi level. We propose this combination of techniques as a promising new route for a precise electronic characterization of complex molecules and other physical properties which have electronic resonances in the tip-sample junction. Scanning tunneling microscopy is a powerful tool for determining the electronic structure of surface adsorbates. Here, carbon monoxide functionalized tips enable more accurate probing of the molecular states of graphene nanorings adsorbed on a gold surface.

new route for a precise electronic characterization of complex molecules and other physical properties which have electronic resonances in the tip-sample junction.

INTRODUCTION
Molecular adsorbates are characterized by electronic resonances that can be spectroscopically detected 1 .Scanning tunneling microscopy (STM) with its sub-molecular imaging resolution and its ability to probe the local density of states (LDOS) by recording the differential conductance dI/dV spectra yields spatially resolved images of molecular states 2,3 .Since scanning tunneling spectroscopy (STS) can access occupied as well as unoccupied molecular states, it has become a valuable tool to unveil the electronic structure of molecular adsorbates 4,5 .
So far, the spatial mapping of molecular resonances has been conducted mostly in the following way: The tunneling bias of the junction is set to the energy value at which the resonance occurs and, while the molecule is scanned by the STM tip, the dI/dV signal is recorded.Such scanning procedure can be performed in two distinct modes: constant-height or constant-current 6 .The obtained dI/dV map is commonly interpreted as a two-dimensional visualization of the molecular state.This approach has been successfully employed to molecular systems on which electronic resonances are well-separated from each other and in many cases decoupled from the metallic substrate by ultrathin insulators 3,7,8 , revealing a clear correspondence between calculated molecular orbitals and the measured electronic resonances.However, for the cases when the molecular resonance spectrum is denser and the detected electronic resonances have contributions from different molecular states 9 , the interpretation of the dI/dV data becomes a difficult problem which has in general not been solved yet.
In this work, we use a graphene nanoring, the cycloarene C108 adsorbed on Au(111), as a model system.Cycloarenes are considered as ideal test cases for studying electron delocalization in extended aromatic systems 10 .The C108 possesses a rich energy spectrum in which several molecular states are only separated in energy by a few meV 11 , making it particularly interesting for disentangling the contributions from the different molecular orbitals in the STM images.
To access, map and have a complete understanding of the C108 electronic structure, we use an approach that aims to overcome the limitation of common dI/dV imaging.For that, we combine CO functionalized tips, shown to provide sub-molecular resolution particularly on sp 2 -hybridized carbon nanostructures [12][13][14] , with feature-detection scanning tunneling spectroscopy (FD-STS) 15 , a statistical analysis tool that allows the unbiased detection of electronic resonances.We find that this combination enables us to classify, resolve and separate orbitals lying close in energy, providing experimental information on the electronic structure of a cycloarene in unprecedented detail.

STM and STS Characterization.
We image C108 at a base temperature of 10 K using a tungsten tip functionalized with a single CO molecule (CO tip), enabling sub-molecular resolution of the geometrical structure (Fig. 1d).Using the same tip, we also perform dI/dV spectroscopy revealing several electronic resonances in the energy range between -2eV and +2eV around the Fermi energy EF (Fig. 1e).The CO tip impacts the measured images and spectra compared to a bare metallic swave tip in two ways: first, at biases close to EF it allows to achieve structural sub-molecular resolution, similar to previous reports 13,[16][17][18] and, second, it enhances significantly the intensity of the electronic resonances.We identify the four main features as positive (PIR) and negative ion resonances (NIR) 11 (orange dots in Fig. 1e).Compared to previously published results, with the increased resolution of the CO tip, we additionally identify four secondary features shifted by approximately 150 mV from the main features to higher absolute bias (green dots in Fig. 1e).Density Functional Calculations.To understand the origin of the resonances, we perform density functional theory (DFT) calculations that use the B3LYP hybrid functional to calculate the density of states (DOS).From the gas phase DFT calculation (Fig. 2a,b) and by comparing with a representative dI/dV spectrum of C108 (red curve in Fig. 2b), we find that the highest occupied molecular orbital (HOMO, H0) at -1.094 eV, along with the orbitals HOMO-1 and HOMO-2 (H-1, H-2), situated both at -1.100 eV, are nearly degenerate and contribute to the measured resonance PIR1.The degenerate orbitals HOMO-3 and HOMO-4 (H-3, H-4), located both at -1.503 eV, contribute to the resonance PIR2 and HOMO-5, HOMO-6 (H-5, H-6) at -1.964 eV and -2.019 eV contribute to the resonance PIR3.On the unoccupied side (V > 0), we observe experimentally a broad NIR peak spanning through the energy position of three different calculated orbitals: the lowest unoccupied molecular orbital (LUMO, L0) at +1.830 eV and the degenerate LUMO+1 and LUMO+2 at +1.911 eV (L+1, L+2).
The significant intensity variations between the NIR and the PIR1-3, which are not found in this extent in the calculated DOS, stem from the strongly varying tunneling probabilities for tunneling into or out of these resonances.At positive V, the tunneling barrier height is effectively reduced, leading to an increase of the NIR, while for increasingly negative V, the barrier for tunneling from these resonances into the tip is successively increased 6 .Molecular Orbital Imaging.To investigate the spatial distribution of the molecular orbitals, we start by performing constant-current dI/dV maps at the bias voltage corresponding to the center of the PIR2, PIR1 and NIR resonances with a CO tip at relatively low tunneling resistances R = V/I ≈ 40 − 500MΩ (Fig. 3).To compare theory and experiment, we must consider the electronic broadening due to the bias modulation Vmod = 20 mV (peak to peak) used for the lockin detection, which limits the energy resolution.Therefore, we compute |∇rΨ(r)| 2 with ∇r as the in-plane gradient and Ψ(r) as the DFT calculated wavefunction and sum up those molecular orbitals which lie within the broadening range of the detection method (Fig. 3ac).The computed p-wave images and the measured constant current dI/dV maps match very well.
The good agreement between theory (Fig. 3a-c) and experiment (Fig. 3d-f) stems from a low hybridization between the C108 and the metallic surface and the energetically well-separated resonances, so that the constant current dI/dV maps do not capture contributions from other states, avoiding intermixing.
Remarkably, we find that the obtained constant current dI/dV maps of PIR1 and NIR (Fig. 3e, f) look very different compared to the constant height dI/dV maps obtained by Fan et al. at presumably significantly larger tip-sample separation. 11These differences are due to the CO tip's p-wave character at low tunneling resistances, i.e. small tip-sample separations, which is sensitive to the nodal planes of the molecular wave function 19 .The metallic s-wave character of the CO tip can be recovered by increasing the tip-sample distance and then the tunneling resistance on which the s-wave character predominates over the faster decaying p-wave character 20 (Supp.Fig. 1).We continue by investigating PIR3 which has contributions from H-5 and H-6.Because in this case H-5 and H-6 are separated by an energy difference of 55 meV, slightly higher than the modulation voltage, we compute p-wave images of the molecular orbitals for H-5 and H-6 separately (Fig. 4a, b) as well as for the sum of H-5 and H-6 (Fig. 4c).Comparing with the constant current dI/dV image, we find that it resembles more the sum of the two orbitals (Fig. 4d).Apparently, we see here the limit of constant current dI/dV map which do not allow to discern between two energetically neighboring molecular orbitals.Feature Detection Scanning Tunneling Spectroscopy Method.Until here, we have only explored the traditional way of how electronic properties of molecules have been characterized by STM.Now, we would like to utilize a method which gives a more complete picture of the molecular electronic properties even at larger tip-sample separations where the high resolution of the CO tip is no longer predominant.
For this purpose, we measure dI/dV spectra on a grid of 64 × 64 points covering an area of 2.4 × 2.4 nm 2 centered on an isolated C108 molecule with the CO tip.At each grid point, we regulate the tip height prior to opening the feedback loop to a tunneling current of I = 80 pA at a bias of V = 200 mV.The bias lies inside the energy gap between PIR and NIR to avoid the influence of molecular orbitals on the tip height and thus the resulting spectroscopic intensities.We then analyze the corresponding matrix of spectra using a recently developed FD-STS algorithm 15 .
The feature detection enables one to find and classify spectral features in a large set of individual dI/dV spectra.This method circumvents the necessity of using a model describing the physical processes and its subsequent time-consuming least-square fits.Because we are mainly interested in finding molecular states which are smeared out by hybridization with the substrate and by thermal broadening, the task is to detect broadened peaks on an unknown background (see the methods section for more details of the peak detection procedure).
To each detected peak in the raw dI/dV signal, we assign a weight W by finding the next lying extrema in the derivative of the signal, i.e. the second derivative dI 2 /d 2 V (Figure 5a).Note that even though W is not identical to the peak intensity in dI/dV, it is a measure commonly used, for instance, in Auger and ESR spectroscopy.By adjusting a threshold Wc, we disregard spurious noise peaks as well as peaks clearly originating from the substrate only (Fig. 5b).
The statistical evaluation of the features detected by FD-STS allows us to plot the (weighted) occurrence center energies Ep of the detected peaks in a certain energy range across the grid in an energy distribution histogram (EDH) (Fig. 5c), or to explicitly map the (weighted) occurrence of spectroscopic features in a feature distribution map (FDM).Finally, it is possible to map the energy position of a given feature, yielding an energy distribution map (EDM).Energy of the molecular states.We start our analysis by normalizing each individual spectrum to the intensity sum of the 3 PIRs that are the result of occupied molecular orbitals in the C108.
This normalization compensates for the different tunneling probabilities occurring at different locations on the C108 (red and gray curves in Fig. 2b) assuming an approximately constant total electron density within the spatial extension of the molecule.By applying the FD-STS algorithm, we are able to detect the center of the peaks in each spectrum.The EDH in Fig. 6a reflects the PIR and NIR main energies and intensities (red in Fig. 6a) as well as a series of side peaks (cyan in Fig. 6a) separated by bias voltages ranging from 140 to 210 meV.
As we will show, FDMs and EDMs are particularly useful for the classification and visualization of molecular orbitals that are very close in energy and are present in different regions of the molecule.In combination with the enhanced spatial resolution of the CO tip's p-wave nature 19,21 , this results in better understanding of the electronic structure than conventional dI/dV maps.
We start by selecting the energy range of the EDH that will be used for the FDM and EDMs (yellow boxes in Fig. 6b-d).The sharpness of the selected peaks suggests that only a single resonance is contributing, despite originating from the contributions of more than one molecular orbital.The produced EDMs use a linear scale to visualize the energy distribution of the selected resonance over the previously defined energy range (Fig. 6e-g).Additionally, the EDMs confirm that the resonances are energetically monodisperse varying approximately only by the energy of the modulation broadening within the molecule.In other words, we do not observe different resonances merged to a larger resonance.Unlike in previous work 15 , for the FDMs we plot the intensity at the center of each detected peak of the normalized dI/dV (Fig. 6h-j).We notice that even by adjusting Wc to zero, i.e. no filtering any peak attributed to spurious noise, or by using different window widths we obtain qualitatively the same EDH and FDM with the caveat of being noisier (Supp.Fig. 2, 3).
Next, we proceed to relate the spatial distribution of the detected features with FD-STS to molecular orbitals.In contrast to Fig. 3 and 4, we consider here an equal contribution of s and a pwave orbitals from the CO tip due to the increased tip-sample distance at the higher tunneling setpoint R = 2.5GΩ 20 .The FDMs reproduce the main features of the calculated s + p wave images (Fig. 6).For instance, the FDM obtained for PIR1 (Fig. 6i) shows lobes situated at the corners of the C108 molecule, in agreement with the corresponding DFT calculations (Fig. 6l).
We note that the molecular resonances decay in their intensity as the tip moves towards the inner part of the ring either in the constant-current dI/dV maps or in the normalized spectra used for the FDM.This effect can be compensated by decreasing the tip-sample distance as shown in Fig. 3 and 4.However, often this is not feasible due to the fragility of adsorbed molecules.This reveals the limits of not only our approach but also to similar experiments conducted on many different molecular systems probed both with metallic s-wave 8,22 and CO tips 23 , in which the molecular states are observed extending beyond the rim of the molecular structure, while the interior of the molecule tends to be electronically transparent.We observe in Fig. 6 that at larger absolute bias voltages (i.e.larger tunneling barriers) the FDMs tend to concentrate the features closer to the periphery of the molecule.
Remarkably, and contrary to the DFT calculation, the FDM as well as constant-current dI/dV maps (Fig. 3, 4 and Fig. 6) show an increased intensity in the central hole of the C108 nanoring.In this position, the tip is placed above the Au(111) surface, which is why we attribute the detection of these electronic resonances to a closer tip-sample distance caused by the directional sharper CO tip that is able to pick up the tails of electronic contributions from the molecular orbitals.On the contrary, such intensity in the center of the C108 molecule is not observed when the same set of measurements is performed with a metallic s-wave tip (Supp.Fig. 4) due to its relative bluntness compared to a CO tip.Vibronic replicas.We now focus on the side peaks separated by bias voltages ranging from 140 to 210 meV from the main peaks.We attribute these peaks to vibronic replica [24][25][26] of the C108 electronic states.This conjecture is supported by a calculation of the vibrational spectrum based on a finite-difference approach at DFT-PBE level; the observed vibronic replica stem most likely from a family of vibrational modes found at frequencies ranging from 1100 to 1600 cm -1 (Supp. Fig. 5), equivalent to 140 to 200 meV.
A close inspection of the vibrational replica (Supp.Fig. 6) reveals a substantial broadening of the EDH compared to the main peak, understood in terms of the symmetry dependence of the vibration assisted tunneling 27 .Additionally, by analyzing the resulting FDMs in Supp.Fig. 6, we observe that the maps reproduce qualitatively the main resonance to which it belongs, confirming that these peaks are the result of vibrational replicas.
Disentanglement of nearby molecular orbitals.Finally, we focus on the PIR3.In contrast to the maps of Fig. 6e-g, the EDM shows a resonance that is composed of two separable regions with slightly different energies: Yellow/green areas with higher peak energy and red/orange areas with lower peak energy (Fig. 7a).A closer look at the EDH reveals that a peak which appears to be monodisperse at first glance is a composition of two separate peaks whose energetic centers are To separate these two orbitals, we now fit this part of the dI/dV spectra with a sum of two Lorenzians with fixed energies at −1.95 and −2.00 eV as determined by the EDH.The resulting intensity maps (Fig. 7c, d) are in remarkable agreement with the DFT calculations (Fig. 7e,f), showing the potential of FD-STS to identify energetically overlapping molecular orbitals that can be separated to obtain a reliable image of the molecular orbitals in contrast to constant-current dI/dV maps showing contributions of both molecular orbitals (Fig. 4d).
Since also NIR contains two orbitals, we follow the same procedure.In this case, the resonance is the contribution of L0 and the degenerate L+1 and L+2 with an energy difference of 90 meV (Fig. 2b).We fix Lorentzian centers at an energy of 1.84 and 1.98 V, corresponding to the positions on which the central peak and the peak shoulders are found (Fig. 1b).In contrast to Fig. 7, however, we do not observe any significant difference in their intensities (Supp.Fig. 7) in agreement with the EDH (Fig. 6d), that only shows a single peak.We relate this to the lack of differences in the wavefunction intensities between the L0 and the degenerate L+1 and L+2.

DISCUSSION
We have shown that the combination of a CO tip and FD-STS is a powerful tool to identify, classify and map the molecular states of the cycloarene C108 in good agreement with theoretical calculations.With this combination, we reach a level of detail in the understanding of the molecule's electronic properties that is neither achievable by constant-current dI/dV maps nor by using a metallic tip and applying FD-STS alone.In particular, our method shows its potential by experimentally identifying and disentangling the close-lying molecular states at -2 V only separated by an energy difference of 50 mV.We note that this method could also be applied to systems on which different types of electronic excitations are observed, e.g., superconductivity, 28 spin excitations 29 or to molecular systems such as those incorporating transition metals or rare earth atoms 30 , on which theoretical calculations are limited due to their complexity.Furthermore, we expect that the methodology presented here will make the studies of the complex electronic structure of larger nanostructures more efficient.In the case of C108, the unambiguous identification of molecular resonances may help in future experiments exploring the intriguing effects of high magnetic fields such as persistent ring currents 31,32 or the Aharonov-Bohm effect 33,34 .

METHODS
Feature Detection Algorithm.We implement the feature detection method in the following way: The spectra containing the bare signal S are in our case dI/dV spectra obtained by modulating the swept bias voltage V with a small sinusoidal modulation voltage Vmod = 20 mV (peak to peak) and recording the first derivative dI/dV of the tunneling current I(V) with the help of a lock-in detection scheme.Because these spectra are equally spaced with respect to V, we can smooth out their intrinsic noise by using the well-established Savitzky-Golay method, which does a polynomial fit on a sliding window with very little computational effort 35 .Since we are mainly interested in finding peaks, a second order Savitzky-Golay filter, which uses local quadratic fits of the form S(x) = a + bx + cx 2 on each point of all spectra is sufficient.The Savitzky-Golay filter provides us not only with the smoothed S = a, but also directly with the smoothed first and second derivative (S' = b and S'' = c) of the signal S. To avoid filtering out meaningful signal components we use a sliding window of approximately 100 mV width, less than the full-width half-maximum of the expected peak structures.Furthermore, we check that the residuals, R = S-a, of the spectra do not significantly exceed the standard deviation σ of all residuals.
To find peaks (dips) in the spectra, we search for zero crossing in S' by evaluating ( ′ ()) − ( ′ ( + 1)) < 0(0), with () =  || ⁄ as the signum function and n as the bin of the linearly increasing bias V.Note that the method can be straightforwardly extended to also detect steps as well as "shoulders": masked peaks which due to some nearby larger structures do not develop a local maximum.
Density Functional Theory calculations.We performed Density Functional Theory (DFT) calculations using the Fritz-Haber-Institute ab initio molecular simulation package (FHI-aims) 36 .
The geometry of the graphene nanoring was optimized using Becke's three parameter hybrid exchange functional and the Lee Yang Parr correlation functional (B3LYP) [37][38][39][40] .Applying the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, we relaxed the structure via two consecutive steps; first by expanding the Kohn-Sham wavefunction with the default numerical "light" basis sets 41 , applied to all atomic species.Once finished, all basis sets were replaced with the default numerical "tight" basis sets 7 and the relaxation was resumed.We continued the relaxation until the maximum force on each atom, in either setup, was less than 10 -3 eV/Å.To aid the convergence in either relaxation step, we applied a 0.05 eV broadening to all states, using a Gaussian occupation smearing scheme 42 .It was later reduced to 0.01 eV for precisely calculating the molecular orbitals and density of states.Pair-wise intramolecular dispersion effects were included using the Tkatchenko-Scheffler (TS) approach 43 and relativistic effects were accounted for via scaled zeroth-order regular approximation (ZORA) 44 To obtain a well converged electronic description of the systems, a threshold of 10 -6 eV for the total energy, 10 -4 eV for the sum of eigenvalues and 10 -6 e/Å 3 for the charge density was applied during all self-consistent field (SCF) cycles.
Synthetic details of the precursor for C108.
Scheme 1 shows the procedure for the formation of Precursor 4. In the first step of the three-step synthesis, 1,3,5-tribromobenzene (1) was reacted with one equivalent of [1,1'-biphenyl]-2ylboronic acid to form terphenyl 2 under Suzuki-Miyaura cross coupling conditions.The corresponding boronic acid pinacol ester 3 was synthesized form 2 by Miyaura borylation reaction.The final reaction step for the formation of precursor 4 was accomplished by reacting 3 and 2 in 2:1 ratio under Suzuki Miyaura cross coupling conditions.

General Information
All reactions were carried out under inert atmosphere (nitrogen) using Schlenk techniques if not mentioned otherwise.All reagents were purchased from commercial sources if not mentioned otherwise and were used without further purification.For thin-layer chromatography, TLC plates from Merck KGaA with silica gel 60 on aluminum with fluorescence-quenching F254 at room temperature were used.All solvents were dried and/or purified according to standard procedures and stored over 3 Å or 4 Å molecular sieves.NMR spectra were recorded in automation or by the service department (faculty of Chemistry, Philipps University Marburg) with a Bruker Avance 300 or 500 spectrometer at 298 K using CD2Cl2 or CDCl3 as solvent and for calibration (residual proton signals).HR-APCI mass spectra were acquired with a LTQ-FT Ultra mass spectrometer (Thermo Fischer Scientific).The resolution was set to 100.000.HR-EI mass spectra were ac-quired with a AccuTOF GCv 4G (JEOL) Time of Flight (TOF) mass spectrometer.An internal or external standard was used for drift time correction.The LIFDI ion source and FD-emitters were purchased from Linden ChroMasSpec GmbH (Bremen, Germany).IR spectra are recorded with a Bruker Alpha FT-IR spectrometer with Platinum ATR sampling.

Figure 2 .
Figure 2. DFT-calculated DOS and molecular orbital energy levels of the cycloarene C108.(a) Simulated DOS of C108 in the gas phase (black curve).(b) Energetic positions of the gas-phase orbitals as found in the DFT calculation (black dots), overlaid with typical experimental dI/dV spectra, red and gray, taken at the positions of the red and gray dots, respectively, in Fig. 1d.Spectroscopy parameters prior to opening the feedback: V = 200 mV, I = 80 pA.Vmod = 20 mV.In both graphs, the region between -0.8 V and 1.2 V is omitted for a better visualization of the resonances.

Figure 4 .
Figure 4. Calculated gas-phase p-wave molecular orbitals and constant-current dI/dV maps of C108 obtained with a CO tip of PIR3.(a-b) Visualization of the molecular wavefunctions as calculated by DFT using a p-wave tip of H-5 = -1.964eV and H-6 = -2.019and the sum of H-5 and H-6 c).d) Constant-current dI/dV map performed with CO tip at the bias voltage VPIR3 = -1.95V and current setpoint of IPIR3 = 40 nA, with Vmod = 10 mV.Image sizes are 2.5 × 2.5 nm 2 .

Figure 5 .
Figure 5. Schematic representation of FD-STS algorithm.(a) dI/dV and d 2 I/dV 2 spectra for different

Figure 6 .
Figure 6.FD-STS with a CO tip and calculated gas-phase molecular orbitals of C108.(a) C108/Au(111) histogram of detected peaks (energy distribution histogram EDH, red) and vibronic replicas (cyan) produced by statistical analysis of a 64 × 64 grid of dI/dV spectra measured over the molecule.(bd) Energy distribution histograms (EDH) of the same resonances.The yellow boxes indicate the selected energy ranges chosen for the FDMs and EDMs in e-g and h-j, respectively.(e-g) Energy distribution maps (EDM) of PIR2, PIR1 and NIR (from left to right), overlaid with the chemical structure of C108.(h-j)

located at − 2 .
00 and −1.95 V (Fig 7b).The strong additional maximum in between the two peaks is an artifact of the FD-STS algorithm: if none of the two resonances is dominant, the superposition of both will result in an apparent single peak close to the intermediate energy, artificially creating an extra peak in the histogram.Remarkably, our DFT calculations predict two non-degenerate molecular orbitals (H-5 and H-6) at approximately −2 eV with a separation of only 55 meV (Fig. 7e, f).