# Graphene optical nonlinearity: From the third-order to the non-perturbative
electrodynamic regime

Alexandros Pitilakis<EMAIL_ADDRESS>School of Electrical and Computer
Engineering, Aristotle University of Thessaloniki  Emmanouil E. Kriezis
School of Electrical and Computer Engineering, Aristotle University of
Thessaloniki

###### Abstract

A non-perturbative model for graphene optical nonlinearity is developed for
the study of ultrafast pulse propagation along a monolayer, as in the case of
graphene-comprising nanophotonic integrated waveguides. This graphene ‘hot
electron’ model (GHEM) builds upon earlier work, based on the Fermi-Dirac
framework for 2D semiconductors, which was aimed mainly at steady-state
absorptive response of a monolayer under free-space laser-beam illumination.
Our extension adapts the GHEM to in-plane light-matter interaction along
graphene monolayers under intense ps-pulse excitation that leads to the
carrier-density saturation regime. We first provide a quantitative overview of
the ‘classic’ perturbative third-order nonlinear response and then study the
static and transient response of graphene as a function of the GHEM
parameters, with focus on the monolayer quality and voltage-tunability. These
results are compared to phenomenological models for saturable absorption and
nonlinear refraction, showing good agreement with recent experimental works.
In conclusion, the GHEM unifies the experimentally observed absorptive and
refractive optical nonlinearities in a single multi-parametric framework
therefore enabling the evaluation of the voltage-tunable light-matter
interaction in structures such as diffraction-limited nanophotonic waveguides.
This formalism can be readily employed to THz frequencies, adapted to multi-
channel (e.g. pump-probe) nonlinear effects, or developed to include carrier
and lattice-temperature diffusion to extend its validity threshold.

††preprint: APS/123-QED

###### Contents

  1. I Introduction
  2. II Linear and Perturbative Third-Order Nonlinear Regime
    1. II.1 Light-Graphene Interaction
    2. II.2 Linear Response
    3. II.3 Perturbative Third-Order Response
  3. III Nonperturbative Regime: Theory
    1. III.1 Mikhailov’s Hot Electron Model
    2. III.2 Extensions and Modifications
      1. III.2.1 Effective Absorbed Intensity
      2. III.2.2 Intraband Rate Equation
      3. III.2.3 Interband Rate Equation
  4. IV Nonperturbative Regime: Results
    1. IV.1 Quasi-Perturbative Regime
    2. IV.2 Nonperturbative Static Response
    3. IV.3 Nonperturbative Transient Response
  5. V Discussion
    1. V.1 Implications for Photonic Waveguides
    2. V.2 Qualitative Trends and Comparison to Experiments
    3. V.3 Outlook
  6. VI Conclusive Remarks
  7. A Fermi-Dirac Framework
  8. B Surface Conductivity Calculation
    1. B.1 Intraband Conductivity
    2. B.2 Interband Conductivity
  9. C Numerical Solution of Equations System
  10. D Pulse Propagation Along Nonlinear Waveguide

##  I Introduction

Graphene monolayers [1] exhibit a broadband nonlinear response [2], from the
MIR/THz to the NIR/VIS spectral region, with diversified absorptive and
refractive features, which are moreover controllable thanks to the material’s
high electro-optic tunability [3]. More than a decade ago intense research on
the subject ignited mostly due to experimental evidence in [4], which
corroborated the first theoretical predictions in [5, 6, 7]. Various turns of
events have marked the way, including orders-of-magnitude shifts in the
absolute values reported, sign changes, and wide deviations in underlying
model parameters. Despite this turbulent course, all high-power evidence [8,
9, 10, 11] indicates that graphene possesses a large magnitude, broadband and
tunable nonlinearity, with plenty of untapped potential both in theory and in
practical implementation. Naturally, interest shifted from free-space optics
to nanophotonics [12] and integrated photonics [13, 14], a platform where
light-graphene interaction can be maximized, owing to the diffraction-limited
confinement that guided modes can routinely provide [15, 16]. Special interest
is given to on-chip pulsed laser platforms [17, 18], where cavity-buildup can
be elevated to its extreme and the interplay between refractive and absorptive
nonlinearities is ever more crucial to the system response.

Nonlinear response of 2D semiconductors like graphene can be distinguished in
two regimes, firstly the perturbative (or parametric), where nonlinear
polarization quasi instantaneously follows the applied field, and secondly the
non-perturbative (or electrodynamic) regime, which is related to free-carrier
effects as in bulk semiconductors. The carrier-related effects in graphene
have been extensively studied in focused laser beam experiments [19, 20, 21,
22, 23, 24, 25], where NIR illumination is mostly considered as a pump
mechanism, exciting nonlinear effects which are probed by THz pulse. Back to
the perturbative regime, more relevant to integrated photonics, complicated
formulas have been theoretically extracted for the associated third-order
surface conductivity in [26, 27, 28], considering monolayers in free-space,
and experiments [29, 30] have validated and slightly twisted the predictions,
especially regarding the sign of the refractive nonlinearities. Theoretical
studies [31, 16, 32] and experiments [33, 34, 35, 36, 37] in integrated
nonlinear waveguide devices followed, which paved the way toward voltage-
tunable ultrafast absorptive and refractive nonlinearity in photonics. While
the interested audience awaits for more integrated nonlinear devices and
experiments, work progresses in designing new components [38, 39, 40, 41, 42,
43, 44, 45, 46, 47] and extending the theoretical modeling horizons [48, 49,
50, 51, 52, 53, 54, 55].

This work presents a comprehensive overview and a comparative assessment of
graphene nonlinearity across diverse operation regimes, and elucidates issues
related to the exploitation of the Kerr-effect and saturable absorption along
integrated NIR photonic waveguides. For this reason, we mainly focus on the
effect of the voltage-tunable chemical potential and the fabrication quality
of graphene monolayers. We first briefly visit the classic linear and third-
order nonlinear regimes, collectively referred to as perturbative response. We
then move to the strongly nonlinear non-perturbative regime within the
framework of an electrodynamic model obeying Fermi-Dirac statistics which we
refer to as graphene ‘hot electron’ model (GHEM). This model is based on
foundational work in [52], which was devised for free-space laser-beam
transmission through a graphene monolayer, Fig. 1(a), mainly studying
absorption change in the steady-state/static (continuous wave, CW) excitation.
In this work, we adapt and modify the GHEM for transient analysis, i.e., for
quantifying the nonlinear response to picosecond pulses propagating along
graphene-comprising waveguides, Fig. 1(b). Moreover, we consider intensities
corresponding to the photogenerated carrier-density saturation regime and
focus on the coupled refractive and absorptive nonlinear effects. Our results
are in good qualitative agreement with recent experimental results, in both
free-space and guided-wave configurations, and provide a range of parameters
that can be fit to measurements for quantitative evaluation. This GHEM extends
the work of [52], to the study of ultrafast pulse propagation in highly
nonlinear integrated waveguides, providing physical insight to the analysis
and design guidelines. Extension of this formalism to other nonlinear 2D
materials [1, 56, 2] can also be envisaged.

Figure 1: High-power pulse (a) normally impinging on an infinite graphene
monolayer lying on a thick dielectric substrate, (b) propagating along an
integrated nanophotonic waveguide whose mode spatially overlaps with a
graphene monolayer.

The rest of the paper is organized as follows: In Section II we set out from
the classical linear and perturbative third-order nonlinear regime,
highlighting their pertinent aspects for pulse propagation in graphene-
comprising waveguides beyond quasi-equilibrium. Section III briefly outlines
the GHEM model theoretical framework and IV presents a comprehensive study of
the parametrically evaluated static (CW) and transient (pulsed) response.
Section V discusses our results’ implications for photonic waveguides,
assesses qualitative trends, compares to recent experiments, and outline
possible extensions of the model. Finally, following the concluding remarks of
Section VI, a set of Appendices contain the bulk of the model formulas and
elaboration on their less intuitive aspects, together with auxiliary proofs
and numerical implementation details.

##  II Linear and Perturbative Third-Order Nonlinear Regime

Graphene’s EM response to wavelengths below UV is well described by surface
conductivity tensors $\tilde{\sigma}^{(n)}$, denoting a 2D/sheet material,
rather than with susceptibility tensors $\tilde{\chi}^{(n)}$, which
characterise 3D/bulk materials; this is due to the 0.35 nm thickness of the
graphene monolayer.

In the absence of strong magnetic biasing, graphene’s surface conductivity is
isotropic and, thus, it can be fully characterized by a single scalar complex
value, $\sigma^{(n)}$. As a 2D/sheet medium, graphene only interacts with
electric fields tangential to their surface, giving rise to a surface current
density according to Ohm’s law,
$\mathbf{J}_{s}=\sigma^{(1)}_{T}\mathbf{E}_{\parallel}$, in units of A/m. In
this expression,
$\sigma^{(1)}_{T}=\sigma^{(1)}_{\mathrm{lin}}+\sigma^{(1)}_{\mathrm{NL}}$ is
the total surface conductivity including both linear and nonlinear
contributions. The nonlinear contribution depends on the local tangential
E-field magnitude, $|\mathbf{E}_{\parallel}|$, and other material-depended
parameters, and can fall in perturbative and non-perturbative regimes; the
former is the subject of this section while the latter is addressed by the
nonlinear GHEM developed in the following sections.

###  II.1 Light-Graphene Interaction

When interacting with EM waves of small photon energies, graphene behaves like
a zero-bandgap 2D semiconductor with a linear energy-momentum band structure,
as opposed to bulk/3D semiconductors which have nonzero bandgap and parabolic
band structure. As a semiconductor, graphene’s electro-optical response can be
deduced from Fermi-Dirac statistics, governed by the Pauli exclusion
principle, in terms of a common ‘hot’ carrier temperature ($T$) and the quasi
Fermi-levels for the conduction ($\mu_{e}$) and valence ($\mu_{h}$) band
occupation, corresponding to electron and hole plasmas. Under thermal
equilibrium, i.e., in the absence of strong optical illumination, it holds
that $T\equiv T_{0}$, where $T_{0}$ is the lattice temperature, and
$\mu_{e}\equiv\mu_{h}=\mu_{c}$, where $\mu_{c}$ is the chemical potential of
graphene. At zero temperature $\mu_{c}=\mathcal{E}_{F}$, i.e., the chemical
potential is equal to the Fermi energy of graphene (where the Fermi-Dirac
carrier distribution function equals 0.5), whereas at room temperature it is
only slightly smaller; consequently, quantities $\mu_{c}$ and
$\mathcal{E}_{F}$ are sometimes used interchangeably. Negative or positive
$\mu_{c}$ denotes the holes or the electrons as majority carriers,
respectively. However, the two types of carriers are symmetric in graphene
(‘ambipolar’ electron and hole transport), meaning that $\mu_{c}=\pm|\mu_{0}|$
will theoretically have identical optical properties, which is the case
throughout this work. Nevertheless, graphene fabrication techniques (e.g., CVD
growing or exfoliation) and its environment (e.g., attachment to dielectric or
semiconductor substrates) lead to either positive or negative charge
densities. For more details on the Fermi-Dirac framework, refer to Appendix A.

Graphene electro-optical properties can be abstracted in Dirac cone diagrams
such as the ones in Fig. 2, depicting band filling, Eq. (15), at thermal
equilibrium and for a few characteristic $\\{T,\mu_{c}\\}$ combinations. It
can be inferred that graphene is a semi-metal, i.e., it can act either as a
conductor or a dielectric, if its $|\mu_{c}|$ is below or above the half-
photon energy, $\hbar\omega/2$, respectively. In its dielectric regime
($|\mu_{c}|\gg\hbar\omega/2$), graphene is highly transparent in the NIR and
above. But, in its metallic regime ($|\mu_{c}|\ll\hbar\omega/2$), it can
exhibit remarkably high conductivity for its sub-nm thickness, from the THz up
to the visible, and even plasmonic behaviour (i.e., equivalent permittivity
with $\mathrm{Re}\\{\varepsilon_{\mathrm{r,eq.}}\\}\ll-1$) in the THz.
Evidently, control over graphene’s $\mu_{c}$ (or equivalently its carrier
density, Appendix A) exerts control over its optical response, i.e., over the
surface conductivity that an electromagnetic (EM) field experiences when in
proximity to the material. There are various schemes for tuning $\mu_{c}$,
either at fabrication or electrically, with biasing or gating, using
appropriate electrode configurations [35].

Figure 2: Energy-momentum graphene cones depicting band-filling at equilibrium
lattice temperature $T_{0}$. The colormap denotes the distribution function of
carriers in the two bands and the horizontal dashed lines denote the chemical
potential, $\mu_{c}$. (a) Pristine graphene at zero temperature. (b) High
electron density at 100 K. (c) Low hole density at 300 K. (d) Quasi-unbiased
graphene super heated at 3000 K. Interaction with 0.8 eV photons gives rise to
both intraband and interband surface conductivity in all cases, except case
(b), where Pauli blocking cancels the interband mechanism.

Figure 2 additionally provides the groundwork for the nonlinear GHEM, i.e.,
depicts the two mechanisms with which graphene interacts with optical photons:
Intraband absorption shifts carriers away from the Dirac point within the same
band whereas interband absorption of photons with
$\mathcal{E}_{\mathrm{ph}}=\hbar\omega>2\mu_{c}$ (above Pauli blocking) moves
carriers between bands, depicted as small colored dots in Fig. 2. Evidently,
the intraband mechanism is equivalent to a carrier temperature shift (heating
and cooling) whereas the interband mechanism is responsible for the
photogeneration of electron-hole pairs (satisfying electro-neutrality
condition, in the absence on carrier injection) which eventually recombine. At
low-power excitation, these two mechanisms act perturbatively meaning that the
graphene’s quantum state, i.e., its $\\{T,\mu\\}$ properties, remain unchanged
(thermal equilibrium is instantaneously reinstated) and the response is
assumed instantaneous and fully described by a limited set of $\sigma^{(n)}$
conductivities in a Taylor-like expansion [57, 58]. In contrast, at
sufficiently high excitation, graphene’s quantum state is changed, i.e., its
$\\{T,\mu\\}$ properties are non-perturbatively affected, which leads to a
nonlinear response, strongly related to the free-carrier (plasma)
spatiotemporal dynamics.

In this work, we are concerned primarily with graphene’s behaviour in the NIR,
$\hbar\omega\approx 0.8$ eV ($\lambda_{0}\approx 1550$ nm), and for
$|\mu_{c}|$ close to and below the half-photon energy; this is interesting for
three reasons: (i) due to the switching of graphene’s qualitative response
manifested around that threshold, (ii) because of the non-trivial physics
required to quantitatively evaluate its nonlinear properties, and (iii) as
most state-of-the-art fabricated graphene-comprising waveguide devices for the
NIR have $|\mu_{c}|$ in that vicinity, when unbiased, and a few-Volt potential
can shift it to cover the whole meaningful range $\mu_{c}=\pm\hbar\omega$.

###  II.2 Linear Response

The linear response of graphene to a low intensity EM wave at harmonic
frequency $\omega_{0}$ can be extracted by a Kubo-like semi-classical
formalism, i.e., tight-binding in the single-atom approximation. There are
multiple pathways that lead to the same formulas [59, 26, 27, 28], which
produce the total linear surface conductivity of graphene,
$\sigma^{(1)}=\sigma^{(1)}_{i}+\sigma^{(1)}_{e}$, as a function of parameters
$\\{\omega_{0},\mu_{e},\mu_{h},T\\}$ and the momentum relaxation rates
$\Gamma_{i,e}=\hbar/\tau_{i,e}$. Subscripts $\\{i,e\\}$ denote $\\{$intraband,
interband$\\}$ throughout the document; the rates ($\Gamma$) are usually
expressed in meV units and the corresponding lifetimes ($\tau$) in fs, where
the $\tau_{\mathrm{[fs]}}\Gamma_{\mathrm{[meV]}}\approx 658$ relation is handy
for conversion.

At finite (nonzero) temperature, the formulas for the $i$\- and
$e$-conductivities entail integration over the energy spectrum, accounting
also for possible energy dependence of the momentum relaxation rates,
$\Gamma_{i,e}=\Gamma_{i,e}(\mathcal{E})$. This is required to model and
balance the effect of various scattering mechanisms that mediate the carrier
relaxation, namely scattering of carriers from impurities, phonons, and other
lattice imperfections. This dependence is of particular importance for the
intraband term and for high temperatures, where most experiments evidence a
dominance of ‘long-range’ scattering with charged impurities (rather than hot
phonons) [25]; specifically, the relaxation lifetime is proportional to the
energy when $\mathcal{E}>\mathcal{E}_{\mathrm{imp}}$, where
$\mathcal{E}_{\mathrm{imp}}$ is a Coulomb energy that increases with the
square root of the density of impurities [52]. More details, together with the
full and simplified expressions for the numerical computation of
$\sigma^{(1)}_{i,e}(\omega_{0},\mu_{e},\mu_{h},T,\Gamma_{i,e,})$, can be found
in the Appendix B. Note that $\mathrm{Re}\\{\sigma^{(1)}\\}>0$ holds for
absorptive part of the total surface conductivity, whereas the refractive
part, $\mathrm{Im}\\{\sigma^{(1)}\\}$, can be positive or negative. In the
NIR, the effect of graphene sheets is mostly perturbative, especially for the
refractive part; however, in the THz band, $\mathrm{Im}\\{\sigma^{(1)}\\}\gg
1$ which gives rise to a plasmonic behaviour,
$\mathrm{Re}\\{\varepsilon_{\mathrm{r,eq.}}\\}\ll-1$.

Assuming thermal equilibrium, $\mu_{e}=\mu_{h}=\mu_{c}$ and $T=T_{0}$, and
using Eq. (23) for the intraband and (30) for the interband conductivity, we
present in Fig. 3 a brief parametric investigation of the total surface
conductivity against the parameters pertinent to this study: chemical
potential ($\mu_{c}$), temperature ($T$), and intraband momentum relaxation
rate ($\Gamma_{i}$). $\sigma^{(1)}$ is normalized to the universal surface
conductivity, $\sigma_{0}=q^{2}/4\hbar\approx 61~{}\mu$S, where $q$ is the
electron charge and $\hbar$ the reduced Plank constant. We note the high
absorption when $\mu_{c}<\hbar\omega/2=0.4$ eV (interband absorption allowed),
the smoothing of the spectral features as the carrier temperature increases,
and the residual absorption slope for $\mu_{c}>0.4$ eV which gets larger as
the graphene quality drops ($\Gamma$ increases). This residual absorption is
caused by the intraband mechanism at high carrier densities and is quantified
either by a high (constant) $\Gamma_{i}$ or by a high impurity energy
$\mathcal{E}_{\mathrm{imp}}$; more details can be found in the Appendix
subsection B.1. Now, concerning the imaginary (refractive) part of
$\sigma^{(1)}$, we note that it is only affected by temperature near the half-
photon energy, it exhibits a sign-flip coinciding with the absorption drop,
and it scales linearly with $|\mu_{c}|$ implying a proportionality to
$\sqrt{n_{T}}$, where $n_{T}=n_{e}+n_{h}$ is the total carrier density,
Eq.(17). Finally, we note that the curves of Fig. 3 are even-symmetric around
the Dirac point, $\mu_{c}=0$, and that sometimes the horizontal axis is
normalized, i.e., $(\mu_{c}/\hbar\omega)^{\pm 1}$; in this sense the term
‘$\mu_{c}$ spectra’ is used in this work as $\omega$ is fixed.

Figure 3: (a) Real and (b) imaginary part of the normalized total
$\sigma^{(1)}/\sigma_{0}$ vs. $\mu_{c}$ at $\lambda_{0}=1550$ nm (0.8 eV). The
interband momentum relaxation rate is fixed at $\Gamma_{e}=0.5$ meV
($\tau_{e}=1.3$ ps) and the intraband $\Gamma_{i}$ is either constant (dotted
line) or energy dependent according to the model of Eq. (27).

The useful conclusions from this analysis are: (i) that good quality graphene,
such as $\mathcal{E}_{\mathrm{imp}}=30$ meV corresponding to 7000 cm2/Vs
mobility, can be switched between the opaque and transparent regime by tuning
$|\mu_{c}|$ from $0.3$ to $~{}0.5$ eV at room temperature, and (ii) that the
refractive part of the conductivity scales with the square-root of carrier
density and is virtually unaffected by sample quality.

###  II.3 Perturbative Third-Order Response

The first attempts to theoretically evaluate the nonlinear response of
graphene naturally focused on its third-order surface conductivity,
$\sigma^{(3)}$, because in the absence of strong magnetic biasing it holds
that $\sigma^{(2)}=0$, owing to graphene’s centrosymmetric crystal lattice
[57]. Moreover, the generalized fourth-rank tensor
$\tilde{\sigma}^{(3)}_{abcd}(\omega_{0};\omega_{1},\omega_{2},\omega_{3})$,
with $\\{a,b,c,d\\}=\\{x,y,z\\}$ and
$\omega_{0}=\omega_{1}+\omega_{2}+\omega_{3}$, is isotropic, meaning that only
one of its elements, $\sigma^{(3)}_{xxxx}=\sigma_{3}$ by convention, suffices
for population of the full tensor.

It must be noted that $\sigma_{3}$ is complex: the sign of its real part
corresponds to induced absorption ($+$) or induced transparency/absorption
saturation ($-$), while the sign of the imaginary part corresponds to focusing
($-$) or defocusing ($+$) refractive Kerr-type nonlinearity, under the
$e^{-i\omega t}$ phase convention. Here, we must stress that a positive or
negative $\mathrm{Im}\\{\sigma_{3}\\}$ only approximately corresponds to
defocusing or focusing nonlinear refraction, respectively. The proportionality
holds only for large $|\mathrm{Im}\\{\sigma_{3}\\}|$ and low-loss media (low
$\mathrm{Re}\\{\sigma^{(1)}\\}$ in the case of graphene), while the exact
value of the real part of the equivalent nonlinear index
($\mathrm{Re}\\{n_{2}\\}$) and its sign-transition is governed by more
complicated expressions that also account for complex-valued $\sigma^{(1)}$,
similar to the case of bulk nonlinear media with losses [60]. The same
approximation holds for the sign of $\mathrm{Re}\\{\sigma_{3}\\}$, roughly
corresponding to the sign of $\mathrm{Im}\\{n_{2}\\}$.

Early estimates for $\sigma_{3}$ were triggered by remarkably high
experimentally measured four-wave mixing [4], followed by a set of thorough
quasi-classical single-atom derivations (perturbative solution of
semiconductor Bloch equations in 2D), starting from zero temperature and zero
relaxation rate [26] and later moving to finite temperature and relaxation
with distinct $\Gamma_{i}\neq\Gamma_{e}$ [27]. These results were immediately
corroborated by a different quantum electrodynamic model [28], applicable only
for $\Gamma_{i}\equiv\Gamma_{e}$ but avoiding the divergent behaviour as
$\mu_{c}\rightarrow 0$; this was initially a zero-temperature model that was
later extended to finite temperature [61] with a special focus on the
Kerr/two-photon absorption (TPA) effect. For identical inputs,
$\\{\mu_{c},T,\Gamma_{i}\equiv\Gamma_{e}$}, the two models produce equivalent
results for $\mathrm{Re}\\{\sigma_{3}\\}$ but their convergence in the
refractive part, $\mathrm{Im}\\{\sigma_{3}\\}$, is limited to high $\Gamma>20$
meV or above one-photon energy, $\mu_{c}>\hbar\omega$. In any case, both
models can be cast to the study of various third-order phenomena, such as
Kerr/TPA, third-harmonic generation
($\omega_{0}/3=\omega_{1}=\omega_{2}=\omega_{3}$), or parametric frequency
conversion ($\omega_{1}\neq\omega_{2}=\omega_{3}$ so that an idler at
$\omega_{0}=2\omega_{\mathrm{pump}}-\omega_{\mathrm{sig}}$ is generated by a
weak and strong signal at different wavelengths). It should be noted that
despite the isotropic nature of the $\tilde{\sigma}^{(3)}$ tensor, extracting
the complex scalar $\sigma_{3}=\sigma_{3}(\omega_{0},\mu_{c},T,\Gamma)$ in its
customary symmetrized form requires complicated calculations of multiple
unsymmetrized tensor components [27, 28].

In this work, we are interested in single-channel self-acting nonlinear
phenomena, so we will solely study nonlinear conductivity related to the Kerr
effect, $\sigma^{(3)}_{xxxx}(\omega;-\omega,+\omega,+\omega)$, and evaluate
its dependence on graphene parameters $\\{\mu_{c},T,\Gamma\\}$ at
$\lambda_{0}=1550$ nm (0.8 eV). Using the model of [28] we study the $\mu_{c}$
spectra of $\sigma_{3}$, equal to $\sigma^{(3)}_{\mathrm{Kerr}}$ normalized
with respect to 1 S(m/V)2, for a few pertinent parameter combinations; note
that, similarly to $\sigma^{(1)}$, the $\sigma^{(3)}$ spectra are also even-
symmetric for negative $\mu_{c}$. In Fig. 4 we study the effect of chemical
potential and scattering rates $\Gamma_{i}=\Gamma_{e}$, at room temperature
(300 K): For $\mu_{c}<\hbar\omega/2$, in all cases we get SA and defocusing
refraction (DFR), both ‘flat’ and vanishing as $\Gamma_{i,e}$ increase, i.e.,
as the quality of the graphene sample decreases. Above half-photon energy, the
magnitude of nonlinearity decreases and experiences sign changes.

Figure 4: Effect of relaxation rate $\Gamma_{i}=\Gamma_{e}$ on $\mu_{c}$
spectra, at $\lambda_{0}=1550$ nm (0.8 eV) and $T=300$ K. (a) Real part of
$\sigma_{3}$, with its sign denoting saturable absorption (SA) or induced
absorption (IA). (b) Imaginary part of $\sigma_{3}$, with its sign denoting
defocusing or focusing refractive nonlinearity.

In Fig. 5 we study the effect of chemical potential and temperature, for a
typical value of $\Gamma_{i}=33$ meV ($\tau_{i}=20$ fs). We notice that the
increasing temperature effect is almost exclusively near and above half-photon
energy, extending the SA and DFR regime to higher chemical potentials; note
the peaking of $\mathrm{Im}\\{\sigma_{3}\\}$ on half-photon energy, whose
resonance increases as temperature decreases.

Figure 5: Effect of lattice temperature on $\mu_{c}$ spectra at
$\lambda_{0}=1550$ nm (0.8 eV) and $\Gamma_{i}=\Gamma_{e}=33$ meV, i.e.,
$\tau=20$ fs. (a) Real part of $\sigma_{3}$, with its sign denoting saturable
absorption (SA) or induced absorption (IA). (b) Imaginary part of
$\sigma_{3}$, with its sign denoting defocusing or focusing refractive
nonlinearity.

Summarizing this subsection, for the typical case of room temperature and
$\mu_{c}\approx 0.2$-$0.3$ eV, one expects SA and DFR, from the real and
imaginary parts of $\sigma_{3}$, respectively, with low sensitivity on
$\mu_{c}$. Moreover, the magnitude of perturbative nonlinearity depends
strongly on the momentum relaxation lifetimes: As the sample quality improves
($\Gamma_{i,e}\rightarrow 0$) both SA and DFR increase by orders or magnitude.
This is the reason for the logarithmic scaling used in the vertical axes of
Fig. 4 and 5.

##  III Nonperturbative Regime: Theory

The theory and governing equations of the GHEM will be outlined in this
section. We set off from the comprehensive model developed by S. A. Mikhailov
in [52], and adapt it to the study of ultrafast pulse propagation in graphene-
comprising nonlinear waveguides. In all this analysis, we have adopted the
Fermi-Dirac framework described in Appendix A for modeling the carrier plasma
in graphene under intense optical excitation that pushes the medium out of
thermal equilibrium.

###  III.1 Mikhailov’s Hot Electron Model

The GHEM developed in [52] relies firstly on instantaneous absorption, i.e.,
as if the carrier-carrier scattering time is zero, and secondly on the
splitting of the effects arising from the intraband and interband absorption
contributions. This is based on the assumption that the timescales of the two
effects are different so that they can be assumed to act independently on the
carrier plasma. Specifically, the intraband absorption, proportional to
$\mathrm{Re}\\{\sigma^{(1)}_{i}\\}$ from Eq. (23), instantaneously heats the
carrier plasma which subsequently cools down in a $\tau_{\mathcal{E}}$
timescale, typically sub-ps and higher than the intraband momentum relaxation
lifetime, $\tau_{i}=\hbar/\Gamma_{i}$. The interband absorption, proportional
to $\mathrm{Re}\\{\sigma^{(1)}_{e}\\}$ from Eq. (24), instantaneously
generates electron-hole pairs which subsequently recombine in a longer few-ps
timescale, i.e., $\tau_{\mathrm{rec}}\gg\tau_{\mathcal{E}}$; note that
$\tau_{\mathrm{rec}}$ has the meaning of recombination time only at very low
excitation intensities. In [52], it was found that the ratio between these two
phenomenological lifetimes governs the response; in this work we set
$\tau_{\mathrm{rec}}=10\tau_{\mathcal{E}}$, unless otherwise stated. The two
absorption mechanisms can thereby be considered to act separately on graphene,
and form two rate equations: one for the total energy density in the hot
plasma, $\mathcal{E}^{D}_{T}$, and one for the photogenerated carrier density,
$n_{\mathrm{PG}}$. Finally, as a transitional state between these effects, the
GHEM assumes that a quasiequilibrium (QE) state is instated shortly after
illumination, where: (i) the plasma is characterised by QE chemical
potentials, $\mu_{e_{0}}\neq\mu_{h_{0}}$, which are moreover different from
graphene’s effective hot $\mu_{e,h}$, (ii) the photogenerated carriers have
not recombined, and (iii) the plasma has cooled down to the lattice
temperature, $T\approx T_{0}$.

The GHEM in [52] has six unknowns:
$\\{\mathcal{E}^{D}_{T},n_{\mathrm{PG}},\mu_{e},\mu_{h},\mu_{e_{0}},\mu_{h_{0}}\\}$,
and so requires six equations. The first two are rate equations related to the
aforementioned intraband and interband absorption, as

$\frac{\partial\mathcal{E}^{D}_{T}}{\partial t}=\langle
I_{\mathrm{abs}}\rangle_{i}-\frac{\mathcal{E}^{D}_{T}-\mathcal{E}^{D}_{\mathrm{QE}}}{\tau_{\mathcal{E}}},$
(1) $\frac{\partial n_{\mathrm{PG}}}{\partial t}=\frac{\langle
I_{\mathrm{abs}}\rangle_{e}}{\hbar\omega}-\frac{n_{\mathrm{PG}}}{\tau_{\mathrm{rec}}}\left(1+\frac{n_{\mathrm{PG}}}{n_{T_{0}}}\right),$
(2)

respectively. In Eq. (1),
$\mathcal{E}^{D}_{T}=\mathcal{E}^{D}_{e}(\mu_{e},T)+\mathcal{E}^{D}_{h}(\mu_{h},T)$
is the total plasma energy density in the out-of-equilibrium ‘hot’ state,
$\mathcal{E}^{D}_{\mathrm{QE}}=\mathcal{E}^{D}_{\mathrm{QE}}(\mu_{e_{0}},\mu_{h_{0}},T_{0})$
is the corresponding energy density at the QE state, and
$\tau_{\mathcal{E}}\leq 1$ ps is the phenomenological energy-relaxation
lifetime. In Eq. (2),
$n_{T_{0}}=n_{e}^{0}+n_{h}^{0}=n_{e}(\mu_{c},\mu_{c},T_{0})+n_{h}(\mu_{c},\mu_{c},T_{0})$
is the total carrier density at thermal equilibrium in the absence of
illumination (where $\mu_{e}\equiv\mu_{h}=\mu_{c}$, the sample’s chemical
potential, and $T=T_{0}$, the lattice temperature), and
$\tau_{\mathrm{rec}}\geq 10\tau_{\mathcal{E}}$ is the effective recombination
(or, more abstractly, decay-rate) lifetime for the photogenerated carriers.
The second term, $n_{\mathrm{PG}}/n_{T_{0}}$, in the parentheses in the RHS of
Eq. (2) arises from the assumption that recombination rate is proportional to
the difference of the product of carrier densities between the hot and the QE
state, $n_{e}n_{h}-n_{e}^{0}n_{h}^{0}$, typical in semiconductors. Calculation
formulas for the energy and carrier densities can be found in the Appendix A.

Evidently, the evolution of graphene’s photoconductive response depends on an
effective absorbed intensity, $\langle I_{\mathrm{abs}}\rangle_{i,e}$,
measured in W/m2. In [52], where we have the zero-dimensional assumption,
i.e., a plane wave of sufficiently wide and uniform focal-spot as in Fig.
1(a), we simply define

$\langle I_{\mathrm{abs}}\rangle_{i,e}=A_{i,e}I_{\mathrm{inc}}.$ (3)

In this spectral domain expression, $I_{\mathrm{inc}}$ is the incident plane-
wave intensity (in W/m2) and $A_{i,e}$ is a (unitless) power absorption
coefficient, that can be deduced from the ‘Tinkham formula’ [62], for the
perturbative case of an ultrathin conductive sheet in the interface between
air and a semi-infinite dielectric substrate, as follows:

$A_{i,e}=A_{i,e}(\omega_{0},\mu_{e},\mu_{h},T)=\frac{4Z_{0}\mathrm{Re}\\{\sigma^{(1)}_{i,e}\\}}{|n_{\mathrm{sub}}+1+\sigma^{(1)}Z_{0}|^{2}},$
(4)

where $n_{\mathrm{sub}}$ is the substrate refractive index at harmonic
frequency $\omega_{0}$ and $Z_{0}\approx 377$ Ohm is the free-space impedance.
Note the full (complex) graphene conductivity,
$\sigma^{(1)}=\sigma^{(1)}_{i}+\sigma^{(1)}_{e}$, in the denominator, as
contrasted to the real part of either intra- or interband term in the
nominator. These surface conductivity spectra depend on the hot (out-of-
equilibrium) graphene plasma parameters,
$\sigma^{(1)}_{i,e}=\sigma^{(1)}_{i,e}(\omega_{0},\mu_{e},\mu_{h},T)$, with
the formulas of Appendix B. The bandwidth of the nonlinear response in this
model depends solely on $\sigma^{(1)}$, whose response is rather broadband as
evidenced by Fig. 3, especially as temperature increases.

The remaining four equations of the GHEM, that complement the two rate
equations, are algebraic and stem from electroneutrality (the number of
photogenerated holes is equal to that of photogenerated electrons) in the hot
and QE states:

$\displaystyle\mu_{e}$ $\displaystyle=+(k_{B}T)F_{1}^{-1}\left[\frac{\pi(\hbar
v_{F})^{2}}{(k_{B}T)^{2}}(n_{e}^{0}+n_{\mathrm{PG}})\right],$ (5a)
$\displaystyle\mu_{h}$ $\displaystyle=-(k_{B}T)F_{1}^{-1}\left[\frac{\pi(\hbar
v_{F})^{2}}{(k_{B}T)^{2}}(n_{h}^{0}+n_{\mathrm{PG}})\right],$ (5b)
$\displaystyle\mu_{e_{0}}$
$\displaystyle=+(k_{B}T_{0})F_{1}^{-1}\left[\frac{\pi(\hbar
v_{F})^{2}}{(k_{B}T_{0})^{2}}(n_{e}^{0}+n_{\mathrm{PG}})\right],$ (5c)
$\displaystyle\mu_{h_{0}}$
$\displaystyle=-(k_{B}T_{0})F_{1}^{-1}\left[\frac{\pi(\hbar
v_{F})^{2}}{(k_{B}T_{0})^{2}}(n_{h}^{0}+n_{\mathrm{PG}})\right],$ (5d)

where $F_{1}^{-1}$ is the inverse of the Fermi-Dirac integral of first order,
defined in Eq. (20).

For the numerical solution of the GHEM, in static (CW) or transient (pulsed)
mode, refer to Appendix C.

###  III.2 Extensions and Modifications

The contribution in this work includes three modifications to the GHEM of
[52]:

####  III.2.1 Effective Absorbed Intensity

For the application of the model on pulse propagation along a graphene sheet,
as in Fig. 1(b), and not impinging on a sheet, as in Fig. 1(a), we first need
to modify the effective absorbed intensity, $\langle I_{\mathrm{abs}}\rangle$.
In nanophotonic waveguides the incident intensity does not have a large and
smooth spot, as in focused laser beams, but it has a mode distribution across
the graphene sheet, close to the diffraction limit. In this case, we first
define a local absorbed intensity for narrowband quasi-harmonic fields [35,
43],

$I_{\mathrm{abs},i/e}(\mathbf{r};t)=\frac{1}{2}\mathrm{Re}\\{\sigma^{(1)}_{i/e}(\mathbf{r};t)\\}|\mathbf{E}_{\parallel}(\mathbf{r};t)|^{2},$
(6)

where $t$ is a retarded timeframe (not the optical cycle), e.g., the pulse-
envelope modulating a NIR carrier frequency $\omega$. Equation (6) can be
derived for time-harmonic versions of Ohm’s law,
$\tilde{\mathbf{J}}=\sigma_{\mathrm{bulk}}\tilde{\mathbf{E}}$, and the power-
density resulting from Joule-heating in a volume $V$, $\partial P_{J}/\partial
V=\frac{1}{2}\tilde{\mathbf{E}}\cdot\tilde{\mathbf{J}}^{*}$, adapted to sheet
conductivity [31, 43]. From Eq. (6), we can subsequently abstract the on-chip
pulse power and waveguide geometry, through the mode profile
$\tilde{\mathbf{e}}(x,y)$ in the cross-section, to produce the effective
$\langle I_{\mathrm{abs}}\rangle_{i,e}(t)$ or $\langle
I_{\mathrm{abs}}\rangle_{i,e}(z,t)$ required when applying the pulsed GHEM to
a free-space sheet or through a $z$-segment of a waveguide, respectively. The
angled brackets in $\langle I_{\mathrm{abs}}\rangle$ denote an averaging in
the transverse directions ($xy$ plane), which translates the vector field
$I_{\mathrm{abs}}(x,y)$ to an effective scalar quantity.

Assuming that the power launched into the waveguide mode is $P$ (in Watt),
then the effective intensity from the mode-graphene overlap is
$I_{\mathrm{eff}}=P/A_{\mathrm{eff}}^{G}$ (in W/m2), where
$A_{\mathrm{eff}}^{G}$ is an effective area defined, here, as

$A_{\mathrm{eff}}^{G}=2Z_{0}\mathcal{P}\frac{\int_{G}|\tilde{\mathbf{e}}_{\parallel}(\omega;x,y)|^{2}d\ell}{\int_{G}|\tilde{\mathbf{e}}_{\parallel}(\omega;x,y)|^{4}d\ell}.$
(7)

In this expression, and with reference to Fig. 1(b) axes, the complex vector
$\tilde{\mathbf{e}}_{\parallel}(\omega;x,y)$ is the E-field component of the
eigenmode that is parallel to graphene at frequency $\omega$, $\int_{G}d\ell$
denotes integration along the trace of the graphene sheet in the waveguide
cross-section ($xy$ plane), $Z_{0}\approx 377~{}\Omega$ is the free-space
impedance, and
$\mathcal{P}=0.5\mathrm{Re}\\{\int\tilde{\mathbf{e}}\times\tilde{\mathbf{h}}^{*}dxdy\\}$
(in Watt) is the eigenmode’s power-normalization constant. The expression of
Eq. (7) is derived from the observation that graphene nonlinearities are
proportional to the local E-field intensity tangential to graphene,
$|\mathbf{E}_{\parallel}(x,y)|^{2}$, and power-weighted in the cross-section.
Note that this $A_{\mathrm{eff}}^{G}$ is not equivalent to the effective area
used in waveguides with bulk nonlinear media, such as silica-core optical
fibers [63] or silicon (semiconductor) integrated structures [64, 65, 66, 15].

The effective area of Eq. (7) is a characteristic of the waveguide mode in the
linear regime and can be assumed unchanged in the nonlinear NIR regime too, as
graphene’s refractive effect (i.e., on the waveguide-mode concentration and
profile) is perturbative across a large bandwidth around $\omega$. Then, using
Eq. (7), (6), and the $I=|E|^{2}/2Z_{0}$ assumption, we can finally define the
spatiotemporal evolution of the effective absorbed power as

$\langle
I_{\mathrm{abs}}\rangle_{i/e}(z,t)=\mathrm{Re}\\{\sigma^{(1)}_{i/e}(z,t)\\}Z_{0}\frac{P(z,t)}{A_{\mathrm{eff}}^{G}}.$
(8)

This expression couples the pulse-envelope power $P(z,t)$ to graphene
transient conductivity $\sigma^{(1)}_{i/e}(z,t)$, via the GHEM [e.g. Eqs. (1),
(2), (5) and the underlying Kubo formulas in the Appendix B] eventually
distorting a high-power pulse as it propagates along the nonlinear waveguide.
We reiterate that in the free-space case studied in [52], which is essentially
a zero-dimensional problem, the absorbed intensity is given as $\langle
I_{\mathrm{abs}}\rangle=AI_{\mathrm{eff,inc}}$, where $I_{\mathrm{eff,inc}}$
is the effective incident optical intensity and $A$ is a power-absorption
coefficient by Eq. (4); note that this is a spectral domain formula, so a
temporal convolution is implied when studying pulsed excitation. Finally, we
stress that diffraction-limited highly-confining nanophotonic waveguides can
have exceedingly small $A_{\mathrm{eff}}^{G}\ll 1~{}\mu$m2 which corresponds
to effective incident intensities $I_{\mathrm{eff,inc}}\gg 1$ MW/cm2, as
inputs to the GHEM. For instance, the TE-polarized mode in the silicon-slot
waveguide of [45] has $A_{\mathrm{eff}}^{G}<0.1~{}\mu$m2, meaning that
$I_{\mathrm{eff}}=1$ GW/cm2 corresponds to a modest on-chip
$P_{\mathrm{peak}}<1$ W.

####  III.2.2 Intraband Rate Equation

In the intraband rate equation, we note that the carrier temperature ($T$,
common for hot electrons and holes) can be used instead of the total plasma
energy. This is more efficient for transient numerical computations, as it
does not require a $\\{\mu_{e},\mu_{h},\mathcal{E}^{D}_{T}\\}\rightarrow T$
mapping [i.e., an inversion of Eq. (18)] to acquire the input parameters
$\\{\mu_{e},\mu_{h},T\\}$ required by the Kubo formulas, Eqs. (23-24), to
extract the nonlinear $\sigma^{(1)}_{i,e}$ response. The rate equation for the
carrier temperature can be easily deduced with the chain rule,

$\frac{\partial T}{\partial t}=\frac{\partial
T}{\partial\mathcal{E}^{D}_{T}}\frac{\partial\mathcal{E}^{D}_{T}}{\partial
t},$ (9)

and properties of the Fermi-Dirac integral derivatives, Eq. (21). The
resulting rate equation for the carrier temperate, which replaces Eq. (1) in
the GHEM equation systems, reads

$\frac{\partial T}{\partial
t}=\frac{T}{3\mathcal{E}^{D}_{T}-\mu_{e}n_{e}+\mu_{h}n_{h}}\bigg{(}\langle
I_{\mathrm{abs}}\rangle_{i}-\frac{\mathcal{E}^{D}_{T}-\mathcal{E}^{D}_{\mathrm{QE}}}{\tau_{\mathcal{E}}}\bigg{)},$
(10)

where the derivation of the first factor in the RHS can be found in the
Appendix A and term in parentheses is actually the RHS term of the plasma
energy density rate equation, Eq. (1). We have verified the equivalence of the
two equation systems.

####  III.2.3 Interband Rate Equation

Recent experiments [37, 36] suggest that the density of the photogenerated
carriers in graphene has an upper limit. This is in contrast to bulk
semiconductors where the carrier density is much higher, even at thermal
equilibrium; moreover, under high-intensity excitation, the recombination rate
increases proportionally to the square of the carrier density, thus making a
saturation density practically unreachable in bulk semiconductors. Now, in
graphene, it is theorized that such a saturation density is indeed reachable
and can be the source of non-saturable losses observed in optical and THz
experiments [67, 68, 69]. The saturation density is expected to depend on the
quality of the graphene sample (carrier mobility) and its equilibrium carrier
density. Now, the rate Eq. (2) does not force a hard limit on
$n_{\mathrm{PG}}$, except for any decrease of
$\mathrm{Re}\\{\sigma^{(1)}_{e}\\}$ predicted by the Kubo formula, Eq. (30),
as the carrier temperature increases. However, this decrease is asymptotic
which means that the GHEM can indeed produce unrealistic values of
$n_{\mathrm{PG}}$ as the illuminating intensity increases. In order to avoid
such a scenario, we heuristically introduce a saturation coefficient in the
generation term of the rate equation, $G_{\mathrm{sat}}$, as follows:

$\frac{\partial n_{\mathrm{PG}}}{\partial t}=G_{\mathrm{sat}}\frac{\langle
I_{\mathrm{abs}}\rangle_{e}}{\hbar\omega}-\frac{n_{\mathrm{PG}}}{\tau_{\mathrm{rec}}}\bigg{(}1+\frac{n_{\mathrm{PG}}}{n_{T_{0}}}\bigg{)}.$
(11)

The phenomenological coefficient $G_{\mathrm{sat}}$ is a monotonically
decreasing function of $n_{\mathrm{PG}}$ and can be used to qualitatively
explain experimental observations such as non-saturable conductivity or fit
the GHEM to such results. More discussion in Section IV (e.g., Fig. 9), but
unless explicitly stated assume $G_{\mathrm{sat}}\rightarrow 1$.

On a less critical note, the quadratic factor in the carrier recombination
rate [the second term in the parentheses in the RHS of Eq. (2) or (11)], which
scales with the minority carrier density, can be dropped for $|\mu_{c}|>0.1$
eV, without deviations in the results. Nevertheless, the numerical burden is
negligible, so we retain the term.

##  IV Nonperturbative Regime: Results

In this section, we will present the results obtained using the GHEM developed
in Section III at photon energy $\hbar\omega=0.8$ eV ($\lambda_{0}=1550$ nm).
This model consists of six unknown graphene variables,
$\\{T,n_{\mathrm{PG}},\mu_{e},\mu_{h},\mu_{e0},\mu_{h0}\\}$, and six
equations: (10), (11) and the set of four electroneutrality Eqs. (5).
Moreover, as we are presently studying only graphene’s nonlinear
photoconductivity, we assume that the instantaneous effective absorbed
intensity is related to the incident E-field amplitude with

$\langle
I_{\mathrm{abs}}\rangle_{i,e}=\frac{1}{2}\mathrm{Re}\\{\sigma^{(1)}_{i,e}\\}|E|_{\mathrm{eff}}^{2}=\mathrm{Re}\\{\sigma^{(1)}_{i,e}\\}Z_{0}I_{\mathrm{eff}},$
(12)

i.e., as if we were in a short (unitary length) waveguide where the modal
electric field that is tangential to graphene has an effective intensity of
$I_{\mathrm{eff}}=|E|_{\mathrm{eff}}^{2}/(2Z_{0})$, in W/m2. This is not
equivalent to the case studied in [52], where there is additionally the full
conductivity $\sigma^{(1)}$ in the denominator of Eq.(3), but nevertheless the
same trends will qualitatively hold as the generation rates are primarily
proportional to $\mathrm{Re}\\{\sigma^{(1)}_{i,e}\\}$, in both cases.
According to the analysis in Section III.2.1, given the instantaneous power
$P$ (in Watt) launched in the waveguide, we can calculate the effective
E-field intensity that feeds the GHEM as
$I_{\mathrm{eff}}=P/A_{\mathrm{eff}}^{G}$ (in W/m2) using Eq. (7) for
$A_{\mathrm{eff}}^{G}$. The term photoconductivity, as used here, refers to
the optically induced complex-valued change in graphene’s surface
conductivity:

$\Delta\sigma^{(1)}(t)=\sigma^{(1)}_{\mathrm{GHEM}}(t)-\sigma^{(1)}_{\mathrm{lin}},$
(13)

where $\sigma^{(1)}_{\mathrm{lin}}$ is the surface conductivity at the absence
of high-power excitation (linear regime) and $\sigma^{(1)}_{\mathrm{GHEM}}(t)$
is the nonlinear conductivity predicted by the GHEM under high-power
excitation, which can moreover be time-dependent. Equation (13) can also be
normalized by the universal conductivity $\sigma_{0}$. Finally, unless
otherwise stated, we use a momentum relaxation rate with energy dependence for
the intraband $\Gamma_{i}$ (refer to Appendix B.1 for details) and a constant
$\Gamma_{e}=0.5$ meV ($\tau_{e}=1.3$ ps) for the interband mechanism.

###  IV.1 Quasi-Perturbative Regime

In this regime, the illuminating intensity produces a nonlinear response
through the GHEM, i.e., a photoconductivity $\Delta\sigma^{(1)}$, but this can
be considered as a small perturbation compared to the linear (low-intensity)
conductivity $\sigma^{(1)}_{\mathrm{lin}}$. In this case, we can use a Taylor
approximation around the effective illuminating (CW) intensity
$|E|_{\mathrm{eff}}^{2}=2Z_{0}I_{\mathrm{eff}}$ to extract a quasi-
perturbative $\sigma^{(3)}_{xxxx}(-\omega_{0},+\omega_{0},+\omega_{0})$ for
the self-acting Kerr-type nonlinearity [35],

$\sigma^{(3)}_{\mathrm{QP}}=\frac{2}{3}\frac{\partial\sigma^{(1)}_{\mathrm{GHEM}}}{\partial|E|_{\mathrm{eff}}^{2}}.$
(14)

Note that this approximation is applicable up to a threshold value of
$I_{\mathrm{eff}}$, for which the condition
$|\Delta\sigma^{(1)}|\ll|\sigma^{(1)}_{\mathrm{lin}}|$ holds. For the
$|\mu_{c}|<\hbar\omega/2$ cases, the threshold intensity for this regime was
found to be in the order of 1 MW/cm2 (few tens of mW on-chip power), whereas
it was much higher for $|\mu_{c}|>\hbar\omega/2$ where nonlinear response
practically vanishes. Note that $\sigma^{(3)}_{\mathrm{QP}}$ is reduced, in
both real and imaginary part, when that intensity threshold is surpassed.

Figure 6: Effect of phenomenological energy relaxation lifetime
$\tau_{\varepsilon}$ on the quasi-perturbatively calculated
$\sigma^{(3)}_{\mathrm{QP}}$ at $\lambda_{0}=1550$ nm. (a) Real and (b)
imaginary parts of the $\mu_{c}$ spectra in units of $10^{-19}$ S(m/V)2.

In Fig. 6 we present the quasi-perturbative $\mu_{c}$ spectra for
$\sigma^{(3)}_{\mathrm{QP}}$, extracted from Eq. (14), for a few
phenomenological energy relaxation lifetimes $\tau_{\mathcal{E}}$. The rest of
the GHEM parameters take their default values, i.e., $T=300$ K,
$\lambda_{0}=1550$ nm (0.8 eV),
$\Gamma_{i}=\Gamma_{i}(\mathcal{E};\zeta,\mathcal{E}_{\mathrm{imp}})$ with
$\\{\zeta,\mathcal{E}_{\mathrm{imp}}\\}=\\{4,30~{}\mathrm{meV}\\}$,
$\Gamma_{e}=0.5$ meV, and $\tau_{\mathrm{rec}}=10\tau_{\mathcal{E}}$.
Evidently, when $\tau_{\mathcal{E}}$ increases, the nonlinearity is in overall
larger, owing to the prolonged thermalized state of the plasma. Moreover, we
observe a qualitative agreement with the perturbative regime, Fig. 4, i.e., a
peaking of defocusing nonlinearity when the chemical potential is tuned near
the half-photon energy and a saturable absorption regime near and below that
energy. However, there are also some important differences in these spectra:
First and foremost we observe much higher peak values of
$\sigma^{(3)}_{\mathrm{QP}}$, in the order to $10^{-19}$ S(m/V)2 for both the
real and imaginary parts, even for low $\tau_{\mathcal{E}}$ values. Secondly,
$\mathrm{Re}\\{\sigma^{(3)}_{\mathrm{QP}}\\}$ vanishes except near and below
half-photon energy, whereas in Fig. 4 it had a flat value from $\mu_{c}=0$ up
until the half-photon energy. Both these features are in accordance with the
experimental observations in [35] [and Fig. S6(a-b) in its supporting
information], where a simpler GHEM was used, neglecting the interband
mechanism and the photogenerated carrier density. Thirdly, we observe a
peaking of self-focusing refractive nonlinearity
($\mathrm{Im}\\{\sigma^{(3)}_{\mathrm{QP}}\\}<0$) below half-photon energy,
which is moreover three times higher than the subsequent defocusing peak on
half-photon energy. This interesting new feature emerges independently of
other GHEM parameters, i.e., momentum relaxation rates ($\Gamma_{i,e}$) or
operating wavelength, and appeared neither in Fig. 4 nor in [35]; we attribute
it to the interband absorption mechanism and to the non-negligible
$n_{\mathrm{PG}}>0.15n_{T_{0}}$ it generates. In any case, this feature
indicates a high sensitivity (sign flip) of the refractive nonlinearity on
$\mu_{c}$, i.e., on electrical tuning.

###  IV.2 Nonperturbative Static Response

When the illuminating power increases beyond the threshold discussed in the
previous section, graphene photoconductivity transcends the perturbative
regime. In terms of the GHEM parameters, this is reflected in a non-negligible
carrier temperature ($\Delta T/T_{0}>1$%) and density
($n_{\mathrm{PG}}/n_{T_{0}}>1$%) increase, as well as a shift in the
‘dominant’ quasi-Fermi level (e.g. in $\mu_{h}$ if $\mu_{c}<0$ ). In this non-
perturbative regime, we cannot use the third-order effect formalism, and
directly study the photoconductivity $\Delta\sigma^{(1)}$ as a function of
graphene parameters and the effective incident CW intensity.

We focus on the NIR spectral region, and specifically to the telecom C-band
wavelength $\lambda_{0}=1550$ nm ($\hbar\omega\approx 0.8$ eV). In this
region, graphene’s response is mostly absorptive and does not contribute to
waveguiding, i.e., the contribution of $\mathrm{Im}\\{\sigma^{(1)}\\}$ is
negligible compared to index guiding. Moreover, for chemical potential well
above $\hbar\omega/2$ graphene is practically transparent and a mild induced
absorption is expected when the illumination intensity increases [52], so we
restrict our analysis to $\mu_{c}<0.5$ eV, i.e., just above the half-photon
energy.

An important conclusion drawn directly from [52], Eq. (47) therein, is that
the nonlinear threshold (NLT) in the intensity is inversely proportional to
the phenomenological lifetimes $\tau_{\mathcal{E}}$ and $\tau_{\mathrm{rec}}$,
governing the intraband energy relaxation and interband recombination,
respectively. Secondarily, in the absence of a carrier-saturation density
[$G_{\mathrm{sat}}=1$ in Eq. (11)], the nonlinear response depends only on the
ratio $\tau_{\mathcal{E}}/\tau_{\mathrm{rec}}$ which should moreover be
smaller than unity to validate the assumptions made. In this work, we assume
values $\tau_{\mathcal{E}}=1$ ps and $\tau_{\mathrm{rec}}=10$ ps.

Figure 7: Saturation curves for (a) carrier temperature and $n_{\mathrm{PG}}$,
(b) out-of-equilibrium and quasi-equilibrium chemical potentials, and (c)
surface conductivity. The three colors denote three different combinations of
$\\{\mu_{c},\tau_{\mathcal{E}}\\}$. The insets in panel (b) depict band-
filling as predicted by the GHEM at very low and very high intensities, for
two of the combinations. In all cases: $T_{0}=300$ K, $\lambda_{0}=1550$ nm
and $\tau_{\mathrm{rec}}=10\tau_{\mathcal{E}}$.

In Fig. 7, we extract the saturation curves for the six graphene variables,
$\\{T,n_{\mathrm{PG}},\mu_{e},\mu_{h},\mu_{e0},\mu_{h0}\\}$, and the induced
conductivity as a function of effective incident intensity,
$I_{\mathrm{eff}}$. We consider three combinations of
$\\{\mu_{c},\tau_{\mathcal{E}}\\}$, i.e., equilibrium chemical potential and
energy relaxation lifetime, respectively. In panel (a), we observe the
exponential increase in carrier density and temperature with the logarithm of
intensity. Also, the NLT for carrier photogeneration onset is almost a hundred
times lower than the NLT for temperature increase; this is anticipated as for
low $|\mu_{c}|$ the interband mechanism dominates the total conductivity.
Moreover, we observe that decreasing $\tau_{\mathcal{E}}$ (energy relaxation
lifetime) increases the NLT and that the temperature increase is
indistinguishable for $|\mu_{c}|=0.1$ and 0.3 eV (black and blue dashed curves
overlap). In panel (b), we see that the minority quasi-Fermi level (e.g.,
$\mu_{e}$ if $\mu_{c}<0$) has a very low NLT and experiences a sign-flip as
$n_{\mathrm{PG}}$ increases. Adversely, the dominant quasi-Fermi level
magnitude increases thus leading to a large $\Delta\mu_{e-h}$ that approaches
$\hbar\omega$. Finally, the quasi-equilibrium (QE) chemical potentials are
plotted with thin dot/dash-dot curves, which overlap with the corresponding
‘hot’ chemical potentials and only deviate (to higher magnitudes) at very high
intensities. Finally, in panel (c), we present the conductivity clearly
showing deep SA beyond the NLT of $10$ MW/cm2; the NLT in this case is what is
phenomenologically referred to as saturation intensity, $I_{\mathrm{sat}}$,
i.e., the intensity where absorption drops to half of its low-power (linear-
regime) value.

It is worth pointing out that both real and imaginary parts of $\sigma^{(1)}$
in Fig. 7(c) bear a qualitative resemblance to the corresponding
$\mu_{c}$-dependence, Fig. 3: the negative peaking of the imaginary part
coincides with the halving of the real part. This resemblance unveils an
association between the logarithm of intensity in the nonlinear regime and the
chemical potential (or square-root of carrier density) in the linear regime.

Figure 8: (a) Real and (b) imaginary surface conductivity vs. effective
incident intensity, $I_{\mathrm{eff}}$, for a few $\mu_{c}$ values. (c) Real
conductivity and (d) imaginary photoconductivity, ${\Delta\sigma^{(1)}}$, vs.
$\\{I_{\mathrm{eff}},\mu_{c}\\}$. In all cases: $T_{0}=300$ K,
$\lambda_{0}=1550$ nm, $\tau_{\mathcal{E}}=1$ ps, and $\tau_{\mathrm{rec}}=10$
ps.

Figure 8 is devoted to the effect of the voltage-tunable equilibrium chemical
potential and illuminating intensity on GHEM-predicted nonlinear conductivity.
Panels (a)-(b) present saturation curves for a few values of $\mu_{c}$ whereas
panels (c)-(d) contain heatmaps with finer resolution. The conclusions drawn
from these results are: (i) The empirical fit for SA,
$\mathrm{Re}\\{\sigma^{(1)}(I)\\}=\sigma^{(1)}_{ns}+\Delta\sigma^{(1)}_{\mathrm{sat}}/(1+I/I_{\mathrm{sat}})$,
where $\sigma^{(1)}_{ns}$ is the non-saturable conductivity and
$\Delta\sigma^{(1)}_{\mathrm{sat}}=\mathrm{Re}\\{\sigma^{(1)}(0)\\}-\sigma^{(1)}_{ns}$,
agrees with the full-GHEM curves; (ii) The saturation intensity, black curve
labeled $I_{\mathrm{sat}}$ in panel (c), decreases over a hundredfold as
$\mu_{c}=0\rightarrow\hbar\omega/2$ and then increases above that; (iii) The
refractive part of the photoconductivity,
$\mathrm{Im}\\{\Delta\sigma^{(1)}\\}$, is negative except near half-photon
energy, where is goes positive.

Figure 9: Effect of the saturation carrier density, $n_{\mathrm{sat}}$, and
the energy of Coulomb impurities, $\mathcal{E}_{\mathrm{imp}}$, on the (a,c)
real and (b,d) imaginary parts of graphene surface conductivity. In panels
(a,b) curve color and line-style denote different $n_{\mathrm{sat}}$ and
$\mathcal{E}_{\mathrm{imp}}$, respectively, while the convention is opposite
in panels (c,d). In all cases: $|\mu_{c}|=0.3$ eV, $T_{0}=300$ K,
$\lambda_{0}=1550$ nm, $\tau_{\mathcal{E}}=1$ ps, and $\tau_{\mathrm{rec}}=10$
ps.

Closing the study of the steady-state solutions of the GHEM, we present the
saturation-carrier density effect on the photoconductivity. We employ the
simple empirical form
$G_{\mathrm{sat}}(n_{\mathrm{PG}})=(1-n_{\mathrm{PG}}/n_{\mathrm{sat}})$ for
the saturation factor in the generation term of Eq. (11), as proposed in [36].
$n_{\mathrm{sat}}$ is the saturation carrier density that $n_{\mathrm{PG}}$
cannot surpass, which would be proportional to the equilibrium carrier
density, i.e., in the order of few $10^{16}$/m2 for the
$\mu_{c}<\hbar\omega/2$ studied in this work. In parallel to
$n_{\mathrm{sat}}$, we also study the effect of the monolayer quality through
the equivalent Coulomb-impurity energy, $\mathcal{E}_{\mathrm{imp}}$, which
governs the energy-dependence of the intraband momentum relaxation lifetime,
$\Gamma_{i}$; for more details refer to Appendix B.1 and Eq. (27).

From Fig. 9(a,b) we see that $n_{\mathrm{sat}}>10^{16}$/m2 results in an
overall quenching of the photoconductivity. The residual
$\mathrm{Re}\\{\sigma^{(1)}\\}$ at high intensities can be fit to non-
saturable losses observed in recent experiments [37] and thus estimate
$n_{\mathrm{sat}}$. Concerning $\mathrm{Im}\\{\sigma^{(1)}\\}$, the carrier-
density saturation seems to quench only the negative feature that peaks near
$I_{\mathrm{sat}}$, leaving the positive up-shooting mostly unimpeded. This
nonlinear refraction trend is also in line with recent experiments [36, 70],
where spectral broadening of negatively chirped pulses was associated with
$\mathrm{Im}\\{\sigma^{(1)}\\}>0$ (under the $e^{-i\omega t}$ phase-
convention) and moreover scaling with the square-root of the carrier density.
Finally, note that some of the curves in Fig. 9 do not extend all the way to
$+40$ dB. We found an unphysical re-rising of $\mathrm{Re}\\{\sigma^{(1)}\\}$
after a threshold intensity, coinciding with a significant deviation of the QE
chemical potentials from their corresponding ‘hot’ counterparts. This anomaly
can be corrected by increasing the lattice temperature ($T_{0}$) after the
real part ceases to decrease, which leads us to believe that photothermal
effects must be accounted for in very high CW intensities, to further extend
the validity range of the GHEM; refer to relevant discussion in Section V.3.
Finally, we investigate the effect of monolayer sample quality through
$\mathcal{E}_{\mathrm{imp}}$, whose lower values correspond to higher carrier
mobility, hence higher quality. From Fig. 9(c), we observe that its effect is
meaningful only for very low $n_{\mathrm{sat}}$, where lower quality
corresponds to lower $I_{\mathrm{sat}}$; note that for
$\mathcal{E}_{\mathrm{imp}}\rightarrow 0$, the SA vanishes. The effect of
sample quality is particularly pronounced on the refractive part of the
surface conductivity, Fig. 9(b,d): An increase in $\mathcal{E}_{\mathrm{imp}}$
(quality decrease) contributes to orders of magnitude decrease in the
intensity threshold where $\mathrm{Im}\\{\Delta\sigma^{(1)}\\}$ crosses from
negative to positive. Evidently, this graphene parameter can have significant
impact on applications relying on refractive nonlinearities, e.g., pulse
shaping, spectral broadening, etc.

###  IV.3 Nonperturbative Transient Response

We now proceed to the transient solution of the GHEM equations (refer to
Appendix C for implementation details) in the retarded time-frame of a
Gaussian-enveloped ps-pulse, with effective intensity
$I_{\mathrm{eff}}(t)=\exp[-2c_{SG}(t/\Delta t_{\mathrm{FHWM}})^{2N_{SG}}]$,
where $c_{SG}=2^{(2N_{SG}-1)}\ln{2}$ is a duration normalization constant
(equal to 1.3863 for regular Gaussian pulse of order $N_{SG}=1$) and $\Delta
t_{\mathrm{FWHM}}$ is the full-width at half-maximum of the pulse power. We
hereby investigate the effect of the controllable GHEM parameters: In Fig. 10
panels (a)-(c), the effect of pulse peak intensity $I_{\mathrm{eff,peak}}$, in
panels (d)-(f), the effect of $\Delta t_{\mathrm{FWHM}}$, and, in panels
(g)-(i), the effect of chemical potential at the absence of illumination,
$\mu_{c}$. When not varied, default values are $I_{\mathrm{eff,peak}}=+30$ dB
vs. $1$ MW/cm2, $\Delta t_{\mathrm{FWHM}}=1$ ps, and $\mu_{c}=-0.3$ eV. The
columns-of-panels of Fig. 10 correspond to the ones in Fig. 7: panels (a,d,g)
depict $\\{n_{GP},\Delta T\\}$ vs. time, panels (b,e,h) depict $\mu_{e,h}$ vs.
time, and panels (c,f,i) depict $\mathrm{Re}\\{\sigma^{(1)}\\}$ and
$\mathrm{Im}\\{\Delta\sigma^{(1)}\\}$ vs. time. In all cases $T_{0}=300$ K,
$\lambda_{0}=1550$ nm, $\tau_{\mathcal{E}}=1$ ps, and $\tau_{\mathrm{rec}}=10$
ps, to which the time axis is normalized.

Figure 10: Transient GHEM response of (a,d,g) photogenerated carrier density
and temperature shift, (b,e,h) quasi Fermi levels, (c,f,i) surface
conductivity. We vary (a,b,c) the impinging intensity, (d,e,f) the pulse
duration, and (g,h,i) the equilibrium chemical potential; default values are
$I_{\mathrm{eff,peak}}=+30$ dB vs. $1$ MW/cm2, $\Delta t_{\mathrm{FWHM}}=1$
ps, and $\mu_{c}=-0.3$ eV, respectively. In all cases $T_{0}=300$ K,
$\lambda_{0}=1550$ nm, $\tau_{\mathcal{E}}=1$ ps, and $\tau_{\mathrm{rec}}=10$
ps, to which the time axis is normalized.

Analyzing the results in Fig. 10, some interesting conclusions can be drawn,
valid for all cases, i.e., when varying either $I_{\mathrm{eff,peak}}$,
$\Delta t_{\mathrm{FWHM}}$ or $\mu_{c}$: (i) The nonlinear photoconductivity
correlates best with $n_{\mathrm{PG}}$, rather than $T$ or $\mu$, implying
that this GHEM parameter is the most important to qualitatively capture the
temporal dynamics. (ii) The $\mathrm{Im}\\{\Delta\sigma^{(1)}\\}$ shows some
rippling around the pulse peak, as expected by the non-monotonic CW curves as
intensity increases, Fig. 7(c); this rippling would impart a chirp on a pulse
that propagates along the monolayer. (iii) The carrier temperature decays much
slower than $n_{\mathrm{PG}}$, despite
$\tau_{\mathcal{E}}=0.1\tau_{\mathrm{rec}}$; this apparent contradiction is
due to the factor $\partial T/\partial\mathcal{E}^{D}_{T}$ in the RHS of Eq.
(10) [and Eq. (22)] that increases the lifetime. We performed a consistency
check on the quantity $\mathcal{E}^{D}_{T}(t)-\mathcal{E}^{D}_{\mathrm{QE}}$,
as in Eq. (1) to which the faster rate $\tau_{\mathcal{E}}^{-1}$ is applied,
which indeed decays faster than $n_{\mathrm{PG}}$.

Some more interesting observations can be extracted from Fig. 10: In panel
(d), we observe that peak $\Delta T$ is proportional to pulse duration, but
peak $n_{\mathrm{PG}}$ is not affected in this range of pulse duration,
$\Delta t_{\mathrm{FWHM}}=0.5$-$2$ ps. From panel (g), we observe that peak
absolute $n_{\mathrm{PG}}$ increases with decreasing $|\mu_{c}|$, which
consequently means that the difference in the normalized
$n_{\mathrm{PG}}/n_{T_{0}}$ would be even higher; this apparent contradiction
can be explained from Fig. 7(a), comparing the $n_{\mathrm{PG}}$ curves for
$\mu_{c}=-0.3$ and $+0.1$ eV at the $+30$ dB power level considered here.
Finally, in panel (i), we note the only case of
$\mathrm{Im}\\{\Delta\sigma^{(1)}\\}>0$, which happens for
$|\mu_{c}|=\hbar\omega/2=0.4$ eV as expected from Fig. 8(d).

We now proceed to the study of the effect of the saturation carrier density,
$n_{\mathrm{sat}}$, on the transient GHEM variables. We now consider a lower
quality monolayer, $\mathcal{E}_{\mathrm{imp}}=100$ meV corresponding to
carrier mobility of about 630 cm2/(Vs) on an air/oxide interface, and a longer
Gaussian pulse duration of $\Delta t_{\mathrm{FWHM}}=3$ ps. Figure 11(a) shows
that $n_{\mathrm{PG}}\rightarrow n_{\mathrm{sat}}$ in all cases on pulse-peak
and that the intraband absorption caused by the high impurity density leads to
a large increase in carrier temperature, over 1000 K. In panel (b), we see
that the rippling of the quasi Fermi levels on pulse-peak translates to a
similar rippling in the refractive (imaginary) part of the surface
conductivity, panel (c), where we also observe that low $n_{\mathrm{sat}}$
values dramatically quench the attainable SA. In panel (b) we also plot the
quasi Fermi levels at the QE (quasi-equilibrium) state, $\mu_{e_{0},h_{0}}$,
with thin dot/dash-dot lines; these differ visibly from the corresponding hot
potentials and have a slower response.

Figure 11: Effect of saturation carrier density $n_{\mathrm{sat}}$ on (a)
$n_{\mathrm{PG}}$ and $\Delta T$, (b) $\mu_{e,h}$ and $\mu_{e_{0},h_{0}}$,
with thick and thin curves, respectively, and (c) $\sigma^{(1)}/\sigma_{0}$.
Default GHEM parameters, i.e., as in Fig. 10, except for $\Delta
t_{\mathrm{FWHM}}=3$ ps and $\mathcal{E}_{\mathrm{imp}}=100$ meV, which
corresponds to lower quality monolayer.

Finally, in Fig. 12, we investigate the effect of the monolayer impurity
energy ($\mathcal{E}_{\mathrm{imp}}\propto\sqrt{N_{\mathrm{imp}}}$, where
$N_{\mathrm{imp}}$ is the effective impurity density) on the
photoconductivity, when the equilibrium chemical potential is biased, panel
(a), exactly on half-photon energy or, panel (b), just above it;
$n_{\mathrm{sat}}=3\times 10^{16}$/m2 and the rest of the GHEM parameters are
the same as in Fig. 11. As expected from the static response curves, Fig. 9,
for $\mu_{c}=0.4$ eV we get $\mathrm{Im}\\{\Delta\sigma^{(1)}\\}>0$, whose
magnitude moreover increases with $\mathcal{E}_{\mathrm{imp}}$, i.e., as the
quality of the monolayer decreases; this feature is predominantly due to the
intraband mechanism and the high thermalization that the carrier-density
saturation brings. For the $\mu_{c}>\hbar\omega/2$ case, the overall
photoconductivity (nonlinearity) is lower: the refractive part exhibits a
positive peak followed by a sign-change with a long trailing negative
$\mathrm{Im}\\{\Delta\sigma^{(1)}\\}$, whereas the real part shows induced
absorption for high impurity densities.

Figure 12: Effect of impurity density $\mathcal{E}_{\mathrm{imp}}$ on surface
photoconductivity, for monolayers biased at (a) $\mu_{c}=\hbar\omega/2$ and
(b) $\mu_{c}=\hbar\omega/1.6$. Default GHEM parameters, i.e., as in Fig. 11,
and $n_{\mathrm{sat}}=3\times 10^{16}$/m2. The shaded area denotes the full
width at 10 dB ($I_{\mathrm{eff}}(t)>100$ MW/cm2).

##  V Discussion

###  V.1 Implications for Photonic Waveguides

Having extensively analyzed the GHEM response per unit length of the waveguide
one naturally wonders how the transient response evolves as the pulse travels
down a waveguide. This effect can be modeled in terms of the nonlinear
Schrödinger equation (NLSE) framework, and numerically solved using the split-
step Fourier method (SSFM); brief introduction for both can be found in
Appendix D.

With these remarks in mind, summarizing the GHEM photoconductivity results
from Section IV, i.e., for $\mu_{c}$ near and slightly below $\hbar\omega/2$,
we conclude that graphene exhibits saturable absorption (SA),
$\mathrm{Re}\\{\Delta\sigma^{(1)}\\}<0$, and positive self-focusing
refraction, $\mathrm{Im}\\{\Delta\sigma^{(1)}\\}<0$, for low and medium
powers. Then, above an intensity threshold, the refraction switches to
negative (defocusing), which moreover coincides with deep SA. The increasing
quality (carrier mobility) of graphene was found to increase the SA depth,
together with the intensity threshold for refraction sign-flipping. Forcing an
upper boundary to photogenerated carrier density decreases the magnitude of
nonlinearity (SA and refraction), without affecting the threshold intensities.
Finally, the model predicts a strong coupling between the absorptive and
refractive nonlinear regimes, with nontrivial boundaries between them.

Consequently, depending on the nonlinear application targeted (i.e., whether
it relies on the absorptive or refractive part of $\sigma^{(1)}$ – or to
both), the graphene-related parameters must be judiciously chosen; these
mainly include the graphene quality, its equilibrium chemical potential, and
the pulse peak power. The dependence of graphene’s nonlinear response on these
parameters is non-trivial and often non-monotonic, meaning that optimal
conditions might arise in banded windows within the parameter space. For
instance: absorption modulation (e.g. for extinction-ratio improvement of an
on/off-keying modulated bitstream) requires deep loss contrast, i.e., deep SA
with a steep slope in the $\mathrm{Re}\\{\sigma^{(1)}\\}$-intensity curve;
spectral broadening requires low losses (or deep SA) together with high
$|\mathrm{Im}\\{\sigma^{(1)}\\}|$; bright-soliton formation requires low
losses (or deep SA) together with
$\mathrm{Re}\\{n_{2}\\}\propto\mathrm{Im}\\{\sigma^{(1)}\\}$ of opposite sign
to the group velocity dispersion coefficient ($\beta_{2}$). In most of these
applications, the optimal chemical potential would be close to half-photon
energy, while the actual value (above or below it) would stem from specific
compromises between its coupled refractive and absorptive parts.

###  V.2 Qualitative Trends and Comparison to Experiments

The dominant nonlinearity in the $\mathrm{Re}\\{\Delta\sigma^{(1)}\\}$ is SA,
with a minimal $I_{\mathrm{sat}}<10$ MW/cm2, found slightly below half-photon
energy. This observation qualitatively agrees both with free-space experiments
and with pulsed waveguide measurements [37]. In the latter, a relatively
shallow SA was measured, corresponding to high non-saturable (residual)
conductivity, with negligible refractive nonlinearity. This behaviour can be
explained by a combination of low carrier mobility and/or low carrier-
saturation density in the framework of this GHEM. Moreover, fitting of the
model parameters to experimentally characterized samples could also improve
quantitative accuracy when optimizing such devices.

The refractive part of the nonlinear photoconductivity,
$\mathrm{Im}\\{\Delta\sigma^{(1)}\\}$, constitutes a rather interesting and
relatively unexplored regime. Our GHEM qualitatively agrees with recent
experiments in silicon and silicon nitride graphene-comprising waveguides [33,
36], i.e., it predicts a positive $\mathrm{Im}\\{\sigma^{(1)}\\}$ regime,
which increases with carrier density. Concerning the quantitative comparison
to the phenomenological model proposed in [36] and theoretically elaborated in
[70], our model predicts a larger threshold stemming from a different
dependence on the carrier density. Specifically, [36] proposes a direct
proportionality $\mathrm{Im}\\{\Delta\sigma^{(1)}\\}\propto\sqrt{n_{T}}$ in
order to explain the superlinear (exponential) rising in spectral-broadening
that was measured, as the waveguide length is increased. Our GHEM predicts a
spectral broadening, but with a different slope (sublinear, e.g., logarithmic)
or at a higher intensity level.

Finally, our GHEM shows that both real and imaginary parts of transient
$\Delta\sigma^{(1)}(t)$ correlate mostly with the $n_{\mathrm{PG}}(t)$. This
means that the carrier density is what predominantly defines the nonlinear
response, and validates the phenomenological models proposed in [36, 37], that
only track the carrier-density, and not the carrier temperature and chemical
potentials. Nevertheless, we stress that using only the interband rate
equation, e.g., Eq. (2), coupled to the pulse or in static regime, produces
results that do not match the full six-equation GHEM response nor the
experimental trends. In Fig. 13 we compare the full transient GHEM response,
for the three pulse-widths of Fig. 10(f), with a ‘look-up’ (LU) approach using
the static saturation curves, e.g., Fig. 8. The dotted and dashed curves in
Fig. 13 is the photoconductivity when the look-up is based on the pulse’s
$I_{\mathrm{eff}}(t)$ or the resulting $n_{\mathrm{PG}}(t)$, respectively. The
dotted curves fail to capture the edges of the pulse (particularly the
trailing edge), while the dashed curves attain a near-perfect fit, with slight
difference only on pulse-peak.

Figure 13:  Comparison of full transient photoconductivity response against
two ‘look-up’ (LU) strategies based on static photoconductivity. Panels
(a,b,c) correspond to pulse duration $\Delta t_{\mathrm{FWHM}}=(0.5,1,2)$ ps.
Thin solid curves are for the full transient GHEM; thick dotted curves are for
look-up using directly the pulse intensity; thick dashed curves are for look-
up using the GHEM-computed photogenerated carrier density.

###  V.3 Outlook

The quantitative divergences between the model and experiments discussed in
the previous subsection hints that more physics should complement the GHEM. In
this way, its validity range can be extended to higher effective intensities,
such as the ones delivered by high-power fiber lasers into highly-confining
waveguides. Photothermal effects, i.e., lattice heating and heat diffusion,
provide a preliminary correction to abnormal induced absorption observed at
very high intensities, particularly in the carrier saturation regime. Thus,
the lattice temperature, $T_{0}$, can be a problem variable complemented in
the GHEM by a differential equation related to Joule heat and its diffusion.
Additionally, we theorize that diffusion of carriers [43], with a possible
dependence of the diffusion coefficient on the carrier temperature according
to the Einstein relation, $D_{\mathrm{diff}}=\mu_{\mathrm{mob}}k_{B}T$, could
also improve agreement with experimental observations.

This model can be directly used in the MIR/THz spectral region, where the
intraband Drude-like contribution dominates the photoconductivity via thermal
effects. This regime can be readily implemented within this GHEM, by forcing
an appropriately low carrier-density saturation.

Finally, spectral bandwidth studies and/or multi-channel (e.g., pump-probe)
effects can be readily incorporated in this GHEM, with the Kubo formulas
accounting for the spectra of the model parameters.

##  VI Conclusive Remarks

In summary, we have developed an electrodynamic model for the study of the
ultrafast absorptive and refractive transient response of graphene monolayers
placed along integrated photonic waveguides. This strongly non-perturbative
model agrees with quasi-classical perturbative derivations of third-order
nonlinearity in the area of interest, i.e., near and slightly below the half-
photon energy, but only at low intensities. For higher intensities, the
qualitative agreement with experimental observations, in both static and
transient regimes, is good and the model’s phenomenological parameters can be
fit to measurements for quantitative analysis and design. Apart from the
illumination parameters, such as the pulse duration and peak intensity,
emphasis was given to the response sensitivity on voltage-tunability (via the
equilibrium chemical potential), on the sample’s quality (via the impurity
density), and on the existence of a carrier-density saturation regime. For
most photonic applications, the optimal chemical potential would be close to
half-photon energy, while the actual compromise between the coupled absorptive
and refractive nonlinearity shift the value slightly above or below it.

##  Appendix A Fermi-Dirac Framework

The graphene hot-electron model (GHEM) adopted from [52] relies on distinct
quasi-Fermi levels (or chemical potentials) for the electron and hole plasmas,
$\mu_{e}$ and $\mu_{h}$, respectively. In general it holds that
$\mu_{e}\neq\mu_{h}$, and the difference becomes more pronounced as the system
is pushed farther from thermal equilibrium especially by photogeneration
(interband absorption). Both plasmas have the same carrier temperature, $T$,
which can (far) surpass the lattice temperature $T_{0}$, if even for a sub-ps
timespan. In the Fermi-Dirac statistical framework, the distribution function
for electrons/holes is given by

$f_{e/h}(\mathcal{E};\mu_{e/h},T)=\frac{H(\pm\mathcal{E})}{1+\exp\left(\pm\frac{\mathcal{E}-\mu_{e/h}}{k_{B}T}\right)},$
(15)

where $k_{B}$ is the Boltzmann constant and $H(x)$ is the step function ($H=1$
for $x>0$, else $H=0$). The values $f_{e/h}(\mathcal{E})$ are color-coded
along the vertical (energy) direction of the Dirac cones presented in Fig. 2,
where the conduction/valence band corresponds to $\mathcal{E}\gtrless 0$ for
this zero-bandgap semiconductor.

The solid-state parameters on which the GHEM finally relies are the carrier
surface densities ($n_{e,h}$ in 1/m2) and the plasma-energy surface densities
($\mathcal{E}^{D}_{e,h}$ in J/m2). Note that, like the distribution functions,
these are distinct for electrons ($e$-subscript) and holes ($h$-subscript) and
depend solely on the corresponding chemical potential ($\mu_{e,h}$) and the
common carrier temperature. The formulas for the carrier and energy densities
are derived from energy integrals of $f_{e/h}(\mathcal{E})$ times the density
of states which, for graphene near the tip of the Dirac cone, takes the linear
form

$N_{\mathrm{DOS}}(\mathcal{E})=\frac{2}{\pi(v_{F}\hbar)^{2}}|\mathcal{E}|,$
(16)

where $\hbar$ is the reduced Plank constant ($\hbar=h/2\pi$) and $v_{F}\approx
c_{0}/300$ is the Fermi velocity in graphene ($c_{0}$ is the speed of light in
vacuum); formally, $v_{F}=\alpha_{0}\gamma_{0}\sqrt{3}/2\hbar\approx
0.874\times 10^{6}$ m/s, for a lattice constant of $a_{0}=2.46$ Å and a
nearest-neighbour coupling energy of $\gamma_{0}=2.7$ eV [26, 71]. The
quantities $\\{n,\mathcal{E}^{D}\\}$ are given by the following expressions:
$n=\int_{0}^{\infty}N_{\mathrm{DOS}}fd\mathcal{E}$ and
$\mathcal{E}^{D}=\int_{0}^{\infty}N_{\mathrm{DOS}}f\mathcal{E}d\mathcal{E}$,
resulting in:

$n_{e/h}(\mu_{e/h},T)=\frac{2(k_{B}T)^{2}}{\pi(\hbar
v_{F})^{2}}F_{1}\left(\pm\frac{\mu_{e/h}}{k_{B}T}\right),$ (17)
$\mathcal{E}^{D}_{e/h}(\mu_{e/h},T)=\frac{2(k_{B}T)^{3}}{\pi(\hbar
v_{F})^{2}}F_{2}\left(\pm\frac{\mu_{e/h}}{k_{B}T}\right).$ (18)

Note the symmetry in these formulas for electrons and holes, with an attention
to the sign used for electrons ($+$) and holes ($-$). The function $F_{m}(x)$
is the Fermi-Dirac integral (FDI) of order-$m$ defined as

$F_{m}(x)=\int_{0}^{\infty}\frac{u^{m}}{1+\exp\left(u-x\right)}du.$ (19)

The total carrier and energy densities in a given state (e.g., at thermal
equilibrium or quasi-equilibrium) are given by the sum of the electron and
hole terms, $n_{T}=n_{e}+n_{h}$ and
$\mathcal{E}^{D}_{T}=\mathcal{E}^{D}_{e}+\mathcal{E}^{D}_{h}$, respectively.

Finally, we stress that Eq. (17) and (18) are (numerically) invertible, which
means that knowing any two variables in the $\\{\mu,T,n\\}$ or
$\\{\mu,T,\mathcal{E}^{D}\\}$ set, we can calculate the third variable. In the
GHEM developed in this work, we specifically require the inversion of Eq.
(17), for the calculation of the chemical potential when the temperature and
carrier density are known,

$\mu_{e/h}=\pm(k_{B}T)F^{-1}_{1}\left[\frac{\pi(\hbar
v_{F})^{2}}{2(k_{B}T)^{2}}n_{e/h}\right],$ (20)

where $F^{-1}_{1}$ is the inverse function of Eq. (19) for $m=1$. Note that
asymptotic expressions can simplify the numerical inversion in extreme cases,
namely $F_{1}(x)=x^{2}/2$ and $F_{1}(x)=e^{x}$, when $x\gg 1$ and
$x\rightarrow 0$, respectively. Finally, fittings exist for the inverted
calculation, i.e., extracting $\\{\mu,T\\}$ from $\\{n_{e},n_{h}\\}$, such as
the one proposed in [70] and its supporting information, applicable in the
equilibrium case, $\mu_{e}\equiv\mu_{h}=\mu$.

###  Temperature Dependence of Plasma Energy Density

The derivative of the FDI in Eq. 19 is defined as

$\frac{dF_{m}(x)}{dx}=\frac{\Gamma_{F}(m+1)}{\Gamma_{F}(m)}F_{m-1}(x),$ (21)

where $\Gamma_{F}$ is the Gamma function, with $\Gamma_{F}(2)=1$ and
$\Gamma_{F}(3)=2$ for the cases of interest. This FDI property is useful in
evaluating the dependence of the total plasma energy density on carrier
temperature. Defining $b=2/(\pi\hbar^{2}v_{F}^{2})$ and
$x_{e,h}=\mu_{e,h}/k_{B}T$ for shorthand, and using Eq. (17) and (18), we
have:

$\displaystyle\frac{\partial\mathcal{E}^{D}_{T}}{\partial T}$
$\displaystyle=\frac{\partial}{\partial
T}\Big{\\{}b(k_{B}T)^{3}\left[F_{2}(x_{e})+F_{2}(-x_{h})\right]\Big{\\}}$ (22)
$\displaystyle=\frac{3}{T}\mathcal{E}^{D}_{T}+\frac{2b(k_{B}T)^{3}}{k_{B}T^{2}}\Big{[}-\mu_{e}F_{1}(x_{e})+\mu_{h}F_{1}(-x_{h})\Big{]}$
$\displaystyle=\frac{1}{T}\left(3\mathcal{E}^{D}_{T}-\mu_{e}n_{e}+\mu_{h}n_{h}\right).$

##  Appendix B Surface Conductivity Calculation

For finite (nonzero) carrier temperatures, the formulas for the intraband
($i$-subscript) and interband ($e$-subscript) surface conductivities entail
integration over the energy spectrum [72]. The full expressions are given
below, for the quasi-equilibrium case where the quasi-Fermi levels for
electrons and holes can be unequal [52], $\mu_{e}\neq\mu_{h}$, and assuming
arbitrary energy-dependent scattering rates
$\Gamma_{i,e}=\Gamma_{i,e}(\mathcal{E})$, as follows

$\displaystyle\sigma^{(1)}_{i}$
$\displaystyle(\omega,\mu_{e},\mu_{h},T)=\sigma_{0}\frac{i}{\pi
k_{B}T}\int_{0}^{\infty}\frac{\mathcal{E}}{\hbar\omega+i\Gamma_{i}(\mathcal{E})}$
(23)
$\displaystyle\times\bigg{[}{\cosh^{-2}\left(\frac{\mathcal{E}-\mu_{e}}{2k_{B}T}\right)}+{\cosh^{-2}\left(\frac{\mathcal{E}+\mu_{h}}{2k_{B}T}\right)}\bigg{]}d\mathcal{E},$

and

$\displaystyle\sigma^{(1)}_{e}$
$\displaystyle(\omega,\mu_{e},\mu_{h},T)=\sigma_{0}\frac{4i}{\pi}\int_{0}^{\infty}\frac{\hbar\omega+i\Gamma_{e}(\mathcal{E})}{[\hbar\omega+i\Gamma_{e}(\mathcal{E})]^{2}-4\mathcal{E}^{2}}$
(24)
$\displaystyle\times\left[\frac{1}{1+\exp\left(\frac{-\mathcal{E}-\mu_{h}}{k_{B}T}\right)}-\frac{1}{1+\exp\left(\frac{\mathcal{E}-\mu_{e}}{k_{B}T}\right)}\right]d\mathcal{E}.$

In these expressions, the $e^{-i\omega t}$ phase-convention is used,
$\sigma_{0}=q^{2}/4\hbar\approx 61$ $\mu$S is the universal optical
conductivity of graphene, and $q$ is the (absolute) electron charge. Note that
$\sigma_{0}$ is responsible for the $2.3$ % absorption through an air-
suspended pristine graphene monolayer, as $0.023\approx
1-4/(2+\sigma_{0}Z_{0})^{2}$, where $Z_{0}\approx 377~{}\Omega$ is the free
space impedance, according to Eq. (4). For specific information regarding the
calculation of the integrals in the $i$\- and $e$-conductivities, refer to the
dedicated subsections, below.

The integrals in Eqs. (23) and (24) can be cast in simple closed-form
expressions in the zero temperature limit. However, at room temperatures and
above (as is the case in this non-perturbative model), the hyperbolic and
exponential functions in square brackets in the integrands, will spread across
wider energy bands and thus demand for more careful considerations, even for
numerical integration.

###  B.1 Intraband Conductivity

The integral in Eq. (23) must be numerically computed in the case where
$\Gamma_{i}$ is strongly energy dependent and/or at high carrier temperatures
arising when graphene is absorptive ($\mu_{c}\ll\hbar\omega/2$) and
illuminated with intense optical radiation. However, the integral can be
analytically computed when $\Gamma_{i}$ is assumed constant (energy
independent), which is a valid assumption at low carrier temperatures. In this
case, the energy spreading of the Fermi-Dirac-like function in the integrand
of Eq. (23) is low, i.e., only the part around $\mathcal{E}\approx|\mu_{c}|$
contributes to the final value. So, $\Gamma_{i}\approx\Gamma_{i}(|\mu_{c}|)$
can be assumed constant and the resulting closed-form expression is

$\sigma^{(1)}_{i}=\sigma_{0}\frac{4i}{\pi}\frac{k_{B}T}{(\hbar\omega+i\Gamma_{i})}\ln\left[\left(1+e^{\frac{\mu_{e}}{k_{B}T}}\right)\left(1+e^{\frac{-\mu_{h}}{k_{B}T}}\right)\right].$
(25)

Moreover, in thermal equilibrium ($\mu_{e}\equiv\mu_{h}=\mu_{c}$), this
expression can be further reduced to the more plainly-evident Drude-like form
[31]

$\sigma^{(1)}_{i}=\sigma_{0}\frac{4i}{\pi}\frac{\mu_{c}}{(\hbar\omega+i\Gamma_{i})}\mathcal{T}\left(\frac{\mu_{c}}{2k_{B}T}\right),$
(26)

where the temperature-dependent function $\mathcal{T}(x)=x^{-1}\ln[2\cosh(x)]$
tends to $\pm 1$ for $|x|>1$, i.e., for $|\mu_{c}|>2k_{B}T\approx 52$ meV for
room temperature. From Eq. (26) and for the wavelength of interest,
$\hbar\omega\approx 0.8$ eV (NIR), it can be readily seen that both real
(absorptive) and imaginary (refractive) parts of $\sigma^{(1)}_{i}$ are
proportional to $|\mu_{c}|$ while only the real part is significantly affected
by $\Gamma_{i}$, typically smaller than 0.1 eV ($\tau_{i}>6$ fs): As rate
$\Gamma_{i}$ increases ($\tau_{i}$ drops) the carriers’ momentum relaxation is
more efficient thus graphene is less conductive (absorptive), i.e.,
$\mathrm{Re}\\{\sigma^{(1)}_{i}\\}$ diminishes.

In this work, we adopt the closed-form model proposed by S.A. Mikhailov in
[52] (and Appendix therein) for the energy-dependence of the intraband
scattering rate, which captures the behaviour both at low energies and at high
energies, where charged-impurity scattering dominates [73, 74, 75]. The
proposed formula is

$\Gamma_{i}(\mathcal{E})=\frac{|\mathcal{E}|}{\zeta/2-1+\sqrt{1+(\mathcal{E}/\mathcal{E}_{\mathrm{imp}})^{4}}}$
(27)

where $\zeta>2$ is the minimal static conductivity of graphene (at the Dirac
point, i.e., when $\\{T,\omega,\mu_{c}\\}\rightarrow 0$) in units of
$\sigma_{0}$ and $\mathcal{E}_{\mathrm{imp}}$ is a Coulomb energy associated
with the density of impurities in the lattice. $\mathcal{E}_{\mathrm{imp}}$
scales with the square-root of the impurity density (typically in the
$0.1$-$10\times 10^{16}/$m2 range) with a proportionality factor depending on
the dielectric properties of the medium in contact with graphene; in the same
framework it is shown that the low frequency mobility of graphene carriers is
inversely proportional to $\mathcal{E}_{\mathrm{imp}}^{2}$. In any case, both
model parameters, $\zeta$ and $\mathcal{E}_{\mathrm{imp}}$, can be fitted to
measurements of static conductivity vs. electrically tuned carrier density at
low temperatures, or extrapolated from room-temperature measurements. The
energy dependence of rate $\Gamma_{i}$ on $\zeta$ and
$\mathcal{E}_{\mathrm{imp}}$ is depicted in Fig. 14, where the comparatively
stronger effect of $\mathcal{E}_{\mathrm{imp}}$ is revealed together with the
inverse trend between low energies ($\Gamma_{i}\propto\mathcal{E}$) and high
energies ($\Gamma_{i}\propto\mathcal{E}^{-1}$). Note that the quality of a
graphene sample increases as its mobility increases, which is inversely
proportional to $\mathcal{E}_{\mathrm{imp}}$. In this work, unless otherwise
specified, we assume typical values $\zeta=4$ and
$\mathcal{E}_{\mathrm{imp}}=30$ meV, corresponding to good quality graphene,
with a mobility of about 7000 cm2/(Vs) when the graphene monolayer lies on the
interface between air and a $\kappa\approx 4$ dielectric substrate.

Figure 14: Energy dependence of the intraband momentum-relaxation rate
$\Gamma_{i}$, (a) for various minimal static conductivities $\zeta=2$ to 10
and fixed Coulomb impurity energy $\mathcal{E}_{\mathrm{imp}}=30$ meV, and (b)
for various $\mathcal{E}_{\mathrm{imp}}=10$ to 100 meV and fixed $\zeta=4$.
Relaxation lifetime: $\tau_{\mathrm{[fs]}}\approx
658/\Gamma_{\mathrm{[meV]}}$.

###  B.2 Interband Conductivity

Unlike the intraband term, the integral in Eq. (24) cannot be analytically
solved, even though in most practical cases the interband momentum relaxation
rate is energy independent and it can even be neglected,
$\Gamma_{e}\rightarrow 0$, under the condition $\Gamma_{e}\ll
k_{B}T_{0}\approx 26$ meV (at room temperature); in most practical cases
$\Gamma_{e}<1$ meV ($\tau_{e}>1$ ps) [76]. This complication is due to the
singularity exhibited by the integrand at half-photon energy,
$\mathcal{E}=\hbar\omega/2$, which can fortunately be circumvented by a
transformation involving a principal value integral. This procedure has been
outlined in [59] and, for the out-of-equilibrium case studied here,
$\mu_{e}\neq\mu_{h}$, it can be extended as follows: We first define the
auxiliary function

$G(\mathcal{E};\mu_{e},\mu_{h},T)=\frac{\sinh\left(\frac{2\mathcal{E}-\mu_{e}+\mu_{h}}{2k_{B}T}\right)}{\cosh\left(\frac{\mu_{h}+\mu_{e}}{2k_{B}T}\right)+\cosh\left(\frac{2\mathcal{E}-\mu_{e}+\mu_{h}}{2k_{B}T}\right)}$
(28)

which can be used in compactly rewriting Eq. (24) as

$\sigma^{(1)}_{e}=\sigma_{0}\frac{4i\Omega}{\pi}\int_{0}^{\infty}\frac{G(\mathcal{E})}{\Omega^{2}-4\mathcal{E}^{2}}d\mathcal{E},$
(29)

where $\Omega=\hbar\omega+i\Gamma_{e}$ (the scattering rate is assumed energy
independent) and the singularity at $\Omega=2\mathcal{E}$ is evident in the
denominator. Adding and subtracting the term $G(\Omega/2)$ in the nominator of
the integrand in Eq. (29), we get one singularity-free term [from the
$-G(\Omega/2)+G(\mathcal{E})$ terms], that can be straightforwardly
numerically computed, and one term that requires a principal-value integral
(PVI) [from the $+G(\Omega/2)$ term]. Now, as the integrand function
[proportional to $1/(\Omega^{2}-4\mathcal{E}^{2})$] is anti-symmetric around
the singularity, the PVI reduces to the proportionality constant times $i\pi$.
The resulting expression for the numerically integrable interband conductivity
is:

$\sigma^{(1)}_{e}=\sigma_{0}\left[G(\Omega/2)+\frac{4i}{\pi}\int_{0}^{\infty}\frac{G(\mathcal{E})-G(\Omega/2)}{\Omega^{2}-4\mathcal{E}^{2}}d\mathcal{E}\right].$
(30)

Assuming $\Gamma_{e}=0$, Eq. (30) reveals that the real and imaginary parts of
the interband conductivity are solely defined by the first and second terms,
respectively. Note that the real part can acquire negative values (‘gain’)
owing to population inversion in strongly non-equilibrium states [52], i.e.,
when $\Delta\mu_{(e-h)}>\hbar\omega$; nevertheless, we restrict our study the
cases where this regime is not entered, ensuring always that the real part of
the total conductivity is positive.

##  Appendix C Numerical Solution of Equations System

The GHEM developed relies on a set of differential/algebraic equations (DAE).
In its most complicated transient form, the DAE systems consists of two
differential equations, (10) and (11) or (1) and (2), and a set of four
electroneutrality Eqs. (5). The first set of rate equations has
$\\{T,n_{\mathrm{PG}},\mu_{e},\mu_{h},\mu_{e_{0}},\mu_{h_{0}}\\}$ as unknowns
or variables, whereas the second set has
$\\{\mathcal{E}^{D}_{T},n_{\mathrm{PG}},\mu_{e},\mu_{h},\mu_{e_{0}},\mu_{h_{0}}\\}$,
just as in [52].

The number of unknown chemical potentials in the GHEM can be reduced in a
number of ways, (i) by assuming a single Fermi level for electrons and holes
in both the hot (out-of-equilibrium) and the QE state, i.e., $\mu_{h}=\mu_{e}$
and $\mu_{h_{0}}=\mu_{e_{0}}$, or (ii) by assuming the absence of QE state,
i.e., $\mu_{e_{0}}=\mu_{e}$ and $\mu_{h_{0}}=\mu_{h}$, or (iii) by a
combination of the previous two cases. This reduces the total equations to
four, four, or three, for cases (i)-(iii), respectively, while the carrier
temperature (or energy density) and photogenerated carrier density rate
equations cannot be dropped, for NIR photons. Even though the qualitative
behaviour of the GHEM equation system is equivalent in cases (i)-(iii), we
have found quantitative differences in the magnitude of the nonlinear effects
and in the intensity threshold between the various regimes. Throughout this
work we retain the full six-equation system.

Now, when time derivatives can be dropped, $\partial/\partial t\equiv 0$, we
have the static (steady-state) or CW (continuous wave) case, which can be
easily handled by MATLAB’s FSOLVE. In the more general transient case, e.g.,
when graphene is excited by a ps-pulse, the DAE can be solved with MATLAB’s
ODE15S, with an appropriate mass-matrix definition.

Apart from the six unknowns (variables) and possible time-dependence, one must
also define the input parameters (constants) of the system. The most important
is the peak effective illuminating intensity, $I_{\mathrm{eff}}$ in W/m2.
Apart from that, and the operating frequency, there is the various material-
parameters for graphene: chemical potential ($\mu_{c}$), lattice temperature
($T_{0}$), phenomenological intra- and interband lifetimes ($\tau_{\epsilon}$
and $\tau_{\mathrm{rec}}$), sample quality (e.g., impurity ‘energy’
$\mathcal{E}_{\mathrm{imp}}$ and residual conductivity $\zeta$), etc.

Finally, all variables and equations should be normalized so that they take
near-unity values, to aid numerical solution, either in static or transient
mode. The four electroneutrality equations, and the corresponding chemical
potentials ($\mu$) are normalized from Joules to electron-volts (eV). The rate
equations are normalized to unitless in static (CW) mode and in 1/s in
transient (pulsed) mode; in both cases, it is more numerically convenient to
compute a normalized temperature $T/T_{0}$ and photogenerated carrier density
$n_{\mathrm{PG}}/n_{T_{0}}$.

Refer to the Supplemental Material for MATLAB code implementing the GHEM
equation system and its static and transient solution.

##  Appendix D Pulse Propagation Along Nonlinear Waveguide

The initial pulse $B_{\mathrm{in}}(t)=B(0,t)$ is $z$-stepped through a segment
of the waveguide that is sufficiently short ($dz\rightarrow 0$) so that it
perturbatively alters the complex amplitude of the pulse envelope; this
$z$-stepping is essentially a transient solution of the GHEM, which produces a
photoconductivity $\Delta\sigma^{(1)}(z,t)$, assumed constant within
$[z,z+dz]$, which afterwards distorts the envelope at the end of the step; the
procedure is iteratively repeated until the end of the waveguide is reached,
applying linear effects (absorption and dispersion) along the way, until we
calculate $B(L_{z},t)$. This reflectionless step-wise propagation is the
essence of the split-step Fourier method (SSFM) [63]. Note that dispersion and
linear absorption depend on constant parameters derived from the linear regime
eigenmode analysis at $\omega_{0}$, including the effect of
$\sigma^{(1)}_{\mathrm{lin}}$ (primarily on absorption). In contrast,
photoconductivity $\Delta\sigma^{(1)}(z,t)$ is a nonlinear term, i.e., with
temporal dependence, which is moreover updated as the the propagation evolves
along $z$, since it is coupled to the pulse envelope, $B(z,t)$.

The NLSE modeling of linear and third-order nonlinear effects stemming from
graphene is well understood [31]. Now, the GHEM photoconductivity is
essentially a free-carrier effect that can be modeled as an extra term in the
NLSE [64, 15], $\delta_{\mathrm{GNL}}(z,t)$, together with the dispersion and
third-order nonlinearity, whose complex amplitude depends on the
spatiotemporal pulse amplitude $B(z,t)$, in a strongly nonlinear manner. A
basic single-channel NLSE for the optical pulse envelope $B(z,t)$ (in units of
$\sqrt{W}$) modulating a carrier frequency $\omega$, under the $e^{-i\omega
t}$ phase convention, is

$\frac{\partial B}{\partial
z}=\bigg{(}-\frac{\alpha}{2}+D+i\gamma_{\mathrm{NL}}|B|^{2}-\delta_{\mathrm{GNL}}(z,t)\bigg{)}B(z,t).$
(31)

In this compact form, the constants $\alpha$ and $\gamma_{\mathrm{NL}}$ are
the power loss coefficient (real positive) and the complex third-order
nonlinear parameter (including Kerr effect and perturbative SA/TPA),
respectively; self-focusing refraction corresponds to
$\mathrm{Re}\\{\gamma_{\mathrm{NL}}\\}>0$. Note that $\alpha$ includes a
contribution from graphene, through
$\mathrm{Re}\\{\sigma^{(1)}_{\mathrm{lin}}\\}$, but $\gamma_{\mathrm{NL}}$
comes exclusively from bulk/3D nonlinear materials and is totally unrelated to
graphene. $D$ is the linear dispersion operator, for group-velocity dispersion
and higher {e.g., check Eq. (4) in [45]}, and $t$ is the retarded-envelope
time frame, moving with the mode group velocity. Finally, the complex-valued
term $\delta_{\mathrm{GNL}}\propto\Delta\sigma^{(1)}$ includes all nonlinear
refractive and absorptive contributions from graphene. Notice how
$\delta_{\mathrm{GNL}}$ is added with a minus in Eq. (31), so that an
absorption saturation characterized by $\mathrm{Re}\\{\Delta\sigma^{(1)}\\}<0$
introduces a ‘gain’ that counteracts the linear loss factor $\alpha$. In the
same sense, the $\pm$ sign in $\mathrm{Im}\\{\Delta\sigma^{(1)}\\}$
corresponds to a defocusing ($+$) or self-focusing ($-$) refraction.

###### Acknowledgements.

The research work was supported by the Hellenic Foundation for Research and
Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to
support Faculty members and Researchers and the procurement of high-cost
research equipment grant.” (Project Number: HFRI-FM17-2086)

## References

  * Ferrari  _et al._ [2015] A. C. Ferrari, F. Bonaccorso, V. Fal'ko, K. S. Novoselov, S. Roche, P. Bøggild, S. Borini, F. H. L. Koppens, V. Palermo, N. Pugno, J. A. Garrido, R. Sordan, A. Bianco, L. Ballerini, M. Prato, E. Lidorikis, J. Kivioja, C. Marinelli, T. Ryhänen, A. Morpurgo, J. N. Coleman, V. Nicolosi, L. Colombo, A. Fert, M. Garcia-Hernandez, A. Bachtold, G. F. Schneider, F. Guinea, C. Dekker, M. Barbone, Z. Sun, C. Galiotis, A. N. Grigorenko, G. Konstantatos, A. Kis, M. Katsnelson, L. Vandersypen, A. Loiseau, V. Morandi, D. Neumaier, E. Treossi, V. Pellegrini, M. Polini, A. Tredicucci, G. M. Williams, B. H. Hong, J.-H. Ahn, J. M. Kim, H. Zirath, B. J. van Wees, H. van der Zant, L. Occhipinti, A. D. Matteo, I. A. Kinloch, T. Seyller, E. Quesnel, X. Feng, K. Teo, N. Rupesinghe, P. Hakonen, S. R. T. Neil, Q. Tannock, T. Löfwander, and J. Kinaret, Science and technology roadmap for graphene, related two-dimensional crystals, and hybrid systems, Nanoscale 7, 4598 (2015).
  * Cheng  _et al._ [2018] J. L. Cheng, J. E. Sipe, N. Vermeulen, and C. Guo, Nonlinear optics of graphene and other 2d materials in layered structures, Journal of Physics: Photonics 1, 015002 (2018).
  * Cox and de Abajo [2014] J. D. Cox and F. J. G. de Abajo, Electrically tunable nonlinear plasmonics in graphene nanoislands, Nature Communications 5, 10.1038/ncomms6725 (2014).
  * Hendry  _et al._ [2010] E. Hendry, P. J. Hale, J. Moger, A. K. Savchenko, and S. A. Mikhailov, Coherent nonlinear optical response of graphene, Physical Review Letters 105, 10.1103/physrevlett.105.097401 (2010).
  * Mikhailov [2007] S. Mikhailov, Non-linear electromagnetic response of graphene, EPL 79, 10.1209/0295-5075/79/27002 (2007), cited By 363.
  * Mikhailov and Ziegler [2008] S. Mikhailov and K. Ziegler, Nonlinear electromagnetic response of graphene: Frequency multiplication and the self-consistent-field effects, Journal of Physics Condensed Matter 20, 10.1088/0953-8984/20/38/384204 (2008), cited By 361.
  * Ishikawa [2010] K. L. Ishikawa, Nonlinear optical response of graphene in time domain, Physical Review B 82, 10.1103/physrevb.82.201402 (2010).
  * Zhang  _et al._ [2012] H. Zhang, S. Virally, Q. Bao, L. K. Ping, S. Massar, N. Godbout, and P. Kockaert, Z-scan measurement of the nonlinear refractive index of graphene, Optics Letters 37, 1856 (2012).
  * Hong  _et al._ [2013] S.-Y. Hong, J. I. Dadap, N. Petrone, P.-C. Yeh, J. Hone, and R. M. Osgood, Optical third-harmonic generation in graphene, Physical Review X 3, 10.1103/physrevx.3.021014 (2013).
  * Yang  _et al._ [2018] H. Yang, H. Guan, N. Biekert, G. Arefe, D. C. Chang, Y. Sun, P.-C. Yeh, X. Liu, S.-Y. Hong, I. D. Marion, M. Kralj, J. C. Hone, R. M. Osgood, and J. I. Dadap, Layer dependence of third-harmonic generation in thick multilayer graphene, Physical Review Materials 2, 10.1103/physrevmaterials.2.071002 (2018).
  * Lin  _et al._ [2018] Y. Lin, Y. Li, J. T. Sadowski, W. Jin, J. I. Dadap, M. S. Hybertsen, and R. M. Osgood, Excitation and characterization of image potential state electrons on quasi-free-standing graphene, Physical Review B 97, 10.1103/physrevb.97.165413 (2018).
  * de Abajo [2013] F. J. G. de Abajo, Graphene nanophotonics, Science 339, 917 (2013).
  * Driscoll  _et al._ [2012] J. B. Driscoll, N. Ophir, R. R. Grote, J. I. Dadap, N. C. Panoiu, K. Bergman, and R. M. Osgood, Width-modulation of si photonic wires for quasi-phase-matching of four-wave-mixing: experimental and theoretical demonstration, Optics Express 20, 9227 (2012).
  * Driscoll  _et al._ [2015] J. B. Driscoll, R. M. Osgood, R. R. Grote, J. I. Dadap, and N. C. Panoiu, Squeezing light in wires: Fundamental optical properties of si nanowire waveguides, Journal of Lightwave Technology 33, 3116 (2015).
  * Pitilakis and Kriezis [2013] A. Pitilakis and E. E. Kriezis, Highly nonlinear hybrid silicon-plasmonic waveguides: analysis and optimization, J. Opt. Soc. Am. B 30, 1954 (2013).
  * Pitilakis  _et al._ [2016] A. Pitilakis, D. Chatzidimitriou, and E. E. Kriezis, Theoretical and numerical modeling of linear and nonlinear propagation in graphene waveguides, Optical and Quantum Electronics 48, 10.1007/s11082-016-0510-5 (2016).
  * Yao  _et al._ [2018] B. Yao, S.-W. Huang, Y. Liu, A. K. Vinod, C. Choi, M. Hoff, Y. Li, M. Yu, Z. Feng, D.-L. Kwong, Y. Huang, Y. Rao, X. Duan, and C. W. Wong, Gate-tunable frequency combs in graphene–nitride microresonators, Nature 558, 410 (2018).
  * Mock [2017] A. Mock, Modeling passive mode-locking via saturable absorption in graphene using the finite-difference time-domain method, IEEE J. Quantum Electron. 53, 1 (2017).
  * Breusing  _et al._ [2011] M. Breusing, S. Kuehn, T. Winzer, E. Malić, F. Milde, N. Severin, J. P. Rabe, C. Ropers, A. Knorr, and T. Elsaesser, Ultrafast nonequilibrium carrier dynamics in a single graphene layer, Physical Review B 83, 10.1103/physrevb.83.153410 (2011).
  * Winnerl  _et al._ [2011] S. Winnerl, M. Orlita, P. Plochocka, P. Kossacki, M. Potemski, T. Winzer, E. Malic, A. Knorr, M. Sprinkle, C. Berger, W. A. de Heer, H. Schneider, and M. Helm, Carrier relaxation in epitaxial graphene photoexcited near the dirac point, Physical Review Letters 107, 10.1103/physrevlett.107.237401 (2011).
  * Malic  _et al._ [2011] E. Malic, T. Winzer, E. Bobkin, and A. Knorr, Microscopic theory of absorption and ultrafast many-particle kinetics in graphene, Physical Review B 84, 10.1103/physrevb.84.205406 (2011).
  * Cox  _et al._ [2017] J. D. Cox, A. Marini, and F. J. G. de Abajo, Plasmon-assisted high-harmonic generation in graphene, Nature Communications 8, 10.1038/ncomms14380 (2017).
  * Baudisch  _et al._ [2018] M. Baudisch, A. Marini, J. D. Cox, T. Zhu, F. Silva, S. Teichmann, M. Massicotte, F. Koppens, L. S. Levitov, F. J. G. de Abajo, and J. Biegert, Ultrafast nonlinear optical response of dirac fermions in graphene, Nature Communications 9, 10.1038/s41467-018-03413-7 (2018).
  * Hafez  _et al._ [2018] H. A. Hafez, S. Kovalev, J.-C. Deinert, Z. Mics, B. Green, N. Awari, M. Chen, S. Germanskiy, U. Lehnert, J. Teichert, Z. Wang, K.-J. Tielrooij, Z. Liu, Z. Chen, A. Narita, K. Müllen, M. Bonn, M. Gensch, and D. Turchinovich, Extremely efficient terahertz high-harmonic generation in graphene by hot dirac fermions, Nature 561, 507 (2018).
  * Hafez  _et al._ [2019] H. A. Hafez, S. Kovalev, K.-J. Tielrooij, M. Bonn, M. Gensch, and D. Turchinovich, Terahertz nonlinear optics of graphene: From saturable absorption to high-harmonics generation, Advanced Optical Materials 8, 1900771 (2019).
  * Cheng  _et al._ [2014] J. L. Cheng, N. Vermeulen, and J. E. Sipe, Third order optical nonlinearity of graphene, New Journal of Physics 16, 053014 (2014).
  * Cheng  _et al._ [2015] J. L. Cheng, N. Vermeulen, and J. E. Sipe, Third-order nonlinearity of graphene: Effects of phenomenological relaxation and finite temperature, Phys. Rev. B 91, 10.1103/physrevb.91.235320 (2015).
  * Mikhailov [2016] S. A. Mikhailov, Quantum theory of the third-order nonlinear electrodynamic effects of graphene, Phys. Rev. B 93, 10.1103/physrevb.93.085403 (2016).
  * Dremetsika  _et al._ [2016] E. Dremetsika, B. Dlubak, S.-P. Gorza, C. Ciret, M.-B. Martin, S. Hofmann, P. Seneor, D. Dolfi, S. Massar, P. Emplit, and P. Kockaert, Measuring the nonlinear refractive index of graphene using the optical kerr effect method, Optics Letters 41, 3281 (2016).
  * Jiang  _et al._ [2018] T. Jiang, D. Huang, J. Cheng, X. Fan, Z. Zhang, Y. Shan, Y. Yi, Y. Dai, L. Shi, K. Liu, C. Zeng, J. Zi, J. E. Sipe, Y.-R. Shen, W.-T. Liu, and S. Wu, Gate-tunable third-order nonlinear optical response of massless dirac fermions in graphene, Nature Photonics 12, 430 (2018).
  * Chatzidimitriou  _et al._ [2015] D. Chatzidimitriou, A. Pitilakis, and E. E. Kriezis, Rigorous calculation of nonlinear parameters in graphene-comprising waveguides, J. Appl. Phys. 118, 023105 (2015).
  * Ooi  _et al._ [2014] K. J. A. Ooi, L. K. Ang, and D. T. H. Tan, Waveguide engineering of graphene's nonlinearity, Applied Physics Letters 105, 111110 (2014).
  * Vermeulen  _et al._ [2016] N. Vermeulen, D. Castelló-Lurbe, J. Cheng, I. Pasternak, A. Krajewska, T. Ciuk, W. Strupinski, H. Thienpont, and J. V. Erps, Negative kerr nonlinearity of graphene as seen via chirped-pulse-pumped self-phase modulation, Physical Review Applied 6, 10.1103/physrevapplied.6.044006 (2016).
  * Alexander  _et al._ [2017] K. Alexander, N. A. Savostianova, S. A. Mikhailov, B. Kuyken, and D. V. Thourhout, Electrically tunable optical nonlinearities in graphene-covered SiN waveguides characterized by four-wave mixing, ACS Photonics 4, 3039 (2017).
  * Alexander  _et al._ [2018] K. Alexander, N. A. Savostianova, S. A. Mikhailov, D. V. Thourhout, and B. Kuyken, Gate-tunable nonlinear refraction and absorption in graphene-covered silicon nitride waveguides, ACS Photonics 5, 4944 (2018).
  * Vermeulen  _et al._ [2018] N. Vermeulen, D. Castelló-Lurbe, M. Khoder, I. Pasternak, A. Krajewska, T. Ciuk, W. Strupinski, J. Cheng, H. Thienpont, and J. V. Erps, Graphene’s nonlinear-optical physics revealed through exponentially growing self-phase modulation, Nat. Commun. 9, 10.1038/s41467-018-05081-z (2018).
  * Demongodin  _et al._ [2019] P. Demongodin, H. E. Dirani, J. Lhuillier, R. Crochemore, M. Kemiche, T. Wood, S. Callard, P. Rojo-Romeo, C. Sciancalepore, C. Grillet, and C. Monat, Ultrafast saturable absorption dynamics in hybrid graphene/Si3N4 waveguides, APL Photonics 4, 076102 (2019).
  * Wang  _et al._ [2015] J. Wang, Z. Cheng, Z. Chen, J.-B. Xu, H. K. Tsang, and C. Shu, Graphene photodetector integrated on silicon nitride waveguide, Journal of Applied Physics 117, 144504 (2015).
  * Shiramin and Thourhout [2017] L. A. Shiramin and D. V. Thourhout, Graphene modulators and switches integrated on silicon and silicon nitride waveguide, IEEE Journal of Selected Topics in Quantum Electronics 23, 94 (2017).
  * Chatzidimitriou and Kriezis [2018] D. Chatzidimitriou and E. E. Kriezis, Optical switching through graphene-induced exceptional points, J. Opt. Soc. Am. B 35, 1525 (2018).
  * Doukas  _et al._ [2018] S. Doukas, A. Chatzilari, A. Dagkli, A. Papagiannopoulos, and E. Lidorikis, Deep and fast free-space electro-absorption modulation in a mobility-independent graphene-loaded bragg resonator, Applied Physics Letters 113, 011102 (2018).
  * Ono  _et al._ [2019] M. Ono, M. Hata, M. Tsunekawa, K. Nozaki, H. Sumikura, H. Chiba, and M. Notomi, Ultrafast and energy-efficient all-optical switching with graphene-loaded deep-subwavelength plasmonic waveguides, Nature Photonics 14, 37 (2019).
  * Chatzidimitriou and Kriezis [2020] D. Chatzidimitriou and E. E. Kriezis, Light propagation in nanophotonic waveguides considering graphene’s saturable absorption, Phys. Rev. A 102, 053512 (2020).
  * Zhang  _et al._ [2021] Y. Zhang, J. Wu, Y. Qu, L. Jia, B. Jia, and D. J. Moss, Optimizing the kerr nonlinear optical performance of silicon waveguides integrated with 2d graphene oxide films, Journal of Lightwave Technology 39, 4671 (2021).
  * Pitilakis  _et al._ [2021] A. Pitilakis, D. Chatzidimitriou, T. Yioultsis, and E. E. Kriezis, Asymmetric si-slot coupler with nonreciprocal response based on graphene saturable absorption, IEEE Journal of Quantum Electronics , 1 (2021).
  * Chatzidimitriou  _et al._ [2021] D. Chatzidimitriou, A. Pitilakis, T. Yioultsis, and E. E. Kriezis, Breaking reciprocity in a non-hermitian photonic coupler with saturable absorption, Physical Review A 103, 10.1103/physreva.103.053503 (2021).
  * Sahoo  _et al._ [2021] A. Sahoo, A. Marini, and S. Roy, Free-carrier-induced nonlinear dynamics in hybrid graphene-based photonic waveguides, Physical Review A 104, 10.1103/physreva.104.063501 (2021).
  * Semnani  _et al._ [2016] B. Semnani, A. H. Majedi, and S. Safavi-Naeini, Nonlinear quantum optical properties of graphene, Journal of Optics 18, 035402 (2016).
  * Semnani  _et al._ [2017] B. Semnani, S. M. Raeis-Zadeh, A. Rohani, A. H. Majedi, and S. Safavi-Naeini, Graphene-integrated plasmonic structure for optical third harmonic generation, IEEE Journal of Selected Topics in Quantum Electronics 23, 156 (2017).
  * Marini  _et al._ [2017] A. Marini, J. D. Cox, and F. J. G. de Abajo, Theory of graphene saturable absorption, Phys. Rev. B 95, 10.1103/physrevb.95.125408 (2017).
  * Semnani  _et al._ [2019] B. Semnani, R. Jago, S. Safavi-Naeini, H. Majedi, E. Malic, and P. Tassin, Anomalous optical saturation of low-energy dirac states in graphene and its implication for nonlinear optics, 2D Materials 6, 031003 (2019).
  * Mikhailov [2019] S. A. Mikhailov, Theory of the strongly nonlinear electrodynamic response of graphene: A hot electron model, Phys. Rev. B 100, 10.1103/physrevb.100.115416 (2019).
  * Soavi  _et al._ [2019] G. Soavi, G. Wang, H. Rostami, A. Tomadin, O. Balci, I. Paradisanos, E. A. A. Pogna, G. Cerullo, E. Lidorikis, M. Polini, and A. C. Ferrari, Hot electrons modulation of third-harmonic generation in graphene, ACS Photonics 6, 2841 (2019).
  * Cheng  _et al._ [2020] J. L. Cheng, J. E. Sipe, and S. W. Wu, Third-order optical nonlinearity of three-dimensional massless dirac fermions, ACS Photonics 7, 2515 (2020).
  * Mikhailov [2021] S. A. Mikhailov, Nonperturbative quasiclassical theory of graphene photoconductivity, Physical Review B 103, 10.1103/physrevb.103.245406 (2021).
  * You  _et al._ [2018] J. You, S. Bongu, Q. Bao, and N. Panoiu, Nonlinear optical properties and applications of 2d materials: theoretical and experimental aspects, Nanophotonics 8, 63 (2018).
  * Butcher and Cotter [1990] P. Butcher and D. Cotter, eds.,  _The Elements of Nonlinear Optics_ (Cambridge University Press, Great Britain, 1990).
  * Boyd [2020] R. Boyd,  _Nonlinear Optics_ (Elsevier Science, 2020).
  * Falkovsky [2008] L. A. Falkovsky, Optical properties of doped graphene layers, Journal of Experimental and Theoretical Physics 106, 575 (2008).
  * del Coso and Solis [2004] R. del Coso and J. Solis, Relation between nonlinear refractive index and third-order susceptibility in absorbing media, Journal of the Optical Society of America B 21, 640 (2004).
  * Savostianova and Mikhailov [2018] N. A. Savostianova and S. A. Mikhailov, Optical kerr effect in graphene: Theoretical analysis of the optical heterodyne detection technique, Physical Review B 97, 10.1103/physrevb.97.165424 (2018).
  * Glover and Tinkham [1957] R. E. Glover and M. Tinkham, Conductivity of superconducting films for photon energies between 0.3 and40ktc, Physical Review 108, 243 (1957).
  * Agrawal [2012] G. Agrawal,  _Nonlinear Fiber Optics_, Quantum electronics–principles and applications (Elsevier Science, 2012).
  * Lin  _et al._ [2007] Q. Lin, O. J. Painter, and G. P. Agrawal, Nonlinear optical phenomena in silicon waveguides: modeling and applications, Opt. Express 15, 16604 (2007).
  * Afshar Vahid and Monro [2009] S. Afshar Vahid and T. M. Monro, A full vectorial model for pulse propagation in emerging waveguides with subwavelength structures part i: Kerr nonlinearity, Opt. Express 17, 2298 (2009).
  * Daniel and Agrawal [2010] B. A. Daniel and G. P. Agrawal, Vectorial nonlinear propagation in silicon nanowire waveguides: polarization effects, Journal of the Optical Society of America B 27, 956 (2010).
  * Bao  _et al._ [2009] Q. Bao, H. Zhang, Y. Wang, Z. Ni, Y. Yan, Z. X. Shen, K. P. Loh, and D. Y. Tang, Atomic-layer graphene as a saturable absorber for ultrafast pulsed lasers, Adv. Funct. Mater. 19, 3077 (2009).
  * Bao  _et al._ [2010] Q. Bao, H. Zhang, Z. Ni, Y. Wang, L. Polavarapu, Z. Shen, Q.-H. Xu, D. Tang, and K. P. Loh, Monolayer graphene as a saturable absorber in a mode-locked laser, Nano Research 4, 297 (2010).
  * Zhang  _et al._ [2010] H. Zhang, D. Tang, R. J. Knize, L. Zhao, Q. Bao, and K. P. Loh, Graphene mode locked, wavelength-tunable, dissipative soliton fiber laser, Applied Physics Letters 96, 111112 (2010).
  * Castelló-Lurbe  _et al._ [2020] D. Castelló-Lurbe, H. Thienpont, and N. Vermeulen, Predicting graphene's nonlinear-optical refractive response for propagating pulses, Laser Photonics Rev. 14, 1900402 (2020).
  * Neto  _et al._ [2009] A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Reviews of Modern Physics 81, 109 (2009).
  * Falkovsky and Varlamov [2007] L. A. Falkovsky and A. A. Varlamov, Space-time dispersion of graphene conductivity, The European Physical Journal B 56, 281 (2007).
  * Sarma  _et al._ [2011] S. D. Sarma, S. Adam, E. H. Hwang, and E. Rossi, Electronic transport in two-dimensional graphene, Reviews of Modern Physics 83, 407 (2011).
  * Trushin and Schliemann [2007] M. Trushin and J. Schliemann, Minimum electrical and thermal conductivity of graphene: A quasiclassical approach, Physical Review Letters 99, 10.1103/physrevlett.99.216602 (2007).
  * Hwang and Sarma [2009] E. H. Hwang and S. D. Sarma, Screening-induced temperature-dependent transport in two-dimensional graphene, Physical Review B 79, 10.1103/physrevb.79.165404 (2009).
  * Gu  _et al._ [2012] T. Gu, N. Petrone, J. F. McMillan, A. van der Zande, M. Yu, G. Q. Lo, D. L. Kwong, J. Hone, and C. W. Wong, Regenerative oscillation and four-wave mixing in graphene optoelectronics, Nature Photonics 6, 554 (2012).