# Analysis of All-Optical Generation of Graphene Surface Plasmons by a Frequency-Difference Process
## 2 Materials and Methods
### Second-Order Conductivities
The relevant transitions for DF generation are shown in Figure 2, which represent two channels [8]: Channel 1 with transitions between the states  
\[|1\rangle=|-1,\mathbf{k}_{1}\rangle\ \,\ \ \ |2\rangle=|1,\mathbf{k}_{2}\rangle\ \,\ \ \ |3 \rangle=|1,\mathbf{k}_{3}\rangle\, \tag{4}\]  
and the intraband transition in the conduction band, and Channel 2 involving the states  
\[|1^{\prime}\rangle=|-1,\mathbf{k}_{1^{\prime}}\rangle\ \,\ \ \ |2^{\prime} \rangle=|-1,\mathbf{k}_{2^{\prime}}\rangle\ \,\ \ \ |3^{\prime}\rangle=|1,\mathbf{k}_{3^{\prime}}\rangle\, \tag{5}\]  
with the intraband transition in the valence band. For the case of DF generation of the field oscillating with the frequency \(\omega=\omega_{b}-\omega_{a}\), we have \(\mathbf{\Omega}=(\omega,\mathbf{q})\) and \(\mathbf{\Omega}_{a,b}=(\omega_{a,b},\mathbf{q}_{a,b})\), the only relevant component of the conductivity tensor, for the experimental arrangement of Figure 1, is the one with \(i=j=k=x\), and the general expression (2) takes the form:  
\[ \begin{split}\sigma^{xxx}(\mathbf{\Omega},\mathbf{\Omega}_{b},-\mathbf{ \Omega}_{a})&=\frac{e^{3}}{\hbar^{2}\omega_{b}\omega_{a}\omega _{b}}\frac{\sum_{23}^{x_{23}-x_{23}^{5}x_{31}\cdot\mathbf{x}_{22}-\mathbf{\Omega}_{a} }\left(\frac{f(\varepsilon_{1})-f(\varepsilon_{3})}{\omega_{1}-\omega_{b}- \varepsilon_{31}}-\frac{f(\varepsilon_{2})-f(\varepsilon_{1})}{\omega_{1} \omega_{2}+\omega_{a}-\varepsilon_{31}}\right)}\\ &+\frac{e^{3}}{\hbar^{2}\omega_{a}\omega_{b}\omega_{a}}\sum_{ \gamma\gamma\gamma}\frac{v_{1\gamma\gamma\gamma}^{2}-\mathbf{q}^{5}v_{1\gamma\gamma \gamma}\cdot\mathbf{v}_{2\gamma\gamma}-\mathbf{q}}{\omega_{2^{\prime}}-\omega- \varepsilon_{3^{\prime}}\gamma\cdot\mathbf{v}_{2\gamma\gamma}}\left(\frac{f( \varepsilon_{3^{\prime}})-f(\varepsilon_{3^{\prime}})}{\omega_{2\gamma}+ \omega_{a}-\varepsilon_{3^{\prime}}\gamma\gamma}-\frac{f(\varepsilon_{4^{ \prime}})-f(\varepsilon_{3^{\prime}})}{\omega_{3^{\prime}}-\omega_{b}- \varepsilon_{3^{\prime}}\gamma\gamma}\right)}.\end{split}  \tag{6}\]  
where we retained only resonant terms. Evaluation of the matrix elements of the modulated velocity operator (3) can be done in different ways, the most straightforward is considering the velocity operator \(\mathbf{v}=v_{F}\mathbf{\sigma}\) (see SI for details). In principle, in the **A**-gauge one has to calculate them with the full Hamiltonian including the interaction [24; 25]. However, in our case, the interaction operator commutes with \(e^{i\mathbf{r}\cdot\mathbf{q}}\) and we can use the unperturbed  
Figure 2: Transitions and states considered in the calculation of the conductivity for DF generation at frequency \(\omega\): Channel 1 (right) and Channel 2 (left). For the process of generation of probe photons (\(\omega_{b}\)), the green and yellow arrows are inverted. Hamiltonian (without **A**). The matrix elements contain Kronecker \(\delta\)-s that cancel out two summations in Equation (6) and yield the conservation of momentum \(\textbf{q}=\textbf{q}_{b}-\textbf{q}_{a}\). To compute (6), we expand the matrix elements \(v^{j}_{mn,\textbf{q}}\), the energy differences, \(\omega_{mn}\), and the Fermi-Dirac distributions to the first order in the wavevectors \(\textbf{q}_{j}\) and \(\textbf{q}\). For example, for a state \(|n\rangle=|1,\textbf{k}^{\prime}\rangle\) with \(\textbf{k}^{\prime}=\textbf{k}+\textbf{q}_{j}\), then:  
\[f(\varepsilon_{n})=f_{c}(\varepsilon_{\textbf{k}^{\prime}})\approx f_{c}( \varepsilon_{\textbf{k}})+\partial_{\varepsilon_{\textbf{k}}}f(\varepsilon_ {\textbf{k}})\textbf{q}_{j}\cdot\
abla_{\textbf{k}}\varepsilon_{\textbf{k}} \tag{7}\]  
We denote by \(f_{c}(\varepsilon)\) the distribution function for the electrons in the conduction band and by \(f_{v}(\varepsilon)=f_{c}(-\varepsilon)\) the one for electrons in the valence band. Additionally, for simplicity, we set equal all phenomenological damping constants, \(\gamma_{nm}=\gamma\). With these approximations, considering both channels the lowest order term in the conductivity is proportional to \(q\).