Second harmonic generation in graphene dressed by a strong terahertz field
Mikhail Tokman Sergei B. Bodrov Yuri A. Sergeev Alexei I. Korytin Ivan Oladyshkin Yongrui Wang Alexey Belyanin Andrei N. Stepanov

	
	
	
	
	
	

 Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod, 603950, Russia 



 Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod, 603950, Russia 



 Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod, 603950, Russia 



 Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod, 603950, Russia 



 Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod, 603950, Russia 

	

 Department of Physics and Astronomy, Texas A&M University,
College Station, TX, 77843 USA




 Department of Physics and Astronomy, Texas A&M University,
College Station, TX, 77843 USA



 Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod, 603950, Russia 


	
	
	We observe enhanced second-harmonic generation  in monolayer graphene in the presence of an ultra-strong terahertz field pulse with a peak amplitude of 250 kV/cm. 	This is a strongly nonperturbative regime of light-matter interaction in which particles get accelerated to energies exceeding the initial Fermi energy of 0.2 eV over a timescale of a few femtoseconds. The second-harmonic current is generated as electrons drift through the region of momenta corresponding to interband transition resonance at an optical frequency. The resulting strongly asymmetric distortion of carrier distribution in momentum space gives rise to an enhanced electric-dipole nonlinear response at the second harmonic. We develop an approximate analytic theory of this effect which accurately predicts observed intensity and polarization of the second-harmonic signal. 
	
	
	
	
	
	Second harmonic generation in graphene dressed by a strong terahertz field
     Andrei N. Stepanov
    December 30, 2023
==========================================================================



§ INTRODUCTION


There has been a surge of interest in the nonlinear and quantum optics of graphene <cit.>. High speed of carriers in monolayer graphene, v _F≈ 10^8 cm/s , and the resulting large dipole moment of the optical transition ∼ ev_F/ω, give rise to high magnitudes of the nonlinear optical susceptibilities. A particularly intriguing question raised in recent studies is whether monolayer graphene can support second-order nonlinear processes such as second harmonic or difference frequency generation and other three-wave mixing processes. Since graphene has an in-plane inversion symmetry, the second-order nonlinear response is forbidden in the electric dipole approximation. Of course graphene, like any surface, has anisotropy between in-plane and out-of-plane electron motion, but the corresponding nonlinear response is very weak <cit.>. A much stronger second-order nonlinearity originates from nonlocal response beyond the electric dipole approximation, i.e. magnetic-dipole and electric-quadrupole response <cit.>. In <cit.> second-harmonic generation (SHG) was observed in suspended graphene and attributed to the inversion symmetry breaking due to wrinkles and tears in a monolayer. In a bilayer graphene the inversion symmetry can be broken by applying a voltage bias in transverse direction <cit.>.
 
The most obvious method of creating in-plane anisotropy is the anisotropic perturbation of carrier distribution in the k-space by a constant electric field <cit.>. This can be also achieved by applying a low-frequency field, in particular at THz frequencies. This SHG mechanism corresponds effectively to a third-order nonlinearity. However, since the THz frequency is much lower than the optical frequencies, it is natural to introduce an effective second-order nonlinear response which depends on the THz field amplitude E_0(t) as a parameter <cit.>. For a weak low-frequency field the perturbation of an originally isotropic distribution is localized mainly near the Fermi surface. In this case the parameter of anisotropy is expressed through a low-frequency Drude current J_Dr. The resulting nonlinear response obtained by perturbations turns out to be proportional to the magnitude of J_Dr (see <cit.>). The magnitude of the SH current J_2ω scales as

    J_2ω∝ J_Dr E_ω^2 ,
 
where J_Dr=e^2 τ_Dr|μ|/4ħ E_0, τ_Dr is the relaxation time of a DC current, μ is chemical potential, E_ω is the electric field amplitude of an optical pumping at frequency ω. Equation (<ref>) is valid when the perturbation of electron distribution at Fermi surface is small.

The experimental results that we obtained below apparently agree with the scaling of Eq. (<ref>). This is however very surprising for ultra-strong fields used in our experiment. Indeed, in recent experiments including our experiment the THz pulses of several hundreds fs duration and field amplitudes E_0≈ 100 - 300  kV/cm were used; see e.g. <cit.>. Such fields distort the original Fermi distribution far beyond small perturbations: an electron acceleration to energies of the order of 0.2 eV (which is a typical Fermi energy for CVD graphene on a glass substrate) happens over a timescale of only a few fs <cit.>. In such strong fields,  a standard expression for the current J_Dr using the Drude relaxation time of ∼ 1 ps adopted in <cit.> would yield a current amplitude much higher than  the maximum possible value e v_F N, where N is the surface density of carriers! Furthermore, in ultra-strong fields the carrier density gets multiplied by a large factor during the THz pulse due to the electron-hole pair creation <cit.>. 

Therefore, to interpret our experiments we need a theory which is not restricted to a standard perturbative model of a weak deformation of the Fermi surface. Within our analytic model the anisotropic deformation of the particle distribution is formed primarily in the vicinity of the interband transition resonance between the optical pump and particle states dressed by a low-frequency field. We consider the case when the energies of resonant particles are much higher than the Fermi energy. A dramatic change in the properties of a quantum system dressed by a strong field is a universal effect; see e.g. <cit.> and references therein. In our case this effect shows up as a broadening of the interband resonance due to particle acceleration by a strong low-frequency field E_0 in the process of an interband transition intiated by an optical field. The broadening of the resonant region in the k-space, δ k_s, corresponds to the frequency bandwidth δω = v_F δ k_s ∼δ t_s^-1=√(e v_F |E_0|/ħ), which is equal to the inverse time of Schwinger pair creation in graphene by the field of magnitude E_0 <cit.>. Under the action of an ultra-strong field E_0 the region of resonant perturbation of carriers in the k-space turns out to be asymmetric with respect to the resonant frequency given by ω = 2 k v_F. This asymmetry gives rise to an anisotropic electric-dipole nonlinear response. The case of an ultrafast electron scattering with characteristic time shorter than the Schwinger time is treated numerically in the Appendix.

Within our model the typical lifetime of the particle perturbation in the resonant region is determined by drift in a low frequency field: ħδ k_s / e E_0 ∼δ t_S, i.e. it corresponds to the Schwinger time δ t_S. Therefore, an increase in the magnitude of a THz field should broaden the applicability region of the perturbation theory with respect to the optical field. For the above parameters of the THz pulse the time δ t_S is of the order of 5-8 fs, which is smaller than, or of the same order as scattering times in graphene. In this case the saturation of the resonant absorption of the optical pump is weak up to the field amplitudes of order E_ω∼ 2 - 3 MV/cm. This fact ensures the validity of the scaling J_2ω∝ E^2_ω despite the use of high-power femtosecond lasers. 

Section 2 describes the experiment. Section 3 contains the basic set of equations describing SHG. Section 4 describes an approximate analytical model. In Section 5 the theory is compared with experiment. Derivation of certain formulas used in the analytic theory and numerical simulations for ultrafast scattering times can be found in the Appendix.







§ EXPERIMENT


The schematic of the experimental setup is shown in Fig. <ref>. A Ti-Sapphire laser system (Spitfire, Spectra-Physics) generated pulses of energy 0.7 mJ, central wavelength of 0.795 nm, duration 70 fs and repetition rate 700 Hz. Optical radiation was injected through mirror PM3 parallel to the THz field. The diameter of the optical beam was ≈ 220 μm FWHM, the intensity of the optical beam on the sample was 3 GW/cm^2.


THz pulses were generated in the LiNbO_3 crystal as in <cit.> and focused on the sample at an angle of 45^∘. The diameter of the THz beam on the sample was ≈ 500 μm  FWHM with respect to the field amplitude. The maximum THz electric field amplitude was 250 kV/cm. The maximum value of the P-polarized field was two times smaller than that of the S-polarized field.  

The sample was a CVD graphene monolayer on borosilicate glass. Interaction with substrate led to p-doping to the level of Fermi energy E_F = 0.2 eV <cit.>, which is much smaller than the particle energy of 0.75 eV corresponding to the interband resonance with pumping.  

The SH signal for graphene on glass and for glass substrate only is shown in Fig. <ref> together with theoretical results. Here p polarization corresponds to the field in the incidence plane, whereas s polarization is orthogonal to this plane. In the notations sss, pss etc. the first index is the polarization of the optical pump, the second index is the polarization of the THz field and the third index is the polarization of SH photons. 	

Clearly, significant SH signal existed only for SH photons with the same polarization as the THz field, in agreement with theory. 

Fig. <ref> shows the number of pss-polarized SH photons N_2ω as a function of the delay time between the fs laser pulse and THz pulse, superimposed on the temporal profile of the THz pulse. Clearly, the SHG signal follows the THz field squared, N_2ω∝ E_0^2. Since the number of SH photons is proportional to the current j_2ω squared, this corresponds to the scaling Eq. (<ref>). Inset to Fig. 3 shows the dependence of N_2ω from the energy of a laser pulse, which also agrees with  Eq. (<ref>). 


Note that the SH radiation propagated in the direction of mirror reflection of the incident optical pump, which proves the coherent nature of the SH signal. 






§ THEORETICAL MODEL

		



 §.§ Equations for the density matrix and current


Consider monolayer graphene located in the (x,y) plane. Not too far from the Dirac point, the effective Hamiltonian for carriers interacting with an electromagnetic field is <cit.>:

    Ĥ = v_Fσ̂·( p̂ + e  A / c) - e ϕ ,
 
where p̂ = -iħ( x_0∂/∂ x +  y_0∂/∂ y), σ̂= x_0σ̂_x + y_0σ̂_y, σ_x,y are Pauli matrices, x_0,y_0 are unit vectors, A and ϕ are vector and scalar potentials of an external EM field. The basis functions of the Hamiltonian Eq. (<ref>), defined at A=φ =0, are

    |  k,s ⟩ =e^i k· r/√(2)([       s; e^iθ_ k ]), s=± 1.
 
Here eigenstates |  k,s ⟩ are determined on a unit area for periodic boundary conditions, indices s=± 1 numerate electron and hole states, θ_ k is an angle between the quasimomentum ħ k and x-axis. The eigenenergy corresponding to states in Eq. (<ref>) is given by E_s = sħ v_F√(k_x^2 + k_y^2) <cit.>.

In the electric dipole approximation an external electric field can be considered uniform, E(t). It can be defined through either vector or scalar potential. The gauge invariance of observables calculated by solving the Schrödinger or master equation for the Hamiltonian (<ref>) has been proven in <cit.>. Therefore, any EM gauge can be used. However, using a scalar potential leads to a slightly simpler derivation when solving a density matrix equation; see the comparison in <cit.>. Furthermore, according to <cit.>, gauge invariance requires that the relaxation operator in the density matrix equation be made dependent on the vector potential, which would be an additional complication. Therefore, we define an external field through the scalar potential, assuming ϕ = -  r· E(t), A = 0 in Eq. (<ref>). 

In the absence of relaxation (which we will take into account later) the Von Neumann equation in the basis of Eq. (<ref>) is

    ρ̇_ss'kk'+iυ_F(sk-s'k')ρ_ss'kk'=-i/ħ∑_s”k”(V_ss^”kk”ρ_s^”s'k”k'-ρ_ss^”kk”V_s^”s'k”k'),
 
where

    V_ss'kk'=⟨k,s|.--eϕ|.-k',s'⟩ =e/2∑_x,y[-∂θ_k/∂ k_x,yδ_kk'+i(1+ss')∂δ_kk'/∂ k_x,y]E_x,y(t),   
       ∂δ_kk'/∂ k_x,y=-∂δ_kk'/∂ k'_x,y,     ∂θ_k/∂ k_x=-sinθ_k/k,    ∂θ_k/∂ k_y=+cosθ_k/k,     k=√(k_x^2+k_y^2).
 

For a uniform field E(t) the matrix element of the interaction Hamiltonian V_ss'kk' is diagonal with respect to k, k'. Switching to a continuous k-spectrum, Eqs. (<ref>) and (<ref>) yield a closed set of equations for matrix elements ρ_ss'kk. Denoting the quantum coherence as ρ_+1-1kk=ρ_k, population difference as ρ_-1-1kk-ρ_+1+1kk=Δ_k, and the sum of populations as ρ_+1+1kk+ρ_-1-1kk=n_Σk, we obtain

    (∂/∂ t+iω_k-eE/ħ∂/∂k)ρ_k=-iΩ_k/2Δ_k,
 

    (∂/∂ t-eE/ħ∂/∂k)Δ_k=-iΩ_k(ρ_k-ρ^*_k),
 

    (∂/∂ t-eE/ħ∂/∂k)n_Σk=0,
 
where ω_k =2v_F k is the interband transition frequency, Ω_k=e/kħ(sinθ_kE_x-cosθ_kE_y) the Rabi frequency. Equation (<ref>) is separated from the rest of the system and we won't use it anymore.

Eqs. (<ref>) and (<ref>) make a closed set of equations, which is a version of semiconductor Bloch equations <cit.>. Note that the “convective” terms ∼∂/∂k in Eqs. (<ref>)-(<ref>) do not originate from some kind of phenomenological assumptions, e.g. an attempt to make it look like a Boltzmann-type equation or quasiclassical equations of motion with a Berry field. These terms rigorously follow from the matrix elements of the interaction Hamiltonian Eq. (<ref>) in a given gauge. If we defined the same field E(t) through the vector potential, we would get the equations in a different form, see e.g. <cit.>. This would not change the observables of course.  

The current operator is defined by Ĵ=-(i e / ħ) [Ĥ,r]=-ev_Fσ̂. The observable current is J = ∑_s,s',k=k'j_ss'kk'ρ_s'sk'k, where j_ss'kk' are matrix elements of the current operator. Using ρ_ss'kk'∝δ_kk'  and performing integration in k-space, we obtain 

    j_x=ev_Fg/4π^2∫_∞d^2k[cosθ_kΔ_k+isinθ_k(ρ_k-ρ^*_k).,
 

    j_y=ev_Fg/4π^2∫_∞d^2k[sinθ_kΔ_k-icosθ_k(ρ_k-ρ^*_k).,
 
where g = 4 is spin and valley degeneracy. The first term in square brackets in Eqs. (<ref>),(<ref>) is the intraband current and the second term is an interband current. 





 §.§ Equations for slow variables


Now we make the following ansatz which separates the slow THz dynamics (subscript 0) from fast optical frequencies:

    E(t) = E_0(t)+ReE_ω(t)e^-iω t
 

    Ω_k(t) = Ω_0k(t)+ReΩ_ωk(t)e^-iω t
 

    Δ_k = Δ_0k(t)+ReΔ_ωk(t)e^-iω t+ReΔ_2ωk(t)e^-i2ω t
 

    ρ_k=ρ_0k(t)+ρ^(-)_ωk(t)e^-iω t+ρ^(+)_ωk(t)e^+iω t+ρ^(-)_2ωk(t)e^-2iω t+ρ^(+)_2ωk(t)e^+i2ω t
 
It also introduces envelopes of the optical fields at the fundamental and second harmonics.

In addition, we add the relaxation operator in its simplest form to the right-hand side of  Eqs. (<ref>), (<ref>). This leads to the replacement  ω_k⟹ω_k-iγ in Eq. (<ref>) and the population relaxation term -Γ(Δ_k-Δ_kF) in Eq. (<ref>). Here γ , Γ are relaxation rates  of the coherence and population difference and  Δ_kF is the equilibrium Fermi distribution. 

The resulting equations below contain only these slow-varying envelopes. However, we still keep counter-rotating terms, as shown below.

    (∂/∂ t+γ +iω_k-eE_0/ħ∂/∂k)ρ_0k=-i/2Ω_0kΔ_0k-i/4Re(Ω^*_ωkΔ_ωk)+e/2ħ(E^*_ω∂ρ^(-)_ωk/∂k+eE_ω∂ρ^(+)_ωk/∂k)
 

    (∂/∂ t+γ +iω_k-iω -eE_0/ħ∂/∂k)ρ^(-)_ωk=-i/4(Ω_ωkΔ_0k+Ω_0kΔ_ωk)+1/2eE_ω/ħ∂ρ_0k/∂k
 

    (∂/∂ t+γ +iω_k+iω -eE_0/ħ∂/∂k)ρ^(+)_ωk=-i/4(Ω^*_ωkΔ_0k+Ω_0kΔ^*_ωk)+1/2eE^*_ω/ħ∂ρ_0k/∂k
 

    (∂/∂ t-eE_0/ħ∂/∂k)Δ_0k+Γ(Δ_0k-Δ_Fk)=-iΩ_0k(ρ_0k-ρ^*_0k)+e/2ħRe(E^*_ω∂Δ_ωk/∂k)+
 

    +Im(Ω^*_ωkρ^(-)_ωk+Ω_ωkρ^(+)_ωk)
 

    (∂/∂ t-iω -eE_0/ħ∂/∂k)Δ_ωk+ΓΔ_ωk=-iΩ_ωk(ρ_0k-ρ^*_0k)-i2Ω_0k(ρ^(-)_ωk-ρ^(+)*_ωk)+e/ħE_ω∂Δ_0k/∂k
 
In Eqs. (<ref>)-(<ref>) the contribution of the terms oscillating at 2ω is neglected since the SH field is very weak. The perturbation at the second harmonics is described by

    (∂/∂ t+γ +iω_k-i2ω -eE_0/ħ∂/∂k)ρ^(-)_2ωk=1/2eE_ω/ħ∂ρ^(-)_ωk/∂k-i/8Ω_ωkΔ_ωk
 

    (∂/∂ t+γ +iω_k+i2ω -eE_0/ħ∂/∂k)ρ^(+)_2ωk=1/2eE^*_ω/ħ∂ρ^(+)_ωk/∂k-i/8Ω^*_ωkΔ^*_ωk
 

    (∂/∂ t-i2ω -eE_0/ħ∂/∂k)Δ_2ωk+ΓΔ_2ωk=-iΩ_ωk(ρ^(-)_ωk-ρ^(+)*_ωk)+e/2ħE_ω∂Δ_ωk/∂k
 
In numerical modeling the solution to Eqs. (<ref>) - (<ref>) was substituted into Eqs. (<ref>) - (<ref>), whereas the solution to Eqs. (<ref>) - (<ref>) was substituted into Eqs. (<ref>) - (<ref>) to find the SH current.






§ AN APPROXIMATE ANALYTIC SOLUTION




 §.§ Equations in rotating wave approximation


Within the rotating wave approximation (RWA) we can solve for the dynamics of carriers in the vicinity of an interband transition at the fundamental frequency of the optical field. The RWA corresponds to the following inequalities:

    |ω_k-ω|,  γ,  Γ≪ω and |eE_0/ħ∂/∂k|∼1/δ t_s≪ω .

In this case one can neglect counter-rotating terms ρ^(+)_ωke^+iω t and ρ^(+)_2ωke^+2iω t in terms  of <ref>. Furthermore, if the following conditions are satisfied,

    |Ω_0k|, |Ω_ωk|≪(γ +1/δ t_s), (Γ+1/δ t_s)≪ω.
 
One can assume that |Δ_0k|≫|Δ_ωk|, |ρ_0k|≪|ρ^(-)_ωk| in <ref> and obtain approximate equations, 

    (∂/∂ t+γ +iω_k-iω -eE_0/ħ∂/∂k)ρ_ωk=-i/4Ω_ωkΔ_0k,
 

    (∂/∂ t-eE_0/ħ∂/∂k)Δ_0k+Γ(Δ_0k-Δ_Fk)=Im(Ω^*_ωkρ_ωk),
 

    ρ_2ωk≈i/ωeE_ω/2ħ∂ρ_ωk/∂k,   Δ_2ωk≈Ω_ωk/2ωρ_ωk,
 
where ρ_k≈ρ_0k+ρ_ωke^-iω t+ρ_2ωke^-2iω t.





 §.§ The stationary phase solution of density matrix equations 


For simplicity we consider the case when optical and THz fields are polarized along x; other polarizations are considered in Sec. IV D. Using Ω_ωk(t)=eħ^-1sinθ_kE_ω(t)/k=k_yeħ^-1E_ω(t)/k^2_y+k^2_x and eE_0(t)/ħ∂/∂k=eE_0(t)/ħ∂/∂ k_x, the formal solution of Eq. (<ref>) can be written as

    ρ_ωk(t)=-i/4∫^t_-t_0e^iΨ_k(t',t)-γ·(t-t')Ω_ωk(t',t)Δ_0k(t',t)dt',
 
where

    Ψ_k(t',t)=-∫^t_t'ω_k(τ ,t)dτ +ω·(t-t').
 
The functions Ω_ωk(t',t), Δ_0k(t',t) and ω_k(τ ,t) are obtained from functions  Ω_ωk(t), Δ_0k(t) and ω_k=2υ_F√(k^2_y+k^2_x) as

    Ω_ωk(t',t) = k_yeħ^-1E_ω(t')/k^2_y+(k_x+e/ħ∫^t_t'E_0(τ)dτ)^2,
 

    Δ_0k(t',t) = Δ_0k(t',k_y,k_x+e/ħ∫^t_t'E_0(τ)dτ),
 

    ω_k(τ ,t) = 2v_F√(k^2_y+(k_x+e/ħ∫^t_τE_0(τ')dτ')^2).
 

It is easy to check that <ref> is an exact solution to Eq. (<ref>), satisfying the initial condition ρ_ωk(-t_0)=0, where t=-t_0 is the time when the optical field is turned on. 

The integral in Eq. (<ref>) can be evaluated by the method of stationary phase. Within this approach one assumes that the main contribution to the integral comes from the vicinity of a stationary point given by the condition  

    .∂Ψ_k(t',t)/∂ t'|_t'=t_s=.ω_k(t',t)|_t'=t_s-ω =0.
 
Expanding Ψ_k(t',t) near the stationary point t'=t_s(t,k), we obtain

    ρ_ωk(t)≈ -i/4e^iΨ_k(t_s,t)-γ(t-t_s)Ω_ωk(t_s,t)Δ_0k(t_s,t)∫^t-t_s_-(t_0+t_s)e^iα_k(t_s,t)τ^2dτ.
 
where Ψ_k(t_s,t)=.Ψ_k(t',t)|_t'=t_s,

    α_k(t_s,t)=1/2.∂^2Ψ_k(t',t)/∂ t^'2|_t'=t_s=1/2.∂ω_k(t',t)/∂ t'|_t'=t_s.
 

The interval  in Eq. (<ref>) which makes the main contribution to the integral is given by

    Δτ =√(1/|α_k|)∼δ t_S=√(ħ/ev_F|E_0|),
 
The following hierarchy of timescales ensures the validity of the stationary phase method:


    δ t_S≪Δ t_opt,  Δ t_THz,  γ^-1 , where Δ t_opt and Δ t_THz are durations of the optical and THz pulses. Furthermore, we assume that Δ t_opt≪Δ t_THz,  δ t_sω≫ 1,


Eq. (<ref>) can be written in the form (see Appendix A) 

    ρ_ωk=-i/4Ω_ωk(t_s,t)Δ_0k(t_s,t)π/2v_F𝔚(Z,δ k_s),
 
where

    𝔚(Z,δ k_s)=e^± iZ^2 /δ k_s(e^∓ iπ/4/√(π)∓2/π∫^Z_0e^∓ ix^2dx),
 
Z=k-ω/2v_F/δ k_s, δ k_s=√(|cosθ_keE_0/ħ v_F|). The top and bottom signs in Eq. (<ref>) correspond to cosθ_kE_0 > 0 and cosθ_kE_0 < 0 respectively. 

The function 𝔚 is normalized as ∫^∞_-∞𝔚(Z,δ k_s)δ k_sdZ=1. In the limit δ k_s→ 0 (i.e. when E_0→ 0) Eqs. (<ref>), (<ref>) give

    ρ_ωk≈1/4Ω_ωkΔ_0k[𝒫/ω -ω_k -iπδ(ω -ω_k∓ϵ)],
 
where 𝒫 indicates a principal value of the integral, ϵ→ +0.

One can see from the expression in Eq. (<ref>) that the scattering-induced broadening of the resonance δω∼γ is replaced by the nonlinear field-induced broadening: δω∼1/δ t_s∝√(|E_0|). The absorption line described by Eq. (<ref>) is asymmetric, namely it is shifted towards ω >ω_k for cosθ_kE_0>0 and towards ω <ω_k for cosθ_kE_0<0. An asymmetric shape is due to the drift of carriers in the k-space in the presence of a THz field.  

The functions Ω_ωk(t_s,t) and Δ_0k(t_s,t) in Eq. (<ref>) are Lagrangian variables Ω_ωk(k(t),t) and Δ_0k(k(t),t), determined at the current moment of time t shifted with respect to a stationary point t_s given by Eq. (<ref>). The shift amount is t-t_s≈∓Z/δ k_sv_F∼ Zδ t_s, see Appendix. The possibility to use a “local” approximation Ω_ωk(t_s,t)≈Ω_ωk(t) and Δ_0k(t_s,t)≈Δ_0k(t) is discussed in Appendix A.  It turns out that the first and the second inequality in Eqs. (28a) always ensure the validity the local approximation for the Rabi frequency. At the same time, the local approximation for the population difference Δ_0k(t) requires in addition that the perturbation of the initial value Δ_0k=Δ_Fk in the resonant region be small enough. The perturbation of the population difference can be estimated using Eqs. (<ref>), (<ref>). The analysis in Appendix A yields 

    |Δ_0k-Δ_Fk|∼δ t^2_s|Ω_ωk|^2/4.
 

Equation (<ref>) looks like a standard result for a two-level system coupled to an optical field if we replace the population relaxation time with δ t_S in the density matrix equations. This is to be expected, since the lifetime of the perturbation in the presence of a strong low-frequency field is determined by the drift of carriers out of the resonant region: ħδ k_s/eE_0∼δ t_S. Using δ t^2_s|Ω_ωk|^2/4≈|E_ω|^2/|E_0|^21/δ t_s^2ω^2 , ω∼ 2.5 fs^-1 and δ t_s longer than 5 fs (i.e. for E_0 smaller than 300 kV/cm) one can obtain |Δ_0k-1|≈ 0.3 for |E_ω|^2/|E_0|^2≈ 45. Here we also assumed Δ_Fk=1, which is the case for high enough electron energies corresponding to the interband resonance ω_k≈ω.





 §.§ The nonlinear current


For the electric field polarized along x, Eqs. (<ref>) and (<ref>) give the following expression for the complex amplitude of the current at SH: 

    j_x2ω=-e v_F g/8ωπ^2∫_∞(Ω_ωkρ_ωkcosθ_k)d^2k.
 
Substituting Eq. (<ref>) into Eq. (<ref>) and assuming Δ_0k(t_s,t)≈ 1 and Ω_ωk(t_s,t)≈Ω_ωk(t) leads to

    j_x2ω≈ ige^3E^2_ω/16πħ^2ω∫^π_-πdθ_k∫^∞_0dk/k𝔚(Z,δ k_s)sin^2θ_kcosθ_k ,
 
where 𝔐(Z,δ k_s) is defined in Eq. (<ref>). Limiting ourselves to the resonant region, we can expand 1/k≈1/k_res-δ k_sZ/k^2_res+δ k^2_sZ^2/2k^3_res… when evaluating the inegral in Eq. (<ref>) (here k_res=ω/2v_F) and change from ∫^∞_0dk… to infinite integration limits when integrating over dZ=dk /δ k_s:

    j_x2ω≈ ige^3E^2_ω/16πħ^2ω∫^π_-πdθ_k∫^∞_-∞dZ(1/k_res-δ k_sZ/k^2_res+δ k^2_sZ^2/2k^3_res…)δ k_s𝔚(Z,δ k_s)sin^2θ_kcosθ_k .
 
The first two terms in the expansion in Eq. (<ref>) give zero after angle integration ∫^π_-πdθ_k… . For the first term this is obvious due to the normalization condition ∫^∞_-∞𝔐(Z,δ k_s)δ k_sdZ=1; the proof for the second term is in the Appendix. To calculate the integral for the third term in the expansion we take into account that the function 𝔐 is a linear combination of the even function e^± iZ^2 and odd function e^± iZ^2∫^Z_0e^∓ ix^2dx, see Eq. (<ref>). As a result, we obtain

    j_x2ω≈ ige^3v^3_FE^2_ω/4πħ^2ω^4∫^π_-πdθ_ksin^2θ_kcosθ_k|cosθ_k eE_0/ħ v_F|∫^∞_-∞e^± i(Z^2- π/4) /√(π)Z^2dZ ,
 
where ∫^∞_-∞e^± i(Z^2- π/4) /√(π)Z^2dZ=∓ i. Recalling that the upper and lower signs correspond to cosθ_kE_0>0 and cosθ_kE_0<0, we get

    j_x2ω≈ge^4v^2_FE^2_ωE_0/4πħ^3ω^4∫^π_-πdθ_ksin^2θ_kcos^2θ_k=e^4v^2_FE^2_ωE_0/4ħ^3ω^4.
 





 §.§ Polarization dependence



Eq. (<ref>) has been obtained for collinear vectors E_ω∥E_0 ∥j_2ω. Consider now a different orientation when E_ω⊥E_0 ∥j_2ω. Let E_0 ∥x_0 and j_2ω∥x_0, but E_ω∥y_0. The expression for the Rabi frequency Ω_ωk takes the form Ω_ωk=-cosθ_k/keE_ω/ħ. It is straightforward to obtain that the relationship between the SH current j_x2ω and Ω_ωk has the same form as in Eq. (<ref>). Therefore, instead of Eq. (<ref>) we obtain

    j_x2ω≈ge^4v^2_FE^2_ωE_0/4πħ^3ω^4∫^π_-πdθ_kcos^4θ_k=3e^4v^2_FE^2_ωE_0/4ħ^3ω^4.
 
Comparing Eqs. (<ref>) with (<ref>) one can see that the power of the SHG signal at E_ω⊥E_0 is 9 times higher than at E_ω∥E_0 for the same magnitudes of the fields.

In the case E_0 ⊥j_2ω we always get j_2ω=0 for any polarization of the optical field. This is expected, since for this geometry the low-frequency field cannot break the inversion symmetry along the direction of j_2ω see the Appendix. 





§ COMPARISON WITH EXPERIMENT


The field dependence and polarization dependence of the observed SH signal coincides with those predicted by the model. To compare absolute numbers of SH photons, one needs to take into account that 

(i) Eqs. (<ref>), (<ref>) contain the components of the fields tangential to the monolayer. Using Fresnel formulas <cit.>, one can get  

    E^(S)_τ=2E^(S)cosθ/√(ε -sin^2θ)+cosθ,  
    E^(P)_τ=2E^(P)cosθ√(ε -sin^2θ)/√(ε -sin^2θ)+εcosθ.
 
Here E^(S,P) is the amplitude of the transverse field of an incident S- or P-polarized wave, E^(S,P)_τ is its tangential component on the monolayer, ε the dielectric constant of the substrate, θ the incidence angle from air with respect to the normal to the surface of graphene.  

(ii) The surface current j^(S,P)_2ω  at frequency 2ω and the amplitude E^(S,P)_2ω of the S- or P-polarized field radiated by this current are related by

    E^(S)_2ω=-4π/cj^(S)_2ω/√(ε -sin^2θ)+cosθ, E^(P)_2ω=-4π/cj^(P)_2ω√(ε -sin^2θ)/√(ε -sin^2θ)+εcosθ,
 
where θ is the reflection angle.

(iii) The powers of the S- and P-polarized THz radiation in the experiment differ by a factor of 4. 

For numerical estimates we assume ε_ω≈ε_2ω≈ 2.25, ε_0≈ 6.25, θ≈45^^∘.

Next we calculate the number of photons in the PSS polarization configuration. From Eqs. (<ref>) - (<ref>) for the above parameters ε_ω, ε_2ω, ε_0 and θ we can estimate the field of the SH signal as

    E^(S)_2ω≈ 0.34π/ce^4v^2_F(E^(P)_ω)^2 E^(S)_0/8ħ^3ω^4,
 
which is equal to 3.8 ×10^-3 esu. For the pulse duration and beam cross-section used in the experiment, and for a 7% experimental efficiency of the detection system, the SH field in Eq. (<ref>) corresponds to about 850 SH photons per series of 60,000 laser pulses. This agrees with an experimentally measured SH signal within the experimental accuracy. 

Relative SHG efficiency for other polarizations can be obtained from Eqs. (<ref>)-(<ref>):

    N^PSS_2ω/N^SSS_2ω=9[α^(P/S)_ω]^2,    N^PSS_2ω/N^SPP_2ω=4[α^(P/S)_ω]^2/α^(P/S)_0α^(P/S)_2ω,   N^PSS_2ω/N^PPP_2ω=36/α^(P/S)_0α^(P/S)_2ω,
 
where

    α^(P/S)_ω ,2ω ,0=(ε_ω ,2ω ,0-sin^2θ)[√(ε_ω ,2ω ,0-sin^2θ)+cosθ/√(ε_ω ,2ω ,0-sin^2θ)+ε_ω ,2ω ,0cosθ]^2.
 
For the same parameters we get α^(P/S)_ω ,2ω≈ 0,84, α^(P/S)_0≈ 1.2, which gives N^PSS_2ω/N^SSS_2ω≈ 6.4, N^PSS_2ω/N^SPP_2ω≈ 2,8, N^PSS_2ω/N^PPP_2ω≈ 36.

These theoretical values are compared with experiment in Fig. 2. As is clear from the figure, there is a good agreement between theory and experiment. 



 
This work has been supported by the RFBR grant No. 18-29-19091. Y.W. and A.B. acknowledge the support by Air Force Office for Scientific Research
through Grants No. FA9550-17-1-0341 and FA9550-14-1-0376.










§ THE STATIONARY PHASE SOLUTION FOR THE DENSITY MATRIX EQUATION


According to the first inequality in Eq. (28a) the optical pulse is longer than the integration interval in Eq. (<ref>). Therefore the lower limit of integration can be taken as -∞; therefore, Eq. (<ref>) yields

    ρ_ωk(t)≈ -i/4e^iΨ_k(t_s,t)-γ(t-t_s)Ω_ωk(t_s,t)Δ_0k(t_s,t)×(e^± iπ/4/2√(π/|α_k(t_s,t)|)+∫^t-t_s_0e^iα_k(t_s,t)x^2dx) ,

where the upper and lower signs in e^± iπ/4 correspond to α_k>0 and α_k<0 respectively. Furthermore, from Eqs. (<ref>), (<ref>) one can get

    .∂Ψ_k(t',t)/∂ t'|_t'=t_s=ω_k(t_s,t)-ω =2v_F√(k^2_y+(k_x+eħ^-1∫^t_t_sE_0(τ)dτ)^2)-ω =0 .

Since the THz pulse is much longer than the optical pulse, one can write the integral ∫^t_t_sE_0(τ)dτ as 

    ∫^t_t_sE_0(τ)dτ≈ E_0 ·(t-t_s) ,

where E_0 can be treated as a constant during the optical pulse. 

The calculations are greatly simplified if we assume 

    e/ħ|E_0|(t-t_s)≪ k∼ k_res=ω/2v_F ;

we will check its validity later. 

Substituting Eq. (<ref>) into Eq. (<ref>) and using Eq. (<ref>), we obtain:

    t-t_s≈ω -ω_k/2cosθ_kv_FeE_0/ħ ;

Using Eq. (<ref>) in Eqs. (<ref>), (<ref>), gives

    Ψ_k(t_s,t) ≈(ω -ω_k)^2/4cosθ_kv_FeE_0/ħ , 
       α_k(t_s,t) ≈ -cosθ_kv_FeE_0/ħ .

Substituting Eqs. (<ref>), (<ref>) into Eq. (<ref>) results in

    ρ_ωk(t)    ≈ -i/4e^i(ω -ω_k)^2/4cosθ_kev_FE_0/ħ - γω -ω_k/2cosθ_kev_FE_0/ħΩ_ωk(t_s,t)Δ_0k(t_s,t) 
       ×(e^∓ iπ/4√(π/|cosθ_kev_F E_0/ħ|)+∫^(ω -ω_k/2cosθ_kv_F eE_0/ħ)_0e^-icosθ_kev_F E_0/ħτ^2dτ) ,

where the upper and lower signs in e^∓ iπ/4 are for cosθ_kE_0>0 and cosθ_kE_0<0 respectively. Taking into account that the integration interval which makes the main contribution near the resonance is |ω -ω_k|δ t_s∼ 1 and using <ref>, we obtain

    eħ^-1|E_0|(t-t_s)/k_res∼1/δ t_s2v_Fk_res∼1/δ t_sω,
 
which means that the last inequality in <ref> ensures the validity of the approximation <ref>. In the region |ω -ω_k|δ t_s∼ 1 the factor exp(- γω -ω_k/2v_Fcosθ_keħ^-1E_0) in <ref> cannot be greater than exp(- γδ t_s), so if the last inequality in <ref> is satisfied, it can be taken as 1. As a result, we obtain Eqs. (<ref>), (<ref>). We also provide here useful asymptotics of the function 𝔚(Z,δ k_s) in Eq. (<ref>) at |Z|≫ 1:

    𝔐(Z,δ k_s)≈1/iπδ k_s1/Z+{[                         0                          ; 2e^± i(Z^2-π/4)/√(π)δ k_s                           ]. ,

where the upper and lower signs e^± i(Z^2-π/4) are taken for cosθ_kE_0>0 and cosθ_kE_0<0 respectively.





§ APPROXIMATE EXPRESSIONS FOR Ω_ΩK(T_S,T) AND Δ_0K(T_S,T) 


The functions Ω_ωk(t_s,t),  Δ_0k(t_s,t) can be written for the same conditions <ref> and <ref> as

    Ω_ωk(t_s,t) 
       ≈Ω_ωk(t)-ω -ω_k/2v_Fcosθ_keE_0/ħ(∂/∂ t-eE_0/ħ∂/∂ k_x)Ω_ωk
       =sinθ_k/ke/ħ(1-ω -ω_k/2v_Fcosθ_keE_0/ħ∂/∂ t-ω -ω_k/kv_F)E_ω(t)       
    
    Δ_0k(t_s,t)
       ≈Δ_0k(t)-ω -ω_k/2v_Fcosθ_keE_0/ħ(∂/∂ t-eE_0/ħ∂/∂ k_x)Δ_0k(t)

 
Taking into account that |ω -ω_k/2v_Fcosθ_keE_0/ħ|∼δ t_s, and δ t_s≪Δ t_opt, we can drop the terms with time derivative ∂/∂ t in <ref> and <ref>. The factor ω -ω_k/kv_F∼δ k_s/k is small too. At the same time the term eE_0/ħω -ω_k/2v_Fcosθ_keE_0/ħ∂Δ_0k(t)/∂ k_x∼∂Δ_0k(t)/∂ k_x·k - ω/2v_F/cosθ_k∼|Δ_0k-Δ_FK|/δ k_sk - ω/2v_F/cosθ_k in <ref> is small only if the perturbation of the population difference in the resonant region is small too. Therefore, we obtain 

    Ω_ωk(t_s,t)≈Ω_ωk(t)         
    
       Δ_0k(t_s,t)≈Δ_0k(t)-∂Δ_0k(t)/∂ k_x·k - ω/2v_F/cosθ_k





§ QUASISTATIONARY PERTURBATION OF POPULATIONS


At high carrier energies ω_k≈ω one can take Δ_Fk=1 in Eq. (<ref>). Consider a stationary solution of this equation:   

    Δ_0k-1=-ħ/eE_0e^ħΓ/eE_0k_x∫^k_x_Ce^-ħΓ/eE_0k'_xIm[Ω^*_ωk'ρ_ωk']dk'_x ,

where the notation k' means k'= k(k'_x,k_y). The boundary C of the integration limit in <ref> is chosen where the effective source Im[Ω^∗_ωk'ρ_ωk'] approaches zero. This choice depends also on the sign of E_0, which determines the direction of particle drift in the k-space: for E_0 >0 we get (k_x-k'_x)<0, whereas for E_0 <0 we get (k_x-k'_x)>0. Next, we substitute expressions Eqs. (<ref>), (<ref>) into <ref> and take into account that the characteristic “size” of the source Im[Ω^∗_ωk'ρ_ωk'] in the k-space δ k_s∼1/δ t_sv_F. As a result, under the condition ħΓ/eE_0δ k_s≪ 1 we arrive at Eq.(<ref>) for the perturbation of populations. Note that the inequality ħΓ/eE_0δ k_s≪ 1 is equivalent to δ t_sΓ≪1. 




§ THE SECOND TERM IN THE NONLINEAR CURRENT EXPANSION


Here we prove that the second term in the expansion in Eq.(<ref>) is equal to zero. Its expression is given by

    I_2=A ·lim_ϵ⟶ 0∫^π_-πdθ_k∫^∞_-∞dZe^-ϵ Z^2 Zδ k^2_s𝔚(Z,δ k_s)sin^2θ_kcosθ_k ,
 
where ϵ >0 and A is a constant. We introduced the factor e^-ϵ Z^2 which makes the proof easier, and we will take the limit ϵ =0 in the end. Taking into account that the function 𝔐 is a linear combination of the even function e^± iZ^2 and odd function e^± iZ^2∫^Z_0e^∓ ix^2dx (see Eq. (<ref>)), we obtain

    I_2=A ·lim_ϵ⟶ 0∫^π_-πdθ_ksin^2θ_kcosθ_k√(|cosθ_k|)Y ,
 
where

    Y=∓∫^∞_-∞Ze^(± i-ϵ)Z^2 (∫^Z_0e^∓ ix^2dx)dZ=√(π/ϵ)/2(i∓ϵ).
 
In the last expression upper and lower signs are given by the signs of cosθ_kE_0, as usual. Therefore we get 

    I_2=A·lim_ϵ⟶ 0√(π/ϵ)·ϵ/ϵ^2+1∫^π/2_-π/2dθ_ksin^2θ_kcos^3/2θ_k=0.
 




§ POLARIZATION SELECTION RULES FOR THZ FIELD-INDUCED SHG


Consider the orientation E_ω∥y_0, E_0∥y_0 and j_2ω∥x_0. In this case one should take Ω_ωk=-cosθ_k/keE_ω/ħ, eE_0(t)/ħ∂/∂k=eE_0(t)/ħ∂/∂ k_y and ρ_2ωk≈i/ωeE_ω/2ħ∂ρ_ωk/∂ k_y in Eqs. (<ref>)-(<ref>). The dependence of j_x2ω on Ω_ωk has the same form as in Eq. (<ref>), whereas in Eq. (<ref>) one has to replace δ k_s=√(|cosθ_k eE_0/ħ v_F|) with δ k_s=√(|sinθ_k eE_0/ħ v_F|). Upper and lower signs in all coefficients in Eq. (<ref>) correspond to the signs of sinθ_kE_0. As a result, instead of Eq.(<ref>) we get

    j_x2ω≈g e^4v^2_FE^2_ωE_0/8πħ^3ω^4∫^π_-πdθ_ksinθ_kcos^3θ_k=0.
 

If E_ω∥x_0, E_0∥y_0 and j_2ω∥x_0, similar considerations lead to

    j_x2ω≈ge^4v^2_FE^2_ωE_0/8πħ^3ω^4∫^π_-πdθ_ksin^3θ_kcosθ_k=0.
 





§ NUMERICAL SIMULATIONS


We used <ref> to simulate the SHG in  graphene illuminated with a strong THz pulse and an optical field beyond the stationary phase approximation. To derive these equations, the amplitudes of optical-frequency coherences  and population differences Δ_ωk and Δ_2ωk were assumed to be slow-varying. When the optical field is far off-resonance from an interband transition, Rabi oscillations can have a frequency comparable to the frequency detuning of the optical field. However, the slow-varying assumption is still valid if Rabi oscillations are strongly damped by ultrafast  dephasing processes. So, we use ultrafast population relaxation and dephasing times for hot photoexcited electrons, Γ^-1 = 25 fs  and γ^-1 = 2.5 fs, which is consistent in order of magnitude with results from related studies <cit.> 
 and can be attributed to strong Coulomb interaction between carriers in graphene which results in ultrafast carrier-carrier scattering through interband and intraband Auger recombination and impact ionization <cit.>. (Auger recombination could be enhanced due to lattice imperfections in CVD graphene <cit.>.) Then we can still assume that  the optical-frequency populations and coherences  follow the source terms adiabatically, namely, we can put the ∂/∂ t to be zero in all equations except those for Δ_0k. Also, for reasons already discussed above we can assume the optical field to be weak enough to treat it in a perturbative way. 

Another technical difficulty is that Ω_k has a singularity at |k|=0, which can lead to divergence in numerical simulations. To avoid this problem, we replaced k in the denominator of Ω_k by k+ϵ. We also assumed the chemical potential μ_F = 200 meV and electron temperature at equilibrium T_e = 300 K. The THz field is chosen to be polarized in y-direction, and the optical field is polarized in x-direction. 








The simulation shows that the SH current is generated predominantly in y-direction, i.e. along the direction of the THz field. In Fig. <ref> we plot the SH current calculated from the simulation, together with the profile of the THz pulse. The SH current generally follows the THz field, except for some small variations originated from the time evolution of the carrier distribution Δ_0k. These variations are likely beyond the detector resolution in the experiment. 

Dressing of electron states by the THz field is expected to play an important role in the SHG process, even in the presence of ultrafast scattering. As an illustration, we assume that the carrier distribution Δ_0k is in thermal equilibrium, and calculate the SH current for two cases, with and without e E_0/ħ∂/∂k terms included. The latter would be similar to four-wave mixing in a two-level medium.  Figure <ref> shows the dependence of the SH current on the THz field amplituce for the two cases. It indicates that the dressed-state effect can enhance the SH current by one order of magnitude, and therefore the signal intensity by two orders. We can also see that the SH current at the highest THz field of 250 kV/cm is about three times higher in Fig. 5 as compared to Fig. 4 where the thermal equilibrium distribution for Δ_0k was assumed. One could say that roughly 1/3 of the SH current comes from direct parametric interaction between a THz field and an optical field, whereas the nonequilibrium distortion of carrier distribution contributes the remaining 2/3.










54
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Bonaccorso et al.(2010)Bonaccorso, Sun, Hasan, and Ferrari]bonaccorso2010graphene
  author author F. Bonaccorso, author Z. Sun,
  author T. Hasan,  and author A. Ferrari, @noop 
  journal journal Nature photonics volume 4, pages 611 (year
  2010)NoStop
[Kumar et al.(2013)Kumar,
  Kumar, Gerstenkorn, Wang,
  Chiu, Smirl, and Zhao]kumar2013
  author author N. Kumar, author J. Kumar,
  author C. Gerstenkorn, author R. Wang, author
  H.-Y. Chiu, author
  A. L. Smirl,  and author
  H. Zhao, 10.1103/PhysRevB.87.121406 journal journal Phys.
  Rev. B volume 87, pages 121406
  (year 2013)NoStop
[Hong et al.(2013)Hong,
  Dadap, Petrone, Yeh,
  Hone, and Osgood]hong2013
  author author S.-Y. Hong, author J. I. Dadap,
  author N. Petrone, author P.-C. Yeh, author
  J. Hone,  and author
  R. M. Osgood, 10.1103/PhysRevX.3.021014 journal journal Phys.
  Rev. X volume 3, pages 021014
  (year 2013)NoStop
[Gu et al.(2012)Gu,
  Petrone, McMillan, van der
  Zande, Yu, Lo, Kwong,
  Hone, and Wong]gu2012regenerative
  author author T. Gu, author N. Petrone,
  author J. F. McMillan, author A. van der Zande, author
  M. Yu, author G.-Q. Lo, author D.-L.  Kwong, author J. Hone,  and author C. W. Wong, @noop  journal journal Nature
  Photonics volume 6, pages 554
  (year 2012)NoStop
[Sun et al.(2010)Sun,
  Divin, Rioux, Sipe,
  Berger, de Heer, First, and  Norris]sun2010
  author author D. Sun, author C. Divin, author J. Rioux, author
  J. E. Sipe, author
  C. Berger, author W. A. de Heer, author P. N.  First,  and author T. B.  Norris, 10.1021/nl9040737 journal journal Nano Letters volume
  10, pages 1293 (year 2010), note
  pMID: 20210362, http://arxiv.org/abs/https://doi.org/10.1021/nl9040737
  https://doi.org/10.1021/nl9040737 NoStop
[Otsuji et al.(2012)Otsuji,
  Tombet, Satou, Fukidome,
  Suemitsu, Sano, Popov,
  Ryzhii, and Ryzhii]otsuji2012
  author author T. Otsuji, author S. A. B. Tombet, author A. Satou,
  author H. Fukidome, author M. Suemitsu, author
  E. Sano, author V. Popov, author M. Ryzhii,  and author V. Ryzhii, http://stacks.iop.org/0022-3727/45/i=30/a=303001 journal
  journal Journal of Physics D: Applied Physics volume 45, pages 303001 (year
  2012)NoStop
[Yao and Belyanin(2012)]yao2012
  author author X. Yao and author A. Belyanin,  10.1103/PhysRevLett.108.255503 journal
  journal Phys. Rev. Lett. volume 108,  pages 255503 (year 2012)NoStop
[Tokman et al.(2013)Tokman,
  Yao, and Belyanin]tokman2013
  author author M. Tokman, author X. Yao,  and  author A. Belyanin, 10.1103/PhysRevLett.110.077404 journal journal Phys. Rev. Lett. volume 110, pages 077404 (year 2013)NoStop
[Tokman et al.(2016)Tokman,
  Wang, Oladyshkin, Kutayiah,  and Belyanin]tokman2016
  author author M. Tokman, author Y. Wang,
  author I. Oladyshkin, author A. R. Kutayiah,  and author A. Belyanin, 10.1103/PhysRevB.93.235422 journal journal Phys.
  Rev. B volume 93, pages 235422
  (year 2016)NoStop
[Cheng et al.(2014a)Cheng, Vermeulen, and Sipe]cheng2016
  author author J. L. Cheng, author N. Vermeulen,   and author J. E. Sipe, http://stacks.iop.org/1367-2630/16/i=5/a=053014 journal
  journal New Journal of Physics volume
  16, pages 053014 (year
  2014a)NoStop
[Cheng et al.(2014b)Cheng, Vermeulen, and Sipe]cheng2014
  author author J. L. Cheng, author N. Vermeulen,   and author J. E. Sipe, 10.1364/OE.22.015868 journal journal
  Opt. Express volume 22, pages 15868
  (year 2014b)NoStop
[Cheng et al.(2017)Cheng,
  Vermeulen, and Sipe]cheng2017second
  author author J. Cheng, author N. Vermeulen,   and author J. Sipe, @noop
   journal journal Scientific reports volume 7, pages 43843 (year
  2017)NoStop
[An et al.(2013)An,
  Nelson, Lee, and Diebold]an2013
  author author Y. Q. An, author F. Nelson, author J. U. Lee,  and author A. C. Diebold, 10.1021/nl4004514 journal journal Nano Letters  volume 13, pages 2104 (year
  2013), note pMID: 23581964, http://arxiv.org/abs/https://doi.org/10.1021/nl4004514
  https://doi.org/10.1021/nl4004514 NoStop
[An et al.(2014)An,
  Rowe, Dougherty, Lee, and  Diebold]an2014
  author author Y. Q. An, author J. E. Rowe,
  author D. B. Dougherty, author J. U. Lee,  and author A. C. Diebold, 10.1103/PhysRevB.89.115310 journal journal Phys.
  Rev. B volume 89, pages 115310
  (year 2014)NoStop
[Avetissian et al.(2012)Avetissian, Avetissian, Mkrtchian, and  Sedrakian]avetissian2012
  author author H. K. Avetissian, author A. K. Avetissian, author G. F. Mkrtchian,  and author K. V. Sedrakian, @noop  journal journal
  Journal of Nanophotonics volume 6, pages 061702 (year 2012)NoStop
[Entin et al.(2010)Entin,
  Magarill, and Shepelyansky]entin2010
  author author M. V. Entin, author L. I. Magarill,
   and author D. L. Shepelyansky, 10.1103/PhysRevB.81.165441 journal journal Phys. Rev. B volume
  81, pages 165441 (year 2010)NoStop
[Glazov and Ganichev(2014)]glazov2014
  author author M. Glazov and author S. Ganichev, https://doi.org/10.1016/j.physrep.2013.10.003
  journal journal Physics Reports volume 535, pages 101  (year
  2014), note high frequency electric field induced nonlinear
  effects in grapheneNoStop
[Glazov(2011)]glazov2011
  author author M. M. Glazov, 10.1134/S0021364011070046 journal
  journal JETP Letters volume 93,  pages 366 (year 2011)NoStop
[Hendry et al.(2010)Hendry,
  Hale, Moger, Savchenko, and  Mikhailov]hendry2010
  author author E. Hendry, author P. J. Hale,
  author J. Moger, author A. K. Savchenko,  and author S. A. Mikhailov, 10.1103/PhysRevLett.105.097401 journal journal
  Phys. Rev. Lett. volume 105, pages
  097401 (year 2010)NoStop
[Mikhailov and Ziegler(2008)]mikhailov2008
  author author S. A. Mikhailov and author K. Ziegler, http://stacks.iop.org/0953-8984/20/i=38/a=384204
  journal journal Journal of Physics: Condensed
  Matter volume 20, pages 384204
  (year 2008)NoStop
[Mikhailov(2016)]mikhailov2016
  author author S. A. Mikhailov, 10.1103/PhysRevB.93.085403 journal journal Phys. Rev. B volume
  93, pages 085403 (year 2016)NoStop
[Mikhailov(2007)]mikhailov2007
  author author S. A. Mikhailov, http://stacks.iop.org/0295-5075/79/i=2/a=27002
  journal journal EPL (Europhysics Letters)  volume 79, pages 27002 (year
  2007)NoStop
[Mikhailov(2011)]mikhailov2011
  author author S. A. Mikhailov, 10.1103/PhysRevB.84.045432 journal journal Phys. Rev. B volume
  84, pages 045432 (year 2011)NoStop
[Mikhailov(2014)]mikhailov2014-2
  author author S. A. Mikhailov, 10.1103/PhysRevLett.113.027405 journal journal Phys. Rev. Lett. volume 113, pages 027405 (year
  2014)NoStop
[Yao et al.(2014)Yao,
  Tokman, and Belyanin]yao2014
  author author X. Yao, author M. Tokman,  and  author A. Belyanin, 10.1103/PhysRevLett.112.055501 journal journal Phys. Rev. Lett. volume 112, pages 055501 (year 2014)NoStop
[Constant et al.(2016)Constant, Hornett, Chang, and  Hendry]constant2016all
  author author T. J. Constant, author S. M. Hornett, author D. E. Chang,
   and author E. Hendry,  @noop  journal journal Nature Physics  volume 12, pages 124 (year
  2016)NoStop
[Wang et al.(2016)Wang,
  Tokman, and Belyanin]wang2016
  author author Y. Wang, author M. Tokman,  and  author A. Belyanin, 10.1103/PhysRevB.94.195442 journal journal Phys. Rev. B volume 94, pages 195442 (year 2016)NoStop
[Wang et al.(2015)Wang,
  Tokman, and Belyanin]wang2015
  author author Y. Wang, author M. Tokman,  and  author A. Belyanin, 10.1103/PhysRevA.91.033821 journal journal Phys. Rev. A volume 91, pages 033821 (year 2015)NoStop
[Smirnova et al.(2014)Smirnova, Shadrivov, Miroshnichenko,
  Smirnov, and Kivshar]smirnova2014
  author author D. A. Smirnova, author I. V. Shadrivov, author A. E. Miroshnichenko, author A. I. Smirnov,  and author Y. S. Kivshar, 10.1103/PhysRevB.90.035412 journal journal Phys. Rev. B volume
  90, pages 035412 (year 2014)NoStop
[Winzer et al.(2013)Winzer,
  Mali ćć, and Knorr]winzer2013
  author author T. Winzer, author E. Mali ćć,  and author
  A. Knorr, 10.1103/PhysRevB.87.165413 journal journal Phys.
  Rev. B volume 87, pages 165413
  (year 2013)NoStop
[Li et al.(2012)Li,
  Luo, Hupalo, Zhang,
  Tringides, Schmalian, and Wang]li2012
  author author T. Li, author L. Luo, author M. Hupalo, author
  J. Zhang, author M. C. Tringides, author J. Schmalian,  and author J. Wang, 10.1103/PhysRevLett.108.167401
  journal journal Phys. Rev. Lett. volume 108, pages 167401 (year
  2012)NoStop
[König-Otto et al.(2016)König-Otto, Mittendorff, Winzer,
  Kadi, Malic, Knorr,
  Berger, de Heer, Pashkin,
  Schneider, Helm, and Winnerl]konig-otto2016
  author author J. C. König-Otto, author M. Mittendorff, author T. Winzer,
  author F. Kadi, author
  E. Malic, author A. Knorr, author C. Berger, author W. A.  de Heer, author A. Pashkin, author H. Schneider,
  author M. Helm,  and author S. Winnerl, 10.1103/PhysRevLett.117.087401 journal journal
  Phys. Rev. Lett. volume 117, pages
  087401 (year 2016)NoStop
[Dean and van
  Driel(2009)]dean2009
  author author J. J. Dean and author H. M. van
  Driel, 10.1063/1.3275740 journal journal Applied Physics Letters volume 95,  pages 261910 (year 2009), http://arxiv.org/abs/https://doi.org/10.1063/1.3275740
  https://doi.org/10.1063/1.3275740 NoStop
[Dean and van
  Driel(2010)]dean2010
  author author J. J. Dean and author H. M. van
  Driel, 10.1103/PhysRevB.82.125411 journal
  journal Phys. Rev. B volume 82,  pages 125411 (year 2010)NoStop
[Lin et al.(2014)Lin,
  Weng, Lyu, Tsai, and  Su]lin2014
  author author K.-H. Lin, author S.-W. Weng,
  author P.-W. Lyu, author T.-R. Tsai,  and author W.-B. Su, 10.1063/1.4898065 journal journal Applied
  Physics Letters volume 105, pages
  151605 (year 2014), http://arxiv.org/abs/https://doi.org/10.1063/1.4898065
  https://doi.org/10.1063/1.4898065 NoStop
[Brun and Pedersen(2015)]brun2015
  author author S. J. Brun and author T. G. Pedersen, 10.1103/PhysRevB.91.205405 journal journal Phys. Rev. B volume
  91, pages 205405 (year 2015)NoStop
[Bykov et al.(2012)Bykov,
  Murzina, Rybin, and Obraztsova]bykov2012
  author author A. Y. Bykov, author T. V. Murzina,
  author M. G. Rybin,  and  author E. D. Obraztsova, 10.1103/PhysRevB.85.121413 journal journal Phys. Rev. B volume 85, pages 121413 (year 2012)NoStop
[Wu et al.(2012)Wu,
  Mao, Jones, Yao,
  Zhang, and Xu]wu2012
  author author S. Wu, author L. Mao, author A. M. Jones, author
  W. Yao, author C. Zhang,  and author X. Xu, 10.1021/nl300084j journal journal Nano Letters volume
  12, pages 2032 (year 2012), note
  pMID: 22369519, http://arxiv.org/abs/https://doi.org/10.1021/nl300084j
  https://doi.org/10.1021/nl300084j NoStop
[Tani et al.(2012)Tani,
  Blanchard, and Tanaka]tani2012
  author author S. Tani, author F. Blanchard,   and author K. Tanaka, 10.1103/PhysRevLett.109.166603 journal journal Phys. Rev. Lett. volume 109, pages 166603 (year 2012)NoStop
[Oladyshkin et al.(2017)Oladyshkin, Bodrov, Sergeev, Korytin, Tokman, and Stepanov]oladyshkin2017
  author author I. V. Oladyshkin, author S. B. Bodrov, author Y. A. Sergeev,
  author A. I. Korytin, author M. D. Tokman,  and author A. N. Stepanov, 10.1103/PhysRevB.96.155401 journal journal Phys.
  Rev. B volume 96, pages 155401
  (year 2017)NoStop
[Kocharovskaya et al.(1999)Kocharovskaya, Radeonychev, Mandel, and  Scully]kocharovskaya1999
  author author O. Kocharovskaya, author Y. V. Radeonychev, author P. Mandel,
   and author M. O. Scully,  10.1103/PhysRevA.60.3091 journal journal Phys. Rev. A volume 60, pages 3091 (year 1999)NoStop
[Vajna et al.(2015)Vajna,
  Dóra, and Moessner]vajna2015
  author author S. Vajna, author B. Dóra,   and author R. Moessner, 10.1103/PhysRevB.92.085122 journal journal Phys. Rev. B volume 92, pages 085122 (year 2015)NoStop
[Lewkowicz et al.(2011)Lewkowicz, Rosenstein, and Nghiem]lewkowicz2011
  author author M. Lewkowicz, author B. Rosenstein,  and author D. Nghiem, 10.1103/PhysRevB.84.115419 journal journal Phys. Rev. B volume
  84, pages 115419 (year 2011)NoStop
[Fülöp et al.(2010)Fülöp, Pálfalvi, Almási,  and Hebling]fulop2010
  author author J. A. Fülöp, author L. Pálfalvi, author G. Almási,  and author J. Hebling, 10.1364/OE.18.012311 journal
  journal Opt. Express volume 18,  pages 12311 (year 2010)NoStop
[Katsnelson and Kat︠s︡nelʹson(2012)]katsnelson2012graphene
  author author M. Katsnelson and author M. Kat︠s︡nelʹson, https://books.google.com/books?id=FgwW3aVRBT0C title
  Graphene: Carbon in Two Dimensions (publisher Cambridge
  University Press, year 2012)NoStop
[Tokman(2009)]tokman2009
  author author M. D. Tokman, 10.1103/PhysRevA.79.053415 journal journal Phys. Rev. A volume
  79, pages 053415 (year 2009)NoStop
[Zhang et al.(2014)Zhang,
  Arikawa, Kato, Reno,
  Pan, Watson, Manfra,
  Zudov, Tokman, Erukhimova,
  Belyanin, and Kono]zhang2014
  author author Q. Zhang, author T. Arikawa,
  author E. Kato, author
  J. L. Reno, author
  W. Pan, author J. D. Watson, author M. J.  Manfra, author M. A.  Zudov, author M. Tokman, author M. Erukhimova,
  author A. Belyanin,  and  author J. Kono, 10.1103/PhysRevLett.113.047601 journal journal
  Phys. Rev. Lett. volume 113, pages
  047601 (year 2014)NoStop
[Al-Naib et al.(2015)Al-Naib, Sipe, and Dignam]alnaib2015
  author author I. Al-Naib, author J. E. Sipe,   and author M. M. Dignam,  doi:10.1088/1367-2630/17/11/113018 journal
  journal New J. Phys. volume 17,  pages 113018 (year 2015)NoStop
[Malic et al.(2011)Malic,
  Winzer, Bobkin, and Knorr]malic2011
  author author E. Malic, author T. Winzer,
  author E. Bobkin,  and author A. Knorr, 10.1103/PhysRevB.84.205406 journal journal Phys.
  Rev. B volume 84, pages 205406
  (year 2011)NoStop
[Landau and Lifshitz(1984)]landau1984
  author author L. Landau and author E. Lifshitz, @noop  title Electrodynamics of
  continuous media, edition 2nd ed. (publisher
  Pergamon, Oxford, year 1984)NoStop
[Winzer et al.(2010)Winzer,
  Knorr, and Malic]winzer2010
  author author T. Winzer, author A. Knorr,   and author E. Malic, 10.1021/nl1024485 journal journal Nano
  Letters volume 10, pages 4839
  (year 2010), note pMID: 21053963, http://arxiv.org/abs/https://doi.org/10.1021/nl1024485
  https://doi.org/10.1021/nl1024485 NoStop
[Xing et al.(2010)Xing,
  Guo, Zhang, Sum, and  Huan]xing2010
  author author G. Xing, author H. Guo, author X. Zhang, author
  T. C. Sum,  and author
  C. H. A. Huan, 10.1364/OE.18.004564 journal journal Opt.
  Express volume 18, pages 4564
  (year 2010)NoStop
[Tan et al.(2017)Tan,
  Argondizzo, Wang, Cui, and  Petek]Tan2017
  author author S. Tan, author A. Argondizzo,
  author C. Wang, author
  X. Cui,  and author
  H. Petek, 10.1103/PhysRevX.7.011004 journal journal Phys.
  Rev. X volume 7, pages 011004
  (year 2017)NoStop
[Huang et al.(2018)Huang,
  Jiang, Zhang, Shan,
  Fan, Zhang, Dai,
  Shi, Liu, Zeng, Zi, Liu, and Wu]huang2018
  author author D. Huang, author T. Jiang,
  author Y. Zhang, author Y. Shan, author
  X. Fan, author Z. Zhang, author Y. Dai, author L. Shi, author K. Liu, author C. Zeng, author
  J. Zi, author W.-T. Liu,  and author S. Wu, 10.1021/acs.nanolett.8b03967
  journal journal Nano Letters volume 18, pages 7985 (year 2018), note pMID: 30451504, http://arxiv.org/abs/https://doi.org/10.1021/acs.nanolett.8b03967
  https://doi.org/10.1021/acs.nanolett.8b03967 NoStop







