Band gap renormalization of diamondoids: vibrational coupling and excitonic effects

We perform ab initio density functional theory calculations along with the frozen phonon approach to study the renormalization of the electronic structure of diamondoids induced by coupling to nuclear vibrations. We obtain an energy gap reduction between 0.1 and 0.3 eV across nine different diamondoids. The energy renormalization of the highest occupied molecular orbital state varies with the sizes and the symmetry of the diamondoids while that of the lowest unoccupied molecular orbital state is nearly constant. This behavior is explained using the localization of the electronic states. The energy gap renormalization is shown to originate mainly from bending C-H and rotation/shear H-C-H modes with frequencies between 800 and 1400 cm−1. We calculate the band gap renormalization due to excitonic effects using a screened configuration interaction approach and find excitonic binding energies of around 2 eV in excellent agreement with experiment for our largest diamondoids.


Introduction
Diamondoids are ultrasmall diamond nanoparticles composed of a few sp 3 bonded 'cages'. They have generated a wealth of research activity in the fields of drug delivery, bio-labeling, electron emitters, organic electronics, and nanolubrication due to their distinctive physical properties and their biocompatibility [1][2][3][4][5][6][7]. In addition to their current and potential use in real applications, diamondoids also attract attention from a fundamental point of view as they provide an ideal benchmark case for advanced quantum many-body calculations because of their simple and well-defined atomic structures and the availability of experimental results. Many-body studies on the electronic and optical properties of diamondoids were first performed in the middle of the last decade by Drummond et al [1] where the size-dependent optical gaps were calculated. Hereafter, Vörös et al did a series of studies on the electronic and optical properties of pure and modified diamondoids using both time-dependent density functional theory (TDDFT) and the GW approximation in combination with the Bethe-Salpeter equation (BSE) [4,7,8] which represent state-of-the-arts electronic calculations. On the other hand, a large value of zero-point energy gap renormalization (ZPR) induced by the electron-phonon coupling was found in bulk diamond [9][10][11][12], so that the coupling to vibrations can be expected to be non-negligible in carbon based diamondoids. However, obtaining the phonon induced energy gap renormalization has been a challenge in semiconductor physics since the fifties of the last century [9]. Following the pioneering investigation on Fan [13] and Debye-Waller (DW) terms [14] of the self-energy, Allen, Heine, and Cardona (AHC) combined both SE and DW interactions to explain the phonon induced energy gap renormalization of bulk semiconductors via semiempirical pseudopotential calculations [15][16][17]. Recently, theories based on the AHC approach have been implemented in a fully ab initio framework [10][11][12][18][19][20][21], path-integral molecular dynamics simulation [22], and first principles frozen-phonon calculations [12,21,[23][24][25]. A dynamical extension of the AHC approach [11] had been applied to small diamondoids [26]  photoemission experiments. Another approach [6] based on the Franck-Condon idea and a Monte Carlo sampling technique was used to calculate the optical absorption spectra of small diamondoids. Both of these very recent contributions [6,26] highlight the importance of the vibrational coupling in these structures.
In this work, we perform ab initio density-functional theory (DFT) and screened configuration interaction (CI) calculations using the frozen-phonon approach to study the energy gap renormalization as a consequence of zero-point motion in diamondoids C 10 32 . We show that the ZPR of diamondoids is around 0.1 to 0.3eV, which is smaller than the corresponding value of bulk diamond (0.65 eV [12]). We further analyze the ZPR of diamondoids and find that the zero-point energy correction of the lowest unoccupied molecular orbital (LUMO) state is nearly constant, while that of the highest occupied molecular orbital (HOMO) state varies with the size and the geometry. This fact can be explained by the different distributions of the electron and hole wavefunctions. The LUMO states are localized just outside the diamondoids [1,4,6,27,28] and have surface-bound character due to the negative electron affinity (see [28] for details). Furthermore, the ZPR of diamondoids is also heavily dependent on which vibrational mode the electronic states couple to. We identify that the bending C-H and the rotation/shear H-C-H modes with frequencies between 800 and 1200cm −1 dominate the ZPR of diamondoids. Finally, we calculate the excitonic contribution to the band gap renormalization using a screened CI approach and obtain exciton binding energies between 2.0 and 2.5eV in excellent agreement with experiment for the larger diamondoids (more than 18 carbon atoms). However, the approach is shown to quantitatively fail for smaller structures.

Theoretical method
The geometry of C 10 H 16 , C 14 H 20 , C 18 H 24 , C 22 H 28 (with C 2 , C 2h , and C 3v point-group symmetry), C 26 H 30 and C 26 H 32 (with T d and C 2v point-group symmetry) diamondoids, which have been fabricated and studied experimentally [29][30][31], are shown in figure 1. The common names of these diamondoids along with their corresponding sizes and point group symmetries are listed in tables 1 and 2.
It is well known that the calculation of semiconductors with a local functional for exchange and correlation, as in the local density approximation (LDA), bares large systematic errors with an underestimation of the band gaps [32]. However, the DFT-LDA works remarkably well to describe electronic wavefunctions (which are the wavefunctions typically used in non-selfconsistent GW approaches and BSE [32]) and lattice vibrations of semiconductors [33], which are used to calculate the band gap renormalizations. In this work, the ground state geometry, electronic structure, and lattice vibrations of the diamondoids are computed using DFT in the LDA approximation. The DFT-LDA calculations are performed using Trouiller-Martin norm-conserving pseudopotentials with an energy cutoff of 30 Ry. In order to avoid the interaction between periodic images of the structure, a simple cubic supercell with lattice parameter of 20Åis used, where the space between the outermost atoms and the closest boundary is more than 5Åfor all the studied diamondoids. The atomic positions of these diamondoids are relaxed until the forces are reduced to less than 3×10 −6 Ha Bohr -1 under constrained symmetry. Based on the harmonic approximation, the dynamical matrix elements are then calculated via finite difference as implemented in the CPMD code [34]. The lattice vibrational eigenmodes and eigenvectors are thereafter obtained by solving the dynamical matrix [35], where, M I and M J are the masses of atoms I and J, ( ) V R is the potential energy, R represents the atomic positions while n X and w n are the vibrational eigenvectors and frequencies of a ν-mode vibration, respectively. We calculate the ZPR using a frozen-phonon approach, where the atoms I are displaced from their equilibrium position by 'freezing-in' a ν-mode vibration as shown in [23,25], with n X I the three components of the vibrational eigenvectors X belonging to atom I and mode ν. The ZPR of state n induced by one of the 3N-6 ν-mode vibration is given by D n E n and is computed from the average of positive and negative displacements  n u I to eliminate the linear term [23,25]. In order to ensure that the energy shifts extracted from positive and negative displacements correspond to the same state, we project the wave functions with negative displacement on those with positive displacement [25]. The temperature-dependent ZPR D ( ) E T n is therefore calculated by summing over all the vibrational modes, We reduce the computational effort significantly by using the symmetry elements of the point group of the studied diamondoids. In order to compare the energy gap renormalization originating from ZPR to the ones originating from the exciton binding energy E ex , we use a screened CI approach [36] to compute the latter. The idea of the CI framework is to expand the exciton wave function Y a ( ) in terms of different 'configurations' where, α denotes the exciton quantum number and F v c , is the Slater determinant of a certain configuration composed of single-particle wavefunctions obtained via DFT. We include twelve valence (HOMO-n) and  Table 1. Common name, diameters excluding/including hydrogen atoms, percentage of carbon atom with one (C-H), two (C-H 2 ) and zero (C-C) hydrogen passivants.
with the Hamiltonian matrix elements where, E c and E v are the single-particle electronic states obtained from the DFT calculations, and are the screened Coulomb and exchange integrals, respectively. By solving the secular equation given in equation (5), the exciton-binding energy is then obtained as

ZPR of diamondoids
In tables 1 and 2 we summarize the structural information about the considered diamondoids and indicate if the optical dipole transition is allowed by symmetry. For the diamondoids C 14 H 20 and C 22 H 28 in the C 2h point group symmetry, the lowest energy HOMO-LUMO transition is forbidden. In these two cases we will focus on the ZPR of the allowed HOMO-(LUMO+1) transition. Note that in these cases, the transition between LUMO and the excited HOMO-1 state is forbidden and transitions to higher excited occupied molecular orbital states lie at higher energy. Our frozen phonon approach underlying the calculation of ZPR allows to separate the effect of ZPR on the HOMO and the LUMO or (LUMO+1) states. We plot in figure 2 the ZPR of the HOMO (red squares), the LUMO (or LUMO+1) states (green triangles) and the ZPR effect on the gap given as the sum of HOMO and LUMO (or LUMO+1) contributions (black circles). From the DFT-based frozen-phonon calculations, the eigenvalues of the unoccupied states (LUMO or LUMO+1 states) decrease while those of the occupied states (HOMO state) increase with the introduction of ZPR. The ZPR-induced energy shifts of both occupied and unoccupied states lead to a reduction of energy gap, and we chose to label the contributions of both HOMO and LUMO (or LUMO+1) states with negative sign (which is another convention than used in [25]) in figure 2. From this figure along with table 2, we see that the total ZPR varies not only with the number of atoms but also with their point group symmetry. We also see that the non-monotonous size and symmetry dependence of the total ZPR originates mainly from the HOMO state, while the ZPR correction of the LUMO (or LUMO+1) state is nearly a constant (around −0.1 eV) for all diamondoids. To understand this behavior, we plot the electronic wavefunctions of the HOMO and the LUMO (or LUMO+1) states of all the diamondoids in figures 3(a)-(i). We see that the HOMO state wavefunctions, which originate from the σ-orbitals of the valence band p-states, vary significantly from one structure to another. However, the LUMO (LUMO+1) states are surface-bound states as described previously [28], which originate from a surface dipole, a curved geometry, a dielectric mismatch, and a negative electron affinity as described in [28]. These states are localized slightly outside the structures and look very similar for all the diamondoids.
In figure 4 we plot the vibrational density of states (VDOS) of our nine diamondoids along with the induced band gap renormalization DE g contributed by each vibrational mode. The VDOS is color-coded to show the character of the different vibrational modes: black, red and blue for carbon atoms connected to two, one and zero hydrogen atoms, respectively, and the VDOS contributed by hydrogen atoms in cyan. The contribution from carbon atoms connected to one hydrogen (red) is dominant, which can be explained by the fact that there is a large fraction of such type of atoms in these structures. Numerical values for the number of each type of carbon atoms are given in table 1 in percent of the total number of carbon atoms. Note that for all the diamondoids we Table 2. Point group symmetry, irreducible representations of the HOMO, LUMO and LUMO+1 states, character (allowed/forbidden) of the dipole transitions HOMO  LUMO (H-L) and HOMO  LUMO+1 (H-(L+1)).

Symm.
HOMO LUMO LUMO+1 H-L H-(L+1) consider the fraction of bulk-type carbon atoms (the ones connected to four neighboring carbon atoms, labeled as C-C) does not exceed 23%. Hence, a bulk-type vibrational DOS with acoustic and optical branches can hardly be recognized in these structures, in contrast to larger quantum dots [33,37]. The ZPR (DE g ) shows contributions mainly from the frequency range between 800 and 1400 cm −1 , which correspond to bending C-H modes and rotation/shear H-C-H modes. Especially for our largest diamondoids the strongest contributions to DE g seem to be at higher frequency, towards the vibron gap between the carbon and the hydrogen vibrations at around 1400 cm −1 . Contributions from (nearly) pure hydrogen at around 3000 cm −1 are very small. Furthermore, the contributions of carbon atoms connected to one hydrogen (red) along with their corresponding hydrogen atoms (cyan) dominates the ZPR of diamondoids. However, with increased fraction of C-C bonds (larger structures), the contributions of bulk-type carbon atoms also start to play a role on the ZPR so that the observations can merely be related to the bare number of C atoms of different types. It is also interesting to note that for diamondoids with 14, 18, and two of the structures with 22 C atoms the positive contributions to DE g (reducing the band gap) are somewhat counterbalanced by negative contributions. From our limited results we cannot identify a clear rule as to when negative contributions occur. A plain symmetry argument does not hold as we have two diamondoids with C 2v point group symmetry and only one of both exhibits negative contributions. A size dependence is also clearly not recognizable. A further analysis of these results in terms of the diagonal DW terms, the non-diagonal DW terms and and the Fan terms, as performed in [38] for small molecules, may lead to a deeper understanding of this effect.

Energy gap renormalization: vibrational coupling
In figure 5(a) we plot the experimentally measured optical band gaps based on two different experiments, exp.1 [29] and exp.2 [4,29], showing a close agreement between the measurements. From these rather accurate experimental values, we subtract our calculated ZPR, leading to a larger band gap. These curves at higher energy are therefore the 'target' for many-body calculations such as TDDFT or GW+BSE in case vibrational coupling is neglected. The TDDFT [4] and GW+BSE [39] are given as right-and left-triangles, respectively. They both seem to agree well for the smallest two diamondoids, while they seem to underestimate the band gap for larger structures. Note that our aim is not to compute the electronic band gaps, since it was accurately done by others [4,39], but calculate ZPR and exciton binding energies, which are rather accurately given using vibronic and electronic states from standard DFT.
In figure 5(b) we plot the experimentally measured ionization potential. Since our method allows the calculation of the ZPR separately for HOMO and LUMO we can extract of the 'bare' ionization potential, that should be the target to many-body approaches that ignore the vibrational coupling. These values are given as stars in figure 5(b).
The ZPR of adamantane (C 10 H 16 ), diamantane (C 14 H 20 ), and triamantane (C 18 H 24 ) have been calculated in with w i 0, and w s i , the vibrational frequencies with index i in the ground (index 0) and excited singlet state (index s). The ZPR reported by using this very different approach are 0.33eV, 0.27eV, and 0.25eV for the three structures with increasing cluster size, which compares generally very well with our results: 0.33eV, 0.12eV, 0.23eV, with some deviation for the structure C 14 H 20 . This latter structure has significant negative contributions to the ZPR (see figure 4) and has a ZPR value, according to our results, which is lower than other diamondoids. The comparison to experiment in figure 5 also shows, that the GW and TDDFT seem to be generally to low (crosses compared to triangles) for all the diamondoids by a value between

Energy gap renormalization: excitonic effect
We show in figure 6 the exciton-binding energies of C 10 H 16 , C 14 H 20 , C 18 H 24 , C 22 H 28 (C 2h ), and C 26 H 32 (T d ) diamondoids calculated using our screened CI approach at the CI-level (black circles) and at the singleconfiguration (SC)-level (red triangles). The difference between the SC and the CI energy is due to correlation effects [40]. We compare our results to BSE [39] (green squares) and to experimental measurements [29,30] (blue crosses). The calculated exciton-binding energies are given for the lowest energy optically active states not necessarily for the lowest energy state. In the screened CI approach [36,40] we use a model dielectric function [41], is a so-called mask function with w 0 = 0.163 a 0 (a 0 is the bulk lattice constant) [42] and is plotted in the inset of figure 6, and  ( ) r r , dot 1 2 represents the dielectric function of bulk diamond.
Details of the screening function used in the screened CI calculation are given elsewhere [40] and we briefly outline it in the following. In the screened CI calculation, the non-local dielectric function is treated using a phenomenological isotropic and uniform model proposed by Resta  where,  ¥ is the high-frequency dielectric constant of the bulk material as 5.7 [43], q is the Thomas-Fermi wavevector, r ¥ denotes the solution of . For R we use half of the diameter of the carbon cluster defined as the distance between the farthest neighbor carbon atoms, which represents a rough Figure 5. (a) Experimental values of the optical gaps are given as squares and circles and are connected by lines. Experimental data subtracted by our results for the ZPR calculated at T = 300 K are shown as two types of crosses. The GW+BSE [39] and TDDFT [4] results are shown as two types of triangles. The theoretical results (triangles) must be compared to the corrected experimental results (crosses). (b) Experimental values of the ionization potential from [30] along with our 'bare' value explicitly excluding the ZPRs, which should give a target value for TDDFT or GW approaches neglecting vibrational coupling. Figure 6. Exciton-binding energies of diamondoids obtained from our CI (black circles) and SC (red triangles) calculations, BSE calculations (green squares) given in [39], and experimental measurement (blue crosses) given in [29,30]. Insert: mask function of C 26 H 32 (T d ) diamondoid. approximation, especially from small diamondoids which shapes deviate significantly from spheres. With these definitions, when r 1 and r 2 are both close to the center of the diamondoids, the dielectric function is the full, unmodified, bulk diamond function. When both are well outside the structure, the interaction is unscreened. When r 1 is inside and r 2 outside the structure (or vice versa), the dielectric function is scaled down according to the mask function. This approach requires three separate calculations since we perform the integration in reciprocal space, namely a fully unscreened Coulomb and exchange integral calculation . The effective screened Coulomb and exchange integrals are then obtained as = + -