{'fraction_non_alphanumeric': 0.044549608355091384, 'fraction_numerical': 0.01864392950391645, 'mean_word_length': 4.038643371017471, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Second harmonic generation in graphene dressed by a strong terahertz field

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.


I. INTRODUCTION
There has been a surge of interest in the nonlinear and quantum optics of graphene .
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 [33,34]. A much stronger second-order nonlinearity originates from nonlocal response beyond the electric dipole approximation, i.e. magnetic-dipole and electric-quadrupole response [9,12,17,18,23,[25][26][27]29]. In [35] 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 [36].
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 [11,13,14,17,37,38]. 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 The experimental results that we obtained below apparently agree with the scaling of Eq. (1). 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. [39,40]. 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 [40]. In such strong fields, a standard expression for the current J Dr using the Drude relaxation time of ∼ 1 ps adopted in [11] would yield a current amplitude much higher than the maximum possible value ev 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 [39,40].
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. [41] 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 , to the inverse time of Schwinger pair creation in graphene by the field of magnitude E 0 [40,42,43]. 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 ω = 2kv 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 /eE 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.

II. EXPERIMENT
The schematic of the experimental setup is shown in Fig. 1. 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 [44] 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.

Grating
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 [39], 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. 2 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. 3 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 2 0 . Since the number of SH photons is proportional to the current j 2ω squared, this corresponds to the scaling Eq. (1). Inset to Fig. 3 shows the dependence of N 2ω from the energy of a laser pulse, which also agrees with Eq. (1). 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.

A. 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 [12,17,27]: 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. (3) is given by . 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 (2) has been proven in [27]. 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 [27]. Furthermore, according to [46,47], 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. (2).
In the absence of relaxation (which we will take into account later) the Von Neumann equation in the basis of Eq. (3) iṡ where 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. (4) and (5) yield a closed set of equations for matrix elements ρ ss kk . Denoting the quantum coherence as where ω k = 2v F k is the interband transition frequency, Ω k = e k (sin θ k E x − cos θ k E y ) the Rabi frequency. Equation (8) is separated from the rest of the system and we won't use it anymore.
Eqs. (6) and (7) make a closed set of equations, which is a version of semiconductor Bloch equations [11,48]. Note that the "convective" terms ∼ ∂ ∂k in Eqs. (6)-(8) 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. (5) 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. [49]. This would not change the observables of course.
The current operator is defined byĴ = −(ie/ ) Ĥ , 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 where g = 4 is spin and valley degeneracy. The first term in square brackets in Eqs. (9), (10) is the intraband current and the second term is an interband current.

B. Equations for slow variables
Now we make the following ansatz which separates the slow THz dynamics (subscript 0) from fast optical frequencies: 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. (6), (7). This leads to the replacement ω k =⇒ ω k − iγ in Eq. (6) and the population (7). 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.
In Eqs. (11)-(15) 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 In numerical modeling the solution to Eqs. The RWA corresponds to the following inequalities: In this case one can neglect counter-rotating terms ρ (+) ωk e +iωt and ρ (+) 2ωk e +2iωt in terms of Eqs. (11) to (18). Furthermore, if the following conditions are satisfied, One can assume that |∆ 0k | |∆ ωk |, |ρ 0k | ρ (−) ωk in Eqs. (11) to (18) and obtain approximate equations, where ρ k ≈ ρ 0k + ρ ωk e −iωt + ρ 2ωk e −2iωt . where The functions Ω ωk (t , t), ∆ 0k (t , t) and ω k (τ, t) are obtained from functions Ω ωk (t), ∆ 0k (t) and It is easy to check that Eqs. (23) and (24) Expanding Ψ k (t , t) near the stationary point t = t s (t, k), we obtain where Ψ k (t s , t) = Ψ k (t , t)| t =ts , α k (t s , t) = 1 2 The interval in Eq. (26) which makes the main contribution to the integral is given by The following hierarchy of timescales ensures the validity of the stationary phase method: where ∆t opt and ∆t T Hz are durations of the optical and THz pulses. Furthermore, we assume that ∆t opt ∆t T Hz , δt s ω 1, Eq. (26) can be written in the form (see Appendix A) where The top and bottom signs in Eq. (30) correspond to cos θ k E 0 > 0 and cos θ k E 0 < 0 respectively.
One can see from the expression in Eq. (30) 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. (30) is asymmetric, namely it is shifted towards ω > ω k for cos θ k E 0 > 0 and towards ω < ω k for cos θ k E 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. (29) 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. (25). The shift amount is t − t s ≈ ∓ Z δksv 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 = ∆ F k in the resonant region be small enough. The perturbation of the population difference can be estimated using Eqs. (21), (29). The analysis in Appendix A yields Equation (31) 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 s ω 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 ∆ F k = 1, which is the case for high enough electron energies corresponding to the interband resonance ω k ≈ ω.

C. The nonlinear current
For the electric field polarized along x, Eqs. (9) and (22) give the following expression for the complex amplitude of the current at SH: Substituting Eq. (29) into Eq. (32) and assuming ∆ 0k (t s , t) ≈ 1 and Ω ωk (t s , t) ≈ Ω ωk (t) leads to where M (Z, δk s ) is defined in Eq. (30). Limiting ourselves to the resonant region, we can . . when evaluating the inegral in Eq. (33) (here k res = ω 2v F ) and change from ∞ 0 dk . . . to infinite integration limits when integrating over dZ = dk δks : The first two terms in the expansion in Eq. (34) give zero after angle integration π −π dθ k . . . . For the first term this is obvious due to the normalization condition ∞ −∞ M (Z, δk s ) δk s dZ = 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 M is a linear combination of the even function e ±iZ 2 and odd function e ±iZ 2 Z 0 e ∓ix 2 dx, see Eq. (30). As a result, we obtain Recalling that the upper and lower signs correspond to cos θ k E 0 > 0 and cos θ k E 0 < 0, we get

D. Polarization dependence
Eq. (35) 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 k eEω . It is straightforward to obtain that the relationship between the SH current j x2ω and Ω ωk has the same form as in Eq. (32).
Therefore, instead of Eq. (35) we obtain Comparing Eqs. (35) with (36) 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.

V. 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. (35), (36) contain the components of the fields tangential to the monolayer. Using Fresnel formulas [50], one can get 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 Ppolarized field radiated by this current are related by 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.
Next we calculate the number of photons in the PSS polarization configuration. From Eqs. (36) − (38) for the above parameters ε ω , ε 2ω , ε 0 and θ we can estimate the field of the SH signal as 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 where For the same parameters we get α According to the first inequality in Eq. (28a) the optical pulse is longer than the integration interval in Eq. (26). Therefore the lower limit of integration can be taken as −∞; therefore, Eq. (26) yields where the upper and lower signs in e ±iπ/4 correspond to α k > 0 and α k < 0 respectively.
Furthermore, from Eqs. (25), (24) one can get Since the THz pulse is much longer than the optical pulse, one can write the integral where E 0 can be treated as a constant during the optical pulse.
The calculations are greatly simplified if we assume we will check its validity later.
Taking into account that the integration interval which makes the main contribution near the resonance is |ω − ω k | δt s ∼ 1 and using Eq. (A5), we obtain which means that the last inequality in Eq. (28b) ensures the validity of the approximation Eq. (A4). In the region |ω − ω k | δt s ∼ 1 the factor exp − γ cannot be greater than exp (− γδt s ), so if the last inequality in Eq. (28a) is satisfied, it can be taken as 1. As a result, we obtain Eqs. (29), (30). We also provide here useful asymptotics of the function W (Z, δk s ) in Eq. (30) at |Z| 1: where the upper and lower signs e ±i(Z 2 − π 4 ) are taken for cos θ k E 0 > 0 and cos θ k E 0 < 0 respectively.
Appendix B: 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 Eq. (A4) and Eq. (A5) as Taking into account that ω−ω k 2v F cos θ k eE 0 ∼ δt s , and δt s ∆t opt , we can drop the terms with time derivative ∂ ∂t in Eq. (B1) and Eq. (B2). The factor ω−ω k kv F ∼ δks k is small too. At the same time the term eE 0 is small only if the perturbation of the population difference in the resonant region is small too. Therefore, we obtain In the last expression upper and lower signs are given by the signs of cos θ k E 0 , as usual.
Therefore we get Appendix E: 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 we get If E ω x 0 , E 0 y 0 and j 2ω x 0 , similar considerations lead to dθ k sin 3 θ k cos θ k = 0.

Appendix F: Numerical simulations
We used Eqs. (11) to (15) 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 [51][52][53] 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 [39,[51][52][53][54]. (Auger recombination could be enhanced due to lattice imperfections in CVD graphene [39].) 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. 4 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 eE 0 ∂ ∂k terms included. The latter would be similar to four-wave mixing in a two-level medium. Figure 5 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.