Properties of Bilayer Graphene Quantum Dots for Integrated Optics: An Ab Initio Study

Due to their bandgap engineering capabilities for optoelectronics applications, the study of nano-graphene has been a topic of interest to researchers in recent years. Using a first-principles study based on density functional theory (DFT) and thermal DFT, we investigated the electronic structures and optical properties of bilayer graphene quantum dots (GQDs). The dielectric tensors, absorption spectra, and the refractive indexes of the bilayer GQDs were obtained for both in-plane and out-of-plane polarization. In addition, we calculated the absorption spectra via time-dependent DFT (TD-DFT) in the linear response regime. The TDDFT results show that a blue shift occurs in the absorption spectrum, which is consistent with the experimental results. In this investigation, we consider triangular and hexagonal GQDs of various sizes with zigzag and armchair edges. Our simulations show that unlike monolayer GQDs, for which light absorption for out-of-plane polarization occurs in the ultraviolet wavelength range of 85–250 nm, the out-of-plane polarization light absorption peaks in the bilayer GQDs appear in the near-infrared range of 500–1600 nm, similar to those in bilayer graphene sheets. The out-of-plane polarization light absorption peaks in the near-infrared range make bilayer GQDs suitable for integrated optics and optical communication applications.


Introduction
Optical absorption spectra are essential engineering tools for analyzing and designing optoelectronic devices. Although most studies use incident light with in-plane polarization, studying optical properties with out-of-plane polarization can be useful for optical integrated circuits because the light in such circuits propagates parallel to the surface. The need for optical absorption spectra and out-of-plane polarization in optical integrated circuits doubles the importance of computing the electronic and optical properties of bilayer systems with out-of-plane polarization.
The experimental study of the optical properties of graphene began in 2004 with the work of Novoselov et al. [1,2], which was preceded by some theoretical studies. As long as graphite has been studied, so has monolayer graphene. Characterizing the optical behaviors of graphene is essential for developing its applications in fabricating new optoelectronic devices. Knowing the details and locations of the peaks of the optical absorption and distribution spectra, transmission coefficients, reflection, optical absorption, and other linear and nonlinear optical properties has allowed us to create suitable platforms for designing and constructing detectors, modulators, switches, waveguides, light sources, and other optoelectronic devices.
As such, the past 16 years have witnessed the publication of valuable works on characterizing and computing the electronic and optical properties of graphene, such as (i) calculating the optical absorption of graphene and its dielectric coefficient tensor [3], (ii) calculating the transmission coefficient of visible light through graphene [4], (iii) measuring the optical conductivity of graphene [5], (iv) determining how gate voltage affects optical transitions in graphene [6], (v) measuring the saturated optical absorption and microwave absorption in graphene [7,8], (vi) measuring the optical absorption in graphene waveguides [9], (vii) determining the bending and rippling effects on the dielectric coefficient of graphene [10], and (viii) calculating and measuring the optical impact of holey graphene with different nano-mesh geometries [11].
One of the main demands in designing optoelectronic devices is to create an energy gap in graphene, thereby allowing its optical properties to be engineered and the dark current to be minimized. In this regard, nano-graphene began to be considered in research on producing a bandgap in graphene. Since 2006, graphene nanoribbons (GNRs) have been considered [12] as a research topic and have been studied seriously and widely, both theoretically and experimentally. However, compared to GNR, graphene quantum dots (GQDs) [13] have a greater variety of engineered factors and are far more attractive for optoelectronic devices. The optical properties of GQDs can be tuned by changing: (i) the GQD's geometry (circular, hexagonal, or triangular), (ii) the edge type (zigzag or armchair), (iii) the number of layers, (iv) the symmetry of the number of sublattice A and B, and (v) other factors.
Graphene oxide and reduced graphene oxide, like graphene quantum dots, have great potential for various applications that require bandgap engineering [14]. Graphene oxide is a graphene plate that includes oxygen atoms-O-and a functional oxygen group-OH-in addition to carbon atoms [15]. The bonding of these oxygen and functional atoms of the oxygen group with carbon atoms at the edges and in the middle region of graphene is covalent and established by a mixture of sp 2 and sp 3 hybridizations. The high surface to volume ratio in graphene oxide is another feature. Compared to graphene quantum dots, due to differences in its electronic structure, graphene oxide has the property of fluorescence in a wide range of wavelengths [16]. The modification and manipulation of chemical composition is one of bandgap engineering methods used for graphene oxide [17,18]. Apart from the optical detection properties of graphene oxide, the electrochemical applications of this material are of greater importance.
There have been numerous studies on the electronic and optical properties of monolayer GQDs [14,[19][20][21][22][23][24][25][26][27][28][29][30][31] and some studies on bilayer graphene sheets and bilayer nano-graphene with different models. A double-layer graphite lattice was first studied in 1992, and its electronic properties were obtained. Double-layer graphite produces features identical to those of bilayer graphene [32]. The bandgap of bilayer graphene nanoribbons is also a topic of research interest [33], as is the fact that applying a gate voltage can engineer a bandgap [34]. In recent years, the energy bands of graphene quantum dots with AA and AB stacking variants have been calculated with tight-binding models [35]. In the recent work of Mora-Ramos et al. [36] on the electro-optical properties of twisted bilayer GQDs, bilayer GQDs have become a more attractive material in optoelectronics.
In the present work, we investigate how the polarization of incident light affects the optical properties of monolayer and bilayer GQDs via density functional theory (DFT) and thermal DFT (th-DFT). Previous work has been unclear on the optical behavior of monolayer and bilayer GQDs when the electric field is polarized perpendicular to the graphene plane. The relevant aspects of using this type of polarization have not been specified for device designers. Therefore, in this paper, we calculate and compare the absorption spectra for the in-plane and out-of-plane polarizations of different GQDs systems and then investigate how the bilayer structure affects the optical characteristics. The results show that bilayer systems can absorb parallel incident light from visible to infrared (IR), whereas monolayer systems absorb only vertical incident light. Therefore, bilayer systems could be suitable for integrated photodetectors in which the light propagates laterally.

Computational Methods
Standard density functional theory (DFT) calculations were carried out by solving the Kohn-Sham equations using the plane-wave module in the Quantum ESPRESSO software package [37,38] and the local density approximation functional for estimating the exchange correlation [39]. The non-local pseudopotential was used with the norm-conserving specifications [40,41]. The applied pseudopotential is known as the Perdew-Zunger pseudopotential and is generated based on the von Barth-Car method [42]. By examining the previous calculations performed by other researchers using the tight-binding model [43], we selected the desired structures among the GQDs to achieve the best optical absorption in optical communication applications and photonic integrated circuits and then performed DFT calculations on those structures. Selecting a GQD structure involves determining the geometric shape, size, number of constituent atoms, type of side edges, and the symmetry of the number of sublattices of A and B. Figure 1 schematically shows the different atomic structures of monolayer GQDs. The vacuum layer around the unit cell was assumed to be around 10 Å thick and was used to create an isolated system (non-periodic).

Computational Methods
Standard density functional theory (DFT) calculations were carried out by solving the Kohn-Sham equations using the plane-wave module in the Quantum ESPRESSO software package [37,38] and the local density approximation functional for estimating the exchange correlation [39]. The non-local pseudopotential was used with the norm-conserving specifications [40,41]. The applied pseudopotential is known as the Perdew-Zunger pseudopotential and is generated based on the von Barth-Car method [42]. By examining the previous calculations performed by other researchers using the tight-binding model [43], we selected the desired structures among the GQDs to achieve the best optical absorption in optical communication applications and photonic integrated circuits and then performed DFT calculations on those structures. Selecting a GQD structure involves determining the geometric shape, size, number of constituent atoms, type of side edges, and the symmetry of the number of sublattices of A and B. Figure 1 schematically shows the different atomic structures of monolayer GQDs. The vacuum layer around the unit cell was assumed to be around 10 Å thick and was used to create an isolated system (non-periodic). The cut-off energy (Ecut) for determining the range of the central atomic core was calculated to be around 80 Ry. These values were obtained by applying optimization to minimize Etot. Since this research was performed with special attention to applications in photonic integrated circuits and the directions of lateral incident light, we used bilayer structures to achieve a suitable cross-section in a course with horizontal incident light. In bilayer nanostructures, it is also possible to tune the energy band gap by applying an external electric field. Figure 2 shows a bilayer GQD system with an interlayer distance d of 3.334 Å, which agrees with the systems used in other works. The cut-off energy (E cut ) for determining the range of the central atomic core was calculated to be around 80 Ry. These values were obtained by applying optimization to minimize E tot . Since this research was performed with special attention to applications in photonic integrated circuits and the directions of lateral incident light, we used bilayer structures to achieve a suitable cross-section in a course with horizontal incident light. In bilayer nanostructures, it is also possible to tune the energy band gap by applying an external electric field. Figure 2 shows a bilayer GQD system with an interlayer distance d of 3.334 Å, which agrees with the systems used in other works. The convergence threshold for the self-consistent-field (SCF) calculation loop is considered as 10 −8 Ry. Since the energy bands that are obtained for the GQDs are of a "flat" type, we used the gamma-point algorithm to choose the k-point. Algorithms such as the Monkhorst-Pack algorithm [44] for selecting the k-points are not required for quantum dots [45]. The total energy and bandgap of the conduction and valance bands are essential computational issues in our work. Given the importance of the electronic structure results and how they affect calculations of the optical properties and determination of the dielectric coefficient tensor, other estimates are performed using a high-level computing package with thermal DFT (th-DFT) ability (NanoDCAL) [46]. In this way, the obtained results, including the total energy and bandgap, are confirmed by the th-DFT output data. The th-DFT takes into account the electronic temperature of the system via an occupation function using Fermi-Dirac statistics [47]. We investigated the temperature dependence of the bandgap and total energy of the GQD structures using the thermal DFT (th-DFT) computational method. In th-DFT, the Mermin-Kohn-Sham equations are solved integrally [47], and the electron temperature is used with the Fermi-Dirac distribution to obtain the occupation function of the electronic states. In standard DFT calculations, the charge density ) (r ρ is derived from Equation (1): where ) (r i φ is the ith Kohn-Sham orbital. However, for th-DFT, the temperature-dependent charge density ) (r τ ρ is obtained by considering the Fermi-Dirac distribution function as follows: where fi is the Fermi occupation factor, which is defined by where µ is the chemical potential, i τ ε is the total energy of an electron located in the i th orbital with the τ thermal energy component (that is, τ = kBT), and kB and T are the Boltzmann constant for the electron temperature in Kelvin, respectively. Additionally, by pseudopotential functional rewriting for the exchange correlation energy computation in DFT calculations, a more accurate temperature dependence can be considered.
In the optical simulations, we assumed that the incident light is polarized in the transverse electric (TE) mode. To evaluate the optical properties of the GQDs, we obtained their optical absorption spectra for both parallel and perpendicular incidence. Under vertical propagation, the The convergence threshold for the self-consistent-field (SCF) calculation loop is considered as 10 −8 Ry. Since the energy bands that are obtained for the GQDs are of a "flat" type, we used the gamma-point algorithm to choose the k-point. Algorithms such as the Monkhorst-Pack algorithm [44] for selecting the k-points are not required for quantum dots [45]. The total energy and bandgap of the conduction and valance bands are essential computational issues in our work. Given the importance of the electronic structure results and how they affect calculations of the optical properties and determination of the dielectric coefficient tensor, other estimates are performed using a high-level computing package with thermal DFT (th-DFT) ability (NanoDCAL) [46]. In this way, the obtained results, including the total energy and bandgap, are confirmed by the th-DFT output data. The th-DFT takes into account the electronic temperature of the system via an occupation function using Fermi-Dirac statistics [47]. We investigated the temperature dependence of the bandgap and total energy of the GQD structures using the thermal DFT (th-DFT) computational method. In th-DFT, the Mermin-Kohn-Sham equations are solved integrally [47], and the electron temperature is used with the Fermi-Dirac distribution to obtain the occupation function of the electronic states. In standard DFT calculations, the charge density ρ(r) is derived from Equation (1): where φ i (r) is the ith Kohn-Sham orbital. However, for th-DFT, the temperature-dependent charge density ρ τ (r) is obtained by considering the Fermi-Dirac distribution function as follows: where f i is the Fermi occupation factor, which is defined by where µ is the chemical potential, ε τ i is the total energy of an electron located in the ith orbital with the τ thermal energy component (that is, τ = k B T), and k B and T are the Boltzmann constant for the electron temperature in Kelvin, respectively. Additionally, by pseudopotential functional rewriting for the exchange correlation energy computation in DFT calculations, a more accurate temperature dependence can be considered.
In the optical simulations, we assumed that the incident light is polarized in the transverse electric (TE) mode. To evaluate the optical properties of the GQDs, we obtained their optical absorption spectra for both parallel and perpendicular incidence. Under vertical propagation, the chosen polarization places the electric field vector in the GQD plane. In contrast, with horizontal light, the polarization is such that the electric field vector lies perpendicular to the GQD plane. Since the electric field of the TE-polarized beam is perpendicular to the propagation direction, the in-plane light polarization is the same as the perpendicular incidence, and the out-of-plane light polarization is the same as the parallel incidence. Nevertheless, to calculate the absorption of light with perpendicular incidence, the diagonal elements x and y from the dielectric coefficient tensor are multiplied by the in-plane elements of the electric field (E) to generate the electric displacement field (D). Furthermore, with parallel incidence, the diagonal element z from the dielectric coefficient tensor is multiplied by the out-of-plane element of the electric field (E) to generate the electric displacement field (D). Therefore, we evaluated the optical behavior by studying the frequency spectra of the dielectric coefficient tensor (real and imaginary parts). The imaginary part was associated with the loss and absorption of light. We plotted the frequency spectra of n and α for both parallel and perpendicular incidences for both monolayer and bilayer GQDs. In this work, we calculated the dielectric coefficient in the range of 0-15 eV. The number of the frequency points of the spectrum was around 1000, meaning that the frequency spectrum for the dielectric coefficient had a resolution of 15 meV. Independent particle approximation (IPA) was used to determine the diagonal components of the dielectric tensor. The considered IPA model used the random phase estimate to calculate the complex dielectric coefficient function ε(ω) = ε 1 (ω) + jε 2 (ω), which depends on the frequency ω [48]. The numerical results obtained from DFT were used to calculate the integral function of the first-order perturbation theory in the format of Hamiltonian matrix elements based on the single-particle Bloch wave function. In the first step, the imaginary part of the dielectric coefficient was obtained. In this way, the obtained results had a physical basis. Then, perturbation theory was used according to the adiabatic parameter from the broadening function. Using the broadening feature resolved the problems caused by the infinite lifetime of the excited states. This method accounted for interband transitions but ignored intraband electronic transitions. A Gaussian distribution was then used for the broadening function. The broadening parameter (or intersmear) was considered as 0.1 eV, which caused the dielectric coefficient spectrum drawing to be displayed contiguously. Finally, we obtained a continuous frequency spectrum for the imaginary part of the dielectric coefficient. This imaginary part was used to calculate the real part of the dielectric factor using the Kramers-Kronig transformation [49]. Combining the imaginary (ε 2 ) and real (ε 1 ) parts give the final equation for the dielectric coefficient. These two quantities for different frequencies, ω, allowed us to calculate other optical properties [50][51][52], such as the refraction index, and the absorption coefficient, For our analytical computations, we used Python and the pandas library [53]. The open-source Gnuplot software tool was used for plotting the curves in this work [54].

Electronic Properties
For monolayer systems, we began by considering a monolayer of graphene with no lateral quantum confinement. Then, we investigated how the energy-band structure changes because of quantum confinements in the x and y directions to transform the graphene into quantum dots. Figure 3a shows the energy levels of the C 132 H 40 GQD, the curvatures of which disappear entirely because of electron localization, thereby rendering the energy bands flat. The flat band structure from Γ point to Z observed in Figure 3a was due to the electron confinement in GQDs arising from the Heisenberg uncertainty principle. A bandgap was also created, and graphene was transformed from a semi-metal into a semiconductor in the GQD form. Table 1 gives the electronic characteristics of monolayer and bilayer GQDs compared to those of single and bilayer graphene. As shown, increasing the size of a GQD decreased its bandgap energy. The GQD size, its geometry (triangular or hexagonal), and its side edge (zigzag or armchair) enabled bandgap engineering. For example, as shown in Table 1 for the HexC 96 zzH 24 and TriC 168 acH 42 systems, a triangular GQD had a larger bandgap compared to a hexagonal GQD. Unlike TriC 168 acH 42 , which was larger than HexC 96 zzH 24 , the bandgap of the TriC 168 acH 42 system was broader than that of the HexC 96 zzH 24 system. Thus, triangular GQDs had a higher bandgap than hexagonal GQDs of the same size.
The total energy of the system includes the interactions between all particles in the system. When considering electrostatic potential, we multiplied the charges. Thus, for a minus charge (electrons) and a positive charge (nucleus), the result was always negative, indicating attraction (but for two negative or two positive point charges, the result was still positive, indicating repulsion). This is why all the total energies in Table 1 were negative. We also found the total energy for the graphene sheet and the GQD system. Table 1 shows that all GQD systems had a total energy value that exceeded that of the graphene sheet. The situation for the electric carrier in the graphene sheet prefers a zero-energy reference point since the zero-energy reference point is a (fictitious) system with all the particles (electrons and nuclei) at rest infinitely far away from each other. Conversely, the electrons in GQDs are confined and under the electrostatic potential of the cores.
To evaluate the accuracy and validity of the results obtained from standard DFT, Table 1 compares these results with those obtained via th-DFT for monolayer GQD systems. The obtained results indicate that the ground-state DFT calculations were valid for determining the electronic properties of GQDs up to 300 K. The real reason for this result is that the thermal energy component of K B T was small compared to that of other contributions to the total energy of GQDs. We also calculated the electronic properties of bilayer GQDs based on AB (Bernal) stacking. Figure 3b shows how the bilayer structure affected energy levels. The comparison in Figure 3a,b of the energy-band structures in the gamma valley (Γ) for monolayer and bilayer GQDs shows that the creation of new energy levels close to the previous levels was due to π-type bonds between carbon atoms in the bottom and top layers of the bilayer GQD. As Figure 3 shows, the energy bandgap of a bilayer GQD was significantly smaller than that of a monolayer GQD. The last four rows of Table 1 summarize the numerical results obtained from the DFT simulation for bilayer graphene and GQD systems. To examine the electronic property results, we conducted a comparative study using two codes for each bilayer system. Table 1 also lists the comparative results for the electronic properties of bilayer GQD systems.
In Table 1, our calculations were limited to a structure of up to 210 atoms. DFT requires enormous computational resources for GQDs with a large number of atoms, which were scaled in the order of N 3 (O(N 3 )). Here, N is the number of free electrons and relates to the size of the problem. Thus, we could tune the semi-empirical parameters of the tight-binding (TB) model by matching the DFT results for small GQDs, and then generalize this adjusted TB model for large GQDs up to many thousands of atoms to obtain results with acceptable accuracy.
As shown in Table 2, we compared our results with those of other articles and found good agreement between them [29,30,55].
For the triangular graphene quantum dots (TGQD), there is another important issue-the creation of degenerate energy levels near the Fermi level caused by breaking the symmetry of the number of sublattices A and B forming the TGQD; subsequently, the degenerative edge states are created [56]. Such differences in the electronic properties of TGQD produce excellent magnetic properties within them. To investigate these magnetic properties, spin polarized DFT calculations are required, which would require another study.

Optical Properties
The dielectric function is a complex quantity that describes the linear response of the structure to electromagnetic wave radiation. Figure 4 shows the real and imaginary parts of the dielectric coefficients for different arrangements of monolayer GQDs, as well as the in-plane and out-of-plane polarization of the incident light.

Optical Properties
The dielectric function is a complex quantity that describes the linear response of the structure to electromagnetic wave radiation. Figure 4 shows the real and imaginary parts of the dielectric coefficients for different arrangements of monolayer GQDs, as well as the in-plane and out-of-plane polarization of the incident light.  As shown in Figure 4b, for out-of-plane polarization, the imaginary part of the dielectric coefficient was negligible for photon energies of less than around 5 eV; hence, no absorption could occur. However, for in-plane polarization, the dielectric coefficient values in the visible and near-IR regions were significant, making it suitable for light absorption applications. Figure 4 shows that changing the GQD's geometry and size changed its absorption coefficient and wavelength. Clearly, with in-plane polarization, reducing the GQD size shifted the optical absorption peak to higher photon energies. Figure 4b shows hexagonal GQDs with a zigzag edge of C42. Here, the first peak of the absorption spectrum occurred at a point close to 2.4 eV; meanwhile, for a more massive GQD structure such as C168, the first peak of the absorption spectrum under in-plane polarization occurred at an energy close to 1.4 eV. However, for out-of-plane polarization, this rule was reversed, i.e., the first absorption peak shifted to a higher value with an increase in GQD.
To summarize the optical properties of the monolayer GQD structures, the frequency positions of the first and second peaks of the absorption spectra for the in-plane and out-of-plane polarizations are tabulated in Table 3. These peaks are due to interband electronic transitions (transitions between the occupied and unoccupied states). Notably, as shown by the results in Table 3, by changing the incident light polarization from in-plane polarization to out-of-plane polarization electric fields, the first absorption and photodetection peak shifted to high photon energy in monolayer systems.
Stacking, the use of bilayered graphene, and establishing an interlayer bond all considerably change the system energy and electronic band structure of the graphene. These changes in electronic As shown in Figure 4b, for out-of-plane polarization, the imaginary part of the dielectric coefficient was negligible for photon energies of less than around 5 eV; hence, no absorption could occur. However, for in-plane polarization, the dielectric coefficient values in the visible and near-IR regions were significant, making it suitable for light absorption applications. Figure 4 shows that changing the GQD's geometry and size changed its absorption coefficient and wavelength. Clearly, with in-plane polarization, reducing the GQD size shifted the optical absorption peak to higher photon energies. Figure 4b shows hexagonal GQDs with a zigzag edge of C 42 . Here, the first peak of the absorption spectrum occurred at a point close to 2.4 eV; meanwhile, for a more massive GQD structure such as C 168 , the first peak of the absorption spectrum under in-plane polarization occurred at an energy close to 1.4 eV. However, for out-of-plane polarization, this rule was reversed, i.e., the first absorption peak shifted to a higher value with an increase in GQD.
To summarize the optical properties of the monolayer GQD structures, the frequency positions of the first and second peaks of the absorption spectra for the in-plane and out-of-plane polarizations are tabulated in Table 3. These peaks are due to interband electronic transitions (transitions between the occupied and unoccupied states). Notably, as shown by the results in Table 3, by changing the incident light polarization from in-plane polarization to out-of-plane polarization electric fields, the first absorption and photodetection peak shifted to high photon energy in monolayer systems.
Stacking, the use of bilayered graphene, and establishing an interlayer bond all considerably change the system energy and electronic band structure of the graphene. These changes in electronic structure can change the shape of the optical absorption spectrum, especially for out-of-plane polarization. Figure 5 compares the dielectric coefficients of the monolayer and bilayer graphene sheets and the GQDs for different incident light polarizations. The appearance of a new absorption peak at low energy is the most crucial change in the absorption spectrum of bilayer systems compared to monolayer systems. The existence of absorption peaks in the IR region is suitable for IR photodetection applications. As shown in Figure 5, this phenomenon can be observed in both the bilayer graphene and bilayer GQDs. The absorption coefficient α(ω) indicates the attenuation percentage of light intensity per unit distance when a light wave propagates in the given material. The refractive index n(ω) characterizes the velocity of light in different mediums caused by the illuminating beam and electron interactions. Figure 6 shows the refractive index and absorption spectra of the monolayer and bilayer GQDs. With out-of-plane polarization, the electric field component of the incident light was perpendicular to the graphene surface. Thus, this field affected the common interface bonds of the two graphene layers.
Photonics 2020, 7, x FOR PEER REVIEW 9 of 16 structure can change the shape of the optical absorption spectrum, especially for out-of-plane polarization. Figure 5 compares the dielectric coefficients of the monolayer and bilayer graphene sheets and the GQDs for different incident light polarizations. The appearance of a new absorption peak at low energy is the most crucial change in the absorption spectrum of bilayer systems compared to monolayer systems. The existence of absorption peaks in the IR region is suitable for IR photodetection applications. As shown in Figure 5, this phenomenon can be observed in both the bilayer graphene and bilayer GQDs. The absorption coefficient α(ω) indicates the attenuation percentage of light intensity per unit distance when a light wave propagates in the given material. The refractive index n(ω) characterizes the velocity of light in different mediums caused by the illuminating beam and electron interactions. Figure 6 shows the refractive index and absorption spectra of the monolayer and bilayer GQDs. With out-of-plane polarization, the electric field component of the incident light was perpendicular to the graphene surface. Thus, this field affected the common interface bonds of the two graphene layers.  Note that the overall shape of the absorption spectrum in the case of in-plane polarization remained relatively unaffected by bilayering the GQDs, as shown in Figure 6. To ensure the generality of this result (i.e., the changes in the absorption spectra due to bilayering under out-of-plane polarization), we calculated the dielectric coefficient tensors for different sizes of hexagonal bilayer GQDs and confirmed our obtained results. The last four rows in Table 3 summarize the first and second peaks of the absorption spectra for the bilayer graphene and GQD systems under illumination via out-of-plane and in-plane polarization. For incident light with out-of-plane polarization, in bilayer systems, the first absorption peak was observed with energy close to 0.8 eV, while in the monolayer systems, this energy was close to 6 eV. Therefore, unlike monolayer GQD systems, the bilayer ones with parallel incidence (i.e., out-of-plane polarization) could be used for photodetection applications. Therefore, the out-of-plane polarization (horizontal incident light) in bilayer systems could absorb IR light (the communication window). Note that the overall shape of the absorption spectrum in the case of in-plane polarization remained relatively unaffected by bilayering the GQDs, as shown in Figure 6. To ensure the generality of this result (i.e., the changes in the absorption spectra due to bilayering under out-of-plane polarization), we calculated the dielectric coefficient tensors for different sizes of hexagonal bilayer GQDs and confirmed our obtained results. The last four rows in Table 3 summarize the first and second peaks of the absorption spectra for the bilayer graphene and GQD systems under illumination via out-of-plane and in-plane polarization. For incident light with out-of-plane polarization, in bilayer systems, the first absorption peak was observed with energy close to 0.8 eV, while in the monolayer systems, this energy was close to 6 eV. Therefore, unlike monolayer GQD systems, the bilayer ones with parallel incidence (i.e., out-of-plane polarization) could be used for photodetection applications. Therefore, the out-of-plane polarization (horizontal incident light) in bilayer systems could absorb IR light (the communication window). All of the first and second peaks that are introduced in Table 3 are related to an electronic transition E ij from the ith valence band to the jth conduction band. Figure 7 shows this phenomenon for the ABhexC 132 H 40 bilayer GQD via a projected density of states diagram (PDOS). For example, the E 11 and E 22 denoted in Figure 7 had excellent agreement with the contents of Table 3 (8th row).  All of the first and second peaks that are introduced in Table 3 are related to an electronic transition Eij from the i th valence band to the j th conduction band. Figure 7 shows this phenomenon for the ABhexC132H40 bilayer GQD via a projected density of states diagram (PDOS). For example, the E11 and E22 denoted in Figure 7 had excellent agreement with the contents of Table 3 (8th row).

DFT vs. TDDFT for Calculating Optical Properties
In addition to models that use an independent electron model based on DFT results to extract the optical absorption spectrum, there is a more accurate method called time-dependent DFT (TDDFT). Standard DFT has two major limitations: (i) in standard DFT, the ground state is calculated without the excitation, and (ii) the model used for standard DFT does not consider many-body effects. Hence, it assumes that electrons are independent of each other and replaces the exact exchange correlation potential of the interactions of electrons in the materials with their approximated pseudopotential to simplify the model of the electron's many-body interactions. Therefore, here we used the turbo TDDFT code based on the multiparticle interaction model for calculating the absorption coefficient via the linear response regime [57,58]. The interactions of

DFT vs. TDDFT for Calculating Optical Properties
In addition to models that use an independent electron model based on DFT results to extract the optical absorption spectrum, there is a more accurate method called time-dependent DFT (TDDFT). Standard DFT has two major limitations: (i) in standard DFT, the ground state is calculated without the excitation, and (ii) the model used for standard DFT does not consider many-body effects. Hence, it assumes that electrons are independent of each other and replaces the exact exchange correlation potential of the interactions of electrons in the materials with their approximated pseudopotential to simplify the model of the electron's many-body interactions. Therefore, here we used the turbo TDDFT code based on the multiparticle interaction model for calculating the absorption coefficient via the linear response regime [57,58]. The interactions of electrons (Hartree and exchange correlation effects) were taken into account by the ab initio fully self-consistent scheme. In this code, the Liouville-Lanczos equation is solved to derive the dipole polarizability tensor χ ij (ω) in the standard batch representation [59], thereby avoiding the need to multiply or invert large matrices. Along with the polarizability, we obtained the oscillator strength S(ω) as where µ B is a physical constant (Bohr magneton) and e is the electron charge. Comparing Equation (6) and Equation (7), which describes the absorption coefficient [60], it can be seen that S(ω) is equivalent to the optical absorption coefficient: Figure 8 shows the absorption spectra calculated using both models (i.e., α(ω) in the IPA model with standard DFT and S(ω) with TDDFT) for hexagonal GQDs with a zigzag edge that have been passivated by hydrogen atoms (C 24 H 12 ). The effect of the TDDFT modifications on the results of the standard DFT, as shown in Figure 8, was a large blue shift in the peak of the light absorption spectrum. To ensure that the blue shift in the absorption spectrum occurred due to the many-body interactions, calculations were also performed for the C 42 H 18 system with an armchair edge. A blue shift at the peak of absorption was also observed in the other graphene and non-graphene nanostructures. This demonstrates that the results from the independent electron model were not exact, whereas the experimental studies confirmed the results of the electron interaction model [29]. However, the most notable achievement of the present work is its investigation into the physical effects of out-of-plane polarization and the bilayering of GQDs on the light absorption behavior in a comparative framework. Since our intention was not to obtain absolute and precise numerical values for the absorption peaks, the use of independent electron approximation, which requires far fewer computational resources, was sufficient for our purposes.
Photonics 2020, 7, x FOR PEER REVIEW 12 of 16 electrons (Hartree and exchange correlation effects) were taken into account by the ab initio fully self-consistent scheme. In this code, the Liouville-Lanczos equation is solved to derive the dipole polarizability tensor χij(ω) in the standard batch representation [59], thereby avoiding the need to multiply or invert large matrices. Along with the polarizability, we obtained the oscillator strength S(ω) as where µB is a physical constant (Bohr magneton) and e is the electron charge. Comparing Equation (6) and Equation (7), which describes the absorption coefficient [60], it can be seen that S(ω) is equivalent to the optical absorption coefficient:  .
(7) Figure 8 shows the absorption spectra calculated using both models (i.e., α(ω) in the IPA model with standard DFT and S(ω) with TDDFT) for hexagonal GQDs with a zigzag edge that have been passivated by hydrogen atoms (C24H12). The effect of the TDDFT modifications on the results of the standard DFT, as shown in Figure 8, was a large blue shift in the peak of the light absorption spectrum. To ensure that the blue shift in the absorption spectrum occurred due to the many-body interactions, calculations were also performed for the C42H18 system with an armchair edge. A blue shift at the peak of absorption was also observed in the other graphene and non-graphene nanostructures. This demonstrates that the results from the independent electron model were not exact, whereas the experimental studies confirmed the results of the electron interaction model [29]. However, the most notable achievement of the present work is its investigation into the physical effects of out-of-plane polarization and the bilayering of GQDs on the light absorption behavior in a comparative framework. Since our intention was not to obtain absolute and precise numerical values for the absorption peaks, the use of independent electron approximation, which requires far fewer computational resources, was sufficient for our purposes.

Conclusions
In this paper, we extracted the electronic structures of monolayer and bilayer GQDs based on DFT and thermal DFT (th-DFT) calculations. The calculated results show that the corresponding bandgaps of monolayer and bilayer GQDs were in the range of 1.2-2.8 eV and 0.9-3.0 eV, respectively. All these quantum dots were direct bandgap semiconductors and had flat band structure characteristics. We compared our results with those other articles and found good agreement between them. In addition, we calculated the permittivity tensors for each structure. The elements of this tensor show that the graphene and graphene quantum dots (GQDs) were anisotropic. Moreover, by comparing the results of the DFT method (ground state) with those of the TDDFT method (excited states in a linear response regime), a blue shift in the absorption spectrum was observed due to migration from the independent electron model to many-body interactions. The results show that monolayer GQDs covered the range from infrared (IR), to visible, to ultraviolet (UV) light. By increasing the number of carbon atoms involved in the GQDs, the optical absorption spectrum changed from visible to IR for the in-plane polarization of the incident light. In contrast, for the out-of-plane polarization, when the GQD size increased, the absorption spectrum moved from UV toward a deep UV range of 85-250 nm. Furthermore, in the bilayer graphene systems, a new absorption peak was produced at a lower incident photon energy with the out-of-plane polarization. Remarkably, the absorption peak was obtained in an IR range of 500-1600 nm under illumination from out-of-plane polarization. Therefore, bilayer GQDs could be suitable for integrated photodetection applications with laterally propagated light.