# Analysis of All-Optical Generation of Graphene Surface Plasmons by a Frequency-Difference Process
## 2 Materials and Methods
### Second-Order Conductivities
At zero temperature, the following analytical expression is obtained:  
\[\sigma^{xxx}(\textbf{\Omega},\textbf{\Omega}_{b},-\textbf{\Omega}_{a})\approx \delta_{\textbf{q}\textbf{q}_{b}-\textbf{q}_{a}}\frac{e^{3}v_{F}^{2}E_{F}}{h^{ 3}\omega_{b}\omega_{a}\pi^{2}}\frac{2\omega_{F}-\omega_{a}}{\left(2\omega_{F}- \omega_{a}\right)^{2}+\gamma^{2}}qA(\omega,q,\gamma)\,, \tag{8}\]  
where  
\[A(\omega,q,\gamma) = \frac{\pi(\omega+i\gamma)}{(v_{F}q)^{4}}\Bigg{\{}(v_{F}q)^{2}+2( \omega+i\gamma)^{2}\Bigg{[}\sqrt{1-\left(\frac{v_{F}q}{\omega+i\gamma}\right) ^{2}}-1\Bigg{]}\Bigg{\}} \tag{9}\] \[\approx -\frac{\pi}{4}\frac{1}{(\omega+i\gamma)}\]  
is an auxiliary function, which is independent of \(q\) in the limit \(q\to 0\). Figure 3 shows the behavior of the real and imaginary parts of \(\sigma^{xxx}(\textbf{\Omega},\textbf{\Omega}_{b},-\textbf{\Omega}_{a})\) as functions of \(\omega\) and \(q\), for two values of the Fermi energy. The results obtained by the analytical formula (8) are confirmed by direct numerical evaluation of Equation (6). In the left panel, for a fixed value of \(q\), the real part of the conductivity peaks at \(\omega\approx v_{F}q\), and the imaginary part becomes considerable for frequencies \(\omega<v_{F}q\). An important feature of Equation (8) for the second-order conductivity \(\sigma^{xxx}(\textbf{\Omega},\textbf{\Omega}_{b},-\textbf{\Omega}_{a})\) is that it goes to zero as \(q\to 0\), as it should for a centrosymmetric material. This physically meaningful feature was lacking in the results of Refs. [8; 23]. The reason for this probably is the problem with intraband matrix elements for the **r** operator, first pointed out by [PERSON] [26]. The analytical formula derived in Ref. [22], although showing the correct limiting behavior for \(q\to 0\), is somewhat different from ours and we were not able to reproduce it following the steps indicated in the supporting material of that article. Let us emphasize that the points in Figure 3 represent the results of direct numerical evaluation of Equation (6), performed without approximations such as expansions in terms of wavevectors, so we believe that these spectra are correct within the model and approach used. It is the imaginary part of the conductivity that is important to describe the SP generation process, so we have from Figure 3 a \(|\text{Im}\sigma^{xxx}(\textbf{\Omega},\textbf{\Omega}_{b},-\textbf{\Omega}_{ a})|\approx 3\leavevmode\
obreak\ esu\approx 10^{-16}\text{Am}/\text{V}^{2}\) for \(\omega=50\) ps\({}^{-1}\). Even though it is larger than the value calculated in Ref. [22], it is below the experimentally estimated value [6] by some 3 orders of magnitude. We shall discuss the possible cause of this discrepancy in Section 3. Similar calculations can be carried out for the conductivities oscillating with the frequencies \(\omega_{b}\) and \(\omega_{a}\), representing the DF destruction/generation of the pump and probe photons, respectively. Our analytical expressions, obtained with the same approximations as Equation (8), are:  
\[\sigma^{(2)}_{xxx}(\mathbf{\Omega}_{b},\mathbf{\Omega}_{a},\mathbf{\Omega})\approx-\delta_{ \textbf{q}\textbf{q}\textbf{q}}-\textbf{q}\frac{e^{3}p_{F}^{2}E_{F}}{\hbar^{3} \omega_{a}\omega\pi^{2}}\frac{q}{2\omega_{F}-\omega_{a}-i\gamma}A(\omega,q, \gamma)\,; \tag{10}\]  
\[\sigma^{xxx}(\mathbf{\Omega}_{a},\mathbf{\Omega}_{b},-\mathbf{\Omega})\approx\delta_{ \textbf{q}\textbf{q}\textbf{q}}-\textbf{q}\frac{e^{3}p_{F}^{2}E_{F}}{\hbar^{3} \omega_{b}\omega\pi^{2}}\frac{q}{2\omega_{F}-\omega_{a}-i\gamma}A(\omega,q,- \gamma)\,. \tag{11}\]  
We notice that  
\[\frac{1}{\omega_{a}}|\sigma^{xxx}(\mathbf{\Omega}_{a},\mathbf{\Omega}_{b},-\mathbf{\Omega })|\approx\frac{1}{\omega_{b}}|\sigma^{(2)}_{xxx}(\mathbf{\Omega}_{b},\mathbf{\Omega }_{a},\mathbf{\Omega})|\approx\frac{1}{\omega}|\sigma^{xxx}(\mathbf{\Omega},\mathbf{\Omega }_{b},-\mathbf{\Omega}_{a})|\,. \tag{12}\]  
This relation can be used for numerical estimates of the former [6]. We also verified Equations (10) and (11) by direct numerical evaluation.