# Analysis of All-Optical Generation of Graphene Surface Plasmons by a Frequency-Difference Process
## 2 Materials and Methods
### Second-Order Conductivities  
We work within the Dirac cone approximation and use the density matrix formalism to calculate the second-order nonlinear conductivity, as it was done in the previous works (the work [22] used a tight-binding Hamiltonian, which is equivalent to the Dirac-cone approximation in the vicinity of the Dirac point) [8; 23]. The Hamiltonian of the electrons coupled to the electromagnetic field is taken in the so-called velocity (or **A**) gauge [24; 25]:  
\[\hat{H}=v_{F}\boldsymbol{\sigma}\cdot\left(\mathbf{p}+\frac{e}{c}\mathbf{A}(t,\mathbf{r})\right), \tag{1}\]  
where \(v_{F}\) is the Fermi velocity (a material parameter), \(\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y})\), \(\sigma_{i}\) is the \(i\)th Pauli matrix, \(\mathbf{p}\) is the usual momentum operator, and \(\mathbf{A}(t,\mathbf{r})\) is the vector potential (the electric field is given by \(\mathbf{E}=-\frac{1}{c}(\partial\mathbf{A}/\partial t)\). Using the density-matrix perturbation theory [18] with the second term in (1), we find a general expression for the DF conductivity (see SI for details):  
\[ \begin{split}\sigma^{kj\left(\boldsymbol{\Omega},\boldsymbol{ \Omega}_{p},\boldsymbol{\Omega}_{q}\right)}&=\frac{e^{3}}{\hbar^{2} \omega_{q}\omega_{p}\mathcal{S}}\sum_{nml}\frac{v_{mn,-(\mathbf{q}+\mathbf{p} )}^{k}v_{lm,\mathbf{p}}^{j}v_{lm,\mathbf{q}}^{i}}{\omega_{lm}-\omega_{q}-i \gamma_{lm}}\frac{f(\varepsilon_{m})-f(\varepsilon_{l})}{\omega_{nm}-\left( \omega_{p}+\omega_{q}\right)-i\gamma_{nm}}\\ &+\frac{e^{3}}{\hbar^{2}\omega_{q}\omega_{p}\mathcal{S}}\sum_{ nml}\frac{v_{mn,-(\mathbf{q}+\mathbf{p})}^{k}v_{ml,\mathbf{p}}^{j}v_{lm, \mathbf{q}}^{i}}{\omega_{lm}-\omega_{q}-i\gamma_{lm}}\frac{f(\varepsilon_{m})- f(\varepsilon_{l})}{\omega_{nm}-\left(\omega_{p}+\omega_{q}\right)-i\gamma_{nm}}\\ &-\frac{e^{3}}{\hbar^{2}\omega_{q}\omega_{p}\mathcal{S}}\sum_{ nml}\frac{v_{mn,-(\mathbf{q}+\mathbf{p})}^{k}v_{ml,\mathbf{p}}^{j}v_{lm, \mathbf{q}}^{i}}{\omega_{nl}-\omega_{p}-i\gamma_{nl}}\frac{f(\varepsilon_{l})- f(\varepsilon_{n})}{\omega_{nm}-\left(\omega_{p}+\omega_{q}\right)-i\gamma_{nm}}\\ &-\frac{e^{3}}{\hbar^{2}\omega_{q}\omega_{p}\mathcal{S}}\sum_{ nml}\frac{v_{mn,-(\mathbf{q}+\mathbf{p})}^{l}v_{lm,\mathbf{p}}^{j}v_{ml, \mathbf{q}}^{i}}{\omega_{nl}-\omega_{q}-i\gamma_{nl}}\frac{f(\varepsilon_{l})- f(\varepsilon_{n})}{\omega_{nm}-\left(\omega_{p}+\omega_{q}\right)-i\gamma_{nm}} \end{split}  \tag{2}\]  
Figure 1: Schematics of all-optical generation of surface plasmons via a frequency-difference NL process of interaction of two optical beams, \(a\) and \(b\), represented by their electric fields \(\mathbf{E}_{a,b}\), including the incident and reflected ones. \(\mathbf{E}_{\text{SP}}\) represent the surface plasmon field. where \(\mathbf{\Omega}_{p}=(\omega_{p},\mathbf{p})\), \(\mathbf{\Omega}_{q}=(\omega_{q},\mathbf{q})\), and \(\mathbf{\Omega}=(\omega_{p}+\omega_{q},\mathbf{p}+\mathbf{q})\), with \(\omega_{p}\) and \(\omega_{q}\) denoting the frequencies of the fields and \(\mathbf{p}\) and \(\mathbf{q}\) the respective in-plane wavevectors. The sums are over different states near the Dirac point and  
\[v^{j}_{nm,\mathbf{q}}=\langle m|v^{j}e^{i\mathbf{r}\cdot\mathbf{q}}|n\rangle \tag{3}\]  
are the matrix elements of the spatially modulated velocity components. Here \(\omega_{mn}=(\varepsilon_{n}-\varepsilon_{n})/\hbar\) is the energy difference between states \(|m\rangle\) and \(|n\rangle\) and \(f(\varepsilon_{n})\) is the occupation number of state \(|n\rangle\), at equilibrium described by the Fermi-Dirac function. Finally, \(\gamma_{nm}\) are phenomenological damping parameters for transitions between states \(|n\rangle\) and \(|m\rangle\).