MATLAB formulas for computing the nonlinearity of graphene, with emphasis on the optical band (e.g., NIR) and the self-acting (Kerr-like) effect. Currently, the repository only includes a function for the perturbative (instantaneous) third-order nonlinearity, but I mean to add more formulas for the non-perturbative/electrodynamic regime, i.e., when light intensity is higher and free-carrier effects come into play, giving rise to transient response. These formulas can be useful in nonlinear optics and ultrafast science.
MATLAB implementation of the quantum-derived formulas in [1], for the calculation of the complex-valued nonlinear surface conductivity sigma(3) in a perturbative regime. The sigma3 value depends on frequency, chemical potential, finite (nonzero) temperature, and finite relaxation rates; the latter can be different for intraband and interband transitions. The function also computes complex-valued sigma(1) spectra.
The core-fuction is sigma3_Graphene_2015JLC.m and I've also included a test-script (script_Repro_2015JLC_Fig4.m) to showcase the execution and validate its outputs. Executing the core-function (with no arguments) should produce the following figure in several seconds, which reproduces the blue curves in Fig. 3(b) of [1], concerning third-harmonic generation (THG).
Fig.1 : Complex spectra of THG sigma3 (xxxx-component) of graphene, at 300 K (room temperature) and at 33 meV relaxation rates.
Executing the test-script should produce the following repro of all four panels of Fig. 4 in [1], conerning the Kerr effect. This should take somewhat longer, e.g., 1-2 minutes.
Fig.2 : Complex spectra of Kerr sigma3 (xxxx-component) of graphene, for various combinations of temperature (zero and room temperature) and inter/intraband momentum relaxation rates. In all cases, graphene's chemical potential (mu) is tuned to 0.3 eV, so note that sharp features appear for frequencies near 0.6 eV, i.e., when mu is tuned near the 'half photon-energy'.
These functions have been developed and tested in MATLAB R2014a. I'd appreciate any comment or feedback :)
Oh, and if you find these useful, please cite my WIP papers (hopefully they can be published soon)
- [3] "Graphene optical nonlinearity: From the third-order to the non-perturbative electrodynamic regime" (DOI=10.48550/arXiv.2203.09373) and/or
- [4] "Kerr Microcombs in Integrated Waveguide Ring Resonators enabled by Graphene Nonlinearity" (DOI=10.36227/techrxiv.171072961.17885911)
This function uses the JLC formulas in [1], that go one step beyond paper [2] by SAM, allowing nonzero temperature and unequal Gi and Ge (paper [2] treats zero temperature and finite Gi==Ge, i.e., equal intra/inter-band rates. However, I should note that the formulas in [2] are generally more stable (non divergent) and can, of course, be converted to finite (nonzero) temperature.
JLC [1] implementation deals with the full arbitrary case of (f1,f2,f2) frequency mix, for any tensor component. This function computes the SYMMETRIZED sigma3_xxxx tensor component of input parameter "omega" (in [rad/sec]) or, equivalently, f (in [Hz]). Depending on the "s3_Formulas" input argument, one can use the full formulas of [1] or the simplified versions for each case (e.g., Kerr or THG). The case used is controlled by the "s3_Mode" input argument (string):
- 'THG' (third-harmonic generation): The frequency 3f is generated by light (at f) interacting with graphene. This (f,f,f) mix of frequencies is the simplest nonlinearity case.
- 'Kerr': The "self-acting" f<->f effect, i.e., the (-f,f,f) mix. The Kerr effect is quantified by the IMAG part of sigma3, which can be self-focusing or defocusing, depending on the sign. The REAL part of sigma3 gives rise to perturbative SA (saturable absorption) or IA (induced absorption, similar to TPA in bulk semiconductors), depending on the sign.
- 'KerrDetune': See below for details.
- 'PFC' (Parametric Frequency Conversion): from the (-f,fp,fp) mix, where the output is at frequency 2*fp-f, with fp=const.
In the Kerr case, note that the full formulas of JLC [1] diverge EXACTLY on f*(-1,1,1), i.e., the computed spectra are "noisy". For this, I recommend using 'KerrDetune' mode, i.e., a slightly detuned mix of frequencies: f*(e-1,1,1), where e = 1e-4. Using e-->0 (e.g., e=1e-6), creates "noisy" (unstable) spectra; using higher e (e.g., e=1e-2) creates smoothed-down spectra, i.e., any sharp Kerr resonances (useful in NL optics) are blurred out.
Paper [1] by JLC also includes simplified formulas for the "pure" cases of THG and Kerr. The THG formulas are OK (meaning that the spectra they give are close to the ones acquired using the Full formulas), but the Kerr ones have an error (in the "L3" term, Eq. (60) in [1], the units in the RHS don't check out). So, for Kerr effect, I recommend using the full formulas at 'KerrDetune' s3_Mode, with the default e~1e-4, as described above.
Note that in the model assumed for gapless doped (mu>0) graphene, there's only 8 nonzero elements in the 4th-rank sigma3 tensor, and these 8 only depend on 3 independent terms: {xxyy,xyxy,xyyx}. In all cases it holds that xxxx=xxyy+xyxy+xyyx, plus the symmetry in swapping x<->y. Now, the formulas in the JLC [1] and SAM [2] paper are given for the UNsymmetrized sigma3_dabc(f1,f2,f3); to go to the SYMMETRIZED sigma3_dabc(f1,f2,f3), which is conventionally used in nonlinear optics, an "averaging" over the up-to-six concurrent permutations of the "dabc" coordinates and the corresponding f_123 frequencies is required, e.g., Eq. (14) in [1]:
d;abc (1,2,3) = 1/6*( d;abc(1,2,3) + d;bca(2,3,1) + d;cab(3,1,2) + d;acb(1,3,2) + d;cba(3,2,1) + d;bac(2,1,3) )
Now, to assess the magnitude of each NL effect (Kerr, THG, etc.) from light-graphene interaction, we conventionally evaluate the xxxx component of symmetrized sigma3. To compute that in the THG case (f1=f2=f3), we need only one permutation of the "dabc coord's", as it holds that the SYMMETRIZED xxyy==xyxy==xyyx = 1/3*xxxx. For the Kerr case (-f,f,f) and the PFC cases, we need only three permutations.
Finally, a note on the stability of the finite-temperature spectra using the smoothing formulas of Appendix B in [1]: In case of "noisy" or inaccurate spectra (e.g., trying to repro curves in [1]), consider decreasing the "Es_stp" (step) variable below and/or increasing the "Es_max". Problems with noisy spectra appear due to the numerical integration (e.g., check Eqs.(B1),(B2),(B4) in [1]) and they are worse for low Gi and/or Ge, i.e., when Gie < 5 meV; low Ge is generally more noise-prone.
[1] JinLuo Cheng (JLC), N. Vermeulen, J. E. Sipe, "Third-order nonlinearity of graphene: Effects of phenomenological relaxation and finite temperature", PRB, 2015. URL/DOI
[2] Sergey A. Mikhailov (SAM), "Quantum theory of the third-order nonlinear electrodynamic effects of graphene", PRB, 2016. URL/DOI
[3] A. Pitilakis and E.E. Kriezis, "Graphene optical nonlinearity: From the third-order to the non-perturbative electrodynamic regime", arXiv, 2022 URL/DOI
[4] A. Pitilakis and E.E. Kriezis, "Kerr Microcombs in Integrated Waveguide Ring Resonators enabled by Graphene Nonlinearity", techRxiv, 2024 URL/DOI