Resonant optical second harmonic generation in graphene-based heterostructures

TL;DR

This paper investigates how resonant optical second harmonic generation (SHG) in graphene heterostructures reveals symmetry breaking and electronic properties, with magnetic fields enabling distinction between graphene signals and substrate effects.

Contribution

It introduces a model for SHG in graphene on substrates, demonstrating how magnetic fields can distinguish substrate signals from graphene responses.

Findings

SHG signal depends linearly on inversion symmetry breaking

Magnetic fields induce Landau levels, creating multiple SHG resonances

Resonant SHG occurs at energies below substrate band-gap under magnetic field

Abstract

An optical Second-Harmonic Generation (SHG) allows to probe various structural and symmetry-related properties of materials, since it is sensitive to the inversion symmetry breaking in the system. Here, we investigate the SHG response from a single layer of graphene disposed on an insulating hexagonal Boron Nitride (hBN) and Silicon Carbide (SiC) substrates. The considered systems are described by a non-interacting tight-binding model with a mass term, which describes a non-equivalence of two sublattices of graphene when the latter is placed on a substrate. The resulting SHG signal linearly depends on the degree of the inversion symmetry breaking (value of the mass term) and reveals several resonances associated with the band gap, van Hove singularity, and band width. The difficulty in distinguishing between SHG signals coming from the considered heterostrusture and environment (insulating substrate) can be avoided applying a homogeneous magnetic field. The latter creates Landau levels in the energy spectrum and leads to multiple resonances in the SHG spectrum. Position of these resonances explicitly depends on the value of the mass term. We show that at energies below the band-gap of the substrate the SHG signal from the massive graphene becomes resonant at physically relevant values of the applied magnetic field, while the SHG response from the environment stays off-resonant.

Equations

$${\hat{H}}_{\mathbf{k}}=\left(\begin{array}{cc}{f}_{\mathbf{k}}+m& S_{\mathbf{k}}\\ S_{\mathbf{k}}^{\ast }& f_{\mathbf{k}}-m\end{array}\right).$$

$${\hat{H}}_{ij}[A]$$

$$\mathcal{S}[A]=-\mathrm{Tr}\sum _{\mathbf{k}\nu }{\hat{c}}_{\mathbf{k}\nu }^{\ast }\left(1\left(i\nu +\mu \right)-{\hat{H}}_{\mathbf{k}}[A]\right){\hat{c}}_{\mathbf{k}\nu }^{}.$$

$$j_{\alpha \omega }[A]-j_{\alpha \omega }[0]$$

$$+\frac{1}{2}\sum _{\beta \gamma ,{\omega }^{\prime }{\omega }^{\prime \prime }}{\frac{{\delta }^2{j}_{\alpha \omega }[A]}{\delta A_{\beta {\omega }^{\prime }}\delta A_{\gamma {\omega }^{\prime \prime }}}|}_{A=0}{A}_{\beta {\omega }^{\prime }}{A}_{\gamma {\omega }^{\prime \prime }}.$$

$${\frac{{\delta }^2{j}_{\alpha \omega }[A]}{\delta A_{\beta {\omega }^{\prime }}\delta A_{\gamma {\omega }^{\prime \prime }}}|}_{A=0}={\frac{{\delta }^3\mathcal{F}[A]}{\delta A_{\alpha \omega }{A}_{\beta {\omega }^{\prime }}\delta A_{\gamma {\omega }^{\prime \prime }}}|}_{A=0}=2e^3{\Pi }_{\alpha \beta \gamma }^{\omega {\omega }^{\prime }{\omega }^{\prime \prime }}.$$

$${\Pi }_{\alpha \beta \delta }^{(2)}(\omega )$$

$$+2\mathrm{Tr}\sum _{\mathbf{k},\nu }{\hat{\mathrm{v}}}_{\alpha \beta }^{(2)}\hat{G}\left(\mathbf{k},\nu +\omega \right){\hat{\mathrm{v}}}_{\gamma }^{(1)}\hat{G}\left(\mathbf{k},\nu \right),$$

$${\Pi }_{\alpha \beta \gamma }^{(3)}(\omega )$$

$${\hat{\mathrm{v}}}_{\alpha }^{(1)}={\frac{1}{e}\frac{\delta {\hat{H}}_{\mathbf{k}}[A]}{\delta A_{\alpha }}|}_{A=0},{\hat{\mathrm{v}}}_{\alpha \beta }^{(2)}={\frac{1}{e^2}\frac{{\delta }^2{\hat{H}}_{\mathbf{k}}[A]}{\delta A_{\alpha }\delta A_{\beta }}|}_{A=0}.$$

$${\bar{\Pi }}_{\alpha \beta \gamma }(\omega )=\underset{ϵ\to 0}{\lim }\left[{\Pi }_{\alpha \beta \gamma }(-i\omega +ϵ)-{\Pi }_{\alpha \beta \gamma }(0)\right].$$

$$\eta (\omega )=\bar{\Pi }(\omega )/\omega .$$

$${\hat{H}}_{\mathrm{D}}=v\left[\tau k_x{\hat{\sigma }}_x+k_y{\hat{\sigma }}_y\right]+m{\hat{\sigma }}_z.$$

$${\hat{H}}_{\text{TW}}=\lambda \left[2\tau k_x{k}_y{\hat{\sigma }}_y-(k_x^2-k_y^2){\hat{\sigma }}_x\right],$$

$$\eta (\omega )=12im\lambda v\sum _{\mathbf{k}}\frac{\tanh \left(\frac{\beta {\epsilon }_{\mathbf{k}}}{2}\right)}{{\epsilon }_{\mathbf{k}}}\left[\frac{1}{{\omega }^2-{\epsilon }_{\mathbf{k}}^2}+\frac{4}{{\omega }^2-4{\epsilon }_{\mathbf{k}}^2}\right],$$

$$S_{\mathbf{k}}$$

$$f_{\mathbf{k}}$$

$${\mathrm{v}}_{\alpha }^{(1)\sigma {\sigma }^{\prime }}(\mathbf{k})=\frac{1}{e}{\frac{\delta H_{\mathbf{k}}^{\sigma {\sigma }^{\prime }}}{\delta A_{\alpha }}|}_{A=0}={\partial }_{k_{\alpha }}{H}_{\mathbf{k}}^{\sigma {\sigma }^{\prime }}-i(r_{\alpha }^{{\sigma }^{\prime }}-r_{\alpha }^{\sigma })H_{\mathbf{k}}^{\sigma {\sigma }^{\prime }},$$

$${\mathrm{v}}_x^{(1)}(\text{k})$$

$${\mathrm{v}}_y^{(1)}(\text{k})$$

$${\mathrm{v}}_{\alpha \beta }^{(2)\sigma {\sigma }^{\prime }}(\mathbf{k})=\frac{1}{e^2}{\frac{\delta H_{\mathbf{k}}^{\sigma {\sigma }^{\prime }}}{\delta A_{\alpha }\delta A_{\beta }}|}_{A=0}=\left[{\partial }_{k_{\alpha }}{\partial }_{k_{\beta }}+i({\rho }_{\alpha }^{{\sigma }^{\prime }}-{\rho }_{\alpha }^{\sigma }){\partial }_{k_{\beta }}+i({\rho }_{\beta }^{{\sigma }^{\prime }}-{\rho }_{\beta }^{\sigma }){\partial }_{k_{\alpha }}-({\rho }_{\beta }^{{\sigma }^{\prime }}-{\rho }_{\beta }^{\sigma })({\rho }_{\alpha }^{{\sigma }^{\prime }}-{\rho }_{\alpha }^{\sigma })\right]H_{\mathbf{k}}^{\sigma {\sigma }^{\prime }}.$$

$${\mathrm{v}}_{xx}^{(2)}(\mathbf{k})$$

$${\mathrm{v}}_{yy}^{(2)}(\mathbf{k})$$

$$\hat{G}(\text{k},\nu )=\frac{1}{1(i\nu +\mu )-{\hat{H}}_{\mathbf{k}}}=\frac{1(i\nu +\mu -f_{\mathbf{k}})+\hat{\mathbf{\sigma }}\cdot {\mathbf{\xi }}_{\mathbf{k}}}{(i\nu +\mu -f_{\mathbf{k}}^{}{)}^2-{\xi }_{\mathbf{k}}^2},$$

$$\frac{1(i\nu +\mu -f_{\mathbf{k}})+\hat{\mathbf{\sigma }}\cdot {\mathbf{\xi }}_{\mathbf{k}}}{(i\nu +\mu -f_{\mathbf{k}}^{}{)}^2-{\xi }_{\mathbf{k}}^2}=\frac{\hat{A}}{i\nu +\mu -f_{\mathbf{k}}-{\xi }_{\mathbf{k}}}+\frac{\hat{B}}{i\nu +\mu -f_{\mathbf{k}}+{\xi }_{\mathbf{k}}}.$$

$$\hat{G}(\mathbf{k},\nu )=\sum _{s=\pm 1}\frac{{\hat{\Lambda }}^s(\text{k})}{i\nu +\mu -f_{\mathbf{k}}-s{\xi }_{\mathbf{k}}},$$

$${\hat{\Lambda }}^s(\text{k})=\frac{1}{2}\left(1+\frac{s}{{\xi }_{\mathbf{k}}}\hat{\mathbf{\sigma }}\cdot {\mathbf{\xi }}_{\mathbf{k}}\right).$$

$${\Omega }_{\alpha \beta \gamma }^{s{s}^{\prime }{s}^{\prime \prime }}(\text{k})=\mathrm{Tr}\left[{\hat{\mathrm{v}}}_{\alpha }^{(1)}{\hat{\Lambda }}^s(\text{k}){\hat{\mathrm{v}}}_{\beta }^{(1)}{\hat{\Lambda }}^{s^{\prime }}(\text{k}){\hat{\mathrm{v}}}_{\gamma }^{(1)}{\hat{\Lambda }}^{s^{\prime \prime }}(\text{k})\right],$$

$${\Pi }_{\alpha \beta \gamma }^{(3)}(\omega ,\mathbf{k})=\sum _{\nu ,s,s^{\prime },s^{\prime \prime }}\frac{{\Omega }_{\alpha \beta \gamma }^{s{s}^{\prime }{s}^{\prime \prime }}(\text{k})}{[i\nu +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^s][i\nu +i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime }}][i\nu -i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime \prime }}]}.$$

$$\frac{1}{\left(i\nu +i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime }}\right)\left(i\nu -i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime \prime }}\right)}=\frac{A}{i\nu +i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime }}}+\frac{B}{i\nu -i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime \prime }}},$$

$$A=-B=-\frac{1}{2i\omega -\left({\xi }_{\mathbf{k}}^{s^{\prime }}-{\xi }_{\mathbf{k}}^{s^{\prime \prime }}\right)}.$$

$${\Pi }_{\alpha \beta \gamma }^{(3)}(\mathbf{k},\omega )=\sum _{s,s^{\prime },s^{\prime \prime }}\frac{-{\Omega }_{\alpha \beta \gamma }^{s{s}^{\prime }{s}^{\prime \prime }}(\text{k})}{2i\omega -\left({\xi }_{\mathbf{k}}^{s^{\prime }}-{\xi }_{\mathbf{k}}^{s^{\prime \prime }}\right)}\sum _{\nu }\frac{1}{[i\nu +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^s]}\left[\frac{1}{i\nu +i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime }}}-\frac{1}{i\nu -i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime \prime }}}\right].$$

$$\sum _{\nu }\frac{1}{[i\nu +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^s][i\nu +i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime }}]}=\frac{1}{i\omega +\left({\xi }_{\mathbf{k}}^s-{\xi }_{\mathbf{k}}^{s^{\prime }}\right)}\sum _{\nu }\left[\frac{1}{i\nu +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^s}-\frac{1}{i\nu +i\omega +\mu -f_{\mathbf{k}}^{}-{\xi }_{\mathbf{k}}^{s^{\prime }}}\right].$$

Full text

I Introduction

The Second-Harmonic Generation (SHG) has become a very important tool to investigate different properties of materials. The sensitivity of the SHG to inversion symmetry and number of layers is an important aspect to perform experiments or to realize devices based on quasi-two-dimensional (2D) heterostructures Rosa et al. (2018). The fact that the SHG is forbidden in materials where the inversion symmetry is preserved Boyd (2008) can also be exploited for the investigation of layered systems composed from different materials. Recently, such heterostructures attracted a lot of attention from the physical community due to their unusual electronic properties Geim and Grigorieva (2013). The state-of-the art method to study structural properties of these systems is the scanning tunneling microscopy (STM). However, the direct STM measurements can probe only the surface states and are not sensitive to structural changes in multilayered heterostructures. On the other hand a much simpler experiment on the SHG can indirectly capture the differences between various combinations of layers even if it occurs not at the surface.

Unfortunately, an application of this technique to a simplest and most extensively studied 2D material, i.e. a monolayer graphene, turns out to be inefficient. Indeed, a pristine graphene exhibits inversion symmetry, which prevents any SHG. The SHG signal in graphene can be observed either by inducing an asymmetry between two sublattices of graphene placing it on top of a band insulator, or considering the fact that a photon momentum $\mathbf{q}$ of the applied light already works as a source of asymmetry. It has been shown that the response caused by the photon momentum is weak, since it is proportional to the momentum itself Glazov (2011). The problem of the SHG in the case of the hexagonal lattice with the broken inversion symmetry has also been considered previously Säynätjoki et al. (2017); Margulis et al. (2013). It is worth mentioning another physical effect, namely the valley polarization, that allows the SHG in graphene Wehling et al. (2015). Although this mechanism could be very useful in the context of valleytronics, addressing the valley polarization experimentally is still a matter of a research.

In this work we investigate the SHG from the graphene disposed on the insulating hexagonal substrates with different band-gaps as a particular example of quasi-2D heterostructures mentioned above. For this aim we perform calculations using the diagrammatic technique based on the full dispersion of the non-interacting tight-binding model with the mass term. The reason for a yet another theoretical study of the SHG in graphene with the broken inversion symmetry is the fact that previous studies on this subject report features that can hardly be explained within physical intuition. For instance, this concerns a stronger SHG response for smaller mass term Margulis et al. (2013), which is very surprising, since the mass is a consequence of the inversion symmetry breaking. Furthermore, an overwhelming majority of theoretical studies of the SHG in graphene-based heterostructures are focused on frequencies of the incident light around the band-gap, which is far from the experimentally accessible regime, where the energy of the incoming photons is usually around 1.5 eV for red light sources. For this reason, we obtain the full SHG spectrum that is needed for a description of the actual experimental data. This allows us to reveal additional resonances in the optical spectrum that correspond to the van Hove singularity and band-width, which cannot be captured by a simplified Dirac model. The comparison between the SHG response obtained for the full tight-binding spectrum and the one for the approximated Dirac picture allows us to define the limits of applicability of this approximation.

Another experimentally relevant problem that stays undiscussed in all previous works is the difficulty to distinguish between SHG signals from the graphene flake placed on the insulating substrate and the rest of the insulating sample. Indeed, we find that, contrary to the result of the Ref. Margulis et al., 2013, the SHG response is proportional to the band-gap, which in the case of graphene disposed on the substrate is small compared to the band-gap of the clean substrate. Therefore, the SHG signal from graphene can hardly be seen on top of the large SHG signal from the band insulator. We show that this problem can be resolved including the homogeneous magnetic field in the system. The presence of the magnetic field results in the formation of Landau levels in the energy spectrum. Since these levels are sharp, we expect intense resonances associated with transitions between Landau levels in the SHG spectrum. Therefore, the presence of the magnetic field introduces a natural amplification of the SHG that can be tuned adjusting the value of the external magnetic field at a fixed laser frequency. The small value of the mass term of the graphene-based hetrostructure allows to find resonances on Landau levels already at energies below the band-gap of the insulating substrate, while the SHG signal of the environment stays off-resonant. An experimental evidence of the applicability of this technique can be found in Ref. Pavlov et al., 2005.

II SHG response from massive graphene

Here, we study optical second harmonic generation from graphene-based heterostructures using the following tight-binding model that describes a behavior of non-interacting electrons on a hexagonal lattice. The corresponding Hamiltonian matrix written in the sublattice space reads

[TABLE]

Here, the off-diagonal term $S_{\mathbf{k}}$ is a Fourier transform of the nearest-neighbour hopping process $t$, and the diagonal one $f_{\mathbf{k}}$ describes the next-nearest-neighbor hoping of electrons $t^{\prime}$ (see Appendix A). Here we also introduce the mass term $m$ that explicitly breaks the inversion symmetry of the system. The mass term describes the sublattice imbalance, which illustrates the situation when graphene is placed on top of a band insulator, such as the hexagonal Boron Nitride (hBN) Woods et al. (2014); Wallbank et al. (2013); Bokdam et al. (2014); Slotman et al. (2015), or Silicon Carbide (SiC) Zhou et al. (2007), as discussed above. This model is found to be relevant for other materials, such as the monolayer MoS2 Ochoa and Roldán (2013) and germanene under the effect of an electric field Acun et al. (2015). Note that for the case of graphene on hBN the appearance of the gap corresponds to the commensurate phase at small enough misorientation angle between their lattices, whereas for the incommensurate phase the average gap is supposed to be zero Woods et al. (2014); Titov and Katsnelson (2014). It is the average gap that matters in SHG experiments, where the typical laser spot size is much larger than the interatomic distance. In the following when discussing graphene placed on the insulating substrate we will keep in mind only the commensurate case. The dispersion relation $E_{\bf k}^{\pm}=f^{\phantom{|}}_{\bf k}\pm\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{|S^{\phantom{|}}_{\bf k}|^{2}+m^{2}}$ for the case of graphene on top of hBN with $t=-2.8$ eV, $t^{\prime}=-0.1t$ Kretinin et al. (2013), and $m=30$ meV is shown in Fig. 1.

Following the procedure adopted from Ref. Wehling et al., 2015, the effect of the applied probe light is accounted via the Peierls substitution introducing a vector potential ${\bf A}$ that represents an external radiation

[TABLE]

Here, $e$ is the modulus of the electronic charge, $c$ is the speed of light, and $\hat{H}_{ij}$ is the lattice $\{i,j\}$ space representation of the Hamiltonian matrix (1).

The SHG response function we aim to obtain in the current work can be derived using the Feynman diagram technique. For this reason, we stick to the path integral formalism with the corresponding action for our problem

[TABLE]

Here, $\hat{c}^{(*)}_{\mathbf{k}}$ is a two-dimensional spinor of Grassman variables that describe creation (annihilation) of an electron, and $\mathbb{1}$ is the $2\times 2$ identity matrix in the sublattice space; the trace is taken over the same space, $\beta$ is the inverse temperature and $\sum_{\mathbf{k},\nu}$ stands for the infinite summation over the fermionic Matsubara frequencies $\nu_{n}=(2n+1)\pi/\beta$ and momentum integration over the Brillouin Zone ($BZ$). The chemical potential $\mu=-3t^{\prime}$ corresponds to the half-filled case (neutrality point).

The electric current density can be defined as the response of the system on the applied vector potential $j_{\alpha\omega}[A]=\delta{\cal F}[A]/\delta A_{\alpha\omega},$ where ${\cal F}[A]=\ln{\cal Z}[A]$ and ${\cal Z}[A]=\int D[c^{*},c]\exp\left\{-\beta{\cal S}[A]\right\}$ is the is the generating functional, that is, the partition function of the problem written in terms of the action (3). Expanding the electric current up to the second order with respect to the vector potential, one gets the usual relation

[TABLE]

The coefficient in front of the liner term is the usual optical conductivity of the system. The coefficient in front of the square of the vector potential describes the second harmonic generation. The latter can also be expressed via the three-particle correlation function $\Pi_{\alpha\beta\gamma}^{\omega\omega^{\prime}\omega^{\prime\prime}}$ as

[TABLE]

Diagrammatic expressions for this correlation function are shown in Fig. 2. The diagram “$b$” will be ignored in the following, because it represents a constant energy shift. Explicit expressions for the triangular $\Pi^{(3)}(\omega)$ (“$a$”) and nonlinear bubble $\Pi^{(2)}(\omega)$ (“$c$” and “$d$”) diagrams are following

[TABLE]

where $\hat{G}({\bf k},\nu)=\left[\mathbb{1}(i\nu+\mu)-\hat{H}_{\bf k}\right]^{-1}$ is the Green’s function of our problem, and velocity operators can be defined in the same way as in Ref. Wehling et al., 2015

[TABLE]

As can be seen from Eq. 6, expressions for diagrams “$c$” and “$d$” are connected by the following simple relation $\Pi_{\alpha\beta\gamma}^{(2)\,c}(\omega)=\Pi_{\alpha\beta\gamma}^{(2)\,d}(2\omega)$. The coefficient $2$ is not included in the definition of $\Pi_{\alpha\beta\gamma}^{(2)\,d}(2\omega)$. Then, the total result for the nonlinear bubble can be written as $\Pi_{\alpha\beta\gamma}^{(2)}(\omega)=\Pi_{\alpha\beta\gamma}^{(2)\,d}(2\omega)+2\Pi_{\alpha\beta\gamma}^{(2)\,d}(\omega)$, which explicitly connects the behavior of the SHG spectrum at double- and single frequencies of the applied light. It is also worth mentioning that the contribution of the $\Pi_{\alpha\beta\gamma}^{(2)\,d}(\omega)$ diagram is missing in Wehling et al. (2015).

The real frequency dependence of correlation functions can be obtained performing an analytic continuation

[TABLE]

In practical calculations $\epsilon$ is taken to be small, but finite.

In the following we take into account that the experimentally measurable quantity for the SHG is the conversion efficiency. It can be defined in the same way as in Ref. Boyd, 2008 and is proportional to the ratio

[TABLE]

The explicit evaluation of the introduced diagrams is shown in the Appendix B. We find that the contribution of the triangular diagram $\Pi^{(3)}(\omega)$ is zero even when a non-zero next-nearest-neighbor hopping $t^{\prime}$ and chemical potential away from the half-filling are considered. The reason is that the integral over momentum $\mathbf{k}$ in equation (6) for $\Pi^{(3)}(\omega)$ averages to zero in the whole Brillouin zone. This is essentially due to the fact that two valleys $K$ and $K^{\prime}$ of graphene contribute to the integral with opposite signs and hence compensate each other. This result is a generalization of the situation considered in Säynätjoki et al. (2017) for the case of a low-energy Hamiltonian for MoS2 material, where the contribution of the triangular diagram is canceled by symmetry with respect to the inversion of $k_{y}$. A non-zero result for the triangular diagram can be obtained introducing a valley polarization that generates an imbalance between the two valleys, as discussed in Ref. Wehling et al., 2015.

Contrary to the triangular diagram, the contribution from the nonlinear bubble is nonzero, and the $\Pi^{(2)}_{\alpha\beta\gamma}(\omega)$ tensor reveals the reduced symmetry $C_{3}$ instead of $C_{6}$ with respect to rotation (see Appendix B). Thus, we find that the contribution $\Pi^{(2)}_{xxx}(\omega)=0$, whereas the result for $\Pi^{(2)}_{yyy}(\omega)$ is nonzero. It can be shown that the only non-zero components of the tensor are $\Pi^{(2)}(\omega)=\Pi_{yyy}^{(2)}(\omega)=-\Pi_{xxy}^{(2)}(\omega)=-\Pi_{yxx}^{(2)}(\omega)=-\Pi_{xyx}^{(2)}(\omega)$.

Remarkably, the account for the next-nearest-neighbor hopping process $t^{\prime}$ also does not change the result for the nonlinear bubble diagram. This can be explained looking at the dispersion relation in Fig. 1. The inclusion of $t^{\prime}$ equally shifts the upper and the lower band at given momentum $\mathbf{k}$, which does not change the energy difference between them. Since we consider only direct excitations at zero momentum, the SHG spectrum depends on energy difference between two bands and hence does not change with the inclusion of the next-nearest-neighbor hopping.

Fig. 3 shows the absolute value of the conversion efficiency $\eta_{yyy}(\omega)$ for different values of the mass term (half of the band-gap). Here, the results for the hexagonal Boron Nitride (hBN, $t=-2.4$ eV, $m=2.78$ eV), graphene on hBN (Gr/hBN, $t=-2.8$ eV, $m=0.03$ eV), and graphene on a SiC substrate (Gr/SiC, $t=-2.8$ eV, $m=0.13$ eV) Bokdam et al. (2014); Margulis et al. (2013) are compared. The data on the right panel for the Gr/SiC is multiplied by a factor of 5, and the one for the Gr/hBN is multiplied by $5\times(m_{\rm Gr/SiC}\div m_{\rm Gr/hBN})$ for clarity. The real and imaginary parts of the conversion efficiency $\eta_{yyy}(\omega)$ are shown in Appendix B.

The SHG is a virtual process that is allowed even if the frequency of the applied light is smaller than the half of the band-gap. The nonzero increasing value of the conversion efficiency $\eta(\omega)$ at small $\omega<m$ frequencies confirms this statement. This consideration appears to be in agreement with the result of Säynätjoki et al. (2017), and also reproduces the trend observed in Brun and Pedersen (2015). Increasing the frequency of the applied light, the SHG spectrum reveals the first resonance at energies when excited electrons reach the band-gap. This resonance appears as a pair of peaks at frequencies $\omega=m$ and $\omega=2m$ in agreement with the frequency dependence of the nonlinear bubble diagram $\Pi^{(2)}(\omega)$ discussed above. These peaks are labeled as “1” in Fig. 3. The double resonance on the band-gap was reported previously in Margulis et al. (2013), but is missing, for example, in Brun and Pedersen (2015).

The use of the full tight-binding dispersion allows us to capture additional resonances in the SHG spectrum labelled as “2” and “3” in Fig. 3. The peak “2” corresponds to the van Hove singularity and appears in the nonlinear optical spectrum at the frequency $\omega=2.8$ eV (Gr/hBN and Gr/SiC) and $\omega=3.67$ eV (hBN) with its replica at $2\omega$. The resonant peak “3” at the highest energy comes from the bandwidth of the system. Note that the use of the full tight-binding dispersion is crucial for a description of these additional optical resonances depicted in Fig. 1 via red arrows, since the low energy expansion (Dirac picture) does not provide the corresponding features of the energy spectrum.

In addition, we observe that the off-resonant SHG response function linearly depends on the value of the mass term. Indeed, the conversion efficiency $\eta(\omega)$ for the Gr/SiC shown in Fig. 3 (right) is almost indistinguishable from the one of the Gr/hBN multiplied by the factor $(m_{\rm Gr/SiC}/m_{\rm Gr/hBN})$. The fact that a smaller mass term leads to a smaller value of the conversion efficiency can also be seen comparing the off-resonant behavior of $\eta(\omega)$ for the hBN with the one for graphene-based heterostructures. For instance, the SHG response from the hBN for the applied red light ($\omega\simeq 1.5$ eV) is five times larger than the one from the Gr/hBN. However, the direct comparison of these signals is complicated by a predominant resonant behavior of the conversion efficiency of hBN at large frequencies. From the physical point of view, this result can be explained as follows. The mass term is exactly the factor that breaks the inversion symmetry and hence is responsible for the SHG. Since the latter is identically zero in systems with the unbroken inversion symmetry, the larger SHG signal is expected when the symmetry breaking is more severe.

It is worth mentioning that our result is in a disagreement with the one of the Ref. Margulis et al., 2013, where the authors report a larger SHG signal for the material with a smaller mass term.

III Influence of magnetic field on the SHG

Now, let us discuss the effect of a homogeneous magnetic field on the SHG. The inclusion of the magnetic field in a general tight-binding model can be done only numerically, which is less intuitive and is much more complicated than an analytical solution of the problem. However, this issue can be resolved in the framework of the Dirac model. The latter can be obtained performing the low-energy expansion of the Hamiltonian matrix (1) in the vicinity of $K$ and $K^{\prime}$ points of the hexagonal Brillouin Zone of graphene Katsnelson (2012)

[TABLE]

Here, $\tau=\pm 1$ is the valley index, and $v=3at/2$ is the electron speed at conical points $K$ and $K^{\prime}$.

The SHG in the case of Dirac electrons is forbidden by Furry’s theorem Peskin and Schroeder (1995). This is represented by the fact that the triangular diagram $\Pi^{(3)}(\omega)$ is identically zero. The nonlinear bubble diagram is absent in the Dirac approximation, since the corresponding low-energy Hamiltonian does not contain any second order term in momentum $\mathbf{k}$ that is responsible for the existence of the the non-linear vertex ${\rm v}^{(2)}$. Moreover, the SHG response calculated on the basis of the low-energy Hamiltonian does not reproduce all features of the SHG spectrum shown in Fig. 3. Therefore, we need to go beyond the Dirac approximation in order to obtain an experimentally relevant result for the SHG. Here, we can benefit from the fact that the frequency of the red light, which is commonly used in SHG experiments, is smaller than the van Hove singularity. For this reason, the inclusion of already the first order correction in momentum $\mathbf{k}$ to the low energy Hamiltonian (10) will be sufficient to describe the SHG in graphene-based heterostructures, although the resulting model will not hold for energies around the van Hove singularities. We indicate the correction to the Diarc Hamiltonian as $H_{\text{TW}}$, which is the so called trigonal warping term

[TABLE]

where $\lambda=3a^{2}t/8$ is the trigonal warping parameter Katsnelson (2012). The crucial role of the trigonal warping for the SHG in graphene at low energies was pointed out in previous works Golub and Tarasenko (2014); Margulis et al. (2013), although the result for the SHG is obtained there for the regime of very low frequencies, which is not accessible experimentally. The Dirac approximation with the trigonal warping term was also considered in Ref. Säynätjoki et al. (2017) for the case of the SHG in MoS2. From the diagrammatic point of view, the role of the trigonal warping is to introduce a non-linear vertex ${\rm v}^{(2)}$ in the theory, which is responsible for the existence of the non-linear bubble diagram $\Pi^{(2)}(\omega)$. In the following, we make an additional approximation expanding the $\Pi^{(2)}(\omega)$ up to the first order in the trigonal warping parameter $\lambda$. Then, the contribution of the trigonal warping remains only in the vertex function ${\rm v}^{2}$, while the Green’s function stays the same as for massive Dirac electrons. This will allow to account for the effect of Landau levels in the Green’s function analytically without any approximations.

In order to estimate limits of applicability of the derived approximation, let us consider the correction to the the dispersion relation due to the trigonal warping. At small values of the mass $m$, the contribution to the energy from Dirac dispersion is approximatelly equal to $\frac{3t}{2}k$ and the contribution from the trigonal warping term is $\frac{3t}{8}k^{2}$. The latter can be considered as a small correction for $k$ up to $0.4$, which corresponds to the energy of about $1.7$ eV. Fig. 4 shows that this estimation is rather conservative, and the SHG response function of the approximate model is in a good agreement for the one obtained using the full tight-binding spectrum up to energies of about $2$ eV. As expected, the breakdown of the approximation is associated with the presence of the resonance on the van Hove singularity in the SHG spectrum, which can not be reproduced without the full tight-binding dispersion.

After all, the conversion efficiency $\eta(\omega)$ at zero magnetic field can be recast in a very simple form (see Appendix C)

[TABLE]

where $\varepsilon_{\mathbf{k}}=\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{v^{2}\mathbf{k}^{2}+m^{2}}$ is the massive Dirac dispersion of electrons in graphene. Remarkably, Eq. 12 shows that the off-resonant value of $\eta(\omega)$ linearly dependents on the mass term $m$, as expected from above discussions.

The homogeneous magnetic field directed perpendicular to the surface of the sample can be introduced in the system via the Peierls substitution ${\mathbf{k}}\to\mathbf{k}+e{\bf A}_{B}$ with the following vector potential ${\bf A}_{B}=\frac{B}{2}(-y,x,0)$. With this substitution, the energy spectrum of the problem changes dramatically from $\varepsilon_{\mathbf{k}}$ to a discrete set of Landau levels described by the following expression $\varepsilon_{n}=\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{m^{2}+2|eBv^{2}n|}$ with $n\in\mathbb{Z}$, (see Katsnelson (2012) and Gusynin et al. (1995) for the case of zero mass term).

In the presence of the magnetic field, the translational symmetry of the initial problem is explicitly broken. However, the symmetry with respect to inversion in $\mathbf{k}$-space in the Dirac model with the trigonal warping is still preserved. This ensures that the contribution from the triangular diagram $\Pi^{(3)}(\omega)$ to the SHG response remains zero. The explicit calculation of the nonlinear bubble diagram $\Pi^{(2)}(\omega)$ with the above approximations is shown Appendix D. The result for the corresponding conversion efficiency for the case of Gr/SiC in the presence of the magnetic field is shown in Fig. 5. Here, we clearly see the multiple-peak structure of the SHG response function due to excitations between Landau levels. A similar picture has been observed experimentally in Pavlov et al. (2005) for the SHG in GaAs material. It can be shown that a selection rule for the allowed transitions between Landau levels is $\Delta n=\pm 1$. A more precise analysis of Fig. 5 allows to distinguish two types of peaks with different intensities that correspond to the contribution of different diagrams “$c$” and “$d$” shown in Fig. 5 to the SHG.

As we observe above, the SHG signal from the Gr/hBN in the absence of the magnetic field is much smaller than the signal coming from the clean hBN at experimentally accessible frequencies. An application of the homogeneous magnetic field results in the formation of Landau level, which in the case of small value of the effective mass of the graphene-based heterostructures appear already at energies below the band-gap of the insulating substrate. The resonance on the Landau levels drastically enhances the SHG signal from the graphene placed on top of the hBN, while the SHG responce from the hBN remains off-resonant at energies below the bang-gap. Therefore, the inclusion of the magnetic field simplifies the detection of the SHG signal from the graphene flake disposed on top of the insulating substrate.

IV Conclusion

In this work we present a consistent calculation of the SHG response from the graphene-based heterostructures. In our calculations we first start with the full tight-binding dispersion and obtain the nonlinear optical spectrum for experimentally relevant frequencies of the applied light. We find that the conversion efficiency has three pairs of resonances that correspond to optical excitations between the band-gap, van Hove singularity and band width. We also observe that the off-resonant behavior of the SHG response function linearly depends on the mass term, contrary to what has been reported in previous studies. The problem of distinguishing the small signal from the graphene on top of the insulating substrate is proposed to resolve here by the inclusion of the magnetic field. The presence of the latter in the system leads to a formation of Landau levels in the energy spectrum at energies below the band-gap of the insulating substrate. This allows to obtain resonant SHG excitations for the considered heterostructure keeping the SHG signal from hBN off-resonant. The magnetic field is included in the theory analytically in the framework of the simplified Dirac model with the trigonal warping. The limit of applicability of this approximation is carefully discussed. For instance, we find that the derived approximation is valid at energies below 2 eV, which is sufficient for a description of the SHG experiment with experimentally relevant frequency of the applied light.

Acknowledgements.

We would like to thank Clement Dutreix, Sergey Semin, Alexey Kimel, and Franca Manghi for fruitful discussions. M.I.K. acknowledges the supported by NWO via Spinoza Prize. The work of E.A.S. was supported by the Russian Science Foundation, Grant 17-72-20041.

Appendix A Derivation of vertex functions

Explicit relation for the Fourier transform of the nearest-neighbor $t$ and next-nearest-neighbor $t^{\prime}$ hopping processes in the case of a hexagonal lattice is following

[TABLE]

The expression for velocities can be derived following the Ref. Wehling et al., 2015

[TABLE]

where ${\bf A}$ is the vector potential of the applied light introduced by the Peierls substitution, and $\sigma,\sigma^{\prime}$ indicate pseudospin degrees of freedom related to a sublattice space. $r$ indicates the atomic position within the unit cell. The resulting expression for velocities becomes

[TABLE]

where ${\rm v}_{x}^{(1)\,AA}(\mathbf{k})={\rm v}_{x}^{(1)\,BB}(\mathbf{k})=-2\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{3}t^{\prime}a\left(\sin(\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{3}k_{x}a)+\sin(\frac{\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{3}}{2}k_{x}a)\cos(\frac{3}{2}k_{y}a)\right)$. In the same way, recalling the equation for the two-photon velocity

[TABLE]

This leads to the following expressions for the components of the second-order velocities

[TABLE]

where ${{\rm v}_{xx}^{(2)\,AA}}(\mathbf{k})={{\rm v}_{xx}^{(2)\,BB}}(\mathbf{k})=-3t^{\prime}a^{2}\left(\cos\left(\frac{\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{3}k_{x}a}{2}\right)\cos\left(\frac{3k_{y}a}{2}\right)+2\cos\left(\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{3}k_{x}a\right)\right)$, which is the expression for the two-electrons-two-photons vertex.

Appendix B Evaluation of diagrams for the SHG

An important step in our derivation is to take explicitly the summation over internal frequencies in the loop of equation 6. This simplifies considerably the expressions to be computed and gives physical insight into the problem, allowing to show explicitly the connection between formulas and physical transitions. The Green’s function $G({\bf k},\nu)$ of the initial problem (1) is

[TABLE]

where $\mathbb{1}$ is the identity matrix and $\boldsymbol{\hat{\sigma}}=(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})$ is the vector of Pauli Matrices in the sublattice space. We also define $\boldsymbol{\xi}_{\bf k}=(\mathrm{Re}\,S_{\mathbf{k}},\mathrm{Im}\,S_{\mathbf{k}},m)$ and $\xi_{\rm k}=\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{|S^{\phantom{|}}_{\bf k}|^{2}+m^{2}}$.

This representation of the Green’s function is not convenient to look for a compact expression for the diagram (6), so we define the so-called spectral representation of the Green’s function (see e.g.Le Bellac (1996)). It’s an expansion of the denominator in the above expression in simple fractions

[TABLE]

Solving for the matrix coefficients $\hat{A}$ and $\hat{B}$, we find the following expression for the propagator

[TABLE]

where we have introduced projectors over positive and negative energy states

[TABLE]

We now define the third-order factor for the diagram as

[TABLE]

Most importantly, this object does not depend on the frequency $\nu$. Omitting the integration over momenta in (6), one gets

[TABLE]

Now we concentrate on the following part of the denominator and we notice that it can be conveniently manipulated exploiting partial fractions

[TABLE]

where

[TABLE]

The useful fact is that the previously complicated evaluation is reduced to the evaluation of two polarization bubbles. The expression for this diagram at finite chemical potential can be found, for example, in  Le Bellac (1996) (p. 157). In our case one gets

[TABLE]

The first term in the second line of Eq. 30 reads

[TABLE]

Now we see that this expression can be summed introducing a convergence factor $e^{i\nu_{n}\eta}$ with $\eta\rightarrow 0$ in every term (see e.g. Fetter and Walecka (2012), p. 272). Then, the evaluation of the Matsubara sum becomes simply the evaluation of the function $\frac{1}{2}-n_{F}(\xi)$ in correspondence of the poles of the function involved, where $n_{\rm F}(\xi-\mu)=\left(e^{\beta(\xi-\mu)}+1\right)^{-1}$ is the Fermi distribution function. In this case there are two poles: $i\nu=-\mu+f^{\phantom{s}}_{\bf k}+\xi^{s}_{\bf k}$ for the first fraction in the square brackets and $i\nu+\omega=-\mu+f^{\phantom{s}}_{{\bf k}}+\xi^{s^{\prime}}_{{\bf k}}$ for the second one, keeping in mind that $i\omega$ is a bosonic Matsubara’s frequency and the exponential of bosonic frequencies gives just a factor 1. We can then derive that the sum over frequencies as

[TABLE]

The second term in the second line of Eq. 30 is given by a similar expression with the replacement $s^{\prime}\to s^{\prime\prime}$. The full expression for the diagram can be obtained putting these two terms together and reads

[TABLE]

Now we take into account that in our case the $s$ indexes can assume only values $\pm 1$. The result becomes different from 0 only if $s\neq s^{\prime}$ or/and $s\neq s^{\prime\prime}$, otherwise one of the numerators in brackets in Eq. 37 vanishes. Therefore, one gets three contributions

[TABLE]

Rearranging terms in (34) and (35), and changing $s$ to $-s$ in the second one, we get the expression

[TABLE]

The next step is to calculate the contribution coming from the non-linear bubble. The evaluation of the Matsubara summation can be done in the same way as we did for the triangular diagram. The integrand of the non-linear bubble $\Pi^{(2)}_{\alpha\beta\gamma}(\omega)$ becomes:

[TABLE]

where

[TABLE]

We notice that this expression is identical to Eq. 32 except for $-i\omega$ in the first factor. Adapting Eq. 32, we can then easily derive the expression for the bubble to be

[TABLE]

As in the previous case, the only possible contribution is given by $s\neq s^{\prime}$, so we can write

[TABLE]

Results for the Real and Imaginary parts of the conversion efficiency $\eta_{yyy}(\omega)\sim\Pi^{(2)}_{yyy}(\omega)/\omega$ for real frequencies is shown in Fig. 6.

Appendix C SHG in Dirac model with trigonal warping

At low energies and the initial Hamiltonian matrix $\hat{H}_{\mathbf{k}}$ can be expanded around $K$ and $K^{\prime}$ points with respect to momentum small momentum $\mathbf{k}$. Then, we get a Dirac approximation with the trigonal warping for the initial problem

[TABLE]

where we have introduced the electron speed at the conical points $v=\frac{3at}{2}$, valley index $\tau=\pm 1$, and the trigonal warping parameter $\lambda=3a^{2}t/8$. Velocities in this case are defined in the low energy limit (continuous limit) as

[TABLE]

Here, $\hat{\rm v}^{(1)}$, $\hat{\rm v}^{(2)}$ and $\hat{H}$ are $2\times 2$ matrices in the sublattice space.

We use the simplest approximation, that amounts to include the contribution of the trigonal warping only in the non-linear velocity ${\rm v}^{(2)}_{\alpha\beta}$. This is justified by the fact that in the dispersion relation the trigonal warping parameter, which is already smaller than the electronic speed, appears multiplied by the squared momentum. So, the corresponding correction to the position of the poles of the Green’s function is assumed to be small. On the other hand, it is the first non-zero contribution to the non-linear velocity, so it cannot be excluded from the consideration. This amounts to have a bubble diagram with a usual velocity ${\rm v}^{(1)}$ that does not depend on $\lambda$, and a non-linear vertex that is linear in $\lambda$. This is the lowest order in $\lambda$ that describes SHG. In this case the velocities can be obtained differentiating the Hamiltonian with respect to momentum as discussed above. So, they are simply proportional to Pauli Matrices as follows.

[TABLE]

We take a further approximation and consider only first order terms in trigonal warping in the diagram $\Pi^{(2)}(\omega)$. Since the $\hat{{\rm v}}^{(2)}_{yy}(\textbf{k})$ is already proportional to the trigonal warping parameter, we neglect the contribution of the trigonal warping in the Green’s function. As we discuss in the main text (see Fig. 4), the approximated result quantitatively agrees with the one of the tight-binding model. The bubble diagram $\Pi^{(2)}$ in this case can be expressed in particularly simple way. Starting from Eq. (41), we can considerably simplify this expression in our Dirac approximation. Starting from

[TABLE]

where we consider that there are three couples of $K$ and $K^{\prime}$ points in the Brillouin zone, $\tau$ is the valley index and $k_{c}$ is the maximum value of the momentum for which the Dirac approximation is valid. Now we notice that every component where $\tau$ enters only once, is zero when summed over valley index. The only non-zero contributions are therefore those involving ${\rm v}^{(1)}_{x}$ and ${\rm v}^{(2)}_{xy/yx}$ or ${\rm v}^{(1)}_{y}$ and ${\rm v}^{(2)}_{xx/yy}$, as discussed in the main text with symmetry considerations. It is straightforward to show that this considerations are all equal to each other up to minus sign. We can then evaluate the following quantity

[TABLE]

Let us also recall that $\,\,n_{\rm F}(s\varepsilon_{\mathbf{k}}-\mu)-n_{\rm F}(-s\varepsilon_{\mathbf{k}}-\mu)=s\tanh\left(\frac{\beta\varepsilon_{\mathbf{k}}}{2}\right)$. Now we can plug these expressions back to the general expression for $\Pi^{(2)}(\omega)$ and we get

[TABLE]

where we integrated over the angular coordinate since there is no angular dependence and thus the last integral is taken over the modulus $k$ of the momentum only.

Appendix D SHG in the presence of the magnetic field

In order to consider the effect of the magnetic field, we add the a vector potential describing the incident light using the Peierls substitution as

[TABLE]

where ${\bf A}_{B}=\frac{B}{2}(-x,y,0)$ is the vector potential describing the external constant magnetic field and ${\bf A}^{\text{rad}}$ describes the radiation field of the incident light. We started the result given in Ref. Gusynin et al. (1995) for the Green’s function of a Dirac particle in a magnetic field. In that paper it was demonstrated that the Green’s function of the problem can be written as

[TABLE]

where $\Phi({\bf r},{\bf r}^{\prime})=\int\limits_{\bf r}^{\bf r^{\prime}}{\bf A}_{B}({\bf z})\cdot d{\bf z}$, $\Phi_{0}$ is the magnetic flux quantum and $\tilde{G}$ is only a function of the difference between $\bf r$ and $\bf r^{\prime}$. If we consider the bubble diagram between two points in real space $\Pi^{(2)}({\bf r},{\bf r}^{\prime},\omega)\sim G({\bf r},{\bf r}^{\prime},\omega)G({\bf r^{\prime}},{\bf r},-\omega)$, we see that the two phases acquired along the paths are equal but opposite in sign, therefore they cancel each other.

[TABLE]

where the velocities coincide with those calculated in the Dirac case without magnetic field. This means that the resulting bubble depends just on the translationally invariant part of the Green’s function. This allows us to define a momentum dependent Green’s function $G({\bf k},\omega)=\sum\limits_{\bf r}\tilde{G}({\bf r},\omega)e^{i\mathbf{k}\cdot{\bf r}}$ and to calculate the response using equation (5). The expression for $G({\bf k},\omega)$ was derived in Ref. Gusynin et al. (1995) and was found to be

[TABLE]

where:

[TABLE]

$L_{n}^{\alpha}(x)$ are generalized Laguerre Polynomials and $\varepsilon_{n}=\mkern-6.0mu\mathop{}\sqrt[\mkern 8.0mu]{m^{2}+2\hbar^{2}\Omega_{c}^{2}|n|}$ are the discrete Landau levels of the system with cyclotron frequency $\Omega_{c}$. The order of magnitude for the cyclotron frequency in Graphene can be estimated replacing the value of the electron velocity $v=\frac{c}{300}$:

[TABLE]

We can now see that the denominator has two poles for each frequency. We can rewrite everything in simple fractions in order to obtain a summation over functions with a single pole. To do so we have to identify terms that depend on $\omega$ in the numerator. It is then useful to split the numerator in two parts. Explicitly this reads $\hat{D}_{n}=\hat{\Sigma}_{0}(n,\mathbf{k})+i\omega\;\hat{\Sigma}_{1}(n,\mathbf{k})$, where we have collected the quantities

[TABLE]

In this case the electron Green’s Function becomes

[TABLE]

The projectors $\Lambda_{n}$ are obtained solving the equations to reduce Eq. (52) in simple fractions and are

[TABLE]

Now only the quantity at the denominator of the Green’s function depends on the Matsubara frequencies and we can repeat exactly the same steps used to derive Eq. as in the case without magnetic field in order to calculate the bubble diagram. In the Dirac approximation, the triangular diagram $\Pi^{(3)}=0$. The bubble diagram with magnetic field then reads

[TABLE]

where

[TABLE]

and we used Eq. 32 to go from second to third line. The use of this expression reduces the computational effort necessary to compute the polarization that appears in other methods, because the only factor that depends on the momentum is the matrix element. The three-particle correlation function can be written simply as

[TABLE]

where $\Theta_{\alpha\beta\gamma}(n,n^{\prime})=\int_{BZ^{\prime}}d^{2}\mathbf{k}\,\Theta_{\alpha\beta\gamma}(\mathbf{k},n,n^{\prime})$ and the $BZ^{\prime}$ indicates that we are integrating in a small region where the quadratic approximation is enough to describe the bands of Graphene around the valleys with $\tau=\pm 1$. The limitations and proper choice of $BZ’$ is discussed in the main text. In natural units $e=0.0854$ and $694\;\mathrm{eV}^{2}=1\;\mathrm{T}$. Finally $\lambda^{\prime}=0.4\;\mathring{A}\;\cdot v$ and $1\;\mathring{A}=\frac{1}{2000\;eV}$, so that $\lambda=0.06$. Using this method for the calculation of the Green’s function, the very big computational time required for calculations in magnetic fields can be overcome. I also approximated considering just few $k$-points per valley (since the $k$ points are weighted by a Gaussian function and they decay quite fast). Analyzing the structure of the matrix element $\Theta_{\alpha\beta\gamma}(n,n^{\prime})$ it is easy to show that there is a selection rule on $\Delta n=n-n^{\prime}$ for the allowed transitions, that is $\Delta n=\pm 1.$ This becomes evident if we realize that the integral in the matrix element contains an integral over $k$ of Laguerre Functions, that are orthonormal by definition. Expression (61) can be rewritten in less compact, but more physically comprehensible way dividing the landau level number $n$ and the band index $s$ as follows

[TABLE]

The sum over $s$ accounts for the selection rule discussed in the previous paragraph and now the summation is taken just over a single positive Landau level index .