How strong is the Second Harmonic Generation in single-layer monochalcogenides? A response from first-principles real-time simulations

TL;DR

This study uses first-principles real-time simulations to investigate the origin and magnitude of second harmonic generation in single-layer monochalcogenides, revealing that many-body interactions are not responsible for the strong experimental signals.

Contribution

The paper provides the first ab-initio real-time simulation analysis of SHG in single-layer monochalcogenides, clarifying the role of excitonic effects and reconciling experimental and theoretical results.

Findings

Theoretical SHG intensity is an order of magnitude smaller than experimental reports.

Many-body interactions are not the primary cause of the strong SHG observed experimentally.

Using a single-layer model aligns experimental data with theoretical predictions.

Abstract

Second Harmonic Generation (SHG) of single-layer monochalcogenides, such as GaSe and InSe, has been recently reported [2D Mater. 5 (2018) 025019; J. Am. Chem. Soc. 2015, 137, 79947997] to be extremely strong with respect to bulk and multilayer forms. To clarify the origin of this strong SHG signal, we perform first-principles real-time simulations of linear and non-linear optical properties of these two-dimensional semiconducting materials. The simulations, based on ab-initio many-body theory, accurately treat the electron-hole correlation and capture excitonic effects that are deemed important to correctly predict the optical properties of such systems. We find indeed that, as observed for other 2D systems, the SHG intensity is redistributed at excitonic resonances. The obtained theoretical SHG intensity is an order of magnitude smaller than that reported at the experimental level. This result is in substantial agreement with previously published simulations which neglected the electron-hole correlation, demonstrating that many-body interactions are not at the origin of the strong SHG measured. We then show that the experimental data can be reconciled with the theoretical prediction when a single layer model, rather than a bulk one, is used to extract the SHG coefficient from the experimental data.

Figures

Tables

Table 1. Table 1: Static limit of the SHG, χ ( 2 ) ​ ( ω = 0 ) superscript 𝜒 2 𝜔 0 \chi^{(2)}(\omega=0) , obtained from real-time simulations at different level of approximation: independent-particles approximation (IPA) and the random-phase approximation (RPA) on top of the DFT electronic-structure, the IPA ( G 0 ​ W 0 subscript 𝐺 0 subscript 𝑊 0 G_{0}W_{0} ) and the real-time BSE (G W 0 0 subscript subscript 𝑊 0 0 {}_{0}W_{0} +BSE) on top of the G 0 W 0 electronic-structure. Units in pm/V.

	${\chi }^{\left(2\right)}$(IPA)	${\chi }^{\left(2\right)}$ (RPA)	${\chi }^{\left(2\right)}$ ($G_0{W}_0$)	${\chi }^{\left(2\right)}$ ($G_0{W}_0$+BSE)
GaSe	110	71.0	54.0	81.2
InSe	135	86	68.0	122

Equations

$${\text{P}}_i={\chi }_{ij}^{(1)}{\mathcal{E}}_j+{\chi }_{ijk}^{(2)}{\mathcal{E}}_j{\mathcal{E}}_k+{\chi }_{ijkl}^{(3)}{\mathcal{E}}_j{\mathcal{E}}_k{\mathcal{E}}_l+O({\mathcal{E}}^4),$$

$$\frac{{\chi }_{\text{bulk}}^{\text{(2)}}}{{\chi }_{\text{sheet}}^{\text{(2)}}}=32\pi F\frac{n_{\text{2D}}(\omega )\sqrt{n_{\text{2D}}(2\omega )}}{(n_{\text{sub}}+1{)}^3},$$

Full text

I Introduction

Nowadays there is a considerable interest in the excited state properties of 2D materials which can be distinctly different from that of their bulk counterpart.Autere et al. (2018) Strongly bound excitons often characterize the optical spectra of two-dimensional (2D) semiconducting crystals and in several cases, intense nonlinear optical spectra have been observed  Wang et al. (2015); Säynätjoki et al. (2017) such as strong second-harmonic generation Malard et al. (2013); Kumar et al. (2013); Janisch et al. (2014); Ribeiro-Soares et al. (2015); Zhou et al. (2015) (SHG)—typically up to one order of magnitude larger than in conventional nonlinear crystals. For this reason these materials are of potential technological interest as frequency converters in nanophotonic circuits. Majumdar et al. (2015) Furthermore, the SHG is an excellent spectroscopic tool for the imaging and characterization of 2D materials Li et al. (2013); Säynätjoki et al. (2017); Malard et al. (2013); Kumar et al. (2013); Zhou et al. (2015); Hsu et al. (2014); Yin et al. (2014) and understand the origin of the signal from theoretical point of view is of fundamental importance for the interpretation.

Recently, extremely large values for the SHG of single-layer monochalcogenides have been reported: for InSe Zhou et al. (2018) a value of 6.39 nm/V at 1.55 eV and for GaSe Zhou et al. (2015) a value of 2.4 nm/V at 1.02 eV. To put these values into perspective, the SHG of single-layer MoS2 is of the order of 0.1–0.4 nm/V Malard et al. (2013); Clark et al. (2014) at resonance, that of bi- and tri- and multilayer Jie et al. (2015); Tang et al. (2016); Karvonen et al. (2015) GaSe ranges from 9–90 pm/V, and that of bulk GaSe is about 60 pm/V. A similar large SHG of 10 nm/V was reported for WS2 and WSe2 monolayers. Janisch et al. (2014); Ribeiro-Soares et al. (2015) However, the extremely large SHG coefficients for single-layer monochalcogenides are not confirmed by first-principle calculationsLin et al. (2018); Hu et al. (2017) which give consistently values smaller by one order of magnitude.

The discrepancy between the theoretical and experimental results may be due to many-body effects that were neglected in the computational studies of Refs. Lin et al., 2018; Hu et al., 2017. In fact, electron-electron interaction and excitonic effects are known to play a crucial role in the description of optical properties, especially in the case of low-dimensional materials. Indeed the incomplete screening of the electron-electron interaction and the quantum confinement may further enhance these effects with respect to bulk. Specifically for single-layer monochalcogenides, the formation of strong bound excitons have been recently shown in Ref. Antonius et al., 2018. Further, their peculiar optical properties have been attributed to saddle-points excitationsLuo et al. (2018) that originate from the Mexican hat-like band structure near the $\Gamma$ point.Rybkovskiy et al. (2014); Chen et al. (2018)

This work addresses the question of whether the extremely large SHG of single-layer monochalcogenides is due to strongly bound excitons. To this end, by using a first-principle real-time approach based on many-body theory (Sec. II), we calculate the electronic structure, dielectric function (Sec. III.1) and the SHG (Sec. III.2) of InSe and GaSe monolayers and, in order to single out many-body effects, we compare with the results at the independent-particle level. Previous studies on h-BN, MoS2 and other semiconducting 2D materials found a significant enhancement at resonance of the SHG signal due to excitonic effects—up to 2–4 times depending on the system. Grüning and Attaccalite (2014); Attaccalite et al. (2015) For single-layer monochalcogenides, we find that intensity is only slightly redistributed by many-body effects respect to the independent-particle approximation. We conclude (Sec. IV) that many-body effects do not account for the difference between theoretical predictionsLin et al. (2018); Hu et al. (2017) and experimental estimates.Zhou et al. (2018, 2015) We turn then the attention to the experimental results and argue that the extremely large SHG coefficient is an artifact of the model assumed to extract the estimate from the experimental data.

II Theoretical framework and computational details

The structural optimization and electronic structure have been calculated using density-functional theory (DFT). The crystal structure of the GaSe and InSe monolayers (Fig. 1) has been obtained from the corresponding bulk structure using the same in-plane lattice parameters ($a=0.3743$ nm and $a=0.4$ nm respectively for the GaSe Schwarz et al. (2007) and the InSe monolayer De Blasi et al. (1982)). Since experiments on GaSe and InSe were performed on different substrates, the lattice constant of these materials could depend on the substrate interaction. Therefore we decided to keep fixed the lattice constants to the bulk values and optimizing only the atomic positions that are not fixed by symmetries, namely the $z$ coordinates of the atoms. In order to simulate an isolated monolayer, we used a 35 a.u. supercell in the $z$-direction.

DFT calculations have been performed with the Quantum-Espresso codeGiannozzi et al. (2009) using the Perdew-Burke-ErnzerhofPerdew et al. (1996) (PBE) functional and the scalar-relativistic optimized norm-conserving VanderbiltHamann (2013) pseudopotentials from the PseudoDojo repository (v0.4).van Setten et al. (2018) The valence configuration for each pseudo-atoms are: $3d^{10}4s^{2}~{}4p^{1}$ for the Ga, $4d^{10}5s^{2}~{}5p^{1}$ for the In and $4s^{2}4p^{4}$ for Se. We used a shifted $18\times 18\times 1$ k-point sampling for the ground-state, a plane-wave cutoff of $90$ Ry for the structural optimization and $70$ Ry for the band-structure calculations.

Calculations of the quasi-particle energies and optical susceptibilities have been performed using the Yambo code.Sangalli et al. (2019) The quasiparticle corrections to the fundamental band gap have been calculated from the standard $G_{0}W_{0}$ approximationMartin et al. (2016) with the Godby-Needs plasmon-pole modelGodby and Needs (1989) and applied as a rigid shift to all the bands. We used a $24\times 24\times 1$ k-point sampling, a 5 Hartree cutoff for the dielectric function, and 200 total bands for Green’s function expansion. We calculate the gap correction at $\Gamma$ point and then shifted all conduction bands by this amount. In this way the quasiparticle band structure has the same band ordering and band-width of the DFT one, with the only difference that the gap has been corrected within the $G_{0}W_{0}$ approximation. Convergence with the number of conduction bands was accelerated by means of the Bruneval-Gonze terminator.Bruneval and Gonze (2008) Optical absorption spectra have been obtained by the solution of the Bethe-Salpeter equationStrinati (1988) using a basis of electron-hole pairs for which we considered 6 valence and 7 conduction states for both GaSe and InSe monolayers. We verified that increasing the number of bands in the calculations does not change the spectra in the energy range considered here.

Non-linear susceptibilities have been obtained from the real-time evolution of Bloch-electrons in a uniform time-dependent electric field following the approach proposed in Ref. Attaccalite and Grüning, 2013; Grüning and Attaccalite, 2014. The effective Hamiltonian for the Bloch-electrons, derived from many-body theory, contains both the electron-hole attraction (through a screened-Coulomb term) and exchange (through a Hartree term) needed to describe excitons and local-field depolarization effects. The whole framework is based on DFT and quasiparticle corrections are included at the level of $G_{0}W_{0}$. In the linear-response limit, the approach is equivalent to the solution of the BSE on top of the $G_{0}W_{0}$ electronic structure Attaccalite et al. (2011) and we therefore refer to this level of theory as $G_{0}W_{0}+$BSE for both linear and nonlinear susceptibilities. Within this framework, we can exclude specific terms from the effective Hamiltonian to investigate how they affect the SHG. In particular, we consider the independent-particle approximation (IPA)—no electron-hole correlation; the random-phase approximation (RPA)—only electron-hole exchange. Further, we perform the simulations on top of either the DFT or the $G_{0}W_{0}$ electronic structure. For the screened-Coulomb and Hartree terms, we use the same parameters as in the BSE calculations.

After the integration of the equation-of-motion, the nonlinear susceptibilities $\chi^{(n)}$ are extracted using Fourier techniques from the macroscopic polarization P—calculated as a dynamic Berry-phase Souza et al. (2004)—and its expansion in power of the total electric field $\bm{\mathcal{E}}$:

[TABLE]

where the $i,j,k,l$ subscripts refer to Cartesian components of the field and polarization.

Specifically, to obtain a single SHG spectrum we perform a series of simulations of about 75 fs for a monochromatic electric field with frequency ranging from 0.1 eV to 7.0 eV. We integrate the equation-of-motion using the Crack-Nicholson method Crank and Nicolson (1947) with a time-step of 0.01 fs. Since the turning-on of the electric-field excites all eigenfrequencies of the system, we introduce a phenomenological decoherence non-Hermitian operator corresponding to a Lorentzian broadening of the spectrum of 0.2 eV. 111It can be shown in Ref.Wismer and Yakovlev, 2018 that for a simple decoehrence operator, as the one used here, there is an equivalence between the real-time dynamics generated by a non-Hermitian Hamiltonian and the dephasing usually employed in the density matrix formalisms. However, the Hamiltonian formalism is more convenient in our case because the polarization can be calculated directly from the bloch-electrons evolution. The intensity of SHG depends on this phenomenological parameter and we chose a value compatible with experimental broadening.Another possibility would be to calculate the dephasing by including the coupling with the atomic motion explicitely in our simulations as done in Ref. Monserrat, 2018. The contribution of the eigenfrequencies to the signal decays exponentially in time with the decay constant being the inverse of the dephasing time. After five times the dephasing time (about 33 fs) the eigenfrequency signal is negligible with respect to the SHG signal and we then perform the Fourier analysis as detailed in Ref. Attaccalite and Grüning, 2013.

The static limit ($\omega=0$), corresponding to the outmost left point of Figs. 3,4 has been obtained by extrapolation of the values at small frequencies of the $\chi^{(2)}(\omega)$.

The external electric field is polarized along the $y$ direction and the polarization is recorded in the same direction, obtaining the $\chi^{(2)}_{yyy}=-\chi^{(2)}_{xxy}=-\chi^{(2)}_{yyx}=-\chi^{(2)}_{xyx}=-\chi^{(2)}_{aab}$ that is the only nonzero component of the second-order susceptibility tensor for the hexagonal $D_{3h}$ crystal class, Boyd (2008) being $a$ and $b$ the in-plane crystal axes. To obtain a SHG signal independent of the dimension of the supercell, we rescaled the calculated $\chi^{(2)}(\omega)$ by an effective thickness of $0.796$ nm for the GaSe and $0.832$ nm for the InSe—corresponding to half of $c$ lattice parameter of the bulk structures. Boyd (2008)

III Results

III.1 Electronic structure and optical properties

In Fig. 2, we report the electronic band structure of both monochalcogenide monolayers and the position of the lowest excitons obtained from the solution of the BSE. Because of the mirror-plane symmetry ($z\rightarrow-z$), electronic states near the edges are either even or odd with respect to this symmetry operation.Zólyomi et al. (2014) As a result, the lowest-energy electron-‐hole transition—depicted by a white arrow—is optically inactive for in-plane polarized light and only slightly coupled to $z$-polarized light.Antonius et al. (2018); Bandurin et al. (2017) In the same figure, we also indicate the lowest optical active transition with a black arrow.

The band structure is in good agreement with previous results obtained by Hu et al.Hu et al. (2017), while we found small differences respect to Antonius et al.Antonius et al. (2018) who reported a larger energy difference between the two top valence bands at the vicinity of $\Gamma$.

In Fig. 3 and Fig. 4 (bottom panels), we report the imaginary part of the dielectric function of GaSe and InSe for in-plane polarized light. The first excitons, obtained by full diagonalization of the Bethe-Sapeter equationsStrinati (1988); Martin et al. (2016) are reported as well. As one can see for both GaSe and InSe we are in presence of bound excitons (i.e. at energies below the fundamental band-gap).The lowest exciton at 3.39 eV in GaSe is dark for light polarized in the plane due to the symmetry of the bands close to the $\Gamma$ point (while it has a small but non-zero dipole for light polarization in the $z$ direction). The first bright exciton is found at 3.51 eV. Our results are in good agreement with those of Antonius et al.Antonius et al. (2018) For the InSe monolayer (bottom panels of Fig. 4) we obtained similar results. The lowest excitation at $2.47$ eV is dark for in-plane polarized light and the first bright exciton is at $2.63$ eV. This is in agreement with luminescence measurements that found that the lowest-energy transition is not optically active in InSe monolayer.Bandurin et al. (2017) Therefore, the lowest excitons do not contribute to the second-order susceptibility for in-plane polarized light.

III.2 Second Harmonic Generation

In Table 1, we report the static $\chi^{(2)}(0)$ at different level of theory: IPA, RPA, independent-particle on top of the G0W0 band structure (G0W0), and finally the full G0W0+BSE. At the independent-particle level we found a large static second order susceptibility for both systems. The susceptibility is sensitive to the electronic structure and in particular to the band gap: a larger band gap corresponds to a smaller $\chi^{(2)}(0)$ as one can see going from DFT to the G0W0. Also, the band structure of InSe is similar to that of GaSe and the gap at DFT level is 1.82 eV against 2.34 eV of GaSe. As a consequence InSe has a larger $\chi^{(2)}(0)$.

The GaSe result at G0W0 level is compatible with the result of Hu et al. that found a $\chi^{(2)}(0)=69.4$ pm/V starting from a DFT band structure obtained by means of the Heyd-Scuseria-Ernzerhof hybrid functional that has a gap of 3.34 eV, in-between the PBE 2.34 eV and the G0W0 4.03 eV band gap values. For the InSe monolayer, Lin et al.Lin et al. (2018) report a value of $91.3$ pm/V, compatible with our result, considering that they assumed a smaller effective thickness. Lin et al. also investigated other possible structures of the InSe monolayer with a stronger SHG but we are not aware of experimental realization of these polymorphs and therefore we limit our analysis to the hexagonal monolayer.

Within the RPA, local field effects—that originate from charge fluctuations induced by the system inhomogeneity—attenuate the second order susceptibility in both GaSe and InSe monolayers. This reduction has been observed as well for other low-dimensional systems and in bulk.Attaccalite and Grüning (2013); Grüning and Attaccalite (2014) Within the G0W0+BSE, excitonic effects tend to increase the $\chi^{(2)}(0)$ and counteract the reduction from both local-field effects and quasiparticle corrections, in agreement with what observed for bulk semiconductors. Chang et al. (2001) The SHG intensity at $G_{0}W_{0}$+BSE level is for both materials closer to the IPA one.

Figures 3 and 4 [(a) and (b) panels] show the imaginary, real part and absolute value of the SHG for GaSe and InSe monolayers within the IPA (a) and $G_{0}W_{0}$+BSE (b). They are compared with the imaginary part of the dielectric function at $\omega$ and $\omega/2$ (bottom panels) at the same level of theory. From the comparison, the main spectral features in the SHG are easily recognizable as two-photon resonances at the van Hove (IPA) and exciton energies ($G_{0}W_{0}$+BSE). We focus here on the features below the onset of the optical absorption, that is in the transparency region (imaginary part of the dielectric function $\varepsilon_{2}$), which is the spectral region of technological interest. For GaSe, comparing the IPA with the $G_{0}W_{0}$+BSE results, we first notice that the transparency region is extended by 1 eV and the SHG intensity (absolute value) redistributed to higher energies. The lowest energy feature in the SHG intensity at the IPA level is a broad shoulder at about 1.3 eV. In the $G_{0}W_{0}$+BSE instead, the lowest feature in the SHG intensity is a peak at about 1.75 eV which corresponds to a two-photon resonance with the lowest optical active excitons. Overall, the SHG absolute value is lower at the $G_{0}W_{0}$+BSE than at the IPA level. The reduction is mainly due to the quasiparticle shift (a larger band gap corresponds to a lower SHG as discussed above) which is only partially compensated by the excitonic effects. Differently from the GaSe, in the InSe the excitonic effects fully compensate for the reduction of the SHG due to the quasiparticle shift and the local field effects. The transparency region is extended by 0.4 eV and again the SHG intensity is redistributed at higher energies and enhanced at excitonic resonances. In this case, because of the smaller quasiparticle corrections—the gap opens of 1.35 eV—the SHG at the $G_{0}W_{0}$+BSE level is slightly stronger than at the IPA level at resonance. In particular, the lowest energy feature is a broad shoulder at 1 eV within the IPA and a stronger peak at 1.32 eV corresponding to the two photons resonance of the lowest optical active exciton.

IV Discussion and conclusions

The first thing to observe is that the inclusion of many-body effects changes the SHG spectrum: both at the level of the electronic structure (quasiparticle corrections within the $G_{0}W_{0}$) and at the level of the response function (inclusion of excitonic effects). Then, many-body effects may be important when considering technological applications (e.g. accurately determining the transparency region and excitonic resonances) and are essential when investigating the physics of excitons of these 2D systems.

On the other hand, the theoretically predicted intensities—both with and without the inclusion of many-body effects—differ substantially from the experimental values. For GaSe, Ref. Zhou et al., 2015 reports 0.7 nm/V at 0.77 eV, 1.7 nm/V at 0.92 eV and 2.4 nm/V at 1.02 eV whereas we find 0.096, 0.105, 0.111 nm/V respectively, that is a factor 7–20 smaller. For InSe, Ref. Zhou et al., 2018 reports 6.39 nm/V at 1.55 eV whereas we find 0.42 nm/V—that is a factor 15 smaller.

Previously, large differences in the experimental estimate of the SHG coefficients—over a range of three orders of magnitude Malard et al. (2013); Kumar et al. (2013); Yin et al. (2014)—were observed for the MoS2 monolayer. Eventually, Clark and coworkers Clark et al. (2014) explained that such differences depend on the model assumed to extract the value of the SHG from the experimental data. Indeed, experimentally the SHG is obtained from the measured second-harmonic intensity. Extracting the SHG coefficient is not straightforward and implies to assume a model for the system: Clark and coworkers Clark et al. (2014) argued that 2D materials should be modelled as a ‘sheet’, as in Refs. Malard et al., 2013; Yin et al., 2014, rather than as ‘bulk’, as in Ref. Kumar et al., 2013. Further, they provide an expression for the difference between $\chi_{\text{sheet}}^{\text{(2)}}$ and $\chi_{\text{bulk}}^{\text{(2)}}$, i.e. the SHG coefficient estimated either using the ‘sheet’ or the ‘bulk’ model:

[TABLE]

where $n_{\text{2D}}$ and $n_{\text{sub}}$ are the refractive index for the 2D material and of the substrate. $F$ is an “overall scaling factor” that depends on the transmittance and reflectance at the corresponding interfaces air-GaSe (InSe) and GaSe(InSe)-substrate as explained in the Supp. Mat. of Ref. Clark et al., 2014. For MoS2, the ratio in Eq. (2) can be as large as 900 Clark et al. (2014) and thus accounts for the differences in the experimental estimates. Malard et al. (2013); Kumar et al. (2013) It also indicates that the ‘bulk’ model leads to too large values for the SHG coefficient. This conclusion is also comforted by accurate theoretical predictions Grüning and Attaccalite (2014); Trolle et al. (2014) that reported values close to the experimental estimates that assume the ‘sheet’ model.

Coming back to the monochalcogenides, all the experimental values reported for the GaSe and InSe assumed a bulk model. To investigate whether this is the reason for the discrepancy with the theoretical estimates we use Eq. 2 with: $n_{\text{sub}}\approx 1.45$ (fused silica Malitson (1965)); $n_{\text{GaSe}}(\omega)=n_{\text{GaSe}}(2\omega)=2.8$, Zhou et al. (2015) and $n_{\text{InSe}}(\omega)$ = 2.86 and $n_{\text{InSe}}(2\omega)$ = 3.35, Zhou et al. (2018). We set $F=1$, rather than the value of approximately 20 indicated in Ref. Clark et al., 2014, since the experimental estimates already account for various gain/loss factors. We thus obtain ratios as high as $\approx 32$ for GaSe and $\approx 36$ for InSe. These large ratios indicate that, the experimental estimate would be significantly smaller and in substantial agreement with the theoretical predictions if using the ‘sheet’ model. Residual differences between the experimental and theoretical value could be due to the effects of the substrate and of the surface contributions to the SHGDeckoff-Jones et al. (2016), which are neglected in the theoretical calculations. Another factor, on the side of the theoretical calculations, could be the chosen dephasing parameter as discussed in Sec. II. The so-corrected experimental estimate for GaSe is as well closer to the SHG coefficients for the bi- and trilayer,Jie et al. (2015); Tang et al. (2016) which have been obtained by modeling the few-layers as a bulk medium and accounting for interference in the multilayer system.

In summary, using accurate first-principles real-time simulations based on many-body theory, we confirm strong SHG coefficients. In particular, for InSe we obtain a SHG coefficient of $\approx 0.5nm/V$ at 1.67 eV—close to maximum gain and laser efficiency of Ti:Sapphire. In substantial agreement with previous theoretical works, we found that the SHG coefficients are not as strong as claimed in the experimental works of Refs. Zhou et al., 2015, 2018, differing by an order of magnitude. We argue that the experimental values are overestimated due to the assumption of a ‘bulk’ rather than a more appropriate ‘sheet’ model to extract the SHG from the measurements and they need to be reviewed. We have shown indeed that a substantial agreement between theoretical and experimental values is recovered when a ‘sheet’ model is used to extract the SHG from the experimental data.

Note added at resubmission. Recently, we became aware of a work, Ref. Leisgang et al., 2018, which reports experimental measurements of SHG in InSe in good agreement with our theoretical predictions and indeed uses the ‘sheet’ model to extract the SHG coefficient.

acknowledgments

The research leading to these results has received funding from the European Union Seventh Framework Program under grant agreement no. 696656 GrapheneCore1 and no. 785219 Graphene Core2. M.P. ackowledges “Tor Vergata” University for financial support through the Mission Sustainability project 2DUTOPI. E.C. acknowledges support by the Programma per Giovani Ricercatori - 2014 “Rita Levi Montalcini”. M.P., E.C and C.A. acknowledge PRACE for computational resources on Marconi at CINECA. C.A. acknowledges A. Zappelli for the management of the computer cluster Rosa. This publication is based upon work from COST Action TUMIEE CA17126, supported by COST (European Cooperation in Science and Technology).