Second order optical nonlinearity of graphene due to electric quadrupole and magnetic dipole effects

We present a practical scheme to separate the contributions of the electric quadrupole-like and the magnetic dipole-like effects to the forbidden second order optical nonlinear response of graphene, and give analytic expressions for the second order optical conductivities, calculated from the independent particle approximation, with relaxation described in a phenomenological way. We predict strong second order nonlinear effects, including second harmonic generation, photon drag, and difference frequency generation. We discuss in detail the controllability of these effects by tuning the chemical potential, taking advantage of the dominant role played by interband optical transitions in the response.

effects in centrosymmetric atoms or molecules arise 5,11-13 ; (2) at an asymmetric interface between graphene and the substrate the centre-of-inversion symmetry is broken, and second order nonlinearities are allowed 14-18 ; (3) similarly, the symmetry can be locally broken due to natural curvature fluctuations of suspended graphene 19 ; (4) the application of a dc electric field can be used to generate an asymmetric steady state, and a second order nonlinear optical response can then arise through the third order nonlinearity 10,[16][17][18]20,21 . Second order optical responses due to the first effect have been shown to be important for photon drag (dynamic Hall effects) 9,22,23 , second harmonic generation (SHG) 13 , and difference frequency generation (DFG) [24][25][26] . [Note added. After the submission, we became aware of related works in preprints 27,28 . Overlapping results in these papers are in agreement in the absence of relaxation]. However, in most of the studies of these phenomena only intraband transitions were considered 13,22,24 . As with the third order optical response, the contribution of interband transitions to the second order nonlinearity can lead to a rich and tunable nonlinear optical response 23,26 . In this work we present analytic results for the second order conductivities of graphene induced by the electric quadrupole-like and magnetic dipole-like effects, within the independent particle approximation and with relaxation processes described phenomenologically.

Model.
We consider the charge current response of a graphene monolayer to electromagnetic fields E(r, t; z) and B(r, t; z) with = +r x y x y , and focus on the second order conductivity σ ω ω q q ( , ; , ) dab (2); 1 1 2 2 , which is defined perturbatively for the weak fields from the second order current response the Roman superscript letters stand for the Cartesian directions x or y, the repeated superscripts imply a sum over all in-plane components, and σ ω ω q q ( , ; , ) dab (2); 1 1 2 2 can be taken to be symmetric in its components and arguments, σ ω ω σ ω ω = q q q q ( , ; , ) ( , ; , ) dab d ba (2); 1 1 2 2 (2); 2 2 1 1 , without loss of generality. In writing Eq. (1) we neglect any response of the graphene to electric fields in the z direction, in line with the usual models for excitation around the Dirac points; thus the Cartesian components only range over x and y. There is no term involving the magnetic field B(r, t; z) in Eq. (1), but it is not neglected. Below we sketch the outline of the derivation from the minimal coupling Hamiltonian, involving the vector and scalar potentials. Keeping powers of q in the expansion of the vector potential introduces the magnetic field, but the final result can be written in the form of Eq. (1) in agreement with the usual convention in nonlinear optics. However, since we focus on the response at reasonably long wavelengths, i.e.,  ω 2 with the electron Fermi velocity v F , then the conductivity can be expanded as We have used the zero second order response to uniform fields σ ω ω = 0 0 ( , ; , ) 0 dab (2); 1 2 , due to the inversion symmetry of the graphene crystal structure. For the D 6h crystal symmetry of graphene, fourth-order tensors S dabc have only three independent in-plane components S xxyy , S xyyx , and S xyxy , and in total eight nonzero in-plane components S yyxx = S xxyy , S yxxy = S xyyx , S yxyx = S xyxy , and = = + + S S S S S yyyy xxxx xxyy x yxy x yyx . The response coefficients S dabc (ω 1 , ω 2 ) completely characterize the second order optical response in the small |q| limit. To calculate them, we begin by writing the minimal coupling Hamiltonian as , where e = − |e| is the electron charge, Ĥ 0 is the unperturbed Hamiltonian, v is the velocity operator in the absence of an external field, A(r, t; z) and φ(r, t; z) are the vector and scalar potentials, respectively. Due to the linear dispersion relation of graphene around the Dirac points, any higher order terms in A(r, t; z) can be neglected, unlike the situation in usual semiconductors, where the calculation can be more difficult; the Zeeman interaction can also be ignored. The vector and scalar potentials are then Fourier expanded, as in Eq. (2), and we write . Considering an electric field described by a vector potential, we get σ (2); 2 1 2 in Eq. (31) [see method], and by expanding both sides at small q 1 and q 2 we find 2 directly, and then for S Q dabc ; in the relaxation free limit, we have checked that the results of S Q dabc calculated using the vector potential are the same as those using the scalar potential only. The simple relaxation time approximation used here is not gauge invariant, which could be recovered by a postprocessing method 29 . In this work, the different calculations give differences only on the order of the relaxation parameters. We leave the gauge invariant relaxation time approximation for a future work.
For atoms, or molecules with center-of-inversion symmetry, the kind of "forbidden" second order processes we are discussing here can be identified with electric quadrupole and magnetic dipole interactions, as opposed to the usual "electric dipole interactions" that typically govern the first order response. Here, however, in a model where electrons are free to move through the graphene, there is no simple way to clearly identify these two processes. We note that while the expression for the full S dabc (ω 1 , ω 2 ) can be derived solely from considering the vector potential, as mentioned above, only its contribution ω ω S ( , ) Q dabc 1 2 can be identified by considering only the scalar potential. Since quadrupole interactions in atoms and molecules exist if only scalar potentials are introduced, while magnetic dipole interactions require a vector potential for their description, we take this as a motivation for ascribing quadrupole effects (or, more properly, quadrupole-like effects) to ω ω S ( , ) Q dabc 1 2 , and for ascribing magnetic dipole effects (or, more properly, magnetic-dipole-like effects) to , and in terms of them the second order current can be written as We now present a microscopic theory to calculate the tensor components. We describe the low energy electronic states ψ sk (r; z) at band index s = ± and wave vector k by a widely used two-band tight binding model based on the carbon 2p z orbitals 6 . Ignoring all response to the z-component of the electric field, the total Hamiltonian can be written as = ω ω Here ε sk is the electron band energy, a sk (t) is an annihilation operator associated with ψ sk (r; z), the integration is over one Brillouin zone (BZ), and H scat is the scattering Hamiltonian described below phenomenologically. The interaction matrix elements at α = 0 give  = D ; , , is the matrix element of the velocity density, with v k s s 1 2 being the velocity matrix elements in the absence of an external field. The dynamics of the system is described by a single particle density matrix ρ  Here the last term describes the scattering effects phenomenologically with one relaxation energy Γ , and is the initial carrier distribution without any external fields, where is the Fermi-Dirac distribution at the chemical potential μ and t he temp erature T . We fo c us on t he c ur rent resp ons e ; , is given by where each quantity is expressed as a 2 × 2 matrix with abbreviated band index, and   . This is not surprising, because such a limit involves the contributions from all electrons in the "− " band, most of which can not be described by the linear dispersion. Nonetheless, the contributions to the second order response from electrons close to the Dirac points are well described by Eqs (12) and (15). Combined with the fact that the conditions ω ω = S ( , ) 0 Q M dabc / ;0,0 1 2 are verified by using the symmetries of the system, the expression in Eq. (15) can be used to describe the response coefficient for optical transitions occurring around the Dirac points.
At finite temperature, we follow the technique used in our previous work 7 to calculate the conductivities as The main results of this section are given in Eqs (7), (15) and (16). Within the linear dispersion approximation that we have assumed, the results are analytic, and any calculation can be performed directly. In the following we discuss the divergences (poles) of the analytic expressions at zero temperature, and then give a quantitative analysis for different second order optical nonlinear phenomena including SHG, one color dc current generation (including current injection effects, photon drag (or dynamic Hall effect)), and DFG.

Features and limitations of the result.
We begin by considering some special limits of the response following from Eqs (7), (15) and (16). In the limit of a chemical potential much greater than any other energies involved, i.e.,  µ ω  i , Γ , k B T, the dominant contribution to the response is expected to come from the intraband transitions between states around the Fermi surface. We can isolate this contribution by considering the , 1 2 and keeping only the leading term that varies as ∝ x −3 , and then setting x = 1. In this limit we find for the intraband contribution

 
Note that except for a sign the results in Eq. (17) are independent of the chemical potential, which through Eq. (16) leads also to an insensitivity to the temperature. There are two kinds of divergences that appear in Eq. (17). The first involves the w i in the denominator; the w i never vanish at real frequencies, and the divergences that would arise were Γ to vanish can be said to be ameliorated by the phenomenological relaxation introduced. The second involves the ω i in the expression for ω ω S ( , ) M xxyy ;intra 1 2 , and are unameliorated. Even in the presence of relaxation these lead to divergences as ω 1 or ω 2 vanishes, and they appear in the term where the vector potential was used in the calculation. In fact, if we evaluate the terms ω ω S ( , ) Q xxyy ;intra 1 2 and ω ω S ( , ) Q xyxy ;intra 1 2 using the vector potential by Eq. (31) instead of the scalar potential, we also find that they acquire unameliorated divergences. We emphasize that in the limit of no relaxation, where unameliorated divergences appear everywhere, the result in Eqs (7), (15) and (16) is independent of the gauge used in the calculations; our results agree with those obtained by Tokman et al. 26 . It is just that the commonly used phenomenological model we have introduced for relaxation is too simple to respect this gauge invariance. So the divergences in our expression for ω ω S ( , ) M xxyy ;intra 1 2 should not be taken seriously; they are artifacts of the relaxation model, and have a parallel in the same way that unameliorated divergences can arise in the linear response of a metal if such a relaxation model is used in conjunction with the use of a vector potential to describe the electric field. We will turn to a more sophisticated treatment of the relaxation in a future work; in this paper our focus will be on features of the response where ω 1 and ω 2 are greater than Γ /ħ from zero, and thus the lack of amelioration of the vector potential divergences will not be crucial.
Generalizing now beyond just the intraband response, we note that in the (q, ω) dependence of the linear conductivity 30,31 σ (1);da (q, ω), there are divergences that arise in the absence of relaxation when one of the resonant Scientific RepoRts | 7:43843 | DOI: 10.1038/srep43843 The first is associated with intraband transitions, and the second with interband transitions. Similar divergences arise here, some involving combinations of wave vectors and frequencies, and their appearance is evident in the denominator of Eq. (12). The general resonant conditions could be met in photonic structures, where |q| can be much larger than ω/c. For light incident from vacuum where the magnitude of the incident wave vectors ω the resonant conditions for incident fields become ω j = 0 for intraband transitions and ω j = ± 2|μ|/ħ for interband transitions, as shown in the analytic expression for ω ω µ S ( , ) dabc ,0 1 2 . For the generated field at ω 1 + ω 2 and q 1 + q 2 , the general resonant condition may be satisfied due to the arbitrary choice of incident angle 24,25 . Since the response to the intraband transitions in Eq. (17) is weakly dependent on the chemical potential, one must rely on the interband contribution to tune the resonant second order response in graphene. All coefficients .
Note that a current perpendicular to q 0 arises only when components of the electric field both parallel and perpendicular to q 0 are present, while a current in the direction of q 0 arises quite generally. In Fig. 1(a,b) we show the response coefficients S 1 (ω) and S 2 (ω) for relaxation parameters Γ = 0.5 and 33 meV at temperatures T = 0 and 300 K, for the chemical potential μ = 0.3 eV. As ω → 0 the description of relaxation is not valid, as we discussed above, so we focus on the behavior away from ω = 0. For zero temperature and a small relaxation parameter the coefficient S 1 exhibits two peaks, one at ħω = |μ| and the other at ħω = 2|μ|; they follow from the analytic expression in Eq. (15). As discussed above, the peaks are as expected for interband resonances, the first associated with a two-photon resonance and the second with a one-photon resonance. With an increase in the relaxation parameter or in the temperature both peaks are lowered and broadened. Although the thermal energy at 300 K (25.8 meV) is slightly smaller than the relaxation parameter 33 meV, it affects both peaks more effectively, which follows from the form of the dependence on the chemical potential in Eqs (15) and (16). The intraband contributions from Eq. (17) are plotted as black dashed curves, which fit with the fully calculated results very well for photon energies ħω < 0.1 eV for the chosen parameters. The coefficient S 2 exhibits similar peaks at ħω = |μ| and ħω = 2|μ|, but there is also a dip at around ħω = |μ|/2. This is also apparent from Eq. (15), and is due to a cancellation of contributions from S Q xxyy and S M xxyy ; in fact, at zero temperature, S 2 (μ/2) ∝ Γ and so it would vanish in the limit of no relaxation.
For the extreme case of q 0 = ω 0 /c, corresponding to light propagating parallel to the graphene sheet, it is natural to introduce an effective SHG susceptibility. Identifying a nominal thickness d gr = 3.3 Å for the graphene sheet, the effective susceptibility can be taken as χ ω ω ω ε 0 is simply proportional to S j (ω 0 ), but its introduction makes it easy to compare the strength of the second order response of graphene with that of materials with an allowed second-order response. The values of the χ ω ( ) j (2) 0 are shown on the right y-axis of Fig. 1. At ħω 0 ~ 0.6 eV, the values for χ ω ( ) 1 (2) 0 vary from 10 5 pm/V at T = 0, Γ = 0.5 meV, to 84 pm/V at T = 300, Γ = 0.5 meV, 154 pm/V at T = 0, Γ = 33 meV, and 63 pm/V at T = 300, Γ = 33 meV. The temperature and relaxation greatly reduce its magnitude. The contribution from the intraband transition is about 51 pm/V, which is insensitive to the temperature and relaxation parameters at this photon energy. At optical frequencies the values obtained here are orders of magnitude smaller than those predicted for the current-induced SHG 10 of graphene and the SHG 32 in a gapped graphene. This is not surprising; effects dependent on the finite size of the wave vector of light are typically weak. However, the values we find are still larger than for most SHG materials 20 , where the process is allowed, which indicates the strong second order optical response of graphene, despite the fact that it must rely on the small wave vector of light.
In Fig. 1(c) we show the chemical potential dependence of |S 1 (ω)| for a fixed photon energy ħω = 0.3 eV. We see that control of the chemical potential can be used to change the size of the SHG coefficients, especially at low temperature.
The second order polarizability constitutes part of the second order response, and has been investigated by Mikhailov 13 . The connection between the nonlinear conductivity discussed here and that polarizability follows from the continuity equation Jr n t z tz ( , ; ) ( , ; ) 0 r t . For the field in Eq. (18), the induced second order charge density is identified as , which is in agreement with Mikhailov's calculation 13 . This also confirms that his expression contains only intraband contributions, and as expected there is no contribution to the second order polarizability from magnetic-dipole-like terms.
Photon drag and one color current injection. For the single-mode incident field in Eq. (18), besides the SHG, the other second order current is a dc one For these coefficients, the poles exist at ħω = 0 or ħω = 2|μ|. Depending on the electric field polarization, the dc current can be in the direction of either q 0 or both q 0 and × z q 0 . The latter can only exist when the electric field has nonzero components along both q 0 and × z q 0 . We check the limit Γ → 0 at zero temperature. The terms in Eq. (22) can be approximated as  At ħω = ± 2μ the coefficients S 3 (ω) and S 4 (ω) diverge as Γ −2 for small enough Γ . So for sufficiently small Γ the small q expansion applied to Eq. (12) becomes suspicious, and the q dependence in the denominator of that equation should be considered explicitly. For small frequencies, both coefficients also diverge as (ħω) −1 . These divergences are associated with the resonant photon drag effect, as discussed by Entin et al. 23 . For S 5 (ω) there is an additional Γ −1 divergence that arises for all photon energy. It only contributes to J dc when ⊥ E 0 and E 0 have different phases, which requires elliptically polarized light. This divergence shows that the dc induced current described by S 5 (ω) behaves as a one-color injected current, similar to that observed in semiconductors without inversion Scientific RepoRts | 7:43843 | DOI: 10.1038/srep43843 symmetry 33 ; here the interference between the two transition amplitudes that can lead to an injected current is associated with the ⊥ E 0 and E 0 components of the electric field. In Fig. 2 we show the response coefficients |S 3 (ω)|, |S 4 (ω)|, and |S 5 (ω)| for relaxation parameters Γ = 0.5 and 33 meV at temperatures T = 0 and 300 K, and chemical potential μ = 0.3 eV. The peaks appearing at ħω = 0 and ħω = 2|μ| are obvious. Similar to the behavior of S 2 (ω) in Fig. 1(b), |S 4 (ω)| in Fig. 2(b) also shows a dip at ħω = |μ| at zero temperature. At finite temperature, the frequency of that dip changes. At zero temperature, S 3 (ω) and S 4 (ω) show a very weak dependence on the relaxation parameters for ħω away from the resonances, while S 5 (ω) shows a significant dependence, and indicates the injection process. The effect of increasing temperature on S 3 (ω) and S 4 (ω) is significant for most of the frequencies studied.

Difference frequency generation.
In the presence of a strong pump field ω E q p p at q p and ω p , the injected signal field