Quantum mechanical analysis of nonlinear optical response of interacting graphene nanoflakes

We propose a distant-neighbor quantum-mechanical (DNQM) approach to study the linear and nonlinear optical properties of graphene nanoflakes (GNFs). In contrast to the widely used tight-binding description of the electronic states that considers only the nearest-neighbor coupling between the atoms, our approach is more accurate and general, as it captures the electron-core interactions between all atoms in the structure. Therefore, as we demonstrate, the DNQM approach enables the investigation of the optical coupling between two closely separated, but chemically unbound GNFs. We also find that the optical response of GNFs depends crucially on their shape, size, and symmetry properties. Specifically, increasing the size of nanoflakes is found to shift their accommodated quantum plasmon oscillations to lower frequency. Importantly, we show that by embedding a cavity into GNFs, one can change their symmetry properties, tune their optical properties, or enable otherwise forbidden second-harmonic generation processes.


I. INTRODUCTION
Owing to their ability to confine and guide light down to nanometer scale, the collective oscillations of electrons in conducting materials, known as plasmons, have generated great expectations for applications ranging from metamaterials, 1-3 quantum optics, 4,5 and photovoltaics 6 to photodetectors 7 and biological sensing. 8,9 Plasmons are commonly observed in noble metallic nanostructures, appearing as pronounced spectral resonances in their optical absorption and scattering spectra. However, noble metals suffer from relatively large ohmic losses, resulting in a limited plasmon lifetime. 10,11 In addition, metal plasmons can hardly be tuned, unless one uses metallic nanoparticles with different shapes, 12 thus severely limiting the operating spectral domain.
Graphene, a monolayer of carbon atoms arranged in a hexagonal lattice, 13 has emerged as a promising alternative to noble metals for nanoplasmonic applications due to its ability to support plasmons with unique properties, including large tunability, 14,15 long lifetime, 15,16 and high degree of optical confinement. 17,18 Strong intrinsic optical nonlinearity has also been observed in graphene and that can be further enhanced by plasmons. By using bottom-up chemical synthesis methods 19 or top-down electron-beam techniques, 20,21 nanometer-sized graphene nanoflakes (GNFs) with various shapes and sizes can be manufactured, providing a versatile platform to investigate plasmonic phenomena at the quantum level. In addition to inheriting remarkable physical properties from the extended sheet, the nanometer-sized graphene supports even more confined plasmons. Moreover, the plasmonic response of nanometer-sized graphene can reach the visible-light spectrum region a Author to whom correspondence should be addressed: Fangweiye@sjtu.edu.cn that is beyond the range of the plasmons in extended graphene sheets, thus extending our ability to manipulate visible light at deep-subwavelength scale. Intense research efforts have recently advanced our understanding of linear and nonlinear optical properties of nanostructured graphene. [22][23][24][25] For extended graphene, the nonlinear response is commonly described using a classical nonlinear conductivity derived from the Boltzmann transport equation, 26 assuming that intraband transitions dominate the optical response. Such classical electrodynamic description is invalid for nanometer-sized graphene, as the optical properties of nanostructured graphene are strongly influenced by nonlocal and finite-size effects. More recently, ab initio methods and the tight-binding (TB) description of the electronic states have been applied to investigate the optical response of the nanometer-size graphene. 27,28 The TB method is not accurate enough in general, since only the interaction of nearest-neighbor atoms is considered. Moreover, the TB method is not suitable for non-tightly bound structures, such as two interacting GNFs. Although ab initio techniques take into account the many-body interactions and therefore are more accurate and general than the TB approach, quantitatively accurate predictions of optical properties are computationally expensive, and therefore their application is limited to fewer atoms than those employing the TB method.
In our study, we use an all-atom coupling approach to calculate the linear and nonlinear polarizabilities of GNFs, which are the finite-size analog of linear and nonlinear optical susceptibilities of bulk optical media. In our approach, we assume that the π-orbitals of each atom are coupled to the core potential of all atoms. The electronic structure is then calculated, and a perturbative approach 29 is subsequently used to evaluate the linear and nonlinear polarizabilities. Importantly, our method inherently accounts for the symmetry properties of GNFs. For example, our method predicts that no second-harmonic is generated in GNF configurations invariant to inversion symmetry transformations as second-harmonic generation (SHG) is forbidden in such centrosymmetric structures. In order to illustrate the flexibility and generality of our method, we apply it to GNFs of different shapes and sizes, cavities in GNFs, and dimers of GNFs of different shapes.

A. Optical response calculation
Under illumination with an optical beam with frequency tuned to a frequency at which plasmons can be excited in the GNF, these graphene structures generate significant nonlinear optical response, which manifests as enhanced second-harmonic generation (SHG) and third-harmonic generation (THG). In this work, we study the nonlinear optical response of GNF with various shapes and topologies, namely, the linear and nonlinear polarizabilities of such structures. We model GNFs as a planar hexagonal distribution of carbon atoms with lattice constant a = 1.42 Å. For example, as illustrated in Fig. 1(a), a triangular GNF consisting of 141 carbon atoms has a side length of 2.2 nm. Considering that the localized electrons in the carbon-carbon σ-bond do not contribute significantly to the low energy optical response, only the π-bond forming p z orbitals are included in the quantummechanical calculation. That is, each carbon atom is represented by a single p z orbital oriented perpendicularly to the graphene plane, as per Fig. 1 In order to investigate non-tightly bound structures, we considered interactions between the p z electrons and all cores in the structure. The Hamiltonian operator for a single electron in the GNF is expressed asĤ where Z eff is the effective core charge, the first term describes the kinetic energy, and the second term represents the potential arising from the Coulomb interactions between the p z electron and all the other nuclei in the GNF. The electronic structure is computed using the Schrödinger equation where E and ψ( r) are the eigenenergy and eigenfunction, respectively. Then, we follow a perturbative approach to evaluate the linear and nonlinear optical response of GNFs (see supplementary material for more details about the quantum mechanical calculations). The structures investigated in this paper are a series of triangular and hexagonal GNFs [ Fig. 1(c)], and dimers of triangular and hexagonal GNFs [ Fig. 1(d)] to demonstrate inter-nanoflake optical coupling. We also study the optical response of hexagonal GNFs with an embedded cavity, as shown in Fig. 1(e). The cavity is located either at the center of the cavity or along the ΓM or ΓK symmetry axes of the cavity.

B. Shape and size effects in optical response of graphene nanoflakes
We assume that the incident electric field is polarized in a direction parallel to one of the sides of the triangle or hexagon, indicated as the x axis in Fig. 1(a), with all the structures considered in this work being assumed to lie in the x-y plane. We first consider hexagonal GNFs. Because they are centrosymmetric structures, SHG is forbidden in hexagonal GNFs, and thus we focus on their first-(linear) and third-order optical response.
In Fig. 2, we show the frequency dependence of linear polarizability and THG nonlinear polarizability of hexagonal GNFs, with the size of the GNF increasing from the top to the bottom panel, as illustrated in the middle panels. Figures 2(a)-2(f) present the real and imaginary parts of linear polarizability, α xx (ω), of hexagonal GNFs. The polarizability spectrum clearly shows that the peak of the imaginary part of the polarizability occurs at the spectral position where its real part becomes exactly zero, which is the condition that defines the existence quantum plasmons. For the smallest size, which is a single benzene, the linear polarizability shows a resonance at 6 eV [see Fig. 2(a)], which corresponds to the transition from the highest-occupied molecular orbital (HOMO) to the lowest-unoccupied molecular orbital (LUMO). With the structural size increasing, some additional resonances expectedly appear as more energy levels are involved into the optical transitions. As expected, our calculations predict that the magnitude of the polarizability increases when the size of the GNF increases, as the density of states increases with the number of atoms in the GNF. We also note that the peak energy of the plasmon resonance is continuously red-shifted with the increase in the structural size of the GNF, which is a consequence of finite-size effects of these GNF that is not included in their classical description but captured by our DNQM description. Importantly, as the size of the GNF increases, additional plasmon resonances can be seen in the spectra, which is similar to what can be observed in the classical regime.
Figures 2(g)-2(l) show the real and imaginary parts of the third-order polarizability, γ xxxx (ω), (also called second hyperpolarizability when referring to molecules) associated with the THG. Note that the nonlinear response of GNFs is also highly sensitive to the structural size, and with the increase in the size of the GNF, the magnitude and the resonance frequency of the nonlinear plasmon oscillations exhibit the same trend as in the linear case. Also note that both linear and nonlinear resonances can reach the visible-range spectrum as long as the size of the GNFs is small enough.
The results for triangular GNFs are presented in Fig. 3. The top panels show the structures of triangular GNFs under consideration with their size increasing from left to right. Unlike the centrosymmetric hexagonal structures, triangular GNFs are non-centrosymmetric, and consequently the second-order polarizability tensor, β, (also called hyperpolarizability) has non-zero components, meaning that the SHG is allowed. In particular, considering the symmetry properties of the triangle, β yyy 0, and thus this is the component we considered in our study.
The incident-frequency dependence of the linear polarizability α yy (ω) on triangular GNFs, calculated for GNFs of different sizes, is presented in Fig. 3(a). Note that increasing the structural size induces a red-shift for the plasmon frequency, similar to the case of hexagonal GNFs. The imaginary and real parts of nonlinear polarizabilities corresponding to SHG ( β yyy ), and THG (γ yyyy ) for triangular GNFs are shown in Figs. 3(b) and 3(c), respectively. Similarly to the case of hexagonal structures, the most pronounced peak in the spectra of nonlinear polarizabilities is red-shifted when the GNF size increases. However, the particular shape of the corresponding linear and nonlinear spectra of the triangular GNFs is markedly different from the analog ones calculated for the hexagonal GNFs, implying that the optical response of GNFs is strongly dependent on their geometrical configuration.
In Figs. 3(d)-3(f), we present an overview of the dependence of the first resonance peak energy for the first-, second-, and third-order polarizabilities, respectively, on the number of carbon atoms of a triangular GNF. Expectedly, the magnitude of the linear and nonlinear polarizabilities increases with the number of carbon atoms in the GNF. Additionally, the resonance peaks display an obvious red-shift with the increase in the number of carbon atoms, the steepest variation occurring in the linear response. In particular, when the number of carbon atoms increases from 13 to 198, the first resonance peak energy of the linear, second-order, and third-order polarizabilities varies by 2.45 eV, 1.25 eV, and 0.825 eV, respectively.

C. Graphene nanoflakes containing a quantum cavity
We next investigate the linear and nonlinear optical response of GNFs containing a quantum cavity. The cavity can be created at desired locations by removing one or more carbon atoms. As expected, the optical response of GNFs can be significantly altered by the presence of the cavity. In particular, the cavity can break the inversion symmetry of the GNF so that, as is the case with SHG, forbidden nonlinear optical interactions become possible. We illustrate this and other cavity effects in Fig. 4, by comparing the linear and nonlinear polarizabilities of two hexagonal GNFs (H3 and H6) that contain a cavity at different locations (see the cavity location in the inset). For a better comparison, the linear and nonlinear spectra of full GNFs are also presented. These calculations reveal several important conclusions. Thus, compared to the GNF of larger sizes (H6), the presence of a cavity into a smaller GNF (H3) produces more dramatic modifications of the spectra of polarizabilities. Interestingly, we note that a cavity introduced at the center of the H3 nanoflake induces a blue-shift of the linear and nonlinear plasmon frequencies, whereas a cavity located off center on the ΓK symmetry axis induces a red-shift [see the top panels in Figs. 4(a) and 4(b)]. Equally important, we find that the H3 flake has resonances in the visible domain, even when it contains a cavity.
For the larger GNF (H6), the spectral changes produced when a cavity is inserted in the nanoflake are relatively small; however, in the presence of the cavity, the second-order polarizability is no longer identical to zero. To be more specific, in the case of full hexagonal structures, SHG is strictly forbidden due to inversion symmetry, whereas in the presence of the cavity along the ΓK or ΓM symmetry axes, inversion symmetry breaking induces intense SHG when the fundamental frequency is near the plasmon resonances, as shown in Fig. 4(c). Of course, when the cavity is at the center of the GNF, the inversion symmetry of the structure is preserved, and thus all even-order polarizabilities of hexagonal GNFs identically vanish. This effect can be used as an efficient approach to engineer the nonlinear optical response of graphene structures. It can also have important practical implications, e.g., to molecular sensing. Thus, similar to the effect of the cavity, a molecule adsorbed by the GNF can drastically alter its symmetry properties and significantly change the nature of generated nonlinear optical signal.

D. Optical response of GNF dimers
Finally, we explore the effects arising from optically coupling two closely spaced GNFs. As two specific examples, we investigate the optical coupling between two hexagonal and triangular GNFs, as schematically shown in Fig. 5. Here we want to emphasize that the widely used tight-binding model does not describe the electronic states for GNFs dimers, as it does not allow coupling between atoms beyond the nearest neighbors, and accordingly is not suitable for the study of the interaction between two GNFs. By contrast, the DNQM approach we proposed captures interactions among all atoms of the structure, and thus it is well suited for the investigation of composite configurations, such as GNF dimers. Figure 5 presents the variation of linear and third-order nonlinear polarizabilities of hexagonal and triangular GNF dimers as a function of the spacing between the two comprising nanoflakes, D.  Note that the second-order nonlinear processes are forbidden in these two dimer configurations due to their inversion symmetry. As Fig. 5 shows, for large D values, the coupling between two flakes is negligible, and consequently the spectrum of the "dimers" is essentially identical with that of a single GNF. Only when D becomes small enough, the optical properties of dimerized GNFs begin to deviate from those of a single GNF. Remarkably, as can be seen from the linear polarizability spectra shown in Figs. 5(a) and 5(b), whereas only one pronounced peak exits for a single-GNF as seen for larger D, for small D, the interaction between the two plasmons of the GNFs induces an energy split that increases as D decreases. In addition, with the increase of D, the additional lowest-energy peak is blue-shifted for the hexagonal GNF dimer [ Fig. 5(a)] and is red-shifted for the triangular GNF dimer [ Fig. 5(b)]. The size of the energy splitting is an indication of the strength of the interaction between the GNFs. Using this measure, one can conclude that hexagonal GNFs interact much more strongly as compared to the triangular ones, a possible explanation of this finding being that more carbon atoms are in close proximity of each other than in the case of triangular GNF dimers. It should also be noted that the coupling between two GNFs can significantly alter the third-order nonlinear polarizabilities, too, as per Figs. 5(c) and 5(d).

III. CONCLUSIONS
In summary, we have computed the linear and nonlinear polarizabilities of GNFs using a distantneighbor quantum mechanical approach and revealed that the optical response of graphene nanoflakes depends significantly on their shape, size, and symmetry properties. In particular, the peak energy of the plasmon resonance is red-shifted, and the magnitude of the linear and nonlinear polarizabilities increases as the size of the graphene nanoflake increases. Significant changes in the optical response were demonstrated by introducing a quantum cavity in the graphene nanoflake. In particular, strong second-harmonic generation is enabled in hexagonal graphene nanoflakes by using a cavity to break the structural inversion symmetry. We have also explored the optical response of graphene nanoflake dimers and found that the strong coupling between two closely spaced graphene nanoflake leads to the energy splitting of the plasmon band, the magnitude of this splitting depending on the shape of the interacting nanoflakes. Importantly, the theoretical techniques and ideas developed in this study are not specific to graphene but can be extended to other two-dimensional (2D) materials, such as MoS 2 30,31 and black phosphorous 32,33 (one only needs to update in our method the corresponding Hamiltonian operator and basis functions of each atom for these new materials). Potential applications of our findings to molecular sensors have also been discussed.

SUPPLEMENTARY MATERIAL
See supplementary material for the supporting content.