Abstract
This paper presents a tight binding and ab initio study of finite graphene nanostructures. The attention is focused on three types of regular convex polygons: triangles, rhombuses, and hexagons, which are the most simple high-symmetry convex structures that can be ideally cut out of a graphene layer. Three different behaviors are evidenced for these three classes of compounds: closed-shells for hexagons; low-spin open-shells for rhombuses; high-spin open-shells for triangles.
1. Introduction
Extended layers of conjugated carbon atoms disposed in a 2-dimensional (2D) honeycomb lattice are the constituent of the common graphite and were regarded as a sort of intellectual curiosity constituting the model of low dimensional materials. Indeed the attitude toward such systems dramatically changed at the beginning of this century with the works of Novoselov et al. [1, 2]. As reported in their seminal paper they were able to produce and characterize a new allotropic form of carbon, constituted by the 2-dimensional sheet that is graphene. Soon after this work, an impressive number of research papers appeared dealing with the structural and electronic properties of graphene. Indeed, because of the robustness of the graphene skeleton, this allotropic carbon form happens to be one of the strongest materials ever produced; consequently its use as reinforcer in novel high-performance composite materials becomes straightforward. Notably its remarkable electron transport properties and the fact that it is a zero-gap semiconductor are making it one of the materials of choice for future applications in molecular electronic devices [3]. Technical applications of graphene can be related to extremely diverse technological areas and, for instance, can be related to single molecule detection, to field effects transistors of even to quantum information processing. From a more fundamental science point of view it is remarkable that graphene allowed the predictions of quantum electrodynamic to be tested in a solid state-system because of its unusual linear electron excitations. Remarkably enough, electronic excitation close to the Fermi level strongly reminds us of the one exhibited by massless Dirac fermions. The 2D structure of graphene can also be exploited as a support to depose atoms, molecules, or aggregates; in some cases it has been shown that the interaction with graphene was able to significantly change their properties. These few and simple considerations certainly can justify the extremely high interest paid nowadays to graphene chemistry and physics. But unusual and remarkable properties are not only the domain of extended and formally infinite graphene sheets. Indeed, finite nanostructures can exhibit novel and fascinating properties due to the confining effects and the finite size. First of all graphene nanostructures, ribbons or island, show an accumulation of electronic density at the border, giving raise to the so-called edge states [4]. The presence of such states can strongly modify the physical and chemical behavior of such materials and has been demonstrated by different research groups, both theoretically and computationally. Moreover a very strong correlation between the geometry of the graphene nanostructure and the properties of their edge-state exists, and the global shape of the island indeed can modulate its electronic properties and structure in a very impressive way. For instance, the very strong difference in the properties of armchair and zig-zag edges structures has been the subject of quite a number of interesting publications. It is useful to cite here that the perspective of achieving a high control on the precise shape of graphene nanostructures and also of the type and nature of their edges not only is a theoretical fascinating suggestion but is becoming a sounding reality ready to be exploited. Indeed in the early years of the, yet still recent, graphene era, nanoscale systems were built by a top-down approach basically consisting in stripping out graphite sheets. Obviously such a procedure offered very limited, or even inexistent, control on the shape of the product. Today, on the contrary, we may assist to the emergence of much more sounding so-called bottom-up techniques in which nanoislands are built by precise deposition on metallic surfaces. These techniques generally allow a very good control on the shape of the resulting product that can be built with atomically defined edges. Quite recently, as suggested by Fürst et al. [5], the possibility of creating periodically perforated graphene structures (the so-called graphene antidot lattice) has been evoked to obtain even much more enhanced and diverse electronic and magnetic properties. Indeed complex graphene structures, having hexagonal, rhomboidal, or triangular shape, can exhibit extremely complex, yet predictable, electronic, and magnetic properties that make them ideal candidates, to rationally design organic giant nanomagnets, also exploiting the possibility offered by the different connectivity of different graphene units [6, 7]. These kinds of structure have recently attracted both experimental and theoretical consideration. For instance, Lin and coworkers [8] proposed a new route of synthesis of hexagonal and triangular (and also rectangular) graphene nanoflakes from carbon nanotubes and characterized them using photoluminescence. Pelloni and Lazzeretti [9] looked at the density response to an external magnetic field of triangular, rhombus, hexagonal, and rectangular structures using the polygonal current model [10]. Jamaati and Mehri [11] studied the influence of the edge length on the conductance of rhombuses nanoflakes. Mansilla Wettstein and coworkers [12] used a DFTB approach to simulate absorption spectra of triangular and hexagonal graphene nanostructures and evidenced the impact of the shape of the edifice. Many of these interesting properties belonging to finite graphene nanostructure arise from the presence and nature of their edge states [13]. Indeed the electronic density accumulation to the borders gives raise to a high density of electronic states close to the Fermi level. This in turn induces a high multireference character that may results in the presence of a very rich low energy electronic spectrum. Different states of different multiplicities can indeed happen to be the ground states, while even the lowest multiplicity ones can be described as open-shell states, with the electron close to the Fermi level being unpaired, and coupled ferromagnetically or antiferromagnetically to the others. The spin multiplicity of the actual ground state and hence the ferromagnetic or antiferromagnetic nature of the nanoisland will be dictated by the shape of the nanostructure. An easy way to formalize their magnetic properties is to recognize that the graphene lattices can be decomposed in two sublattices; therefore one will have to deal with and carbon atoms. Note that centers will only have atoms as first neighbors and vice versa. This implies that, for instance, zig-zag edges are all composed of the same type of atoms, while armchairs edges alternate between and atoms. Subsequently one should recall Lieb’s theorem [14] that states that the spin multiplicity of a given structure will be given by the balance between the number of atoms ( or ) belonging to the two sublattices Note that in case one uses the simple tight binding approach or single-reference treatment with a minimal basis set, the difference will also be equal to the numbers of degenerate eigenvalues lying at the Fermi level, that is, whose energy is exactly zero. Note also that Lieb’s theorem implies that two-carbon center will be ferromagnetically coupled if they belong to the same lattice and antiferromagnetically coupled if they belong to the opposite. Finally it has to be noted that magnetization arising from edge states in case of large islands and because of the high number of states close to the Fermi level may give rise to spin instability. Recently some of us have shown that coherently with Lieb’s theorem triangular structures have high-spin ground states and are ferromagnetic, while hexagonal and rhomboidal ones are open-shell low-spin structures of very high multireference character. Usually graphene nanostructures are studied at semiempirical level, or when at ab initio level using DFT methodologies, also because of the size and computational cost required to achieve large cluster sizes. On the contrary the presence of high-correlated structures giving rise to a very rich and subtle low energy spectrum, with an impressive magnetic behavior, would call for a wavefunction based multireference treatment. In this contribution triangular, T_n, hexagonal, H_n, and rhomboidal, R_n, structures, where is the number of hexagons per side, will be studied at ab initio level using multireference perturbation theory formalism, a level of theory that represents the best balance between accuracy and computational cost in treating both static and dynamic correlation. On the other hand the energy band structure will be reported also at semiempirical (Hückel) level. This will allow reaching a very high cluster size and, therefore, achieving the limit of infinite structures. The different character of the density of states close to the Fermi level will be considered also to easily rationalize the emergence of magnetic properties. In particular the presence of a gap at Fermi level, whose value strongly diminishes with the system size, will be underlined in the case of low-spin systems, while a continuous band structure is shown for the other structure. Moreover the evolution of the magnetic coupling with the island shape and size will also be particularly taken into account. In Figures 1, 2, and 3, the largest structures considered in this work, for each of the three different types (i.e., H_6, R_7, and T_7), are illustrated. and sublattices are highlighted in different colors. It appears that the hexagons (H) structures have , to which a singlet ground state corresponds. Also the rhombus (R) structures are associated with a zero value of , but this happens by compensation between the two different equilateral-triangle halves of the rhombus. This explains the singlet open-shell character of the corresponding ground states. Finally, the triangles (T) have a nonzero value of and a high-spin open-shell ground state.
2. Computational Details
We will consider three different topologies of graphene nanoislands built with polygon’s shape that can be described as follows: triangles (T_n), hexagons (H_n), and rhombuses (R_n), being the number of benzene units by sides of the nanoisland. We chose these structures since they are the only flat convex islands having all sides of the same length. Those structures are derived as fragments of an infinite and graphene sheet and are therefore composed by a honeycomb network of carbon atoms, while the dangling bond at the border has been saturated by hydrogen atoms (see Figures 1, 2, and 3). Due to their simple polygonal shapes, all the structures considered in this work have a high degree of symmetry, with a space group of for triangles, for hexagons, and for rhombuses, respectively.
The electronic structure of those compounds has been investigated also by ab initio calculations. Due to the possible multireference character of these highly conjugated systems, static correlation is taken into account using Complete Active Space Self Consistent Field (CAS-SCF) level of theory, and effects of dynamic correlation are included using multireference second-order perturbation theory and in particular N-Electron Valence Perturbation Theory (NEVPT2) method. To enlarge the size of the cluster toward nanometric size and thus reaching the limit of extended systems, we used a simple semiempirical approach based on the tight binding formalism. At ab initio level structures with going from 2 to 7 are considered for triangles and rhombuses, while for hexagons, only and are taken into account, due to computational burdens. The computational methods used as well as the strategies employed along the study will be described in the following subsections.
2.1. Tight Binding
Tight binding calculations have been performed using a dedicated code written by the authors which has already been described in a previous paper [4]. The Hamiltonian matrix, based on the Hückel formalism, has been constructed from the connectivity of every carbon atom in each nanoisland. As usual in Hückel type calculations only the subsystem has been included in the Hamiltonian while the electrons have been disregarded.
The Hamiltonian matrix () elements have been constructed following the usual assumptions in particular:where represents the self-interaction and has been set to zero without loss of generality, while represents a topological matrix that gives the interaction between different sites. The interaction matrix elements have been set equal to zero if the sites and are not connected, while they have been set equal to 1 for directly connected sites. This means that only interactions between first neighbors are considered. Since hydrogen atoms only participate through the network they are not considered at tight binding level. In all the different nanoislands considered in this work, the optimized C-C bond length happens to be rather similar, regardless of the actual shape or size (see the geometries in Supplementary Material). For this reason, we assumed the same value of for all the connected pairs of atoms. In particular, without loss of generality, the parameter was fixed at in all our Hückel calculations.
2.2. Ab Initio
The smallest systems were treated at CAS-SCF [15–17] and multireference NEVPT2 [18–20] level; this represents one of the first ab initio highly correlated multireference treatment of significant size graphene nanostructures.
The minimal STO-3G [21–23] basis set for both carbon and hydrogen was mainly used in this work. Of course, one can not expect to obtain quantitative results using this minimal basis set paragon, but it has been widely used to get qualitative trends for many families of molecules, both organic and inorganic. In particular, its ability to give the right tendencies concerning the behavior of low dimensional systems has been evidenced by some of us [4, 24]. In order to further assess this statement, we compared the Hartree-Fock and NEVPT2 energies (at the STO-3G geometries) for the H_2, R_3, and T_3 nanoflakes using STO-3G and different extended ANO basis sets optimized by Widmark and coworkers [25]. Thus, we tested basis sets of VDZ (), VTZ (), and VQZ () quality. As can be seen in Table 9 of Supplementary Material, the absolute STO-3G energies are of course significantly different from the large basis sets ones, but the relative energies of the three systems only differ by a few percents. Also, the STO-3G correlation energies are systematically underestimated by 0.10 to 0.15 hartree from the smallest to the largest extended basis set. However, the global picture remains unchanged, and the reduced size of the STO-3G allows the investigation at a correlated ab initio level systems whose size is relatively large. Furthermore, in order to evaluate the basis size effect on the quality of our results, the two smallest structures in each series were also studied using the VDZ ANO basis set. Indeed in the case of T_2, R_3, and H_2 a contraction for C and for H has been used. The small differences (a few mhartree) that can be seen by comparing the corresponding Hartree-Fock results of Tables 4 and 10 with those of Supplementary Material reflect very tiny changes in the geometries and further support the use of the STO-3G for the present study.
For all the different clusters we used the following protocol: as a first step a restricted open-shell Hartree-Fock (ROHF) calculation was performed to optimize the orbitals for the high-spin wavefunction that can be described by only one spin-adapted Slater determinant. The high-spin wavefunction was constructed by singly occupying all the orbitals that at Hückel level have been found to be quasi degenerate and close to the Fermi level. These same orbitals and electrons also constitute the active space for the following multiconfigurational treatment. The procedure is trivial for the H_n clusters, since all the considered hexagons have a closed-shell wavefunction. It is also straightforward for the T clusters, where the tight binding calculations are able to give information on the number of open-shell orbitals. Afterwards, a calculation on an ionic system, having as many electrons removed as the number of singly occupied orbitals, gives a set of quasi degenerated low-lying orbitals that describe precisely the required open-shell manifold. Things are less trivial for rhombuses, since the definition of the open-shell manifold is ambiguous. In this case, we performed a series of high-spin ROHF calculations using no symmetry constraints, in order to find the spin corresponding to a minimum of the energy, and we chose these orbitals and multiplicity to define the open-shell orbitals manifold. The energies of the different ROHF calculations, in the case of the STO-3G basis set, are reported in Table 1.
Our experience with this kind of systems is that one obtains, in this way, orbitals that are close to the CAS-SCF ground-state orbitals, at a much lower computational cost [4, 26]. Such a procedure is often used in order to find optimum orbitals in magnetic systems, and it gives MO that are of CAS-SCF quality if the occupation numbers in the open-shell manifold are close to one. Once the orbitals for the ground state have been determined, we optimized the geometries at the R(O)HF level using the corresponding module [27] of the MOLPRO quantum-chemistry package [28], by choosing as reference wavefunction these high-spin Slater determinants. Due to restrictions of the MOLPRO package all the calculations were done using abelian point groups but imposing , , and symmetries on the hexagons, rhombuses, and triangles coordinates, respectively.
Subsequently, this high-spin wavefunction was used as a guess for the CAS-SCF treatment to obtain the solutions for the other low-lying states of different multiplicity. Note that in the post-HF calculation we did not optimize again the molecular orbital coefficients; therefore we can say that our results are actually of CAS-CI type. This choice can be advantageous when in presence of strictly degenerate orbitals, since it avoids CAS-SCF instabilities that can burden the iterative procedure, making the convergence quite hard to be reached. After the static correlation recovered with the CAS-CI step using the multimodule of MOLPRO [16, 17], the effect of dynamical correlation has been added on each state using the partially contracted NEVPT2 code [20] as available in the MOLPRO package. Note that these energies were systematically very close to the strongly contracted ones, differing by at most a few hartree. In the case of the H_2, R_4, and T_3 systems, we also checked that enlarging the active space by adding two of four orbitals changes the correlation energy by at most five mhartrees (see last table of Supplementary Material).
3. Results and Discussions
We present and discuss here separately the tight binding (i.e., Hückel) and ab initio results. In both cases, the behavior of the three types of structures is clearly defined. The hexagons H_n are closed shells, with a large gap between the Highest-Occupied Molecular Orbital (HOMO) and the Lowest-Unoccupied Molecular Orbital (LUMO) energies. The rhombuses R_n have a set of quasi degenerated (QD) orbitals, with a very small HOMO-LUMO gap. At ab initio level, these orbitals have an occupation that is very close to one. However, their ground state is a singlet, in accordance with the “Ovchinnikov rule” [29]. Finally, the triangles T_n have also a set of degenerated or quasi degenerated orbitals at the Fermi level, but they have a high-spin ground state at ab initio level.
3.1. Tight Binding
In Figures 1, 2, and 3, the largest structures for each type we considered are shown. The two different networks are evidenced. In general, 60-degree corners in the nanoisland are associated with an unbalanced number of one type over the other one, the type of the apical carbon being the less-abundant one. This explains the high-spin nature of the wavefunctions for triangles. In the case of corners having 120-degree angles, on the other hand, we have overall the same number of atoms of each network type. However, the two networks exchange their role on the two edges shared by the corner. Therefore, the wavefunctions of hexagons are closed shells. Rhombuses, on the other hand, have a zero-excess of one type over the other one as a result of compensation between the two triangles that form the rhombus, and their wavefunction is an open-shell singlet. The orbital spectrum for the H_5, R_8, and T_10 systems is plotted in Figures 4, 5, and 6. Obviously these systems present a different number of orbitals; thus the orbital numbers have been normalized in the segment for an easy comparison. As expected from Lieb’s theorem, the number of zero energy orbitals strongly differs from one type of system to the other. For the hexagons, being closed-shell systems, there is a clear gap around the Fermi level. For the triangles, one finds exactly such orbitals which correspond to the difference of type (75) and (66) carbon atoms. For the rhombuses, we have two quasi degenerate orbitals resulting from local excess of type and atoms at the 60° edges. For large enough edifices (as it is the case here), the interaction between these two orbitals vanishes and they become degenerate at the Fermi level. This is very similar to the behavior of the open rhombuses studied in our previous paper [4]. It can be seen that for H_5 the coincidence with the extrapolated infinite graphene case (red curve on the plots, obtained from the analytical solution for a graphene sheet with periodic boundary conditions of size 100 × 100 hexagons) is almost perfect. For the other structure, some larger discrepancies arise in the Fermi region. In Figure 7, the tight binding energy gaps are reported for the three types of considered structures. Notice that for the triangles T_n the gaps are defined by discarding the strictly zero set of magnetic Singly Occupied Molecular Orbitals (SOMO), or the gaps would be strictly zero. Since all these clusters converge to an infinite graphene layer (a zero-gap structure) in the limit , it is not surprising that the gaps of the three structures go to zero when is increased. However, for small values of , the difference among the three types of systems is quite remarkable. In Figure 8, the same gaps as a function of are reported using a logarithmic energy scale. The three behaviors are linear to a very good accuracy, suggesting for the energy gaps an exponential law of the type The values for and for the three structures have been optimized through a nonlinear fitting procedure using the GNUPLOT package [30], and the results are reported in Table 2. Although the three laws are exponential, the values of are extremely different, and this explains the qualitative different behaviors among the different types of systems.
3.2. Ab Initio
We present in this section separately the results for the three types of structures. All structures are very rigid, with C-C and C-H bonds varying only a few percent as a function of the bond position in the cluster. In fact, the C-C bond lengths go from a minimum of 2.512 to a maximum of 2.811 bohr (see Supplementary Material for the detailed geometries), while the corresponding C-H values are comprised between 2.046 and 2.049 bohr.
3.2.1. Hexagons
Because of the small size of the hexagon structures, the open-shell character on large systems resulting from local magnetization on the edges [31] does not appear in our results. In Table 3, the HOMO and LUMO orbital energies of the structures ranging from H_2 to H_4 are reported.
The closed-shell nature of these structures is clear, and this can be confirmed by CAS-SCF calculations (not reported in the present work). The frontier orbitals tend to be concentrated on the cluster edges, as can be seen in Figures 9, 10, 11, and 12 for the case of H_4. However, the effect is not very pronounced for these values of . It is only for much larger structures that the frontier orbitals become localized on the edges of the cluster, giving rise to the open-shell character of these systems as shown by Hubbard-type and DFT calculations [31]. In Table 4 the optimized total energies are reported. At HF level, the energies are very well fitted by the straightforward expression where and are the total number of carbon and hydrogen atoms, respectively. A best fit on the HF energies gives optimal values and , with an error of about one mhartree only on the total energies. In this expression, all the carbon atoms are treated on the same foot, as covalent-bonded atoms, and this fact confirms the closed-shell character of the wavefunction in the considered range of values of .
3.2.2. Rhombuses
In Table 5, the orbital energies and occupancies of the open-shell manifold are reported. The energies are comprised in an extremely narrow range of values, showing a picture that is essentially the same as the one obtained at tight binding level. In Figures 13, 14, 15, and 16, the open-shell orbitals for the R_6 system are shown. The orbitals are concentrated in the system edges, with a preference for the regions of the 60-degree corners. In Table 6, the ROHF energies on the different spin multiplicity with respect to the lowest one are reported. As explained in the Section 2.2, these calculations were done in order to determine how many and which orbitals are to be put in the active space. Once the active space orbitals have been chosen, calculations on all spin multiplicity and all symmetries were performed. The results are reported in Table 7 (absolute energies) and Table 8 (energy differences). The singlet is systematically found to be the ground-state spin multiplicity, in accordance with the Ovchinnikov rule, at least as far as the CAS-CI calculations are concerned. In the case of NEVPT2 results, it appears that the “magnetic manifold” (the neutral states that differ for different spin couplings) is usually well described by this formalism. Things are different for the ionic (as revealed by the CAS-SCF wavefunction) states that in graphene nanoislands lie systematically above the magnetic manifold. The NEVPT2 formalism can, in fact, strongly overcorrect their energies, and in many case they can even be placed below the magnetic manifold. This seems to be an artifact of perturbation theory, and in our opinion the NEVPT2 results for the ionic states should be taken very cautiously.
3.2.3. Triangles
In the case of triangles, the magnetic manifold is easier to be found, because it is unambiguously defined at tight binding level. The orbital energies, computed at RHF level, of triangular clusters are given in Table 9, while in Figures 17, 18, 19, 20, 21, and 22 these magnetic orbitals are shown for the case T_7. It is important to stress the fact that, contrary to what happens in the case of rhombuses, these magnetic orbitals have a real physical meaning, since they correspond to a net spin density of the system ground state. They are located on the edges of the systems, and therefore it is in this same region that the ground-state spin excess will be located. The orbital spectra for the case T_2 to T_7 are given in Figure 23.
The magnetic states are comprised in a very narrow range, less than ten kcal/mol, as can be seen in Table 10. Table 11 clearly shows that he high-spin states are systematically the lowest ones, again in agreement with the Ovchinnikov rule. However, the total energies, reported in Table 10, do not seem to follow a simple Heisenberg picture. Work is in progress in order to associate a magnetic Hamiltonian to these patterns of energy levels.
4. Conclusions
We presented in this paper a tight binding and ab initio (CAS-CI and NEVPT2) study to evaluate the electronic structure of graphene nanoislands. The calculations were carried out on three types of nanoislands of different size: triangles, rhombuses, and hexagons. They represent three broad classes of possible behaviors. In fact, at the size we explored in this work, hexagon islands have a closed-shell ground state and a rather large HOMO-LUMO gap at one-electron level. Notice, however, that things change for much larger hexagonal islands, where the ground state becomes a singlet open-shell. Rhombuses and triangles have an open-shell ground state, with several energy levels close to the Fermi level. However, rhombuses have a low-spin singlet ground state, while triangles have a high-spin ground state whose multiplicity grows linearly with the linear size of the system. These results that hold systematically at tight binding and ab initio level are in perfect accord with the Ovchinnikov rule.
The geometries of the systems were relaxed and optimized, imposing however the symmetry constraint for hexagons, for rhombuses, and for triangles. The bond lengths for the C-C and the C-H bonds lie in very narrow ranges, implying that these structures are characterized by a rather rigid graphene skeleton.
We are currently investigating how these structures can be used to build up more complex structures having remarkable and potentially interesting electronic properties. From this point of view, triangles and rhombuses are the most promising ones, because of their quasi degenerate manifold of low-lying states.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors thank the University of Toulouse and the French CNRS for financial support. The ANR-DFG action ANR-11-INTB-1009 (Project MITLOW) and the Erasmus Mundus program (Master TCCM) are also gratefully acknowledged. This work was granted access to the HPC resources of CALMIP under the allocation 2011-p1048. Thierry Leininger and Stefano Evangelisti thank the support of the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement no. 642294.
Supplementary Materials
The supplementary material file contains the optimized geometries for the hexagonal (Tables 1–3), rhombuses (Tables 4–6), and triangular (Tables 7 and 8) structures. Table 9 shows the basis set effects on the energies at the STO-3G geometries for the H_2, R_3, and T_3 systems. The last table presents the active space size effects on the NEVPT2 energies for the R_4 and T_3 nanoislands. Table 1: C=C and C-H distances and angles in hexagonal nanoislands. Table 2: C=C and C-H distances and angles in hexagonal nanoislands, Part 2. Table 3: C=C and C-H distances and angles in hexagonal nanoislands, Part 3. Table 4: C=C and C-H distances and angles in rhombus nanoislands. Table 5: C=C and C-H distances and angles in rhombus nanoislands, Part 2. Table 6: C=C and C-H distances and angles in rhombus nanoislands, Part 3. Table 7: C=C and C-H distances and angles in triangle nanoislands. Table 8: C=C and C-H distances and angles in triangle nanoislands, Part 2. Table 9: basis set effects for the smallest edifices (energies in hartree at the STO-3G optimized geometry). Table 10: active space size effect on NEVPT2 energies (in hartree). (Supplementary Materials)