Optical Properties of Graphene Nanoflakes: Shape Matters

In recent years there has been significant debate on whether the edge type of graphene nanoflakes (GNF) or graphene quantum dots (GQD) are relevant for their electronic structure, thermal stability and optical properties. Using computer simulations, we have proven that there is a fundamental difference in the calculated absorption spectra between samples of the same shape, similar size but different edge type, namely, armchair or zigzag edges. These can be explained by the presence of electronic structures near the Fermi level which are localized on the edges. These features are also evident from the dependence of band gap on the GNF size, which shows three very distinct trends for different shapes and edge geometries.


I. INTRODUCTION
Graphene, a 2D sheet of carbon atoms arranged in a honeycomb lattice, has arguably been the most promising and investigated material since it was first isolated in 2004 1 . Due to its remarkable properties such as unusually high current capacity 2 , mechanical strength 3 and thermal conductivity 4 , several applications have been conceived for this material. However, its zero band gap and semi-metal character represent a limitation that needs to be sorted out for carbon-based electronics to be feasible. Some ways to open up a gap are nanostructuring graphene into graphene nanoribbons (GNR) and nanoflakes (GNF) or carbon nanotubes 5,6 as well as chemical functionalization 7 . Fortunately these compounds have been extensively studied from the experimental point of view, both for GNR [8][9][10][11] and armchair GNF of triangular 12 , hexagonal 13,14 , and other shapes 15,16 .
Quasi-zero-dimensional nanostructures derived from graphene such as GNF (also known as called graphene quantum dots) are good candidates for the aforementioned applications. The quantum confinement caused by reduced dimensionality of these structures opens a tunable bandgap 7,17 . Several studies have characterized the electronic properties of GNF 5,18-39 and some have shed light on their optical features [20][21][22][23][24]27,33,34,[40][41][42][43][44] . Yamijala et al. 42 performed a systematic study of the linear and nonlinear optical properties for ∼400 graphene quantum dots at the ZINDO/S semiempirical level. However, there is not, to the best of our knowledge, a study relating the optical features with structural and electronic aspects, for various shapes, sizes and types of edges.
Actually, much of what is known for GNR is related to the properties of GNF: in armchair GNR there is always a band gap, which can be attributed to edge effects, while for zigzag nanoribbons there are always zero-energy states mainly confined along the edges, causing a peak of the density of states at the Fermi level 19 . These edge states have been experimentally observed 45 . Some of these features also hold for GNF: armchair flakes show a band gap,although in zigzag edges gap states occur only in triangular flakes 18,28 . Moreover, the electronic 46,47 , magnetic 48 and optical properties 49-52 of GNR are well understood from the theoretical point of view.
There seems not to be a general agreement on whether edge and corner terminations affect the electronic properties of GNF. On the one hand, it was shown that the edges and corners in GNF have little influence on the distribution of the highest occupied molecular orbital (HOMO) 26 , although the Fermi level of hexagonal flakes is independent on the type of edges but rather sensitive to corner reconstructions 53 . On the other hand, zigzag edges lower the band gap with respect to armchair edges in GNF 38 and the electronic density of states is modified by edge and corner geometries even in GNF with the same shape and similar size, both for hexagonal 53 and triangular structures 31 .
In this work we show calculated optical absorption spectra of triangular and hexagonal GNF for a range of sizes and edge structures, and give an explanation based on the electronic structure for the different trends. Both edge and corner terminations effects on the excitation energy spectrum are analyzed. Results suggest that previously characterized edge states are responsible for the differences observed in the calculated absorption spectra.

II. COMPUTATIONAL METHOD
In this work we have studied four classes of hydrogen passivated GNF with different shape, edge and size, namely hexagonal and triangular flakes with armchair and zig zag edge. A representative structure for each of these systems is shown at the bottom of Figure 1. Throughout the paper we have defined N as the total number of atoms, and adopted the following abbreviation for the structures: TAC for triangular armchair, TZZ for triangular zigzag, HAC for hexagonal armchair and HZZ for hexagonal zigzag. The electronic structure of the system was described using the self-consistent charge density functional tight binding (SCC-DFTB) method, which is based on a second-order expansion of the Kohn-Sham energy around a reference density of the neutral atomic species 54 . This method has been extensively used for the study of graphene and graphene-based structures 26,31,38,42,53,55 and has been benchmarked with respect to DFT for several graphene based systems with defects, with bond length discrepancies around 2% and formation energies matching DFT values within 1.5% 56 . We have used the DFTB+ implementation of DFTB 57 , with the mio-1-1 set of parameters 58 , to obtain the hamiltonian, overlap matrix and the initial ground state (GS) single-electron density matrix. All purely electronic properties such us the (total and projected) density of states and HOMO-LUMO energy gaps are calculated from these results.
The use of the DFTB hamiltonian allows for the calculation of the optical properties of very large systems (results for up to 1400 atoms are shown in this work) that would be unfeasible using ab-initio methods. The DFTB model provides a well tested approximation to the self consistent electronic structure and its response to external electric fields which can however deal with systems containing thousands of atoms.
It was been shown that structures with zigzag edge are magnetic 47,48,59,60 . Therefore, spin-polarized calculations were carried out for all the TZZ systems and the two largest HZZ systems using the DFTB implementation for collinear spin polarization 54 . The rest of the them present no net polarization in the GS. The total spin polarization obtained agree to those for which there exists data in the literature 59 . Structural optimization was carried out by minimizing the force restricting the z coordinate to keep a planar configuration.
The methodology applied for the calculation of optical properties is based on the real-time propagation of the density matrix of the system, and has been successfully used in our group to study a range of molecular and nanostructured systems [61][62][63][64][65][66][67][68][69] is applied to the initial GS density matrix, which then evolves according to the Liouville-von Neumann equation of motion in the non-orthogonal basis (1) which is numerically integrated. HereĤ denotes the hamiltonian matrix, S the overlap matrix andρ the density matrix.
For spin polarized systems both up and down spin density matrices are propagated within the description provided by the DFTB model. The time step for the integration was set to 0.005f s. 20000 steps of dynamics were done, simulating a total time of 100 fs. The intensity of the electric field was E 0 = 0.001 VÅ −1 . In the linear response regime the dipole moment can be calculated as shown in equation (2), where E(τ ) = E 0 δ(τ − t 0 ) and α(t − τ ) is the polarizability along the axis over which E(t) is applied.
The absorption cross-section is proportional to the imaginary part of the frequency dependent polarizability, which is obtained from the Fourier transform of eq. (2), leading to eq. (3) We have damped the signal with a damping factor of 0.1 fs −1 to obtain uniform broadening of the peaks. The absorption spectra are calculated averaging the polarizabilities over the three cartesian directions.
Triplet excitations where obtained by the same method but applying a magnetic instead of electric initial pulse.
The spin-unpolarized, spin-polarized singlet and spinpolarized triplet spectra are shown in Figure S5 for the TZZ flake of 264 atoms to show the importance of including the spin polarization 70 . Throughout the rest of this work, spectra of spin-polarized systems only consider singlet transitions, since triplet transitions have negligible intensity.

III. RESULTS AND DISCUSSION
A. Optical absorption spectra The optical absorption spectra of zigzag and armchair, triangular and hexagonal GNF, with sizes ranging from 2.9 nm to 7.3 nm, are shown in Figure 1. It has been reported that the lowest-energy excitation usually involves delocalized π electrons 19 . Supplemental Figure S1 shows a comparison of the total density of states (DOS) and projected DOS on p z orbitals for the states involves in these transitions, supporting the π character 70 . Further evidence of collective excitation can be obtained from the electrostatic potential that arises when a resonant field is applied resonant to the lowest-energy signal, plotted for the maximum and minimum dipole moment in Figure 2, for one structure of each of the four types. Section S2 (Supporting Information) contains the details of how this calculation is done 70 . As it is evident, charge polarization is spread in the surface for HAC, TAC and TZZ structures, while it is only delocalized in the edges for the HZZ flake.
A comparison of the main spectral features with results reported elsewhere can be enlightening. Cocchi et al. also calculated absorption spectra for GNF of the HAC family, albeit for smaller structures, using the ZINDO/S model based on configuration interaction including single excitations 44  eV and report two dark transitions that are nevertheless detected in experiments. We observe transitions at energies close to where these dark transitions should appear. A redshift of the lowest-energy excitation occurs as the GNF size increases, since the HOMO and LUMO approach as the electronic structure tends to that of bulk graphene. However, different trends in the peak intensity are observed: while for triangular ( Fig. 1 a and b) and armchair hexagonal (Fig. 1 c) flakes the intensity increases, for zigzag hexagonal structures (Fig. 1 d) it decreases. The first and second peaks for HZZ structures can be told apart given the fact that they follow different trends, as it can be drawn from Figure S2 70 . The intensity of the first peak for the two largest structures is so weak that it could not be detected, and hence could not be included in the forthcoming analysis. Figure 3 offers further insight on the observed trends. First of all, the relation between the excitation energy and the inverse square root of the number of atoms N is linear, which resembles the same type of relation for collective excitations in metallic nanoparticles with N −1/3 , namely, the surface/volume ratio 68 . Hence, for GNF the lowest excitation energies depend linearly on the perimeter/surface ratio, a typical quantum confinement feature. Secondly, a similar relation holds for the intensity as a function of the square root of the number of atoms, as a consequence of the increasing polarizability linked to a larger number of electrons involved in the excitation for large structures. However, this is not the case of HZZ flakes. The reason is described in the following paragraph.
The edge character of the excitation observed in Figure  2 for HZZ flakes suggests that only edge states participate in the transition; by inspecting the projected density of states per atom (see section III B of a description), the number of such states was found to be 12 one-particle states, half of them filled at 0K, for all the structures. A visual representation of such states can be found in Figure S3 70 . As a consequence, as the GNF size increases the dipole density decreases, since the a constant number of electrons is spread in an increasing edge perimeter. Visual evidence of this phenomenon is obtained by plot-ting the transition dipole density (TDD). To calculate the TDD, we consider only the two orbitals that contribute the most on this transition, φ initial and φ f inal , which are obtained as described in Section S5 and Figure S4 70 . A real space picture of the TDD is calculated as where r = (x, y, z). It is essentially a spatially weighted inner product between initial and final orbitals. In Figure 4, a comparison of ρ for two HZZ structures is done, where the same density isosurface value has been considered. The decrease in the transition dipole moment for the larger structure is evident, which in turn explains the decreasing absorption cross section.

B. Electronic structure
The main features that arise from time-dependent calculations are consistent with what is expected from electronic structure analysis. The density of states (centered at the Fermi level E F ) shown in Figure 5 was plotted with a Gaussian envelope at each eigenenergy ǫ i according to where ǫ denotes the energy and Γ is the standard deviation of the distribution (a value of Γ = 0.05 eV was used in all cases). A HOMO-LUMO gap around E F is found in armchair, caused by the same reason than in armchair GNR 19 . Zigzag triangular flakes show several "gap states" close to E F , which in spin-unpolarized calculations appear as zero-energy states whose multiplicity increases with size 18 . However it has been proven that these nonbonding states arise due to topological frustration (impossibility for all π bonds to be satisfied simultaneously). When the spin polarization is included, the asymmetry of spin breaks the degeneracy opening a small gap at the Fermi level 18,59 . Furthermore, the total spin ground state can be very large since it scales linearly with size 59 , as it is shown in Figure S6 70 .
Zigzag hexagonal flakes, on the other hand, show a continuous distribution of states with no clear energy gap. In the latter, the density of states close to the Fermi level increases as the flake becomes larger (Figure 5 d), giving rise to very distinct DOS profiles: a small gap is present for smaller N = 252 flakes (inset), but a continuous DOS is seen around E F for larger ones (N = 936). This phenomenon is visible for structures larger than ≈ 500 atoms, which are also magnetic. Transitions among these states are responsible for low-energy, low-intensity transitions in the absorption spectra.
In general, the calculated DOS are consistent with previous calculations using a simple nearest-neighbor, oneorbital Hückel model 22,29,30,39 and DFT 28,33 . For structures comparable to the ones analyzed here, usually a wide range of energies is considered in the literature, which have very similar DOS profiles among the different structures 53 . On the contrary, we observe a remarkable difference in the lowest electronic states for different type of edges. This implies that a close inspection of this range of energies is crucial for a correct understanding of the practically relevant optical features. The trends shown in Figure 3 are directly related to the dependence of the HOMO-LUMO gap with structural properties. In Figure 6 (a) the correlation between the gap energy and the inverse square root of the number of atoms is depicted. The general observation is a decrease in the gap energy with increasing the size due to quantum confinement 71 , also found in previous calculations 22,42,72 . This explains the red-shift in the absorption spectra. The higher excitation energies values seen in armchair with respect to zigzag flakes also holds when HOMO-LUMO gaps are considered instead of excitations (Figure 6), and agrees with previous calculations 38 . There is a linear relation between the HOMO-LUMO gaps and the excitation energies (see Figure 6 (b)), showing that the selfinteraction term of the excitation density is small, due to the spatial delocalization of these excitations 73 . A natural classification in three groups arises by observing this plot, suggesting that although armchair flakes have very similar electronic properties, zigzag flakes have intrinsic shape-dependent features. To explore the nature of low-energy states in zigzag GNF the projected density of states (PDOS) contribution per atom was calculated. The PDOS of orbital α in a non-orthogonal basis set is calculated according to where S and C are the overlap and the molecular orbital coefficient matrices respectively. The PDOS is analogous to the DOS but weighted by the factor |(SC) iα | 2 indicating the contribution of localized orbital α to molecular orbital i. We observe that carbon p z orbitals are the only ones that contribute to the total DOS close to the Fermi level 70 (see Figure S1), which explains the success of oneorbital tight-binding or Hückel models for the study of these systems 22,29,30 . Figure 7 shows the sum of the contributions of all orbitals of each atom to the state at (or close to) E F for zigzag triangular (a-b) and hexagonal (c-d) GNF. It is evident that such states are mainly localized on the edge atoms, and are hence edge states. For TZZ flakes, the lowest-energy absorption peak involves transitions from edge states to delocalized states, which give rise to a delocalized electrostatic potential profile as shown in Figure 2 (b). For HZZ structures, such transition involves only edge orbitals, as evidenced in Figure 2 (d).
Moreover, there is a strong difference in the magnitude of the contribution of each atom to the total DOS between zigzag triangular and hexagonal GNF. In the former, the PDOS per atom increases by two orders of magnitude from N = 264 to N = 936, while in latter the PDOS only doubles from N = 252 to N = 936. This difference in the hexagonal structures gives rise to the different profiles in the DOS plots shown in Figure 5 d. The origin of these states has been identified in GNR 74 and explained for GNF 18,22 . Due to the gap around E F in armchair GNF, this analysis cannot be done for such structures. It is clear that the frontier orbitals have large edge dependence in zigzag flakes. This is key for the understanding of the optical properties of such structures. Barnard and Snook 53 examine the possible reconstructions that corner atoms undergo in two corner terminations of armchair hexagonal GNF, analyzing the forces and the density of states in a wide range of energies. However, low-energy electronic aspects needs to be considered in order to explain optical properties. In Figure  8 (a) the spectra of two armchair hexagonal nanoflakes (similar in size but not equal) with two possible (zigzag) corner terminations, ending in one benzene ring (HAC-I) or two benzene rings (HAC-II) are shown. The spectra are similar up to the lowest-energy peak, but show different profiles for higher energies. The comparison between the size dependence of the excitation intensity of the first and third peaks is done in Figure 8 (b). Although they show similar intensity for the first peak, when the third peak is considered the intensities differ considerably. A larger absorption intensity is caused by a large transition dipole moment. To gain insight on the structural dependence of this stronger absorption, we calculated the transition dipole density for these strcutures. The same isovalue of ρ(x, y, z) for the orbitals involved in the third excitation is shown in Figure 9 for both structures analyzed previously. It is evident from the larger volume enclosed by the isosurfaces, that the transition density is highly polarized in the HAC-II structure (Figure 9 (b)), giving rise to the stronger optical absorption see in Figure  8.

IV. CONCLUSIONS
In this work we demonstrate that shape and edge geometry do matter when optical properties of GNF are considered. In fact, GNF with different geometries could be identified by their optical absorption spectra which show very distinct profiles and excitation energies. A classification in three families (A,B,C) arises from this analysis: armchair GNF belong to group A, zigzag triangular GNF to group B and zigzag hexagonal GNF to group C. Group A GNF show first excitation peaks with linearly increasing intensity and decreasing energy with respect to size, and higher excitation energies are expected since the HOMO-LUMO gap is larger. In group B the trends in intensity and energy are the same than group A, but the excitation energies are smaller because of gap states. In group C intensity decreases with size because of the edge dependence and constant number of involved states. Different corner terminations do not affect this classification since they become relevant when higher excitation energies of the spectrum are analyzed.
Given the structural polydispersity that is naturally achieved in prepared experimental samples of these nanomaterials at present days, studies that contribute to the understanding of relevant properties in different structural configurations are imperative. The design of either robust materials against disperse structural features, or samples with uniform properties or small property dispersity, are benefited from studies like this one. We show that, by performing a statistic analysis over the optical properties of a mixture of GNF structures, it could be possible to measure the proportions of GNFs with a certain shape and edge type, and tell apart two corner terminations. Further studies should be performed, nevertheless, to extend the characterization and provide useful insight for the design of novel applications.