Interband transitions in closed-shell vacancy containing graphene quantum dots complexed with heavy metals †

High-performance optical detection of toxic heavy metals by using graphene quantum dots (GQDs) requires a strong interaction between the metals and GQDs, which can be reached through a functionalization/ immobilization procedure or doping eﬀect. However, commonly used surface activation approaches induce toxicity into the analysis system and, therefore, are ineligible from the environmental point of view. Here, we show that artificial creation of vacancy-type defects in GQDs can be a helpful means of intentional control of the active sites available for reaction with cadmium (Cd), mercury (Hg) and lead (Pb). Using restricted density functional theory (DFT) and time-dependent DFT (TD-DFT) methods, we predict the eﬀect of vacancy complexes not previously studied to describe the binding ability of GQDs towards metal adsorbates. We also show that the interband absorption in closed-shell GQDs complexed with Cd, Hg and Pb is strongly dependent on the vacancy type and can be eﬃciently tuned to attain the desired coloration of the analysis system. The results suggest that the vacancy defects play an important role in governing the hybridization between locally-excited (LE) and charge-transfer (CT) states of the GQDs. Based on the molecular orbital analysis and in-depth knowledge of excited states, the mechanisms underlying the interband absorption are discussed.


Introduction
Toxic heavy metals (HMs) as a factor of environmental pollution are still one of the main problems of the present day. 1,2 Their high toxicity at relatively low concentrations 3 and their exceptionally large bioaccumulation ability 4 have a harmful effect not only on individual living creatures, 5,6 but also on the planetary ecosystem. 7,8 In response to this global challenge, the World Health Organization (WHO) has implemented environmental regulations to protect humankind from the most dangerous heavy metal pollutants. 9 In accordance with the WHO requirements, first of all, it is necessary to control the content of the three most toxic heavy metals, namely cadmium (Cd), 10 mercury (Hg), 11 and lead (Pb), 12 in food, soil, drinking water, outdoor and indoor air, and biological fluids. At present, the permissible levels of these hazardous elements in each of the categories have been determined. [10][11][12] To meet the WHO guidelines for minimal acceptable heavy metal concentrations in water, air and food, it is vitally important to realize fast methods for real-time detection of HMs and their removal/filtration. During the last decades, extensive research efforts have mostly been focused on development/improvement of effective ''bulky'' detection methods (chromatography, mass spectrometry, etc.), [13][14][15] design of sensitive materials and devices (electrodes, hybrid structures, field emission transistors, and Schottky diodes), [16][17][18][19][20][21][22][23][24][25] and fabrication of potent adsorbents (filters, and membranes). [26][27][28][29][30][31][32][33][34][35] Considering the drawbacks of the traditional ''bulky'' time-consuming techniques, numerous research groups, nowadays, investigate graphene-family materials (graphene quantum dots, reduced graphene oxide, epitaxial graphene, etc.) as possible candidates for conceptualization of both sensors and adsorbents. [36][37][38][39][40][41] Being a two-dimensional (2D) monoatomic sheet of sp 2conjugated carbon atoms, 42 graphene possesses unprecedentedly high sensitivity towards detection of individual atoms and molecules adsorbed onto its planar surface. 43 The basic physics underlying such detection is related to the possibility of local modulation of the charge-transfer reactions. It is reasonable to assume that providing as high a rate as possible of the electrontransfer reaction between graphene and adsorbates is an imperative requirement to consider graphene as an appropriate sensitive and adsorbing material. In other words, the most important parameters of a sensing device and an adsorbent (namely, the response time and the adsorption capacity) strongly depend on the interaction strength between graphene and the adsorbates we are interested in. The optimization of these parameters requires a deep understanding of the interaction mechanism. With respect to heavy metals, there are two fundamental processes governing this mechanism: (i) chemisorption of the HM divalent ions and (ii) physisorption of the HM neutral atoms. An ideal defect-free basal plane of graphene offers several high-symmetry adsorption sites (so-called bridge, hollow and on-top sites) for reaction with HMs. While the divalent Cd, Hg and Pb ions are found to be strongly bonded preferentially to these sites (providing large charge transfer conditions), the elemental HMs are differently physisorbed and bonded to graphene through weak dispersion (van der Waals) forces. 44 Furthermore, in the case of charged ions and neutral atoms of HMs immersed in liquids the interaction strength is very sensitive to the solvent effect and is significantly weakened due to the so-called solvent-involved cation-p interaction. 45 A major disadvantage of the basal plane of graphene towards HM adsorption is the poor sensitivity in those cases when the interaction mechanism is driven predominantly by the physisorption of neutral metal adatoms. For example, the most frequently used electrochemical techniques for HM detection (anodic stripping voltammetry, cyclic voltammetry and chronoamperometry) imply that the output signal of the graphene working electrode is related to HM-involved spontaneous redox reactions. 46 As the result of such reactions, the subsequent signal is proportional to the number of metal adatoms physisorped onto the graphene surface. Additionally, optical detection methods, namely UV-vis absorption spectroscopy and colorimetric detection, are also inappropriate approaches for recognition of elemental heavy metals because of insufficient quenching of the optical signal and negligible color changes, which remain invisible to the unaided eye. 46 Therefore, the effective adsorption of HMs requires that the graphene-family materials are appreciably functionalized/modified to provide the necessary reactive surface area.
Various materials were suggested for the activation of graphene surfaces to improve the sensitivity towards Cd, Hg and Pb. [47][48][49][50][51][52][53][54][55][56][57] Among them, special emphasis has been given to hybrid composites based on graphene-family nanomaterials [47][48][49][50][51][52] and surface-functionalized graphene with chemical moieties and immobilized recognition elements. [53][54][55][56][57] Unfortunately, in most cases the entire preparation procedure of the working electrode, sensitive solution electrolyte and active recognition element in sensing devices is still very complicated as it involves surface functionalization and multistage chemical reactions. Furthermore, as has been demonstrated for other carbon materials, especially for multi-wall carbon nanotubes, the surface functionalization with carbonyl, carboxyl, and hydroxyl groups can increase their toxicity. 58 In the light of global trends towards development of ''green'' fabrication technologies for new less harmful materials and chemicals, 59 utilizing highly-toxic functionalized graphene-based materials for detection and removal of heavy metals becomes less safe and less desirable. There is, therefore, a strong demand for development of an easy-to-use versatile platform (including only the key benefits of graphene) that allows both high-performance detection and effective removal of heavy metals. In this regard, introducing defects (mono-vacancies and vacancy complexes) into graphene-family materials can be an effective tool to enhance the adsorption ability of heavy metals, thereby avoiding surface functionalization/immobilization steps. Recently, it has been demonstrated that the defect-engineering approach in the case of graphene can promote the increase of the binding energy of alkaline-earth metals and hydrogen, 60 gas molecules and water, 61,62 Li,63 Pt, 64 Cd, 65 Co 66 and 3d transition metals (from Sc to Zn). 67 Nevertheless, there is still a limited amount of data in the literature on the role of vacancy-related defects in the complexation between elemental toxic heavy metals and graphene-like structures.
In this article, we report a first-principles study of vacancyassisted interactions between elemental heavy metals (Cd, Hg and Pb) and graphene quantum dots. Since the formation of vacancy-type defects may be accompanied by the appearance of unpaired electrons, we realize that heavy metal adsorption may occur on GQDs in the singlet and triplet states, i.e. on closedshell and opened-shell structures. Nevertheless, the main aim of this work is to understand how the metal adsorbates will affect structural, electronic and optical properties of closedshell GQDs having vacancy complexes. In this case, the increase in the local density of states at dangling bonds is expected to be responsible for the tuning of reactivity of GQDs, while the effect of the unpaired electrons (in the triplet state) on the binding ability is not considered here. Special focus will be given to understanding of the interband transitions that drive the optical response of GQDs after HM exposure and to the influence of the vacancy type on the excited states of the interacting constituents. We anticipate that this understanding will importantly advance the research on GQD-based sensing devices. The choice of GQDs as a representative of the graphene-family materials is justified by the fact that GQD-related systems are multifunctional platforms, which can be used for production of micron-sized solid sheets (a promising working electrode material for electrochemical detection of heavy metals), 68 catalysts accelerating the reaction with HMs, 69 modified working electrodes, 70,71 aqueous solution electrolytes for optical and colorimetric detection of heavy metals [72][73][74][75][76][77][78][79][80][81][82] and adsorbent agents for removal of heavy metals. [83][84][85]

Theoretical approach
All quantum chemical calculations were carried out with the Gaussian 09 Rev. D.01 program package, using the default convergence criteria. 86 The zigzag-edged C 54 H 18 (circumcoronene) polycyclic aromatic hydrocarbon (PAH) comprising 19 hexagonal benzene-like rings was chosen as a model structure of GQDs. This type of GQDs can absorb light very efficiently in the visible part of the spectrum and, therefore, can be used for colorimetric detection of heavy metals. 45 The creation of a monovacancy and complex vacancies can be simply modelled by removing one carbon atom and the nearest neighbour 2, 3, and 4 carbon atoms from the GQDs. Thus, five sets of geometrical configurations of GQDs were considered -the pristine defect-free GQDs and others with monovacancy, divacancy, trivacancy and tetravacancy defects in closed-shell singlet states (Fig. 1). All structures were terminated with the same number of hydrogen atoms. The equilibrium HM adsorption configurations were found by relaxing the HM@GQD (HM = Cd, Hg and Pb) complexes under gasphase conditions using the PBE1PBE-D3 level of restricted DFT (which includes the empirical dispersion correction) 87 with a 6-311G(d,p) basis set for carbon and hydrogen atoms as well as a basis set developed by the Stuttgart-Dresden-Bonn group (SDD) for the heavy metal atoms. 88 Starting coordinates before geometry optimization were the same for all three metals in GQD@HMs complexes. In other words, we placed the HMs directly on the vacancy-type site of GQDs and then the structures were allowed to fully relax. Therefore, it is reasonable to assume that unique binding patterns for the GQDs and HM atoms only originate from the unique metal properties and the specific interaction between defects and metals. It should be mentioned that the PBE1PBE-D3 method takes explicitly London dispersion forces (or van der Waals forces) into account by using an additional term (which is proportional to R À6 ) and, therefore, this approach is suitable to describe the weak physisorption of elemental HMs onto the sp 2 plane. Since the main aim of the current work is to investigate the possible ways of enhancement of the interaction strength in the ''HM@GQD'' complex via artificial creation of vacancy defects (playing a role of reactive sites), a special focus is placed on the achievement of full control of the binding energy and charge transfer from adsorbates to GQDs. The binding energy of the heavy metal species to the defect-free and vacancy-containing GQDs was calculated by using the following equation: (1) where E GQD+HM is the total energy of the GQD after complexation with the HM, E GQD is the total energy of the isolated individual GQD, and E HM is the total energy of an individual HM. The binding energies have been corrected for the basis set superposition error (BSSE) by means of the counterpoise method. 89 The charge transfer is calculated using Mulliken population analysis (MPA). 90 It should be mentioned that MPA is sensitive to the basis set and can, therefore, cause overestimation of the charge transfer magnitude in the interacting complexes, especially in those cases when diffuse and polarization functions are added to basis sets. 91 These disadvantages mainly originate from the fact that MPA deals with nonorthogonal basis sets. On the other hand, Natural Population Analysis (NPA) can help to overcome this bottleneck, since this model implies the use of orthonormal natural atomic functions. 92,93 In contrast to MPA and NPA, the Hirshfeld model 94 can provide a clear partitioning of the electron density between all atoms belonging to the interacting complex and is regarded as a less basis-set dependent approach to predict atomic charges. To investigate the applicability of the MPA model for calculation of partial atomic charges on elemental heavy metals placed on GQDs, we compare the Mulliken population analysis to the Hirshfeld analysis, popular natural population analysis, and available reference data.
In the next step, in order to obtain the absorption spectra, density transition matrices and oscillator strengths for the most probable excited states, time-dependent density functional theory (TD-DFT) calculations were carried out on the full set of the GQD structures before and after complexation with elemental heavy metals using the same level of DFT theory (PBE1PBE-D3/6-311G(d,p)/SDD). The first 12 excited states were investigated and analyzed. In most cases any changes in the ultraviolet-visible (UV-vis) absorption spectra of GQDs after complexation with heavy metals can lead to changes in coloration of the interacting complexes. To estimate the possible colour changes perceived by the naked eye, we performed ab initio colorimetric calculations under daylight conditions based on the absorption spectra and colour-matching functions. The theoretical approach for simulation of the perceived colour can be found elsewhere. 95 It is obvious that understanding the fundamental physics underlying the colour changes of GQDs after interaction with Cd, Hg and Pb is an important step towards designing a printed and easy-to-use sensing platform for the real-time detection of elemental heavy metals.

Equilibrium structures, reactivity and binding energies
The mean size of vacancy-type defects should unambiguously define the local adsorption of heavy metals at the atomic scale. 96 In this regard, one of the most important parameters which allows us to estimate the role of the vacancies in affecting the interaction energy between elemental heavy metals and graphene quantum dots is the ability of the vacancy defects to trap electric charges. It is reasonable to assume that the changes in the spatial distribution of the wave-function near the carbon vacancies imply an enhancement of the electron density between the closest to the vacancy carbon species. This can be interpreted as the formation of shorter C-C bonds within inner carbon rings and even five-membered pentagonal carbon rings (Fig. 1). The first type of vacancy defect considered in our calculations is the monovacancy (Fig. 1b). In this configuration, after a vacancy is created, three sp 2 chemical bonds are broken and the atoms surrounding the vacancy site displace from their initial lattice positions. Consequently, the electrons that participate in broken bonds tend to localize around the single vacancy. Due to this reason localized vacancy energy levels can appear among the allowed molecular orbitals in GQDs. The second defect structure (the so-called 5-8-5 divacancy defect) is related to the removal of two carbon atoms, and, consequently, to the formation of an octagon at the vacancy site with two carbon pentagons (Fig. 1c). It is interesting to note that Liu et al., 97 in a recent study, demonstrated that the most stable electronic state of the evennumbered vacancy-type graphene nanofragments (with the 5-8-5 divacancy defect) is a closed shell singlet state without unpaired electrons, while the presence of unpaired electrons was observed for the less stable structure with the 555777 divacancy in the triplet state. In contrast to this, Yumura et al. 96 showed that the triplet state of C 96 H 26 with vacancy defects is more stable than singlet states and, thus, the available unpaired electrons improve the reactivity of the graphene structure towards platinum clusters. Although the single vacancy and double vacancy geometries are expected to be the most stable in graphene quantum dots, the formation of high-order vacancycomplexes, namely trivacancy 98 and tetravacancy 99 is also possible ( Fig. 1d and e). A complete description of the binding ability of the GQDs containing vacancy complexes towards heavy metals requires accurate knowledge not only of the vacancy-induced redistribution of the carrier wave-function in GQDs, but also the DOS of a certain structure. For instance, the energetic positions of the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) levels with respect to the metal work function can affect the direction of the charge transfer between the metal adsorbate and the GQDs. 100 Before discussing the DOS for each structure presented in Fig. 1, the molecular electrostatic potential (MEP) will be examined to analyse the reactivity of the defective GQDs towards charge-controlled interaction with metals adsorbates. It should be mentioned that the deep analysis of the MEP surface enables determination of the most reactive sites across the sp 2 -plane of the GQDs, namely electron-sufficient and/or electron-deficient sites. 101 As demonstrated in Fig. 2a, the basal plane of defect-free graphene does not provide efficient reactive sites for the chemical adsorption of elemental heavy metals. Indeed, according to electrochemical investigations of graphene, 102 the basal plane of graphene is inactive. It is demonstrated by the homogenous distribution of the electrostatic potential (Fig. 2a). On the other hand, graphene is highly electrochemically active at the edge plane and near defect sites. 103 In the case of defective GQDs, the MEP exhibits some local inhomogeneities near the vacancy defects ( Fig. 2b-e). We note that we found a significant difference between the structures with odd and even numbers of missing carbon atoms. The spatial regions within monovacancy and trivacancy defects become more negatively charged than those within divacancy and tetravacancy defects. From the point of view of the GQD reactivity, the presence of negatively charged regions means that the vacancy defects with odd numbers of removed carbon atoms have a preferential ability to attract electrons belonging to elemental heavy metals compared to the structures with even numbers of missing carbon atoms (divacancy and tetravacancy). In other words, negatively charged carbon atoms around the monovacancy and trivacancy would allow enhancement of the interaction strength between graphene quantum dots and neutral metal adsorbates. A similar trend was observed for vacancy-containing graphene nanofragments interacting with Pt. 96 To support the aforementioned findings and to elucidate the nature of the enhanced chemical reactivity we further consider the effect of vacancy defects on the HOMO-LUMO energy gap (the difference between the energies of the lowest unoccupied molecular orbital and the highest occupied molecular orbital) of the GQDs. Table 1 summarizes the calculated dipole moment, HOMO (LUMO) energy position and HOMO-LUMO energy gap for each considered structure. We found that the main parameters of the molecular systems that we used to mimic GQDs are strongly correlated with the chemical reactivity. The calculated energy gap of the defect-free GQDs is approximately 3.07 eV, suggesting a pronounced semiconductor behavior of the GQDs. While the energy gap is found to decrease, when a vacancy defect is added. It is generally accepted that the HOMO-LUMO energy gap is related to the chemical hardness of the molecular system. In this regard, the systems with larger HOMO-LUMO gaps are more stable than those with small energy gaps. This is due to the fact that the narrow-band-gap materials have a reduced resistance to polarization of the electron cloud after interaction with external adsorbates. We also noticed that the energy gap of ONVTGQDs decreases much more compared with the even-numbered-vacancy-type GQDs (ENVTGQDs). One can assume that the presence of negatively charged regions in the defective GQDs plays an important role in providing high chemical reactivity. Naturally, the vacancy-induced HOMO-LUMO gap shrinking is explained by lowering of the LUMO level and upward movement of the HOMO level. Since the LUMO energy is the lowest energy orbital that can accept electrons, one can expect that its lowering will promote the charge transfer from the metallic adsorbate to the GQDs. Bearing this in mind, from Table 1 it is clearly seen that the GQDs with monovacancy and trivacancy defects are highly desired materials for strong chemisorption of neutral metal adatoms. The high magnitudes of the induced dipole moment for these systems would also translate into their enhanced ability to interact chemically with adsorbates. Finally, we revealed that the vacancy defects in GQDs play a role of perturbations that break the degeneracy of the electronic states near the HOMO and LUMO levels, thereby causing a lowering of the energy of unoccupied states and raising that of the occupied states.
To find the equilibrium structures and the position of the adsorbed heavy metals we carried out geometry optimization using relaxed structures of the vacancy-containing GQDs presented in Fig. 1. By positioning heavy metal adatoms (Cd, Hg and Pb) above the centre of the vacancy defect, we investigated the favourable binding sites and binding energies. As can be seen from Fig. 3a, f and k, the physisorption of Cd, Hg and Pb on defect-free graphene quantum dots, which is predominantly regulated by van der Waals interactions, suggests no significant structural changes and therefore the GQDs complexed with HMs retain their planar structure. Expectedly, among the three heavy metals, our calculations yielded the lowest total energy for Cd and Hg atoms located at the hollow site (the center of the hexagonal ring), and 3.41 Å and 3.29 Å above the surface, respectively. While the Pb adatom tends to occupy the bridge site (above the centre of the C-C bond) at the equilibrium height 2.81 Å. It is interesting to note that the heavy metal adsorption onto defective GQDs induces a strong buckling into the system, propagating away from the adsorption site (Fig. 3). Such structural modifications accompany strong chemical bonding. We found that Pb adatoms located above the vacancy defects as if they saturate the vacancyinduced dangling bonds form shorter bonds with the closest carbon atoms belonging to GQDs than other metals. Even though all heavy metals for each considered configuration behave as electron donors (Table 2), the electron transfer occurring from Pb to GQDs is always significantly enhanced. The difference in the charge transfer between GQDs and Cd 0 , Hg 0 and Pb 0 adatoms can be understood in terms of the ionization energy of the metal. 104 The ionization potential of Pb 0 is approximately equal to 7.41 eV, 105 whereas Cd and Hg have larger ionization energies of 8.99 eV and 10.43 eV, 105 respectively. Thus, Pb 0 transfers electrons much more easily to the GQDs than Cd 0 and Hg 0 adatoms. This clearly explains why the Pb adatom behaves differently compared to the other atoms, which form weaker chemical bonds with the carbon atoms. As was predicted by the MEP mapping analysis, the charge transfer between GQDs with trivacancy defects and heavy metals reaches the highest value in comparison to the other interacting systems. We note that the aforementioned discussion is mainly related to explanation of the difference in partial atomic charges on heavy metals calculated by Mulliken population analysis. Nevertheless, as was mentioned before the results of MPA analysis are still questionable and the accuracy of this method needs to be addressed. In principle, the accuracy of population analysis is typically estimated by comparing it to data from experimental measurements. 106 To quantify the match between predicted and real data the experimental dipole moments of the molecules of interest can be used. 107 Nevertheless, due to the lack of high accuracy experimental benchmark data, it is a challenge to validate the existing charge analysis models. In our case, the results of supercell calculations using periodic boundary conditions to describe the charge transfer between graphene and elemental heavy metals can be, to some extent, regarded as reference data. According to the results of the Bader charge analysis of the 4 Â 4 graphene supercell interacting with heavy metals, 44 the atomic charges on Cd, Hg and Pb were estimated to be approximately +0.02, +0.05 and +0.27, respectively.  It means that they behave as typical electron-donating adsorbates, providing effective charge transfer from the metal to the graphene. The tendency of charge transfer is consistent with that of ionization potential of the considered metals, implying that Pb can transfer electrons more easily to the graphene than other metals. Table S1 and Fig. S1 (ESI †) describe the calculated atomic charges on elemental heavy metals bound to the vacancy-free and vacancycontaining GQDs, in which the charge transfer directions predicted by Mulliken population analysis for all three metals on vacancy-free GQDs are similar to reference data. The results of Hirshfeld analysis and NPA are quite surprising because they contradict in some way reference data (with some exceptions). The Hirshfeld charge is found to be in good agreement with the Bader charge only for cadmium, while the NPA model well describes only the lead case. In other cases, these two models predicted the negative charges on heavy metals and, consequently, electron-accepting behaviour. Despite the numerical instability (dependence on the basis set) of the Mulliken population analysis, this model of partial atomic charge gives a reasonable picture of charge transfer processes between elemental heavy metals and GQDs and can be effectively used to understand the overall trend of the binding ability of GQDs to HMs. Having determined the relaxed structures, we next estimate the binding energy. All calculated parameters are summarized in Table 2. It is obvious that the artificial creation of vacancies in GQDs causes a significant increase of the binding energy of elemental heavy metals, allowing switching from a weakly physiosorbed state to a strongly chemisorbed state. These changes in the nature of the interaction between metal adsorbates and GQDs can be interpreted in terms of redistribution of the components contributing to the total binding energy. Considering the case of adsorption on defect-free GQDs, it is reasonable to assume that the interaction between Cd and Hg and GQDs occurs predominantly through dispersion forces, while the binding of Pb includes both the charge transfer and van der Waals interactions. On the other hand, the chemisorption of metals on defective GQDs involves also electrostatic interactions. Due to this reason, the binding energies are almost an order of magnitude higher compared to those describing the metal adsorption on the basal plane of defectfree GQDs. A significant point to note, however, is that the charge transfer contribution to the binding energy is more pronounced for Pb than Cd and Hg. Except for Cd and Hg adsorbed on GQDs with trivacancy defects, the binding in all remaining structures involves strong electrostatic interactions. It is not too surprising since Pb, by virtue of its low ionization potential, is able to easily donate electrons to GQDs, independently of the type of vacancy defects. Moving from the adsorbatefree defective GQDs (Table 1) to such with adsorbed HMs ( Table 2) we notice that similarly to the previous case the increase of the removed number of carbon atoms causes significant HOMO-LUMO gap shrinking after Cd (Hg) adsorption. While HOMO-LUMO gap widening takes place for defective GQDs complexed with Pb adatoms. This is mainly due to the appearance of the hybridization states in vacancy-mediated interacting systems (this will be discussed below).

Interbrand absorption
Bearing in mind the knowledge about the adsorption ability of defective GQDs to heavy metals, we are aiming to understand the fundamental correlation between the binding energies of HMs and the optical properties of vacancy-containing GQDs. Such a holistic understanding will promote the correct interpretation of the experimental absorption spectra, thereby facilitating design of high-performance sensing devices. At the first step of our analysis we resolve the absorption spectra into partial contributions corresponding to the most probable excited interband transitions between occupied and unoccupied energy levels. Then we focus on the orbital composition analysis, which provides information about the nature of the molecular levels involved in these transitions.
3.2.1 Excited states in defective GQDs. The resulting spectra of the adsorbate-free defective GQDs are shown in Fig. 4a, while corresponding peak assignments can be seen in Tables S2-S6 (ESI †).
From this figure we conclude that the intensity of the absorption peak drastically decreases with increasing the vacancy defect size and the absorption band line-width tends to become broader. This effect can be attributed to the reduction of the oscillator strength caused by the defect-induced symmetry breaking. In addition, we notice the difference between absorption spectra of structures with odd and even numbers of missing carbon atoms. Particularly, ONVTGQDs exhibit a shift of the absorption peak towards the long-wavelength side and larger quenching of the absorption intensity, as is seen in Fig. 4a and 5(a-e). It is noteworthy that the red-shift of the absorption spectra in the visible range leads to changes in perceived coloration of the GQDs from a mustard green colour (defect-free GQDs) to dark blue (GQDs with monovacancies and trivacancies, respectively) and even black (GQDs with tetravacancies). While the system with the divacancy demonstrates no coloration changes compared to defect-free GQDs. The corresponding simulated colours are displayed in Fig. 4g. To understand the origin of the absorption peaks we performed DOS calculations for all considered structures. Since the HOMOÀ1 and HOMO as well as the LUMO+1 and LUMO are degenerate p and p* orbitals in defect-free GQDs (Fig. 5a), the absorption peak at 412 nm originates from the degenerate multilevel electronic transitions from the ground state to the S3 and S4 singlet excited states. Molecular orbitals which are involved in these transitions are shown in the bottom panel of Fig. 5a. Wavefunctions of the aforementioned molecular orbitals are highly delocalized over the GQDs. Contrary to the doubly degenerate LUMO and HOMO in defect-free GQDs, the lowest orbitals in vacancy-containing GQDs split into nondegenerate orbitals (see Fig. 5b-e). Such a splitting is responsible for the change in the nature of the electronic excitations: from the doubly degenerate excited state in defect-free GQDs to a set of excited states in defective GQDs. The detailed description of the electronic transitions which involve the corresponding excited states is provided in ESI. † To shed more light on the nature of the interband absorption in GQDs we computed the transition density matrix 108 for the lowest excited states in each GQD (Fig. 6a-e). It should be mentioned that in most cases the excited state can have two components: a locally-excited (LE) state and a charge-transfer (CT) state. 109 Detailed analysis of the matrix elements allows us two distinguish between them. Indeed, the diagonal terms of the transition density matrix correspond to the charge variation of the corresponding atoms and, consequently, represent the LE character of the excited state, while the off-diagonal terms are related to the CT component exhibiting electron-hole coherence between the corresponding atoms. 110 The 2D-grid map for the defect-free GQDs (Fig. 6a) indicates the LE character dominates in the S4 state. Additional evidence for the correct determination of the type of electron excitation can be found through comparing the charge-transfer length (Dr), the electron-hole wave-function overlap integral (S) and the distance between centroids of holes and electrons (D) for the corresponding excited state. 111 All these parameters are listed in Dataset 1 (ESI †). Dr and D of the S4 state are negligibly small, as shown in Dataset 1, while this state is characterized by a large overlap integral. It means that the electrons and holes are  localized without obvious charge transfer. This is because the graphene quantum dot is a highly delocalized system and thus a strong coherence between sp 2 -bonded carbon atoms is observed. It is interesting to note that the S0 -S4 transition preferentially involves ''edge'' carbon atoms labeled 25-54, while central carbon atoms (labeled 1-24) belonging to inner hexagonal rings are less involved in excited transitions. In the case of ONVTGQDs (Fig. 6b and d), the contribution of the CT component to the lowest excited states, namely S5 and S11 (for monovacancy and trivacancy), becomes larger compared to that of the defect-free GQDs and the bright regions within the corresponding maps are distributed more homogeneously over the GQD plane. In line with this, S5 and S11 states have larger charge-transfer lengths (0.71 Å and 1.16 Å, respectively), and smaller overlap integrals. Therefore, one can classify these states as LE-CT hybridized states. When two carbon atoms are removed from a GQD, the brightest zones are concentrated across the diagonal of the color-filled map (Fig. 6c), indicating the partial LE character of the S10 state. This implies preferential participation of the edge carbon atoms labeled 41-52 and inner carbon atoms labeled 16-21 (which surround the divacancy complex) in optical transitions. We also noticed that due to the large charge-transfer length (1.935 Å) the S5 excited state in GQDs with tetravacancies (Fig. 6e) can be explained as originating from the formation of a hybridized state with increased contribution from the CT component.
3.2.2 Excited states in defective GQDs interacting with elemental HMs. In comparison to the adsorbate-free defective GQDs, the charge transfer from metal adsorbates to GQDs leads to changes in the DOS spectra of GQDs and, consequently, to a redistribution of the molecular levels involved in the excited states. It clearly manifests itself in many ways, from coloration changes (Fig. 4e) to signal quenching (Fig. 4g). Furthermore, we revealed a direct correlation between the perceived colour and binding energy. Weakly interacting systems (Cd or Hg adatoms with defect-free and divacancy-containing GQDs) do not change their coloration compared to the corresponding adsorbate-free GQDs. Contrary to this, the complexes with high interaction energy exhibit large difference in the perceived colour. Converting our theoretical data into colorimetric parameters, we also noticed that the ONVTGQDs are more desirable systems for colorimetric sensing of Cd and Hg. Since Pb interacts with defect-free and vacancy-containing GQDs in a much stronger manner than Cd and Hg, we observed completely different pictures of the interband absorption with large contributions from hybridized and localized energy levels to excited states. This is consistent with the coloration of GQDs complexed with Pb. Fig. 4g shows that monitoring of the interband absorption in defective GQDs is a sensitive method for metal identification, since f 0 /f increases after interaction with metals in a unique manner depending on the type of the vacancy defect. In this case, the oscillator strength ratio can be related to the quenching efficiency. From a practical point of view, this means that sensor arrays consisting of GQDs with different types and densities of vacancy defects can enable discriminative multimetal analysis. To move toward a more detailed description of the excited transitions in ONVTGQDs and ENVTGQDs after adsorption of elemental heavy metals, in the following we consider the features of the interband absorption for each of the considered metal adsorbates. The main trends we noticed are that the charge-transfer length (Dr) for the most probable excited states increased significantly with increasing the size of vacancy-type defects in GQDs interacting with Cd and Hg (see Dataset 1, ESI †), while Dr was found to decrease for excited states in defective GQDs after lead adsorption. This implies the fundamental role of carbon vacancies in the redistribution of the charge density during the interband transition in strongly interacting complexes. Therefore, it is important to gain deep insights into the nature of the excited states depending on the type of vacancy and kind of heavy metal.
3.2.3 Excited states in defective GQDs interacting with elemental Cd. As can be seen from Fig. 4b, Fig. 7(a-e) and Table 2, the interband absorption in GQDs complexed with neutral Cd adatoms depends on the binding energy. The increase in the number of removed carbon atoms causes a stronger interaction between Cd and GQDs and significant HOMO-LUMO gap shrinking, followed by the red-shift of the absorption peak, except in the case of GQDs with divacancies (where a slight blueshift is observed). It is interesting to note that the adsorbate-free GQDs with monovacancies and trivacancies yielded a larger redshift of the broad excitation band than similar GQDs with adsorbed Cd (Fig. 7). This is mainly due to the vacancy-assisted strengthening of the interaction between Cd and GQDs. For the defect-free GQDs complexed with Cd through weak van der Waals interactions (Fig. 7a), the absorption spectrum is dominated by an intense peak at 413 nm with an oscillator strength of 0.9509. Like the pristine GQDs, two degenerate excited states contribute to this absorption band (S5 ' S0 and S6 ' S0). The corresponding electronic transitions which are involved in these excited states are listed in Table S7 (ESI †). Wave functions for the molecular orbitals are highly delocalized over the defect-free GQDs (Fig. 8a). Fig. 9a indicates that the electron-hole coherence for the S5 state (dominated by p-p* transitions) is quite strong Fig. 6 Two-dimensional colour-filled representation of the TDM for the most probable S4 ' S0, S5 ' S0, S10 ' S0, S11 ' S0, and S5 ' S0 excited states in GQDs without interaction with elemental heavy metals: (a) no defects, (b) monovacancies, (c) divacancies, (d) trivacancies and (e) tetravacancies. A colour scale bar is shown on the right of each contour plot. A brighter colour indicates a higher probability of finding the generated electron-hole pair.
within the inner carbon atoms labelled 1-6. The electrons and holes are all localized in the bottom left corner of the map, and therefore the S5 state can be classified as a LE state. The absorption spectra of the GQDs with monovacancies and divacancies after interaction with elemental Cd originate from the S11 excited states observed at 444 nm and 410 nm (Tables S7 and  S8, ESI †), respectively. The 444 nm-peak is associated with a combination of electronic transitions: HÀ6 -LUMO (15%), HÀ2 -L+1 (25%), HÀ1 -L+2 (29%), and HÀ1 -L+3 (26%). We noticed the hybrid nature of the HÀ2 level, with partial contribution of Cd, as low as 9% (Fig. 8b), while the strongly hybridized HOMO level is not involved in these transitions. According to the TDM map for the S11 excited state (Fig. 9b), the electrons still strongly cohere with holes within inner carbon atoms labelled 1-6, but the contribution of the intermolecular charge transfer component to the excited state increases. The nature of the S11 state in GQDs with divacancies can be interpreted as a combination of the following electronic transitions: HÀ5 -LUMO (10%), HÀ1 -L+2 (50%), and HOMO -L+1 (24%), while the second less intense feature at 426 nm (the S10 state) can be assigned to HÀ7 -LUMO (23%), HÀ1 -L+1 (47%) and HOMO -L+2 (23%) transitions (Table S9, ESI †). The orbital composition analysis provides evidence that the hybridized HÀ1 (Fig. 8c) and LUMO molecular orbitals are involved in these transitions. The contribution of Cd to these levels is found to be 41% and 5%, respectively. Since the Cd-related occupied HÀ1 molecular level plays an important role in the S11 excited state, the transition density matrix should reflect this role. Indeed, as can be seen from Fig. 9c intermolecular charge transfer dominates over the LE component. It is obvious that electronhole coherence occurs at both the right and top edges of the map. In other words, S11 is an intermolecular electric charge transfer state, in which electrons and holes cohere between cadmium and outer carbon atoms belonging to GQDs. Analysis of the absorption spectrum of the GQDs with divacancy defects indicates that the S11 excited state has a larger electronic transition probability compared to other excited states ( Fig. 7d and Table S10, ESI †). The most notable result of the complexation of trivacancy-containing GQDs with elemental Cd is the appearance of the hybridized CT-LE state (Fig. 9d). Since the hybridized LUMO level participates in electronic transitions, we observe that electrons and holes cohere between cadmium and edge carbon atoms. In addition, the bright regions in the left bottom corner of the map suggest both electrons and holes are localized on the inner carbon atoms. The interband absorption in GQDs with tetravacancies ( Fig. 7e and Table S11, ESI †) is characterized by the appearance of two dominating features at 507 nm (with an oscillator strength of 0.1892) and 632 nm (0.1097). The band at 507 nm is assigned to the HOMO -L+2 transition, while the most likely transitions contributing to the 632 nm absorption are HÀ2 -LUMO (23%), and HÀ1 -L+1 (61%). The full assignments of the corresponding transitions are summarized for each transition in Table S11 (ESI †). The orbital compositional analysis suggests that among all aforementioned molecular orbitals only HÀ2 is highly hybridized (Cd contribution, 74%), whereas the wave-functions of the remaining molecular orbitals are delocalized over the defective GQDs. Since the Cd-related orbitals are not involved in the S7 excited state (Fig. 8e), both electrons and holes are localized on the inner carbon atoms, thereby indicating the dominating role of the LE state.
3.2.4 Excited states in defective GQDs interacting with elemental Hg. Fig. 7(f-j) illustrates the transitions from the occupied electronic states to the empty levels, which contribute to the intensive peak in the absorption curves demonstrated in Fig. 4c. Based on the analysis of the existing transitions in defect-free GQDs interacting with Hg (Fig. 7f), one can observe that the dominant contribution to the S4 excited state at 413 nm arises from the combination of the electronic transitions: HÀ1 -LUMO (25%), HÀ1 -L+1 (22%), HOMO -LUMO (22%), and HOMO -L+1 (25%). The wave-functions for all these molecular orbitals are delocalized over the GQD (Fig. 8f). Like the adsorbate-free GQD case, the 2D-grid map for defectfree GQDs interacting with Hg (Fig. 9f) indicates the LE character dominates in the S4 state. Nevertheless, electron-hole coherence occurs mainly within the inner carbon atoms in the central part of the GQDs. Three intense features can be distinguished in the absorption spectrum of the monovacancy-containing GQDs complexed with Hg (Table S12, ESI †). The dominant feature at 416 nm (with an oscillator strength of 0.360) can be ascribed to HÀ1 -L+1 electronic transitions. According to the orbital compositional analysis, neither HÀ1 nor L+1 have contributions from Hg (Fig. 8g). Only GQD-related levels are involved in these transitions. Since Hg-related orbitals do not participate in the electronic transitions, the LE character is prevailing in the S8 state (Fig. 9g). In the case of the divacancy the absorption peak is blue-shifted by 16 nm (Fig. 7h), and the nature of the electronic transitions becomes more complicated (Table S14, ESI †). This is due to the fact that the hybrid HÀ6 (GQD: 89%, Hg: 11%) orbital starts to participate in the excited transitions. The most intense peak at 400 nm is due to the combination of two close transitions, such as HÀ1 -L+2 (54%), and HOMO -L+1 (14%). These frontier molecular orbitals are completely delocalized over the GQD (Fig. 8h). The mercury adatom does not participate in the optical excitation S11 ' S0 (Fig. 9h). In this case, the nature of excited transitions can be understood in terms of a strongly hybridized CT-LE state. On one hand, we observe a large electron transfer (off-diagonal coherence) from the edge carbon atoms to the inner carbon atoms and vice versa. On the other hand, one can see the bright regions in the bottom left corner of the map and the top right corner of the map (diagonal coherence), indicating partial LE character of the S11 excited state. A further increase in the number of the removed carbon atoms (from 2 to 3) causes localization of the HOMO and the LUMO level mainly on the Hg adsorbate (Fig. 8i), followed by a significant drop in the absorption intensity and a red-shift of the absorption band from 400 to 479 nm (Fig. 7i). The adsorption spectrum of tri-vacancy containing GQDs is dominated by the single excited state S12 with an oscillator strength of 0.2915. The nature of this state can be attributed to the following transitions: HÀ2 -L+1 (23%), HÀ1 -L+1 (14%), and HOMO -L+2 (38%). Because of the appearance of the hybridized HOMO (the Hg contribution is about 86%) level, we revealed off-diagonal coherences between Hg and edge carbon atoms, which implies that the CT character dominates in the S12 state. Being complexed with Hg tetravacancy-containing GQDs demonstrates an observable decrease in the oscillator strength of the major transitions. We can ascribe this decrease to the non-symmetric nature of the GQD, which manifest itself in large splitting of the lowest occupied and highest unoccupied levels belonging to the GQD into a set of well-distinguished energy levels (Fig. 7j). As an important consequence, the transitions from HÀ2 to the LUMO (23%) and from the HOMO to L+2 (45%) become allowed and generate the most intense absorption feature at 511 nm (with an oscillator strength of 0.1219). The second peak with an oscillator strength of 0.1168 appears at 617 nm and involves the transitions from HÀ1 to L+1 (83%) and from the HOMO to L+2 (11%). Among these molecular levels only HÀ2 has a small contribution from the Hg adsorbate (less than 19%). Analysis of the TDM map indicates that the S5 excited state is dominating the matrix elements (the bottom left corner of the map), implying a dominant contribution of the LE component (Fig. 9j).
3.2.5 Excited states in defective GQDs interacting with elemental Pb. Moving our attention to the absorption spectra of GQDs complexed with Pb, in the following we gain insights into the origin of the electronic transitions in the presence of vacancy defects. The corresponding absorption spectra are illustrated in Fig. 4c, while peak assignments are presented in Tables S16-S20 (ESI †). In the absorption spectrum of defectfree GQDs interacting with Pb three weak bands at 539, 584 and 672 nm can be observed (Table S17, ESI †). The most intense peak at 584 nm can be identified as an overlapped multi-level transition, including contributions from HOMO -L+7 (14%), HOMO -L+9 (16%) and HOMO -L+10 (59%). In this case, hybrid HOMO and L+10 levels are involved in the S10 excited state (Fig. 8k). As can be seen for the TDM map, the nature of the S10 excited state can be interpreted as a hybridized state with large contribution from the CT component. As can be seen from the Dataset (ESI †) this state is characterized by the large value of Dr (2.13 Å) and small value of the overlap integral (0.29). The dualistic nature of this state is related to the fact that the electrons and holes are localized on the inner carbon atoms labelled 1-21 (see the bright region in the bottom left corner of the map) and the electron-hole coherence is strong because of the electron transfer between Pb and GQDs (see the bright regions across the right edge and top edge of the map). Like in the vacancy-mediated interband absorption in GQDs interacting with Cd and Hg, the vacancies also play an important role in the redistribution of the excited states in GQDs after complexation with Pb. We explain these phenomena by vacancy-induced strengthening of the charge transfer and electrostatic interaction between Pb and GQDs. A clear blueshift of the most prominent absorption peaks can be seen in Fig. 7(k-o). Passing from the pristine GQDs to monovacancycontaining GQDs, we noticed that the absorption spectrum is mainly dominated by a peak at 467 nm with an oscillator strength of 0.1896. A detailed analysis of the calculated results indicates that three transitions contribute to the S10 excited state, including HÀ7 -LUMO (26%), HÀ6 -LUMO (26%) and HOMO -L+1 (27%). These transitions contribute almost equally to the absorption peak. The contribution of Pb to the LUMO, HÀ6, and HÀ7 is found to be negligibly small (Table S18, ESI †). Fig. 8l illustrates the contributing orbitals which are involved in the corresponding transitions. In this case, we noticed the presence of the CT-LE state with increased contribution from the LE state (Fig. 9l). This is also evidenced by a simultaneous increase in the overlap integral and a decrease in the charge-transfer length (see Dataset 1, ESI †). Further increase of the number of removed carbon atoms from 1 to 2 and 3 leads to a subsequent shift of the absorption peak of Pb 2+ @GQD towards the short-wavelength side and an enhancement of the absorption intensity, as seen in Fig. 4c. Considering the absorption spectrum of GQDs with divacancies one can see that two bands, at 421 and 441 nm, can be distinguished. The dominating peak at 441 nm is related to the following transitions: HÀ1 -L+1 (60%) and HÀ1 -L+3 (15%). The matrix elements associated with the S8 state (Fig. 9m) are rather homogeneously distributed across the map, suggesting a large degree of hybridization between CT and LE modes. The absorption spectrum of the trivacancycontaining GQDs complexed with Pb mainly originates from two excited states -S6 at 482 nm and S11 at 409 nm. The S6 state is dominated by two transitions, from the HOMO to L+2 (61%) and from the HOMO to L+3 (15%). The Pb is contributing to both L+2 (10%) and L+3 (32%). The last feature at 409 nm is dominated by three transitions from the GQD-related occupied levels to unoccupied levels: HÀ3 -LUMO (14%), HÀ1 -L+2 (36%) and HOMO -L+4 (20%). For L+4, the contribution of Pb is approximately 30%. The transition density matrix describing the S11 state shows strong mixing between LE and CT states (Fig. 9n) implying a high degree of hybridization. It is evidenced by the simultaneous presence of electron-hole coherence in both diagonal and offdiagonal directions (bright regions in the bottom corner of the map as well as off-diagonal matrix elements corresponding to the charge transfer between inner and outer carbon atoms). Finally, due to significant structural modifications, which appear in tetravacancy-containing GQDs after adsorption of elemental Pb adatoms, the interband absorption is characterized by small oscillator strengths for the main transitions (see Table S21, ESI †). The major feature at 476 nm corresponding to the S10 excited state (see Fig. 7o) is related to a mixture of three transitions: HÀ2 -L+1 (11%), HÀ1 -L+2 (61%) and HOMO -L+1 (17%). The wavefunctions for all these orbitals are shared between the Pb adsorbate and the defective GQD (Fig. 8o). The matrix elements related to the S10 state are distributed in both the diagonal and off-diagonal directions (Fig. 9o), indicating the hybridized nature of this state. The brightest region is located in the bottom left corner of the map, representing the dominating role of LE character, which includes inner carbon atoms labelled 1-15.

Conclusions
We have employed DFT and time dependent (TD) DFT methods to theoretically investigate the structure, electronic and optical properties of vacancy-containing closed-shell GQDs (monovacancy, divacancy, trivacancy and tetravacancy) complexed with elemental heavy metals (Hg, Cd and Pb). We have provided evidence that the vacancy engineering can promote a strong interaction between graphene quantum dots and metal adsorbates, which allows for better optical detection of Cd, Hg and Pb. We found that the monovacancy and trivacancy defects are very active for the reaction with HMs compared to the defect-free and even-numbered vacancy type GQDs, which originates from the presence of negatively charged regions facilitating the charge transfer between metal adsorbates and GQDs. Since the interaction mechanism between Pb and GQDs is mainly governed by charge transfer, the adsorbed Pb atoms tend to bind more strongly near vacancy defects compared to Cd and Hg. An important finding is that the vacancy-assisted interaction with HMs causes atomic-scale buckling of the initially flat sp 2 plane of GQDs, which becomes larger with an increase in the size of the vacancy defect. Therefore, a trade-off between the stability of GQDs and their reactivity must be reached to provide effective local modulation of the density of active surface sites. We revealed the fundamental correlation between the coloration of the GQDs and interaction strength of the HM@GQD complexes implying that the color changes can be distinguished only for strongly interacting systems. Furthermore, we find that there are two different mechanisms of interband absorption in GQDs interacting with HMs depending on the presence of the vacancy defects: the locally-excited character dominates the lowest excited singlet transitions in the defect-free GQDs complexed with Cd and Hg, while the hybrid coupled states with chargetransfer and locally excited components are more favorable for vacancy-containing GQDs after adsorption of Cd and Hg. Through the analysis of the excited states in GQDs interacting with Pb, we found strong charge-transfer channels between defect-free GQDs and Pb as well as a decrease of the chargetransfer length in the presence of vacancies. The transformation of the excited state character from the pure locally-excited (LE) state and/or pure CT state to the hybridized state is found to be accompanied by redistribution of the molecular energy orbitals which are involved in the allowed electronic transitions as well as participation of the HM-related orbitals affecting the electronhole coherence under excitation conditions. Our work provides deep insights into the nature of the mechanisms underlying the vacancy-mediated interaction between graphene quantum dots and elemental heavy metals. The output from this advanced understanding will promote the fabrication of colorimetric sensor arrays based on vacancy-containing GQDs.

Conflicts of interest
There are no conflicts to declare. Dr I. Shtepliuk acknowledges the support from Wallenberg foundation and Ångpanneföreningens Forskningsstiftelse (Grant 16-541). All calculations were performed using the supercomputer resources of the Swedish National Infrastructure for Computing (SNIC), National Supercomputing Center (NSC).