Origins of All-Optical Generation of Plasmons in Graphene

Graphene, despite its centrosymmetric structure, is predicted to have a substantial second order nonlinearity, arising from non-local effects. However, there is disagreement between several published theories and experimental data. Here we derive an expression for the second order conductivity of graphene in the non-local regime using perturbation theory, concentrating on the difference frequency mixing process, and compare our results with those already published. We find a second-order conductivity (σ(2) ≈ 10−17 AmV−2) that is approximately three orders of magnitude less than that estimated from recent experimental results. This indicates that nonlinear optical coupling to plasmons in graphene cannot be described perturbatively through the electronic nonlinearity, as previously thought. We also show that this discrepancy cannot be attributed to the bulk optical nonlinearity of the substrate. As a possible alternative, we present a simple theoretical model of how a non-linearity can arise from photothermal effects, which generates a field at least two orders of magnitude larger than that found from perturbation theory.


Results and Discussion
perturbation theory. Formally, the interaction of an electron with a vector potential → A can be incorporated into a Hamiltonian via the substitution ˆ→ → → + → p p eA , where → p is the canonical momentum and e is the elemental charge. Absent the vector potential, an electron in a periodic crystal potential V r ( ) → , with Hamiltonian where m is the bare mass of an electron, can be formally diagonalized to produce a band structure. In the case of graphene, the Hamiltonian is typically taken to be of a tight-binding form. The Peierls substitution 15 formally enables one to incorporate the vector potential into such a postulated model, avoiding the need to actually solve for the eigenstates of the Hamiltonian with the replacement → → → + → p p eÂˆ. The tight-binding Hamiltonian of graphene with → A thus reads: operators for the two sublattices in graphene, with → R i , i = 1, …, N denoting the sublattice sites, and τ l , l = 1, 2, 3 denoting the vectors from a lattice site to its three nearest neighbors. The current density operator can then be obtained by ˆĵ R H A ( ) / i → → = ∂ ∂ → 16 .
In weak electromagnetic (EM) fields, both Ĥ and j→ can be expanded in terms of A → . Ĥ can be broken into a non-interacting part H 0 and an interacting part H Î , with = + +ˆˆ H H H I I I (1) ( 2) and the superscripts indicating the order of (1) ( 2) → = → + → + → + ⋅ ⋅ ⋅: Crucially, the Peierls substitution yields terms in Eqs (4), (6) and (7), which cannot be obtained by replacing → pˆ with → + → p eA in the Dirac Hamiltonian (as done in previous works 8,9,11 ). In fact, it is the term in Eq. (6) that cancels out a term in Eq. (5) that could otherwise cause a divergence in the linear conductivity 11  and solving the density matrix n n ( ) ρ ρ = ∑ˆ perturbatively in powers of Ĥ I . We leave the detailed derivations of the spinor formalism and the currents to the Methods, and only quote the result for the nonlinear current here. We consider the response to EM fields described by a potential www.nature.com/scientificreports www.nature.com/scientificreports/→ , where the electric field components parallel to the graphene layer are related to the potentials by E A t / mx m = −∂ ∂ ; q m and ω m are the wavevectors and angular frequencies, and x is a unit vector along the x-direction. Following relevant experiment 14 the fields for m = 1, 2 are called "pump" and "probe" respectively, and we illustrate their configurations in Fig. (1). Here we are interested in the case of DFG and look for the nonlinear current at difference frequency ω 3 = ω 2 − ω 1 and wavevector q 3 = q 2 − q 1 , which can be formally written as 1 2 , where σ (2) defines the second-order conductivity. Under the relevant experimental conditions 7,14 where in ω F and v F are the Fermi angular frequency and velocity.
Comparison of perturbation theory with experiment. X. Yao 6 were the first to derive the nonlinear conductivity in graphene relevant to DFG. In Fig. 2, we plot this derived second order nonlinear response (Eq. (5) of X. Yao 6 converted to σ (2) ) along with experimental results from Constant 14 , recent perturbative calculations 11 and our own Eq. (8). It is important to note first that the results of Wang 11 and Cheng 9 strongly agree with our   www.nature.com/scientificreports www.nature.com/scientificreports/ own suggesting this is the correct prediction from pertubation theory (in addition Tokman 8 and Rostami 17 achieve the same result except a factor of 2 which may be due to definitions we could not clearly identify). However there is a rather large discrepancy between Eq. (8) and the model derived by X. Yao 6 , which were derived for identical conditions using perturbation theory. Moreover the conductivity derived by X. Yao 6 has a non-physical divergence for → q 0 3 . While it is not clear from where this unphysical behavior arises, in a centrosymmetric material such as graphene this behavior is paradoxical. Meanwhile, the conductivity from Eq. (8) tends to zero as q 0 3 → , as it must in graphene. We note that, depending on the value of q 3 in Fig. 2, the magnitude of σ (2) predicted by Eq. (8) is at least 4 orders lower than that found by X. Yao 6 .
The first experimental signatures attributed to DFG of plasmons were found by Constant 14 . In this experiment, by illuminating the graphene with two tunable, femtosecond laser pulses ("pump" and "probe") with well-defined angles of incidence but different frequencies, Constant 14 was able to phase-match to the plasmon. The geometry of the experiment is the same as that chosen for our theoretical calculation defined in Fig. 1. The graphene supports tightly guided plasmons with a dispersion relation ω pl (k). The differential reflectance of the probe pulse ΔR was seen to change significantly whenever the difference frequency and wavevector were aligned to the plasmon dispersion relation, suggesting efficient plasmon excitation via DFG. In practice, a range of difference frequencies and wavevectors were scanned by continuously varying the pump wavelength, and by choosing different discrete incident beam angles.
Constant investigated three experimental geometries (noted in Table 1) with different angles of incidence, θ. We examine one of the resonant conditions for each of the three experimental geometries, as defined in Table 1. Assuming the differential reflection signals arise from DFG, one can use the model introduced in the Supplementary Information of Constant 14 (briefly reviewed in the Methods section) to estimate a value for σ (2) for each measurement. The results of this analysis, i.e. values of σ (2) which describe the experimental signals, are also shown in Table 1. Figure 2 compares the experimental values of σ (2) from Constant 14 and theoretical predictions from Eq. (8), Wang 11 and X. Yao 6 . Firstly, the q 3 dependence of Yao 6 clearly differs greatly from that of both the experiment and the near-linear predictions of other theoretical derivations. The experimental magnitudes of σ (2) are also significantly lower than the prediction of X. Yao 6 , and several orders higher than those from recent perturbative works 8,9,11 . As found by 17 , it is only possible to find agreement between the experiments 14 and perturbative second order calculations if one invokes an unphysically low decay rate for the plasmon (resulting in extraordinarily narrow resonances). More recently, a similar experiment has been carried out by B. Yao 7 in a waveguiding geometry, and the theory from X. Yao 6 was used to model the experimental signals. While the geometries of B. Yao 7 and Constant 14 are significantly different, similar signals were observed in each experiment. Therefore, ignoring the unphysical results in X. Yao 6 , the large discrepancy between both experiments and the theoretical consensus points to a second order response that is not purely perturbative, as originally interpreted. In the remainder of this paper, we therefore discuss other possible contributions which might account for the discrepancy between perturbation theory and the experiments of B. Yao 7 and Constant 14 . substrate Response. In this section we consider contribution of the second order nonlinearity of the quartz substrate used in experiment 14 . The analysis significantly simplifies if the nonlinear polarization is generated far from a phase-matching condition of the bulk, and depletion can be ignored, as should be the situation for Constant 14 . In this case, the pump and probe fields generate a polarization in quartz are given by where χ (2) is the second order susceptibility of the substrate, k Ti → , t i , and E Ii denote the wavevector on the transmitted (substrate) side, transmission coefficient, and the incident field amplitude of the pump (i = 1) and the probe (i = 2) fields respectively. The transmission coefficients t i are given in the Methods. The subscript i = 3 indicates quantities corresponding to the difference frequency signal at ω 3 = ω 1 − ω 2 . As charge density waves in graphene are driven by an electric field, we must relate the nonlinear polarization to the field generated in the quartz, which satisfies the wave equation Here, the subscript "s" denotes that this is an effective source field that will later drive a response in the graphene (distinct from the resulting plasmon field). Also, n( ) indicates the permittivity of quartz evaluated at the difference frequency, with the model of the frequency dependent n(ω) given in the Appendix. Due to the plane-wave nature of P 3 , E 3s takes on the same spatial and frequency dependence. In our regime of interest, the spatial derivative of the field, , is significantly larger than the time derivative. This is because the pump and probe fields are chosen to phase-match with surface plasmons in graphene (thus the associated wavevectors are much larger than free-space fields at the difference frequency). Thus the field amplitude created by the nonlinear polarization is well approximated by s T T In particular, it should be noted that a large wavevector mismatch results in strong suppression of the field. To simplify the discussion, we will assume the scenario which produces the highest field, i.e. in which E 3s is www.nature.com/scientificreports www.nature.com/scientificreports/ completely polarized along x (parallel to the graphene sheet) so that it maximally drives a charge density wave in graphene. As we see below, even in this best case scenario, the generated field is rather small.
Since the nonlinear response is considered here to be completely within the substrate, which provides an effective source field E 3s , the remaining part of the calculation is completely linear in its nature. Using the conventions in Fig. 1, we take "reflected" and "transmitted" field components of unknown amplitude, which correspond to the plasmon field on the vacuum and substrate sides. The wavevector along x for these fields is equal to q 3 = q T1 − q T2 , where q T1 and q T2 are the in-plane components of → k T1 and → k T2 , while the perpendicular components must satisfy the respective dispersion relations for each side of the interface, e.g., 2 . Similar to the procedures to solve the pump (probe) field laid out in the methods, the two unknown field amplitudes can be readily solved by taking E 3s to be the incident field on the substrate side, and enforcing electromagnetic boundary conditions at the vacuum-graphene-quartz interface, which yields the following parallel-field component on the substrate side, evaluated at the graphene layer (z = 0), Here, σ (1) (ω 3 ) is the linear conductivity of graphene evaluated at frequency ω 3 . Specifically we can numerically evaluate E pl for geometry (b) in 14 . Taking a value of 0 3 pmV (2) 1 χ = . − for quartz 18  . The modeling in 14 predicts a considerably larger value for the inferred plasmon field in experiment of 8  ≈ × − . We therefore do not believe that the substrate nonlinearity contributes significantly to the signals observed in 14 . However, we note that our substrate model does not consider any surface enhanced nonlinearity. Its theoretical modeling would require experimental measurements of surface nonlinear coefficients relevant to our system, which we were unable to find in literature.

Photothermal Effect.
Here, we present an alternative mechanism by which a plasmon field at the difference frequency and wavevector can be generated. Fundamentally, the effect discussed below arises from the linear Seebeck effect, so that the total current is described by The first term on the right hand side describes the normal linear relationship between the current and field, while the second describes the conventional second order electronic nonlinearity. The third term, and most important here, arises due to photothermal effects, and accounts for the Seebeck current emerging due to a temperature gradient ( ) dT dx in a material described by Seebeck coefficient S. As we discuss below, this term can give rise difference frequency currents even in the absence of a nonlinear conductivity (i.e. even when σ (2) = 0).
It is known that excitation of graphene carriers by intense femtosecond pulses (Φ ≈ . − 0 1 mJcm 2 with pulse width ~100 fs in 14 ) is not perturbative in nature. The electron temperature is raised by several thousand kelvin under such excitation 19 and is not in equilibrium with the phonon temperature. Furthermore, under geometries similar to those used in B. Yao 7 and Constant 14 , heating due to optical illumination is not homogeneous. When two or more light sources of similar frequency are incident on an interface at oblique angles, the result will be a near stationary interference pattern such as that shown in Fig. 3. When the frequencies are slightly different, the pattern will propagate in the plane with a velocity equal to ω 3 /q 3 . The spatially dependent interference pattern will give rise to a temperature modulation that propagates along the graphene sheet, which, via the Seebeck effect, can also generate a current and thereby drive a plasmon excitation when phase-matching conditions are satisfied. www.nature.com/scientificreports www.nature.com/scientificreports/ Below, we present a theoretical model of this effect, and show that it can generate a difference frequency field that is orders of magnitude larger than that calculated from the perturbation theory.
The intensity pattern imprinted on the graphene sheet can be expressed as Here E 1x and E 2x are the in-plane components of the pump and probe fields respectively, and we have assumed the probe field is much weaker than the pump field. The ratio of the in-plane components can be calculated as 2 1 , with I 1 and I 2 being the incident intensities of the pump and probe beams respectively, and the linear transmission is given by Eq. (38) considering only the linear optical conductivity of graphene Eq. (40). Using the parameters in the experiment of Constant 14 , we obtain ≈ .

E E
/ 01 . The intensity pattern acts as a heat source for the temperature distribution which, in linear response theory, satisfies a diffusion equation: where α is the diffusivity, β is the heating rate due to the intensity pattern, and y is the relaxation rate back to the equilibrium temperature T 0 . Due to the linearity of the equation, it can be readily solved in the Fourier domain, in which when taking into account Eq. (14) the solution takes the form Here, T dc is the (large) position-and time-independent temperature increase arising from the incident lasers, while T ac represents a position-and time-varying temperature oscillation that must necessarily be generated in the presence of moving intensity interference pattern. Ψ denotes a phase offset between the intensity and temperature modulations, whose specific form is not relevant here. Substituting T(x, t) into Eq. (15), in the regime of interest y ≪ ω 3 , one obtains |T ac |/T dc ≈ 0.2y/ω 3 . As expected, the temperature modulation T ac is reduced significantly as the oscillation frequency ω 3 increases with respect to the damping rate. It is known that intense, femtosecond pulses similar to those in 14 lead to T dc of approximately 2000 K 19 . The relaxation rate y is due to electron-phonon scattering, and we take a value of y ≈ 1/(100 fs) 19 . At ω 3 = 2π × 10 THz these parameters give a temperature modulation T ac = 60 K.
The Seebeck effect enables the generation of a source current in the presence of a temperature gradient, and this can be described by j q qS T ( , ) ( ) , where S is called the Seebeck coefficient. In principle, the Seebeck coefficient could be frequency and wavevector dependent. However, this dependence has not been measured carefully in literature, nor is it straightforward to calculate from first principles. There have been several measurements of the Seebeck effect in graphene, both in DC experiments (S ≈ 5 × 10 −5 V/K 20 , S ≈ 8 × 10 −5 V/K 21 ) and under illumination from 100 fs pulses (S ≈ 10 −4 V/K 22 ). Whilst it is hard to predict how the Seebeck effect behaves on 10 fs timescales relevant here (corresponding to peak to peak propagation time of the intensity pattern in Fig. 3), it is likely that photothermal effects will be higher on ballistic timescales, as with other materials 23 . Here we use a conservative value of S ≈ 10 −4 V/K reported in 22 . Now using the standard EM boundary conditions at the graphene layer (see Fig. 1), with the aid of the charge continuity equation, one can find the relation between the electric fields and the surface current density at the difference frequency ω 3 . Note now the Seebeck effect contribution needs to be added to the surface current density: Then solving the equations of the boundary conditions (see Methods for more details), we obtain for the electric field at the difference frequency: ( ) is the linear transmission coefficient at frequency and wavevector ω 3 , q 3 . For geometry (b) of Constant 14 we find a magnitude of E 3x , when on plasmon resonances, of E 2 3 10 pl 3 ≈ . × V/m. Just as the pump and probe fields can generate a plasmon field through the Seebeck effect, a back-action effect (involving Seebeck mixing of the plasmon and pump fields) results in a change of the probe differential reflectance. In principle, this could be rigorously calculated in a manner similar to above, but this would require knowledge of the Seebeck coefficient at optical frequencies, which has never been measured or calculated. However, we can nonetheless obtain an approximate value for the differential reflection, by exploiting conservation of energy. In particular, in steady state, the number of plasmons dissipated per unit time must equal the rate of photons removed from (added to) the pump (probe) beam. There are two contributions to the energy dissipation at the difference frequency: both the graphene layer and the substrate will exhibit absorption. For graphene, the power loss per unit area from absorption can be found by In the experiment, the substrate www.nature.com/scientificreports www.nature.com/scientificreports/ itself (quartz) can provide a non-negligible loss, through coupling with phonons. The corresponding power loss per unit area can be calculated as 0  3  2  3  3 2 , where our model for the frequency-dependent refractive index n(ω) is provided in the Methods. Then, the number of photons, at the difference frequency, absorbed per unit time and area is ω At the individual photon level of a DFG process, an incoming pump photon breaks down to an outgoing probe photon and a plasmon, and therefore the number of plasmons created is also equal to the number of newly generated probe photons that enter either the reflected or transmitted beam. The number of photons per unit time and area in the incident probe beam is sim-  . Thus the order of magnitude of the differential reflectance of the probe beam can be estimated as Δ ∼Γ Γ R Rn / / d in . For configuration (b) in Table 1, we estimate the peak differential reflectance after normalized by the fluence (0.1 mJ cm −2 ) to be ~7 × 10 −7 mJ −1 cm 2 . Although we emphasize that the Seebeck effect and the electronic nonlinearity σ (2) are completely independent effects, nonetheless to facilitate a better comparison, one can ask what hypothetical value of nonlinear conductivity S (2) σ would be required, in order to produce the same current as predicted from the Seebeck effect, i.e.
We extract a value of 3 1 10 AmV While this model predicts a value slightly smaller than experiment ( 7 5 10 AmV , the Seebeck coefficient could be larger on ballistic timescales relevant here (≈10 fs peak to peak propagation time) as is expected for other materials 23 . Nevertheless such an effect is fundamental to the experiments and is significantly larger than the predictions of perturbation theory.
We note that photothermal effects will be prominent in the waveguiding geometry of B. Yao 7 (in such a geometry, even though the absolute field intensities are lower, the considerably larger propagation length can compensate). Interestingly, the power dependence of such a photothermal signal would not necessarily follow that of conventional difference frequency generation (j E E (2) 1 2 ⁎ σ = ), and could explain those observed by Constant 25 . Investigating the intensity dependences of these nonlinear signals could provide great insight into the origins of these effects. We also note that the photo-Dember effect can similarly induce local intensity dependent currents and is surprisingly large in graphene on ultrafast timescales 26 . However, since the photo-Dember effect depends on mobility asymmetry between electrons and holes, it will be sample and substrate specific, making it difficult to estimate.

Conclusions
We have derived a second order conductivity of planar graphene (σ ≈ − − 10 AmV (2) 1 7 2 ) with non-local perturbation theory, addressing the long wavelength divergence in 6,7 and divergent linear current in 11 . However, while our result is in agreement with recent calculations 8,9,11,17 , it is insufficient to explain observations from experiment 7,14 .
We also show that this discrepancy cannot be attributed to the bulk nonlinearity of the substrate.
We also discuss the possibility of photothermal effects in experiments, wherein a spatial intensity pattern resulting from interference of incident beams leads to local inhomogeneous heating of the sample and show that these effects will give rise to frequency mixing currents. We derive a rigorous model for DFG arising from photothermal effects (with the only uncertainties arising from knowledge of material properties such as the Seebeck coefficient), and conservatively estimate a DFG current which is at least two orders of magnitude larger than that found from perturbation theory, significantly closer to experimental estimates from 14 . Microscopic modeling of such local photothermal effects (and other non-equilibrium processes) presents a considerable challenge, and it would be interesting to develop theoretical techniques to do so. We believe that such efforts would shed further light on discrepancies between recent experiments 7,14 and theory [6][7][8][9][10][11] for all-optical plasmon generation processes in graphene, and enable the strengths of these nonlinear processes to be optimized for future nonlinear optical applications.

Methods spinor formalism. We can perform a Fourier expansion on the operators
l in the tight-binding Hamiltonian in terms of the operators in the reciprocal lattice space: , and 1 , n which Ω B denotes the first Brillouin zone, k → is the electron momentum and N is the number of sites in one sublattice. One can then substitute these expansions into the expressions of the Hamiltonian Eqs (2-4) and current density operators Eqs (5-7). As usual, near the two Dirac points K → and K − → 27 , the operators can be expanded One can then derive equivalent spinor forms for Eqs (2)(3)(4)(5)(6)(7) in the first quantization picture. If only the terms to lowest order of → k are kept, then , with a being the lattice constant of the underlying sublattices, and σ σ σ → ≡ + x y x ŷˆˆˆˆ with x y , σ representing the Pauli spin matrices. Meanwhile, for Ĥ I : (1) www.nature.com/scientificreports www.nature.com/scientificreports/ with "+" at → K and "−" at K − → , and we have for simplicity assumed → A is along the x-axis. The second order component is meanwhile given by where q m is the in-plane component of the momentum and only kept the terms that give rise to a perturbation at ω = ω 1 − ω 2 and q = q 1 − q 2 (DFG). For H 0 we find that the single-particle eigenenergies E sv k Here s = ±1 is the band index, ξ is the area of the graphene sheet, and θ → k is the polar angle of k → . Similarly, the spinor forms of the current densities can be found aŝˆ⁎ ( , ) at , and 32 " "at and " "at (25) x It should be noted that in addition to the "typical" contribution to the current density, Eq. (24), which is associated with the Bloch momentum, there is also the "diamagnetic" term of ĵ x (1) in Eq. (25). We will see in the following that this term, and in fact a higher correction to this term of order k → , cancels the term that could otherwise cause a divergence in the linear conductivity.
Density matrix method. The time evolution of the density matrix under the interaction with EM fields described by H Î is given by ˆˆˆρ . From perturbation theory, the matrix elements n m nm ρ ρ ≡ˆ of the i th order perturbation of ρˆ satisfy , and γ is a phenomenological dissipation term introduced universally for all matrix elements. We take the density matrix in absence of fields to be the Fermi distribution f n with Fermi energy ω F  at zero temperature, , where Θ is the Heaviside step function. A solution to Eq. (26) can be written as inear conductivity. In previous works calculating graphene conductivities using the vector potential 8,9,11 , the authors typically replace p → by p eA → + → in the Dirac Hamiltonian. However, the linear current thus calculated has a term that diverges when the integration limit of the electronic momenta is taken to be infinity. This issue was fixed in 11 by adding an artificial quadratic term to the Dirac Hamiltonian. In this section we show that this problematic term is actually canceled by a term in Eq. (7) in the Results section, thus no artificial term needs to be introduced to regularize the calculation.
To begin with we consider the current response to an in-plane electric field described by a vector potential where q qx → =ˆ. The current generated at → q and ω is calculated through the expectation value of ˆĵ q e j r ( ) 2/ ( ), The matrix elements of the current density operators can be obtained by using Eqs (23)(24)(25). Now we can substitute these results into Eq. (28). We then replace the summation on states by an integral over Bloch momenta k → , www.nature.com/scientificreports www.nature.com/scientificreports/ introducing an upper bound on the range of integration < k k c (which approximately captures the edge of the Brillouin zone). The first term of Eq. (28)  where k F is the Fermi wavevector, and the second term gives zero. If k c is extended to infinity like in a free-electron gas, Eq. (30) will yield a divergent linear current, as discussed by Wang 11 . We show next how this strong dependence on k c is cancelled by the lowest order non-zero contribution to the second term in Eq. (28).
We expand the linear current density operator Eq. (6) to first order of k → , and obtain an additional term to ĵ x (1) in Eq. (25) which we label as ˆ′ j x (1) : x F x x y y (1) 2 2 2ˆˆˆŵ here "+" or "−" sign corresponds to the Dirac point K → or − → K . Now the second term in Eq. (28) has an additional term which gives a finite contribution: . The summation can be transformed to an integral, which in general needs to be evaluated numerically. However we can also expand the kernel in terms of q 3 , and extract the leading order contributions. Under the experimental conditions of both 7  (2)