Kekulene: On-Surface Synthesis, Orbital Structure, and Aromatic Stabilization

We revisit the question of kekulene’s aromaticity by focusing on the electronic structure of its frontier orbitals as determined by angle-resolved photoemission spectroscopy. To this end, we have developed a specially designed precursor, 1,4,7(2,7)-triphenanthrenacyclononaphane-2,5,8-triene, which allows us to prepare sufficient quantities of kekulene of high purity directly on a Cu(111) surface, as confirmed by scanning tunneling microscopy. Supported by density functional calculations, we determine the orbital structure of kekulene’s highest occupied molecular orbital by photoemission tomography. In agreement with a recent aromaticity assessment of kekulene based solely on C–C bond lengths, we conclude that the π-conjugation of kekulene is better described by the Clar model rather than a superaromatic model. Thus, by exploiting the capabilities of photoemission tomography, we shed light on the question which consequences aromaticity holds for the frontier electronic structure of a π-conjugated molecule.

Single crystals of 5 appropriate for X-ray diffraction study were obtained by slow diffusion of n-pentane in a saturated solution of 5 in chloroform. The molecular structure of 5 derived from X-ray diffraction is shown in Figure 1b (top). In solid state, the three phenanthrene moieties of 5 are rotated out of the molecular plane made up by the double bonds. As a consequence, 5 has a highly distorted nonplanar geometry. The solid state packing of 5 (Figure 1b, bottom) is driven by π-stacking interactions of the phenanthrene moieties with intermolecular distances of 3.34 Å between the carbon atoms of neighbored molecules, which is only little larger than the interlayer distance in graphite (equilibrium interlayer distance is 3.231 Å [S1.3] ).

S1.2. General Information
All reactions were carried out under inert atmosphere of 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), 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. Nuclear magnetic resonance (NMR) spectra were recorded in automation or by the service department (Department of Chemistry, Philipps University Marburg, Germany) with a Bruker Avance 300 or 500 spectrometer at 298 K using CD2Cl2 or CDCl3 as solvent and for calibration (residual proton signals). High-resolution atmospheric pressure chemical ionization (HR-APCI) mass spectra were acquired with a LTQ-FT Ultra mass spectrometer (Thermo Fischer Scientific). The resolution was set to 100000. High-resolution electron ionization (HR-EI) mass spectra were acquired with a AccuTOF GCv 4G (JEOL) time-of-flight (TOF) mass spectrometer. An internal or external standard was used for drift time correction. The liquid injection field desorption ionization (LIFDI) ion source and field desorption (FD) emitters were purchased from Linden ChroMasSpec GmbH (Bremen, Germany). The data collection for the single-crystal structure determination was performed on a Stoe Stadivari diffractometer or a Bruker D8 Quest diffractometer (X-ray Service Department, Philipps University Marburg, Germany). Information concerning the hardware and software used for data collection, cell refinement and data reduction as well as structure solution and refinement can be reviewed in the attached crystallographic information files (CIF). After the solution (Shelxt) [S1.4] and refinement process (Shelxl 2017/1) [S1.5] the data was validated by using Platon. [S1.6] All graphic representations were created with Diamond 4. [S1.7] Infrared (IR) spectra are recorded with a Bruker Alpha Fourier-transform infrared spectrometer (FT-IR) with platinum attenuated total reflectance (ATR) sampling.

S3.2. Density functional theory calculation of kekulene/Cu(111)
In a search of the energetically most favorable adsorption site and configuration of kekulene on Cu(111), geometry relaxations of kekulene/Cu(111) have been performed using two different theoretical methodologies. In both cases we employed density functional theory (DFT).
In the first approach (theoretical setup 1), the Vienna Ab Initio Simulation Package (VASP) [S3.4, S3.5] , which uses a plane wave basis and the projector augmented wave (PAW) method to account for the separation of core and valence electrons, was employed. Exchange and correlation effects were treated by the PBE-GGA [S3.2] with van der Waals corrections according to Tkatchenko and Scheffler (TS) [S3.6, S3.7] . For the structural model of kekulene on Cu(111), we used the optimal lattice parameter of Cu within PBE-GGA+TS (3.55 Å) and modeled the surface with 3 Cu layers in a repeated slab approach with a sufficiently large vacuum layer to separate periodic replica from each other. In order to counterbalance the artificial electric field arising from the asymmetric slab, a planar dipole layer with a self-consistently determined dipole moment was added in the middle of the vacuum region [S3.21] . Four different adsorption sites were tested as the starting point for the relaxations: the top, bridge, hollow fcc and hollow hcp positions. Additionally, two azimuthal orientations of kekulene differing by 30° were considered for which we used a supercell size of 7x7x3. For the Cu substrate all degrees of freedom have been relaxed. The obtained adsorption heights and energies for different structural configurations are provided in Table S3. The energetically most favorable adsorption configuration was found to be the hcp hollow site with kekulene's zigzag edge oriented along the [11 # 0] direction of the copper substrate. To test the robustness of the found optimal adsorption configuration, we performed additional calculations using the experimentally measured lattice parameter of Cu (3.61 Å [S3.8] ) and the fixed positions of the copper atoms, taking into account 6 Cu layers and using the same theoretical setup 1. We found the adsorption energy for the same adsorption configuration (hollow hcp, kekulene's zigzag edge along [11 # 0] direction) slightly reduced (5.925 eV), while the adsorption height remained almost unaffected (3.05 Å).
To cross-check the results obtained using the above described theoretical methodology, the adsorption calculations were independently repeated using the SIESTA code [S3.9,S3.10] (theoretical setup 2). SIESTA uses a localized basis set of numerical orbitals [S3.11, S3.12] S3.7]. To avoid basis set superposition errors (BSSE) from our localized basis in this setup 2, our structural optimization proceeds as follows.
After an initial optimization of the geometry, we calculate the binding energy for several surfacemolecule distances by shifting the rigid molecule. To obtain the binding energy without BSSE, we additionally perform two calculations for each height, where we only include the surface or molecular atoms but use the same full basis set as before. We then perform a new BSSE-corrected geometrical optimization where we keep the average binding distance fixed at the minimal position of the BSSE corrected binding energy curve calculated before. From that BSSE-corrected geometry we then obtain the final binding energy curve. The resulting binding energies and distances for all absorption sites and orientations are shown in the last two columns of Table S3. However, the geometrical optimization in this setup 2 will not converge at the bridge position but the molecule shifts in the hollow hcp adsorption site for both orientations. Note that we calculated the hollow hcp position with kekulene's zigzag edge parallel to the [11 # 0] direction in this setup 2 also with an 8x8x6 geometry for comparison, but obtain nearly the same values for the binding energy and distance. Apparently, the geometries obtained by setup 2 are about 0.15 Å closer to the surface but slightly less bound in terms of adsorption energy, reflecting the uncertainties of computational physisorption still present today. Apart from these details, both theoretical setups yield very similar trends concerning the positioning and orientation of kekulene relative to the Cu(111) surface, with the hollow hcp site and orientation along [11 # 0] being the optimum geometry.
The momentum map of photoemission from the HOMO state of kekulene adsorbed on Cu (111)  The typical ring-like appearance of kekulene as also seen experimentally can be clearly recognized in Figure S4.

S3.3. Quantification of aromaticity with the harmonic oscillator model of aromaticity
In order to assess the aromaticity of kekulene quantitatively, we have employed the harmonic oscillator model of aromaticity (HOMA) [S3.18-S3.22] . The optimal C-C bond length Ropt was estimated by averaging the C-C single and double bonds' lengths of the reference molecule trans-1,3-butadiene. In order to get results best comparable to the kekulene calculations discussed above, the geometry of the reference trans-1,3-butadiene molecule was optimized using the same methodology as for free kekulene (cf. S3.1.   The HOMA values for the DFT-optimized structures of free kekulene and kekulene adsorbed on Cu(111) are given in the two central columns of Table S4. Note a moderate effect of the substrate on the aromaticity of different conjugation paths. The increase of H [18] and HA (13 and 24%, respectively) points at an "enhanced" superaromaticity of kekulene upon adsorption on Cu(111).
The Clar and superaromatic models of the free kekulene were constructed using the bonds labeled as in Figure S5. For this, the lengths of single and double C-C bonds were taken from the reference molecule trans-1,3-butadiene as mentioned above. For the Clar model, Rc = Rd = Rf = Ropt to fulfill the ideal aromaticity of the phenanthrene terminal ring of kekulene marked in Figure S5 with B, i.e., HB = 1. For the superaromatic model, Ra = Rb = Rc = Re = Rf = Ropt to fulfill the ideal aromaticity of two annulene paths -inner ([18]annulene) and outer ([30] (double) bond of the reference trans-1,3-butadiene. The inter-bond angles are adjusted in order to conform the six-fold symmetry of kekulene and to fulfill the constraints mentioned above. Here the angle β inside sextet A was used specifically to adjust the length of bond a.

S3.4. Density functional theory calculation of model kekulene structure
For the model structures introduced in S3.3, i.e., the Clar model and the superaromatic model, we have calculated the electronic structure using the same methodology as in S.3.1, however, with the molecular geometry fixed. In particular, the spatial distributions of the wave functions and the total charge densities of occupied states were calculated. The former were used to simulate the theoretical momentum maps of photoemission distribution following ref [S3.3].

S3.5. Hückel's tight binding model calculations
To analyze how the electronic coupling between inner ([18]) and outer ([30]) annulene paths of kekulene affects the charge density distribution, we have employed a simple Hückel model (HM) [S3.23-S3.26] . The molecular structure of kekulene was fixed according to the superaromatic model introduced in S3.3. The HM assumes that molecular orbitals can be written as a linear combination of atomic pz orbitals. For planar, unsaturated hydrocarbons the linear combination is given by the following Hamiltonian matrix .
The eigenvalues of H are the Hückel molecular orbital energies expressed in terms of ε and t. The value of ε represents the energy of the 2p electron relative to an unbound electron and is given by the on-site integral , while t is given by the hopping integral which represents the energy of an electron localized in a 2p orbital compared to a delocalized p orbital formed by the 2p orbitals of adjacent atoms. In the case of kekulene, the Hückel matrix is a (48x48) matrix, which is depicted in Figure S6.
Using this model, we varied the hopping parameter td of the bond connecting the inner and outer annulene paths of kekulene from td = 0 (no coupling between annulene paths) to td = -3.5 eV (typical coupling between C atoms in graphene). ta, tb, tc, te, and tf were kept constant and equal -4.1 eV and ε was set to -3.35 eV, resulting in a good agreement of the orbital energies compared to the DFT results. Note that in the paper and S3.6. td is referred to as t and ta, tb, tc, te, tf are not mentioned.

S3.6. HOMO and HOMO-1 orbitals, their charge density distributions and momentum maps
The real space distributions of the doubly degenerate states HOMO and HOMO-1, their charge density distributions and corresponding momentum maps for four different model cases of free kekulene are shown in Figure S7, namely for the Clar geometry calculated with DFT (first column) and the superaromatic geometry calculated with DFT (second column) and HM for t = -3.5 eV (third column) and 0 (fourth column).
We note a qualitative similarity of the orbitals and the charge density distributions, as well as the simulated momentum maps, for the first three cases ( Figure S7, first, second and third columns). Only if the coupling between outer [30]annulene and inner [18]annulene is artificially suppressed (t = 0), the electronic state of kekulene (i.e., its HOMO and HOMO-1) and the expected distribution of photoelectrons differ significantly. In the case of t = 0, HOMO and HOMO-1 and their charge densities are restricted, correspondingly, to the outer [30]annulene and inner [18]annulene paths ( Figure  S7d,h and S7p,t, respectively), while for all other cases they are delocalized over the entire molecule. Turning on the coupling between [30]annulene and [18]annulene results in the adiabatic change of the kekulene HOMO toward a lobular structure that comprises benzene-type "HOMOs" at the terminal rings of the phenanthrene units (sextet B). This is illustrated in "inter-annulenecoupling.mp4" video provided electronically. Thus, even a small coupling between the annulene paths favors the formation of Clar sextets. The corresponding modification of the electronic structure, in particular the spatial distribution of charge density, is immediately reflected in the photoemission patters of the related molecular states (cf. simulated momentum maps in Figure S7 and the attached video). Figure S7: (a-d) Degenerate HOMO of the free kekulene molecule, (e-h) charge density distribution, and (i-l) simulated momentum maps for each degenerate component of the HOMO (upper part) and their incoherent sum (lower part) for the Clar geometry of kekulene calculated with DFT (first column) and for the superaromatic geometry of kekulene calculated with DFT (second column) and with the Hückel model for the inter-site hopping parameter t = -3.5 eV (third column) and 0 (fourth column). Analogous data for the HOMO-1 of free kekulene is plotted in (m-p), (q-t), and (u-x). Figure S8 shows the correlation between the C-C bond lengths of free kekulene, calculated with DFT, and a normalized bond order parameter from Hückel theory, calculated taking into account either only the lowest 12 occupied π-orbitals (panel a) or only the uppermost 12 occupied π-orbitals (panel b). The opposite slope of the corresponding linear fits reveals an opposite correlation. The observed tendency of an enhanced effect of the higher occupied π-orbitals on the C-C bond lengths is further illustrated by Figure S9. Here, the bond order is calculated by including subsequent groups of four orbitals into summation. The slope of the linear fit of the corresponding plots ( Figure  S9b) changes sign, thereby pointing at the prevailing influence of the uppermost molecular states.