# Two-photon absorption in semiconductors: a multi-band length-gauge analysis

W.-R. Hannes  M.F. Ciappina  Physics Program, Guangdong Technion–Israel
Institute of Technology, Shantou, Guangdong 515063, China  Technion—Israel
Institute of Technology, Haifa 32000, Israel

###### Abstract

The simplest approach to deal with light excitations in direct-gap
semiconductors is to model them as a two-band system: one conduction and one
valence band. For such models, particularly simple analytical expressions are
known to exist for the optical response such as multi-photon absorption
coefficients. Here we show that generic multi-band models do not require much
more complicated expressions. Our length-gauge analysis is based on the
semiconductors Bloch equations in the absence of all scattering processes. In
the evaluation we focus on two-photon excitation by a pump-probe scheme with
possibly non-degenerate and arbitrarily polarized configurations. The theory
is validated by application to graphene and its bilayer, described by a tight-
binding model, as well as bulk Zincblende semiconductors described by
${\bm{k\cdot}}{\bm{p}}$ theory.

##  I Introduction

A broad variety of theoretical methods exists to investigate certain nonlinear
susceptibilities of crystals. In the past, these methods have been often
applied to simplified band-structure models, in many cases just two-band
models [1, 2]. For certain materials like graphene [3] or hBN [4], such models
might capture the most relevant aspects of the optical response. However, for
more complicated materials including Zincblende semiconductors, they are
typically not sufficient for two reasons: (i) more than two states are
optically coupled for the considered excitation frequency (e.g., due to
valence band degeneracies) and (ii) the influence of further bands (like
higher conduction bands) can be considerable, for instance by leading to
anisotropies in the optical response [5, 6].

The development of approaches to the nonlinear response of bulk semiconductors
in more rigorous band-structure models has a long history, even though the
respective results are not in agreement with each other in all details. For
the case of Zincblende semiconductors, the Kane band-structure model has been
used to evaluate two-photon absorption (2PA) coefficients based on the
transition rate approach. [7, 8, 6] Murayami and Nakayama [9, 10, 11]
evaluated 2PA spectra using an ab initio approach based on the microscopic
response theory originally developed for atoms. [12] A more profound approach,
which takes into account the subtleties of the dipole operator in an extended
crystal, is the length gauge analysis by Aversa and Sipe. [13] This lead to
general expressions for the non-linear susceptibility $\chi^{(3)}$-tensor in a
multiband model. However, to our knowledge the expressions in Ref. [13] have
not yet been explicitly applied to any rigorous band-structure model. We can
suspect two reasons: (i) due to the use of the length gauge they contain
$\mathbf{k}$-derivatives of matrix elements, which are not directly available
from numerical band-structure programs due to the random phase problem; (ii)
The appearance of resonance denominators in these expressions renders a direct
numerical evaluation of the sum over $\mathbf{k}$-values difficult, in
particular in the limit of vanishing relaxation. Recently, Passos et al [14]
presented a similarly general analysis of the nonlinear optical conductivity
in the independent particle approximation, which elucidated its structure as
well as some general characteristics. However, the 2PA coefficients that could
be derived from the presented (to be symmetrized) third-order conductivities
still contain many terms, and the numerical evaluation was limited to rather
simple models.

In this sense we can state that according to common belief, multi-photon
absorption coefficients, like other nonlinear susceptibilities, are bulky
expressions (after expanding commutators and performing differentiations),
which, moreover, still need to be integrated over the relevant region of
reciprocal space. Here we show that rather simple expressions are obtained
with appropriate limitations. Our perturbative analysis is carried out in the
length gauge as the previous works, [13, 14, 2] but is more specialized on
multi (mostly two) -photon absorption coefficients in a pump-probe setting,
where the probe delivers one photon and the pump the remaining photons to the
absorption process. Taking into account arbitrary (under the limitations) pump
and probe frequencies and linear polarization directions (with respect to each
other and the crystallographic axes), such experiments can reveal all relevant
information about the nonlinear response (real part of optical conductivity)
of crystals. Since we consider the limit of absent scattering, the presented
expressions for $\ell$PA coefficients are integrals over the
$\mathbf{k}$-points at which two states are optically coupled, i.e., where the
corresponding interband transition is resonant with the total energy of $\ell$
absorbed photons. This means, that for a 3D crystal, only a 2D contour
integral is to be evaluated over the resonance isosurfaces, which greatly
reduces numerical cost. The integrand contains only three factors: (i) the
square modulus of a function ${\cal G}^{(n)}_{vc\mathbf{k}}$, which depends on
the excitation parameters, on transition energies, transition dipole moments,
as well as their derivatives with respect to $\mathbf{k}$ (along the pulse
directions); (ii) $(f_{v\mathbf{k}}-f_{c\mathbf{k}})$, the difference in band
populations, which thus do not enter in a more complicated way, e.g., via
their $\mathbf{k}$-derivatives as one could expect [14]; (iii) a factor
$\left\|\nabla_{\mathbf{k}}({\omega}_{v\mathbf{k}}-{\omega}_{c\mathbf{k}})\right\|^{-1}$,
the inverse norm of the transition group velocity. In the linear order, ${\cal
G}^{(1)}_{vc\mathbf{k}}$ is simply the transition dipole component along the
probe polarization direction, while for $n>1$ ($n$ being odd), other
transitions enter as well.

A different approach within the same length gauge framework is to model the
response by dynamical simulations, that is by numerically solving the time
evolution for finite-length pulses. While this has been carried out using a
multiband parallel transport gauge in Ref. [15], the computational effort is
immense and increases quadratically with the number of modeled bands. It also
increases higher than linearly with the duration of the exciting pulses, since
longer pulses require higher $\mathbf{k}$-resolution. For this reason,
dynamical simulations are an alternative but not particularly convenient tool
for the description of the nonlinear-optical response of crystals. Since
dynamical aspects do not play a crucial role in absorptive processes, a
steady-state description is sufficient to even analyze the response to ultra-
short pulses in the fs-range, as long as other processes (dephasing and
relaxation) happen on longer time scales, which is typically the case in clean
samples. For spectrally broad pulses, it might be necessary to take all
frequency components of each pulse into account. In this work we include a
comparison between the steady-state description, approximating each pulse by a
single carrier frequency, with dynamical simulations, and find that these are
in excellent agreement with each other.

Our article is organized as follows. In Sec. II we derive the perturbative
equations for the generic multiband model and relate the expectation value of
the perturbative current with the optical conductivities / absorption
coefficients. We also show the transformation of the $\mathbf{k}$-integral to
the contour integral. Sec. III contains the explicit evaluation/simplificaiton
of 1PA and 2PA coefficients in the general multiband model. The obtained
results are applied in Sec. IV to specific materials, namely mono- and bilayer
graphene and Zincblende semiconductors. In the nearest-neighbour tight-binding
model for monolayer graphene we obtain analytical results for both 1PA and
2PA. For GaAs and ZnSe, we study the dependencies of 2PA on polarization and
non-degeneracy parameter, and compare with dynamical simulations. We conclude
in Sec. V.

##  II General formalism

###  II.1 Perturbation theory

For brevity, we derive the perturbative analysis directly from the
semiconductor-Bloch equations (SBE) rather than the von Neumann equation.
Scattering terms are neglected from the beginning. Using the length gauge, the
dynamical equation for the single-particle density matrix $\rho_{\mathbf{k}}$
is

$\displaystyle\frac{\partial}{\partial t}\rho_{\mathbf{k}}$
$\displaystyle=\frac{i}{\hbar}\big{[}\epsilon_{\mathbf{k}}-e_{0}\bm{E}(t)\cdot\bm{\xi}_{\mathbf{k}},\rho_{\mathbf{k}}\big{]}+\frac{1}{\hbar}{e_{0}\bm{E}}(t)\cdot{\nabla}_{\mathbf{k}}\rho_{\mathbf{k}},$
(1)

where $\epsilon_{\mathbf{k}}$ is the diagonal matrix of band energies
$\hbar\omega_{\lambda\mathbf{k}}$, and $\bm{\xi}_{\mathbf{k}}$ is the
nonabelian Berry connection. Since $\epsilon_{\mathbf{k}}$ is diagonal, we
have restricted ourselves to proper gauges which is in contrast to the
multiband gauge applied in Ref. [15]. The reason is that for non-diagonal
(gauge-transformed) $\epsilon_{\mathbf{k}}$ the matrix elements of
$\rho_{\mathbf{k}}^{(n)}$ become coupled within each order, making an
analytical solution rather involved.

The perturbative equations are obtained from the series
$\rho_{\mathbf{k}}=\sum_{i=0}^{\infty}\rho_{\mathbf{k}}^{(i)}(t)$, where
$\rho_{\mathbf{k}}^{(0)}=\rho_{\mathbf{k}}^{0}$ is the initial density matrix
with
$\rho_{\lambda\mu\mathbf{k}}^{0}=\delta_{\lambda\mu}f_{\lambda\mathbf{k}}$.
For any $n>0$ we have

$\displaystyle\frac{\partial}{\partial t}\rho_{\mathbf{k}}^{(n)}$
$\displaystyle=i\omega_{\mathbf{k}}\circ\rho_{\mathbf{k}}^{(n)}-i\frac{e_{0}}{\hbar}{\bm{E}}(t)\cdot\left[\bm{\xi}_{\mathbf{k}},\rho_{\mathbf{k}}^{(n-1)}\right]+\frac{e_{0}}{\hbar}{\bm{E}}(t)\cdot{\nabla}_{\mathbf{k}}\rho_{\mathbf{k}}^{(n-1)},$
(2)

where
$\omega_{\lambda\mu\mathbf{k}}\equiv\omega_{\lambda\mathbf{k}}-\omega_{\mu\mathbf{k}}$,
and $\circ$ is the elementwise product. The electric field is Fourier
transformed as

$\displaystyle\bm{E}(t)=\int\frac{\mathrm{d}\omega}{2\pi}{\bm{E}}_{\omega}e^{-i\omega
t}$ (3)

where the reality condition implies ${\bm{E}}_{\omega}={\bm{E}}_{-\omega}$.
The temporal dependence of the density matrix can also be written in the form
of a Fourier transform

$\displaystyle\rho_{\mathbf{k}}^{(n)}(t)$
$\displaystyle=\left(\frac{e_{0}}{2\pi\hbar}\right)^{n}\\!\\!\int\\!\\!\mathrm{d}\omega_{1}\dots\mathrm{d}\omega_{n}e^{-i\omega_{1}t}\dots
e^{-i\omega_{n}t}E_{\omega_{1}}^{\alpha_{1}}\dots
E_{\omega_{n}}^{\alpha_{n}}{\cal
P}_{\mathbf{k}}^{(n);\alpha_{1}\dots\alpha_{n}},$ (4)

with $\rho_{\mathbf{k}}^{(0)}={\cal P}_{\mathbf{k}}^{(0)}$. Here superscripts
$\beta,\alpha_{i}$ indicate Cartesian components, which are to be summed over
if repeated. The tensor ${\cal
P}_{\mathbf{k}}^{(n);\alpha_{1}\dots\alpha_{n}}$ depends on frequency
arguments $(\omega_{1},\dots,\omega_{n})$, which have been dropped, since they
can be inferred from the Cartesian components. Substitution into Eq. (2)
immediately gives the recursive solutions

$\displaystyle{\cal P}_{\mathbf{k}}^{(n);\alpha_{1}\dots\alpha_{n}}$
$\displaystyle=\frac{1}{\omega_{1}+\dots+\omega_{n}+\omega_{\mathbf{k}}}\circ\left(i\partial^{\alpha_{n}}{\cal
P}_{\mathbf{k}}^{(n-1);\alpha_{1}\dots\alpha_{n-1}}\right.$
$\displaystyle\qquad\qquad\qquad\qquad+\left.\left[\xi_{\mathbf{k}}^{\alpha_{n}},{\cal
P}_{\mathbf{k}}^{(n-1);\alpha_{1}\dots\alpha_{n-1}}\right]\right),$ (5)

where $\partial^{\alpha_{i}}\equiv{\partial}/{\partial k^{\alpha_{i}}}$

Instead of splitting the inhomogeneity in Eq. (1) (and in the following
analysis) into the dipole matrix and the gradient terms, one could also divide
it into an interband dipole matrix (with diagonal elements set to zero) and a
generalized derivative (defined in Eq. (46)) term. This would correspond to a
proper inter/intra-band term splitting; however, it does not seem to be
beneficial for simplifying the analysis presented here.

###  II.2 Macroscopic perturbative response

The macroscopic current density is given by

$\displaystyle\bm{J}=-\frac{e_{0}}{L^{3}}\sum_{\mathbf{k}}\mathrm{Tr}(\bm{v}_{\mathbf{k}}\rho_{\mathbf{k}}),$
(6)

where $L^{3}$ is the sample volume (including both confined and unconfined
dimensions). The velocity matrix is

$\displaystyle\hbar\bm{v}_{\mathbf{k}}$
$\displaystyle={\nabla}_{\mathbf{k}}\epsilon_{\mathbf{k}}+i\big{[}\epsilon_{\mathbf{k}},\bm{\xi}_{\mathbf{k}}\big{]},$
(7)

or

$\displaystyle\bm{v}_{\lambda\mu\mathbf{k}}=\left\\{\begin{array}[]{ll}\nabla_{\mathbf{k}}\omega_{\lambda\mathbf{k}},&\lambda=\mu\\\\[8.0pt]
i\omega_{\lambda\mu\mathbf{k}}\bm{\xi}_{\lambda\mu\mathbf{k}},&\lambda\neq\mu\\\
\end{array}\right.$ (10)

The perturbative response is obtained simply by substituting the perturbative
density matrix. In particular,

$\displaystyle J^{(n)\beta}(t)$
$\displaystyle=-\frac{e_{0}}{L^{3}}\sum_{\mathbf{k}}\mathrm{Tr}\left(v_{\mathbf{k}}^{\beta}\rho_{\mathbf{k}}^{(n)}(t)\right),$
(11)

The relation between the applied fields and the current response is given by a
conductivity tensor,

$\displaystyle J^{(n)\beta}(t)$
$\displaystyle=\int\frac{\mathrm{d}\omega_{1}\dots\mathrm{d}\omega_{n}}{(2\pi)^{n}}e^{-i\omega_{t}t}E_{\omega_{1}}^{\alpha_{1}}\dots
E_{\omega_{n}}^{\alpha_{n}}{\sigma}^{(n)}_{\beta\alpha_{1}\dots\alpha_{n}},$
(12)

where $\omega_{t}=\omega_{1}+\dots+\omega_{n}$. Comparison with (4) and (11)
shows that the conductivity tensor is

${\sigma}^{(n)}_{\beta\alpha_{1}\dots\alpha_{n}}=-\frac{e_{0}}{L^{3}}\left(\frac{e_{0}}{\hbar}\right)^{n}\sum_{\mathbf{k}}\mathrm{Tr}\left(v_{\mathbf{k}}^{\beta}{\cal
P}_{\mathbf{k}}^{(n);\alpha_{1}\dots\alpha_{n}}\right).$ (13)

The arguments $(\omega_{t};\omega_{1},\dots,\omega_{n})$ have been again
dropped for shortness. The symmetrization with respect to the $n$ components
and corresponding frequency arguments is carried out specifically for the
setting described in the next subseciton; a complete symmetrization at this
point is not useful.

###  II.3 Probe absorption induced by the pump

In this work we investigate multi-photon absorption in a pump-probe setting,
where the pump and probe beam are assumed to be nearly co-propagating along
the normal to the sample surface. Following Ref. [6], we denote the probe by
${p}$ and the pump by ${e}$ (for ‘excitation’). The probe, for example, is
characterized by intensity $I_{p}$, frequency $\omega_{p}$, polarization
direction $\hat{\mathbf{p}}$, and propagation direction $\hat{\bm{k}}_{p}$. In
this setting, there are typically multiple absorption processes scaling with
different powers of $I_{p}$ and $I_{e}$; here, when analyzing $\ell$PA with
some $\ell\geq 1$ we exclusively consider the process proportional to
$I_{p}I_{e}^{\ell-1}$, i.e., where the probe delivers one and the pump the
remaining photons of the absorption process. The total intensity change of
$I_{p}$ over the propagation coordinate is made up of all processes; however,
the coefficients for the other processes can be deduced from the one
calculated here. For example, one-beam absorption can be easily deduced as a
special case from the two-beam results presented here, as will become clear
below. Two-beam experiments with linearly-polarized light are sufficient to
completely investigate the third-order conductivity tensor. [6]

Even with the above restriction, there are different processes in the
considered order ($\sim I_{p}I_{e}^{\ell-1}$). In the language of the
resonance-based analysis of Ref. [14], we consider only $\ell=(n+1)/2$-photon
resonant terms when analyzing the $n$-th order response. This means, only
field products leading to a resonance denominator for
$\omega_{p}+(\ell-1)\omega_{e}$ are considered. This is justified mainly by
the fact that, typically, the lower-number photon absorption is off-resonant
and thus vanishing in the assumed limit of absent scattering. For gapless
materials like graphene, the off-resonance condition of lower number photons
cannot be met; however, the other response typically stems from a different
set of $\mathbf{k}$-points and does not directly interact with $\ell$PA. By
this assumption we also exclude from our analysis other nonlinear-optical
effects such as photovoltaic effects, DC field-induced second order response,
and electro-optic effects.

The probe absorption in this $\ell$PA process is determined by the $n$-th
order current response propagating in the probe direction, where $n-1$ out of
$n$ field components belong to the pump. Here and in the following, we
implicitly assume $\ell=(n+1)/2$. For the ‘true’ $\ell$PA as defined above,
the positive frequencies come in the lowest orders ($i=1,\dots,\ell$),
followed by the negative frequencies ($i=\ell+1,\dots,2\ell-1$); the sum is
$\omega_{t}=\omega_{p}$. For example, for 2PA ($n=3$), the combinations are
${\cal E}_{p}{\cal E}_{e}{\cal E}_{e}^{*}$ and ${\cal E}_{e}{\cal E}_{p}{\cal
E}_{e}^{*}$, [6] where the field amplitudes ${\cal E}_{p}$ are defined as in
A. For arbitrary $n$, the effective conductivity determining the absorption
coefficient $\alpha^{(n)}_{p\\{e\\}}$ is thus

$\displaystyle\sigma^{(n)}_{\mathrm{eff}}$ $\displaystyle=\sum_{\cal
P}{\sigma}^{(n)}_{p{\cal P}(p\\{e\\})\\{e\\}}(\omega_{p};{\cal
P}(\omega_{p},\\{\omega_{e}\\}),\\{-\omega_{e}\\}),$ (14)

where ${\cal P}$ stands for permutations, and $\\{\\}$ for $(\ell-1)$-fold
repetition of an element in a list. We have also directly substituted the
pump/probe indices in place of the cartesian indices; for instance,
${\sigma}^{(1)}_{p}\equiv{\sigma}^{(1)}_{\beta}p_{\beta}$. This convention
will be used throughout also for $\mathbf{k}$-dependent quantities such as
$\bm{v}_{\mathbf{k}}$, ${\cal P}_{\mathbf{k}}^{(n);\alpha_{1}\dots\alpha_{n}}$
etc. It will be also convenient to introduce a symmetrized (over relevant
terms) form of ${\cal P}_{\mathbf{k}}^{(n)}$ as $\bar{\cal
P}_{\mathbf{k}}^{(n)}\equiv\sum_{\cal P}{\cal P}_{\mathbf{k}}^{(n);p{\cal
P}(p\\{e\\})\\{e\\}}$. Using the slowly varying envelope approximation, in A
it is shown that the relation between the conductivity tensor and the
absorption coefficient is

$\alpha^{(n)}_{p\\{e\\}}=K_{n}\mathrm{Re}\left[{\sigma}^{(n)}_{\mathrm{eff}}\right],\qquad
K_{n}\equiv\frac{1}{n_{0{p}}n_{0{e}}^{(n-1)/2}}\frac{2}{(2\epsilon_{0}c)^{(n+1)/2}},$
(15)

with background refractive indices, e.g., $n_{0{p}}\equiv n_{0}(\omega_{p})$.
For degenerate and co-polarized pump and probe, the relation between the pump-
probe and the single-beam $\ell$PA coefficient is
$\alpha^{(n)}_{p\\{e\\}}=\ell\alpha^{(n)}$. As mentioned above, the latter can
thus be easily deduced from the more general pump-probe study. Substituting
Eq. (13), we obtain

$\alpha^{(n)}_{p\\{e\\}}=-\tilde{K}_{n}\frac{1}{L^{3}}\sum_{\mathbf{k}}\mathrm{Re}\left[\mathrm{Tr}\left(v_{\mathbf{k}}^{p}\bar{\cal
P}_{\mathbf{k}}^{(n)}\right)\right],$ (16)

with $\tilde{K}_{n}\equiv e_{0}K_{n}\left(e_{0}/\hbar\right)^{n}$. Remember
that $\alpha^{(n)}_{p\\{e\\}}$ depends on the frequencies and the polarization
directions of the probe and pump beams.

As a final note, the absorption coefficient $\alpha^{(n)}_{p\\{e\\}}$ is
related here with an effective conductivity tensor
${\sigma}^{(n)}_{\mathrm{eff}}$, which in turn can be expressed in terms of
usual tensor components ${\sigma}^{(n)}_{\beta\alpha_{1}\dots\alpha_{n}}$. The
common practice would be to first compute the set of relevant tensor
components. However, using $\bar{\cal P}_{\mathbf{k}}^{(n)}$ we have the
option to compute the absorption individually for each polarization setting,
for example, oblique polarization directions with respect to crystal axes.
This strategy emerged mainly because the present study is a follow-up of a
previous work based on dynamical simulations [15]. It also allows us to
validate the polarization dependence with the one following from the tensor
symmetry for the crystal class of the studied material.

###  II.4 $\mathbf{k}$-integration

The $\mathbf{k}$-summation in Eq. (16) is replaced by an integration (of same
dimensionality ${D}$ as the sample),
$L^{-3}\sum_{\mathbf{k}}=\kappa_{D}\int\mathrm{d}^{D}k$, with
$\kappa_{D}\equiv(2\pi)^{-{D}}(L_{\mathrm{c}}^{{D}-3})$. The factor
$L_{\mathrm{c}}^{3-{D}}$ contains the sample extensions in the confined
directions. In this way we obtain a bulk absorption coefficient, with
units/dimensions independent of $D$. Without the factor
$L_{\mathrm{c}}^{3-{D}}$, e.g. $\kappa_{D}\equiv(2\pi)^{-{2}}$ for quasi-2D
materials, a reduced dimensionality absorption coefficient is obtained,
denoted by $\alpha^{(n)}_{D\mathrm{D}}$. In a thin-film approximation, where
the relative intensity change $\mathrm{d}I_{p}/I_{p}$ over the confined sample
dimension, along which the probe propagates, is small,
$\alpha^{(n)}_{2\mathrm{D}}$ corresponds to the absorption integrated over
$z$, which is just $\mathrm{d}I_{p}/I_{p}$.

Using $\mathrm{Re}(c)=\mathrm{Im}(ic)$ (with $c$ a complex number) and
possibly using partial integration in order to take out the resonance
denominator from $\mathbf{k}$-derivatives, we can write Eq. (16) in the form

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$ $\displaystyle=-\tilde{K}_{n}\
\kappa_{D}\sum_{\lambda\mu}\int\mathrm{d}^{D}\mathbf{k}\
\mathrm{Im}\frac{{\cal
F}^{(n)}_{\lambda\mu\mathbf{k}}}{\omega_{\Sigma}-{\omega}_{\mu\lambda\mathbf{k}}},$
(17)

where $\omega_{\Sigma}=\omega_{p}+(\ell-1)\omega_{e}$ is the excitation
frequency sum. In this work, we reserve the band index $\lambda$ for valence
bands (VBs) and the index $\mu$ for conduction bands (CBs). At some places,
the symbols $\lambda$ and $\mu$ are interchanged by symbols $v$ (VB) and $c$
(CB), respectively. When further (VB/CB) band indices are required, we use the
symbols $\eta$ and $\eta^{\prime}$.

The function ${\cal F}^{(n)}_{\lambda\mu\mathbf{k}}$ does not contain further
$\omega_{\Sigma}$-resonance denominators. It does contain other resonance
denominators in general, which may lead to divergences. In the models /
parameters which we consider, this is typically not problematic. For example,
for gapped semiconductors, if both excitation frequencies are below gap, the
other resonance denominators are for inter-valence or inter-conduction band
transitions and never vanish at the same $\mathbf{k}$-point as the 2PA
resonance denominator.

In the case of simple band-structure models, like those which are isotropic
around the band center, it is obvious that the $\mathbf{k}$-integration can be
carried out analytically. This has been demonstrated, e.g., for the low-energy
description of graphene [3, 14] or for a simplistic parabolic model of direct-
gap semiconductors. [14] However, it is in fact generally possible to
transform the ${D}$-dimensional integral into (a sum of) resonance surface
integrals (typically one for each resonant transition) of dimension ${D}-1$.
This greatly reduces the numerical expenses and will be the procedure followed
here.

We first choose a particular direction, denoted by unit vector
$\hat{\mathbf{k}}_{\parallel}$ (not to be confused with a propagation unit
vector), along which we carry out the integration by means of the
Sokhotski–Plemelj theorem. We thus write
$\int\mathrm{d}^{D}\mathbf{k}=\int\mathrm{d}^{{D}-1}\mathbf{k}_{\perp}\int\mathrm{d}k_{\parallel}$
and then we transform the $k_{\parallel}$-integral to frequency space.

To this end we split it at each extremum of ${\omega}_{\mu\lambda\mathbf{k}}$
along the integration line, so that in each segment
${\omega}_{\mu\lambda\mathbf{k}}$ is monotonous. We label the extrema
locations, including the start and end point of the $k_{\parallel}$-integral,
by $\bar{k}_{j}$.

$\displaystyle\int\mathrm{d}k_{\parallel}$
$\displaystyle=\sum_{j}\int_{\bar{k}_{j}}^{\bar{k}_{j+1}}\mathrm{d}k_{\parallel}$
(18)
$\displaystyle=\sum_{j}\int_{\bar{\omega}_{\bar{k}_{j}}}^{\bar{\omega}_{\bar{k}_{j+1}}}\mathrm{d}{\omega}_{\mu\lambda\mathbf{k}}\left(\frac{\partial{\omega}_{\mu\lambda\mathbf{k}}}{\partial
k_{\parallel}}\right)^{\\!\\!-1}$ (19)
$\displaystyle=\sum_{j}\int_{\bar{\omega}_{j,\mathrm{min}}}^{\bar{\omega}_{j,\mathrm{max}}}\mathrm{d}{\omega}_{\mu\lambda\mathbf{k}}\left|\frac{\partial{\omega}_{\mu\lambda\mathbf{k}}}{\partial
k_{\parallel}}\right|^{-1}.$ (20)

This introduces singularities whenever
$\partial{\omega}_{\mu\lambda\mathbf{k}}/\partial k_{\parallel}$ vanishes, but
this problem is relieved below by another transformation. In the last step we
have just reversed the direction of integration in the segments with
$\bar{\omega}_{\bar{k}_{j}}>\bar{\omega}_{\bar{k}_{j+1}}$. Thus,

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$ $\displaystyle=-\tilde{K}_{n}\
\kappa_{D}\sum_{\lambda\mu}\int\mathrm{d}^{{D}-1}\mathbf{k}_{\perp}\sum_{j}\int_{\bar{\omega}_{j,\mathrm{min}}}^{\bar{\omega}_{j,\mathrm{max}}}\mathrm{d}{\omega}_{\mu\lambda\mathbf{k}}\left|\frac{\partial{\omega}_{\mu\lambda\mathbf{k}}}{\partial
k_{\parallel}}\right|^{-1}\\!\\!\\!\mathrm{Im}\frac{{\cal
F}^{(n)}_{\lambda\mu\mathbf{k}}}{\omega_{\Sigma}-{\omega}_{\mu\lambda\mathbf{k}}}$
(21)

Now we carry out the integral by adding a damping term $i\gamma_{vc}$ to the
resonance denominator and taking the limit $\gamma_{vc}\to 0^{+}$, in order to
apply the Sokhotski–Plemelj theorem.

$\lim_{\gamma_{vc}\to 0^{+}}\int
d{\omega}_{\mathbf{k}}\frac{{g}_{\lambda\mu\mathbf{k}}}{\omega_{\Sigma}-{\omega}_{\mathbf{k}}+i\gamma_{vc}}=\left.-i\pi{g}_{\lambda\mu\mathbf{k}}\right|_{\omega_{\Sigma}={\omega}_{\mathbf{k}}}+{\cal
P}\int
d{\omega}_{\mathbf{k}}\frac{{g}_{\lambda\mu\mathbf{k}}}{\omega_{\Sigma}-{\omega}_{\mathbf{k}}},$
(22)

where ${\mathcal{P}}$ denotes the Cauchy principal value. If ${\cal
F}^{(n)}_{\lambda\mu\mathbf{k}}\in\mathds{R}$, then, after taking the
imaginary part, the second term vanishes. We will see that for $n=1$ and $3$
this is indeed the case. The final result becomes

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$ $\displaystyle=\pi\tilde{K}_{n}\
\kappa_{D}\sum_{\lambda\mu}\int\mathrm{d}^{{D}-1}\mathbf{k}_{\perp}\sum_{k_{\parallel,\mathrm{res}}}\left|\frac{\partial{\omega}_{\mu\lambda\mathbf{k}}}{\partial
k_{\parallel}}\right|^{-1}{\cal F}^{(n)}_{\lambda\mu\mathbf{k}}.$ (23)

Here $k_{\parallel,\mathrm{res}}$ are the discrete values of $k_{\parallel}$
(for given $\mathbf{k}_{\perp}$) at which
${\omega}_{\mu\lambda\mathbf{k}}=\omega_{\Sigma}$.

The integral (for ${D}>1$) over the perpendicular $\mathbf{k}$-components,
$\int\mathrm{d}^{{D}-1}\mathbf{k}_{\perp}\sum_{k_{\parallel,\mathrm{res}}}$,
is not in a suitable form for the numerical evaluation, because some of the
lines $\mathbf{k}_{\perp}=const.$ necessarily touch the resonance isosurface
of ${\omega}_{\mu\lambda\mathbf{k}}$ (nearly) tangentially. At these points
the integrand has a singularity due to the factor
$\left|\partial{\omega}_{\mu\lambda\mathbf{k}}/\partial
k_{\parallel}\right|^{-1}$, unless it is compensated by a zero of the function
${\cal F}^{(n)}_{\mathbf{k}}$, but this is not the case except for specific
models. Even though the integral converges, it cannot be numerically evaluated
easily. Fortunately this singularity is removed by transforming the
$\mathbf{k}_{\perp}$-integral to a path/surface integral over the resonance
isosurface. The proof is given in B. With this transformation, the absorption
coefficient (23) becomes

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$
$\displaystyle=\pi\tilde{K}_{n}\kappa_{D}\
\sum_{\lambda\mu}\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}\frac{{\cal
F}^{(n)}_{\lambda\mu\mathbf{k}}}{\left\|\nabla_{\mathbf{k}}{\omega}_{\mu\lambda\mathbf{k}}\right\|}.$
(24)

The denominator does not vanish since we assume to have a well-defined (i.e.,
sharp) resonance contour. The equivalence between Eqs. (23) and (24) is
trivial for ${D}=1$, where $\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}$ becomes
$\sum_{k_{\parallel,\mathrm{res}}}$.

##  III General evaluation

###  III.1 Linear absorption coefficient

To obtain the linear absorption coefficient we need to substitute Eq. (II.1)
for $n=1$,

$\displaystyle{\cal P}_{\mathbf{k}}^{(1);{p}}$
$\displaystyle=\frac{1}{\omega_{p}+\omega_{\mathbf{k}}}\circ\left(i\partial^{{p}}{\cal
P}_{\mathbf{k}}^{(0)}+\left[\xi_{\mathbf{k}}^{{p}},{\cal
P}_{\mathbf{k}}^{(0)}\right]\right),$ (25)

into Eq. (16), also for $n=1$:

$\displaystyle\alpha^{(1)}\equiv\alpha^{(1)}_{p}$
$\displaystyle=-\tilde{K}_{1}\kappa_{D}\int\mathrm{d}^{D}\mathbf{k}\
\mathrm{Im}\left[\mathrm{Tr}\left(iv_{\mathbf{k}}^{p}\bar{\cal
P}_{\mathbf{k}}^{(1)}\right)\right],$ (26)

where $\bar{\cal P}_{\mathbf{k}}^{(1)}={\cal P}_{\mathbf{k}}^{(1);{p}}$. If
$\omega_{p}$ is larger than any inter-valence or inter-conduction band
transitions, which we assume here, then

$\displaystyle\alpha^{(1)}$ $\displaystyle=-\tilde{K}_{1}\
\kappa_{D}\sum_{\lambda\mu}\int\mathrm{d}^{D}\mathbf{k}\
\mathrm{Im}\frac{{\cal
F}^{(1)}_{\lambda\mu\mathbf{k}}}{{\omega}_{\lambda\mu\mathbf{k}}+\omega_{p}},$
(27) $\displaystyle{\cal F}^{(1)}_{\lambda\mu\mathbf{k}}$
$\displaystyle={\omega}_{\mu\lambda\mathbf{k}}\left|\xi_{\lambda\mu\mathbf{k}}^{{p}}\right|^{2}\left(\rho_{\lambda\lambda\mathbf{k}}^{0}-\rho_{\mu\mu\mathbf{k}}^{0}\right)={\omega}_{\mu\lambda\mathbf{k}}\left|\xi_{\lambda\mu\mathbf{k}}^{{p}}\right|^{2}f_{\lambda\mu\mathbf{k}},$
(28)

with $\lambda\in v,\mu\in c$. Since ${\cal F}^{(1)}_{\mathbf{k}}$ is real, we
may use the integration method from Sec. II.4, by simply substituting into Eq.
(24) to obtain

$\displaystyle\alpha^{(1)}$ $\displaystyle=\pi\tilde{K}_{1}\kappa_{D}\
\sum_{\lambda\mu}\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}\frac{{\cal
F}^{(1)}_{\mathbf{k}}}{\left\|\nabla_{\mathbf{k}}{\omega}_{\mu\lambda\mathbf{k}}\right\|}$
(29) $\displaystyle=\pi\omega_{p}\tilde{K}_{1}\kappa_{D}\
\sum_{\lambda\mu}\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}\frac{\left|\xi_{\lambda\mu\mathbf{k}}^{{p}}\right|^{2}f_{\lambda\mu\mathbf{k}}}{\left\|\nabla_{\mathbf{k}}{\omega}_{\mu\lambda\mathbf{k}}\right\|},$
(30)

which agrees with the Fermi’s golden rule.

###  III.2 2PA coefficient

We now turn to 2PA in the order $n=3$. Due to the considered pump-probe
setting we can study, for instance, both degenerate and non-degenerate 2PA, as
well as co- or cross-polarized excitation. Equation (16) for $n=3$ reads

$\displaystyle\beta_{pe}\equiv\alpha^{(3)}_{pe}$
$\displaystyle=-\tilde{K}_{2}\kappa_{D}\int\mathrm{d}^{D}\mathbf{k}\
\mathrm{Im}\left[\mathrm{Tr}\left(iv_{\mathbf{k}}^{p}\bar{\cal
P}_{\mathbf{k}}^{(3)}\right)\right].$ (31)

The symmetrized ${\cal P}$ to be substituted follows from Eq. (II.1) for
$n=3$,

$\displaystyle\bar{\cal P}_{\mathbf{k}}^{(3)}$
$\displaystyle=\frac{1}{\omega_{p}+\omega_{\mathbf{k}}}\circ\left(i\partial^{{e}}\bar{\cal
P}_{\mathbf{k}}^{(2)}+\left[\xi_{\mathbf{k}}^{{e}},\bar{\cal
P}_{\mathbf{k}}^{(2)}\right]\right),$ (32)

with $\bar{\cal P}_{\mathbf{k}}^{(2)}={\cal P}_{\mathbf{k}}^{(2);{pe}}+{\cal
P}_{\mathbf{k}}^{(2);{ep}}$. Substituting Eq. (32) into (31) and carrying out
a partial integration (to take out the 2PA resonance denominator from the
$\mathbf{k}$-derivative) in the first term gives

$\displaystyle\alpha^{(3)}_{pe}$
$\displaystyle=-\tilde{K}_{2}\kappa_{D}\int\mathrm{d}^{D}\mathbf{k}\
\mathrm{Im}\sum_{\lambda\mu}\left[\bar{\cal
P}_{\lambda\mu\mathbf{k}}^{(2)}\partial^{{e}}{\cal
V}_{\lambda\mu\mathbf{k}}^{p}+i\left[\xi_{\mathbf{k}}^{{e}},\bar{\cal
P}_{\mathbf{k}}^{(2)}\right]_{\lambda\mu}{\cal
V}_{\lambda\mu\mathbf{k}}^{p}\right],$ (33) $\displaystyle{\cal
V}_{\lambda\mu\mathbf{k}}^{p}$
$\displaystyle=\frac{v_{\mu\lambda\mathbf{k}}^{p}}{\omega_{p}+\omega_{\lambda\mu\mathbf{k}}}.$
(34)

In the first order we have

$\displaystyle{\cal P}_{\eta\eta^{\prime}\mathbf{k}}^{(1);{p}}$
$\displaystyle=\delta_{\eta\eta^{\prime}}\frac{i\partial^{{p}}f_{\eta\mathbf{k}}}{\omega_{p}}-\frac{\xi_{\eta\eta^{\prime}\mathbf{k}}^{p}f_{\eta\eta^{\prime}\mathbf{k}}}{\omega_{p}+\omega_{\eta\eta^{\prime}\mathbf{k}}}\equiv{\cal
P}_{\eta\eta^{\prime}\mathbf{k}}^{(1);{p}[tra]}+{\cal
P}_{\eta\eta^{\prime}\mathbf{k}}^{(1);{p}[ter]}.$ (35)

In the second order, only the $vc$-elements are relevant,

$\displaystyle{\cal P}_{\lambda\mu\mathbf{k}}^{(2);{pe}}$
$\displaystyle=\left[i\partial^{{e}}{\cal
P}_{\lambda\mu\mathbf{k}}^{(1);{p}[ter]}+\sum_{\eta}\left(\xi_{\lambda\eta\mathbf{k}}^{{e}}{\cal
P}_{\eta\mu\mathbf{k}}^{(1);{p}[tra]}-{\cal
P}_{\lambda\eta\mathbf{k}}^{(1);{p}[tra]}\xi_{\eta\mu\mathbf{k}}^{{e}}\right)+\right.$
$\displaystyle\quad\left.\sum_{\eta\in
v}\xi_{\lambda\eta\mathbf{k}}^{{e}}{\cal
P}_{\eta\mu\mathbf{k}}^{(1);{p}[ter]}-\sum_{\eta\in
c}\xi_{\eta\mu\mathbf{k}}^{{e}}{\cal
P}_{\lambda\eta\mathbf{k}}^{(1);{p}[ter]}\right]\frac{1}{\omega_{p}+\omega_{e}+\omega_{\lambda\mu\mathbf{k}}}.$
(36)

The other term ${\cal P}_{\mathbf{k}}^{(2);{ep}}$ in the symmetric $\bar{\cal
P}_{\mathbf{k}}^{(2)}$ is obtained by interchanging ${p}$ with ${e}$. By
substitution into Eq. (33) and comparison with Eq. (17) we can determine
${\cal F}^{(3)}_{\lambda\mu\mathbf{k}}$ as

$\displaystyle{\cal F}^{(3)}_{\lambda\mu\mathbf{k}}$
$\displaystyle=\left(i\partial^{{e}}{\cal
P}_{\lambda\mu\mathbf{k}}^{(1);{p}[ter]}+\sum_{\eta}\left(\xi_{\lambda\eta\mathbf{k}}^{{e}}{\cal
P}_{\eta\mu\mathbf{k}}^{(1);{p}[tra]}-{\cal
P}_{\lambda\eta\mathbf{k}}^{(1);{p}[tra]}\xi_{\eta\mu\mathbf{k}}^{{e}}\right)+\right.$
$\displaystyle\qquad\left.\sum_{\eta\in
v}\xi_{\lambda\eta\mathbf{k}}^{{e}}{\cal
P}_{\eta\mu\mathbf{k}}^{(1);{p}[ter]}-\sum_{\eta\in
c}\xi_{\eta\mu\mathbf{k}}^{{e}}{\cal
P}_{\lambda\eta\mathbf{k}}^{(1);{p}[ter]}+\\{{p}\leftrightarrow{e}\\}\right)$
$\displaystyle\quad\left(\partial^{{e}}{\cal
V}_{\lambda\mu\mathbf{k}}^{p}+i\sum_{\eta}\left(\xi_{\eta\lambda\mathbf{k}}^{{e}}{\cal
V}_{\eta\mu\mathbf{k}}^{p}-\xi_{\mu\eta\mathbf{k}}^{{e}}{\cal
V}_{\lambda\eta\mathbf{k}}^{p}\right)\right).$ (37)

Note that there are poles in the elements of the matrices ${\cal
V}_{\mathbf{k}}^{p}$ and ${\cal P}_{\mathbf{k}}^{\mathrm{(1)}}$. For the terms
into which their $vc$ elements enter they result in singularities when one
excitation frequency approaches zero. This is to be expected because all
electron scattering processes are neglected. On the other hand, when $vv$ or
$cc$ elements enter, there are singularities when one excitation frequency
matches an inter-valence or inter-conduction band transition. Since ${\cal
F}^{(3)}_{\lambda\mu\mathbf{k}}$ is to be evaluated at 2PA resonant
$\mathbf{k}$-points only, these poles are only problematic if $\omega_{e}$ is
resonant with a $vc$ transition at the same point. In this case the analysis
in Sec. II.4 would not hold due to the presence of multiple poles. However,
for the models studied here, this is not problematic. For the gapped
semiconductors we will consider each frequency to be smaller than the band
gap. For bilayer graphene (Sec. IV.1.2), a gapless multiband model, it can be
problematic, but only at discrete values of excitation frequencies, where the
singularities occur, so it is still safe to use this approach for our
purposes.

We now proceed with the algebraic simplification of ${\cal
F}^{(3)}_{\lambda\mu\mathbf{k}}$. As long as $\omega_{\mu\lambda\mathbf{k}}$
occurs outside of any $\mathbf{k}$-derivative, it may be replaced by
$\omega_{\Sigma}$. The terms in the first factor of ${\cal
F}^{(3)}_{\lambda\mu\mathbf{k}}$ can be evaluated as follows:

$\displaystyle i\partial^{{e}}{\cal
P}_{\lambda\mu\mathbf{k}}^{(1);{p}[ter]}=i\frac{\xi_{\lambda\mu\mathbf{k}}^{{p}}\partial^{{e}}f_{\lambda\mu\mathbf{k}}}{\omega_{e}}+i\frac{f_{\lambda\mu\mathbf{k}}\partial^{{e}}\xi_{\lambda\mu\mathbf{k}}^{{p}}}{\omega_{e}}+i\frac{\xi_{\lambda\mu\mathbf{k}}^{{p}}f_{\lambda\mu\mathbf{k}}\partial^{{e}}\omega_{\lambda\mu\mathbf{k}}}{\omega_{e}^{2}}$
(38) $\displaystyle\sum_{\eta}\left(\xi_{\lambda\eta\mathbf{k}}^{{e}}{\cal
P}_{\eta\mu\mathbf{k}}^{(1);p[tra]}-{\cal
P}_{\lambda\eta\mathbf{k}}^{(1);p[tra]}\xi_{\eta\mu\mathbf{k}}^{{e}}\right)=-i\xi_{\lambda\mu\mathbf{k}}^{{e}}\frac{\partial^{{p}}f_{\lambda\mu\mathbf{k}}}{\omega_{p}}$
(39) $\displaystyle\sum_{\eta\in v}\xi_{\lambda\eta\mathbf{k}}^{{e}}{\cal
P}_{\eta\mu\mathbf{k}}^{(1);p}=-\sum_{\eta\in
v}\frac{\xi_{\lambda\eta\mathbf{k}}^{{e}}\xi_{\eta\mu\mathbf{k}}^{{p}}f_{\eta\mu\mathbf{k}}}{\omega_{p}+\omega_{\eta\mu\mathbf{k}}}$
(40) $\displaystyle-\sum_{\eta\in c}\xi_{\eta\mu\mathbf{k}}^{{e}}{\cal
P}_{\lambda\eta\mathbf{k}}^{(1);p}=\sum_{\eta\in
c}\frac{f_{\lambda\eta\mathbf{k}}\xi_{\lambda\eta\mathbf{k}}^{{p}}\xi_{\eta\mu\mathbf{k}}^{{e}}}{\omega_{p}+\omega_{\lambda\eta\mathbf{k}}}=-\sum_{\eta\in
c}\frac{f_{\lambda\eta\mathbf{k}}\xi_{\lambda\eta\mathbf{k}}^{{p}}\xi_{\eta\mu\mathbf{k}}^{{e}}}{\omega_{e}+\omega_{\eta\mu\mathbf{k}}}.$
(41)

The terms in the second factor of ${\cal F}^{(3)}_{\lambda\mu\mathbf{k}}$ can
be evaluated as follows:

$\displaystyle-i\partial^{{e}}{\cal V}_{\lambda\mu\mathbf{k}}^{p}$
$\displaystyle=\omega_{p}\left[\frac{1}{\omega_{e}^{2}}\xi_{\mu\lambda\mathbf{k}}^{{p}}\partial^{{e}}\omega_{\mu\lambda\mathbf{k}}-\left(\frac{1}{\omega_{p}}+\frac{1}{\omega_{e}}\right)\partial^{{e}}\xi_{\mu\lambda\mathbf{k}}^{{p}}\right]$
(42) $\displaystyle\sum_{\eta=\lambda,\mu}\Big{(}\dots\Big{)}$
$\displaystyle=\omega_{p}\left[\frac{1}{\omega_{p}^{2}}\xi_{\mu\lambda\mathbf{k}}^{{e}}\partial^{{p}}\omega_{\mu\lambda\mathbf{k}}-i\left(\xi_{\lambda\lambda\mathbf{k}}^{{e}}-\xi_{\mu\mu\mathbf{k}}^{{e}}\right)\left(\frac{1}{\omega_{p}}+\frac{1}{\omega_{e}}\right)\xi_{\mu\lambda\mathbf{k}}^{{p}}\right]$
(43) $\displaystyle\sum_{\eta\notin\\{\lambda,\mu\\}}\Big{(}\dots\Big{)}$
$\displaystyle={\xi}_{\mu\lambda\mathbf{k};{e}}^{p}-{\xi}_{\mu\lambda\mathbf{k};{p}}^{e}+i\omega_{p}\sum_{\eta\notin\\{\lambda,\mu\\}}\left(\frac{\xi_{\mu\eta\mathbf{k}}^{{p}}\xi_{\eta\lambda\mathbf{k}}^{{e}}}{\omega_{p}+\omega_{\eta\mu\mathbf{k}}}+\\{{p}\leftrightarrow{e}\\}\right),$
(44)

where we made use of the following sum rule derived by Aversa and Sipe [13]:

${\xi}_{\mu\lambda\mathbf{k};{p}}^{e}-{\xi}_{\mu\lambda\mathbf{k};{e}}^{p}=i\sum_{\eta\notin\\{\lambda,\mu\\}}\left(\xi_{\mu\eta\mathbf{k}}^{{p}}\xi_{\eta\lambda\mathbf{k}}^{e}-\xi_{\mu\eta\mathbf{k}}^{{e}}\xi_{\eta\lambda\mathbf{k}}^{p}\right),$
(45)

along with the generalized derivative notation [13]

$\displaystyle{\xi}_{\mu\lambda\mathbf{k};{\alpha_{1}}}^{\alpha_{2}}$
$\displaystyle=\partial^{\alpha_{1}}{\xi}_{\mu\lambda\mathbf{k}}^{\alpha_{2}}-i\left(\xi_{\mu\mu\mathbf{k}}^{\alpha_{1}}-\xi_{\lambda\lambda\mathbf{k}}^{\alpha_{1}}\right){\xi}_{\mu\lambda\mathbf{k}}^{\alpha_{2}}.$
(46)

Applying these steps we obtain

$\displaystyle{\cal F}^{(3)}_{\lambda\mu\mathbf{k}}$
$\displaystyle=\omega_{p}f_{\lambda\mu\mathbf{k}}\left[\left(\frac{\partial^{{e}}\xi_{\lambda\mu\mathbf{k}}^{{p}}}{\omega_{e}}+\frac{\xi_{\lambda\mu\mathbf{k}}^{{p}}\partial^{{e}}\omega_{\lambda\mu\mathbf{k}}}{\omega_{e}^{2}}+i\sum_{\eta}\frac{\xi_{\lambda\eta\mathbf{k}}^{{e}}\xi_{\eta\mu\mathbf{k}}^{{p}}}{\omega_{p}+\omega_{\eta\mu\mathbf{k}}}\right)+\\{{p}\leftrightarrow{e}\\}\right]$
$\displaystyle\quad\times\left[\left(\frac{\xi_{\mu\lambda\mathbf{k};{e}}^{{p}}}{\omega_{e}}-\frac{\xi_{\mu\lambda\mathbf{k}}^{{p}}\partial^{{e}}\omega_{\mu\lambda\mathbf{k}}}{\omega_{e}^{2}}-i\\!\\!\\!\sum_{\eta\neq\\{\lambda,\mu\\}}\frac{\xi_{\mu\eta\mathbf{k}}^{{p}}\xi_{\eta\lambda\mathbf{k}}^{{e}}}{\omega_{p}+\omega_{\eta\mu\mathbf{k}}}\right)+\\{{p}\leftrightarrow{e}\\}\right].$
(47)

Some further algebra shows that the second bracket is just the conjugate of
the first one, so that (choosing the form of the second bracket)

$\displaystyle{\cal F}^{(3)}_{\lambda\mu\mathbf{k}}$
$\displaystyle=\omega_{p}f_{\lambda\mu\mathbf{k}}\left|\left(\frac{\xi_{\mu\lambda\mathbf{k};{e}}^{{p}}}{\omega_{e}}-\frac{\xi_{\mu\lambda\mathbf{k}}^{{p}}\partial^{{e}}\omega_{\mu\lambda\mathbf{k}}}{\omega_{e}^{2}}-i\\!\\!\\!\sum_{\eta\neq\\{\lambda,\mu\\}}\frac{\xi_{\mu\eta\mathbf{k}}^{{p}}\xi_{\eta\lambda\mathbf{k}}^{{e}}}{\omega_{p}+\omega_{\eta\mu\mathbf{k}}}\right)+\\{{p}\leftrightarrow{e}\\}\right|^{2}$
(48)

For a two-band model, only the first two terms in the bracket are relevant.
The function ${\cal F}^{(3)}_{\lambda\mu\mathbf{k}}$ is in this form directly
seen to be real-valued and thus it can be substituted into the contour
integral formula (24) to obtain the 2PA coefficient

$\displaystyle\beta=\frac{1}{2}\beta_{pe}\equiv\frac{1}{2}\alpha^{(3)}_{pe}$
$\displaystyle=\frac{\pi}{2}\tilde{K}_{3}\kappa_{D}\
\sum_{\lambda\mu}\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}\frac{{\cal
F}^{(3)}_{\lambda\mu\mathbf{k}}}{\left\|\nabla_{\mathbf{k}}{\omega}_{\mu\lambda\mathbf{k}}\right\|}.$
(49)

The dimension of $\beta$ is that of length over power. It is also directly
obvious from these equations that $\beta_{pe}/\hbar{\omega_{p}}$ is invariant
with respect to an exchange of probe and pump indices
($\\{{p}\leftrightarrow{e}\\}$), which is expected, since it corresponds to
the photon number absorption rate and each beam delivers one photon. This
symmetry is to be distinguished from the better known intrinsic permutation
symmetry [12] of susceptibilities/conductivities, since it additionally relies
on the consistency of the band-structure model. The reason is that we are not
symmetrizing over all indices, but exclude permutations in Eq. (14)
contributing to pump absorption, such as
$\tilde{\sigma}^{(3);epep}(\omega_{e};\omega_{p},\omega_{e},-\omega_{p})$.
Interestingly, this kind of setup/investigation thus provides a check of sum
rules, which are often thought to be more relevant in velocity gauge analyses.

###  III.3 Implementation and local gauge fixing

Figure 1:  Example path (dark blue) in the algorithm for determining an
ensemble of $\mathbf{k}$-points on a resonance contour line (red numbered
points connected by dashed lines).  Figure 2:  Local gauge restoration scheme
to determine band-structure properties at each resonant $\mathbf{k}$-point
(red circles). Eigenfunctions are required for all points of the ‘neighbour
grid’ and fixed in gauge as shown by arrows. Connection elements
$\bm{\xi}_{\mathbf{k}}$ are determined at the nearest neighbours (black
circles), in order to obtain $\nabla_{\mathbf{k}}\bm{\xi}_{\mathbf{k}}$ at the
central point.

In the case of a bulk crystal, the resonance integral is further split into a
$k_{z}$-integral (which is regularly discretized) and a resonance line
integral (for each discrete value of $k_{z}$). The task is now to find the 2D
resonance contour and to determine the band-structure properties for each
discrete point on the contour line.

First we determine the band energies on all points of a regular
$\mathbf{k}$-grid (‘interpolation grid’). In order to find the resonance line
for a particular VB-CB transition, we follow a particular scheme illustrated
in Fig. 1, where the grid (step size d$k$) is coarse for demonstration
purposes. A starting point (here $k_{x}=k_{y}=0$) is required that is
encircled by the searched resonance line. From here, we move straight along
the $k_{x}$-axis until a resonance line is crossed. The precise position is
determined by (linear/cubic) interpolation, and we also determine the normal
vector
$\hat{\bm{n}}_{\mathbf{k}}=\nabla_{\mathbf{k}}\omega_{\lambda\mu\mathbf{k}}/||\nabla_{\mathbf{k}}\omega_{\lambda\mu\mathbf{k}}||$
(within the $k_{x},k_{y}$-plane) at the determined resonance point
$\mathbf{k}_{1}$. This is the first point for the resonance line. The
subsequent points $\mathbf{k}_{i+1}$ are iteratively found by moving on a
circle around the point $\mathbf{k}_{i}$. The radius of this circle is about
the interpolation grid step size. The start angle is chosen as the direction
of the tangent vector at the previous point, i.e., the normal of the
determined $\hat{\bm{n}}_{\mathbf{k}}$. The precise angle is again found using
interpolation, here on a regular angle grid. When
$|\bm{k}_{i}-\bm{k}_{1}|<k_{\mathrm{circ}}$ for some $i>2$, the end of the
closed contour line is reached ($N\equiv i$). The first and last point
($\bm{k}_{1}$ and $\bm{k}_{N}$) are weighted by
$(|\bm{k}_{N}-\bm{k}_{1}|+k_{\mathrm{circ}})/2$ instead of $k_{\mathrm{circ}}$
in order to account for the shorter distance.

The thus produced list of resonant $\mathbf{k}$-points (along with weights) is
fed as input to the electronic structure program, this time determining all
the required properties
($\omega_{\lambda\mathbf{k}},\bm{\xi}_{\mathbf{k}},\nabla_{\mathbf{k}}\omega_{\lambda\mathbf{k}},\nabla_{\mathbf{k}}\bm{\xi}_{\mathbf{k}}$).
In order to derive the connection elements and their derivatives using
balanced difference quotients, we use the local neighbour grid points shown in
Fig. 2. Here the gauge of each point is maximally aligned to its neighbour,
using the overlap matrix as in a parallel transport gauge (PTG), [16] as shown
by the arrows, starting from the central (resonant) $\mathbf{k}$-point. If two
bands are nearly degenerate (at the central $\mathbf{k}$-point), they are
treated as a group, within eigenfunctions are mixed. Note that the derived
properties are not necessarily smooth over the ensemble of resonant
$\mathbf{k}$-points, so this gauge does not work for dynamical simulations,
but in this analytical approach each resonant $\mathbf{k}$-point can be
evaluated independently. Alternatively to our choice of neighbour grid, one
could also use a local grid of 4$\times$4 points with the resonant
$\mathbf{k}$-point in the center (not being a grid point itself), and derive
the required properties by bicubic interpolation. In this scheme, a local
radial gauge could also be applied.

##  IV Application to specific materials

###  IV.1 Graphene (Nearest-neighbour tight-binding model)

In this section we apply our approach to models with analytically known
eigensystem, so that no gauge fixing (like the one described in Fig. 2) is
needed. For monolayer graphene, the evaluation of 1PA and 2PA is fully
analytical. The numerical algorithm for finding the resonance contour lines is
applied for bilayer graphene. Some results in this section are possibly more
straightforward to obtain in the velocity gauge. [17]

####  IV.1.1 Graphene monolayer

Here we evaluate 1PA and 2PA in intrinsic undoped graphene over the full
spectral range of the $\pi,\pi^{*}$-bands. We choose the nearest-neighbour
(NN) tight binding model,
$H=-\gamma_{0}\sum_{\mathbf{R}}\sum_{\alpha}\hat{a}_{\mathbf{R}}^{\dagger}\hat{b}_{\mathbf{R}+\bm{\delta}_{\alpha}}+\mathrm{h.c.}$,
where $\hat{a}_{\mathbf{R}}$ ($\hat{b}_{\mathbf{R}}$) is the annihilation
operator for an electron on site $\mathbf{R}$ belonging to sublattice $A$
($B$), $\bm{\delta}_{\alpha}$ is a NN vector, and $-\gamma_{0}$ is the hopping
constant. For simplicity, we neglect the overlap integral [18] between
orbitals on different sites, even though this is only a good approximation in
the low-energy regime. At the same time, the model still has an improved
validity range compared to the effective Hamiltonian valid near the $K$
points, which has been used previously in the evaluation of 1PA and 2PA
coefficients. [17, 19, 3, 20] Nevertheless, we are able to obtain analytical
results for 1PA and 2PA coefficients.

The parameter $\gamma_{0}$ is usually in the range between 2.5 and 3.0 eV;
here we set it to the value $\gamma_{0}=3.0$ eV. The spin index is omitted;
the spin degeneracy is accounted for by a factor of $g_{s}=2$ at the end. The
kernel of the Fourier transformed Hamiltonian is

$\displaystyle{\cal H}_{\mathbf{k}}$
$\displaystyle=-\gamma_{0}\left(\begin{array}[]{*{2}{c}}0&s_{\mathbf{k}}\\\
s_{\mathbf{k}}^{*}&0\end{array}\right),\quad
s_{\bm{k}}=\sum_{\alpha=1}^{3}e^{-i\mathbf{k}\cdot\bm{\delta}_{\alpha}}.$ (52)

The NN vectors are chosen as

$\displaystyle\bm{\delta}_{1}=\delta\begin{pmatrix}0\\\ 1\end{pmatrix},\
\bm{\delta}_{2}=\delta\begin{pmatrix}-\sqrt{3}/2\\\ -1/2\end{pmatrix},\
\bm{\delta}_{3}=\delta\begin{pmatrix}\sqrt{3}/2\\\ -1/2\end{pmatrix},$ (53)

where $\delta=a_{0}/\sqrt{3}$ is the NN distance and $a_{0}=2.46$ Å is the
lattice constant. The eigenvalues and normalized eigenvectors are

$\displaystyle\hbar\omega_{\eta\mathbf{k}}$
$\displaystyle=\sigma_{\eta}\gamma_{0}\left|s_{\mathbf{k}}\right|,\quad\Psi_{\eta\mathbf{k}}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{*{2}{c}}-\sigma_{\eta}s_{\mathbf{k}}/|s_{\mathbf{k}}|\\\
1\end{array}\right),\quad$ (56)

with $\eta=c,v$ (which are just single band indices in this two-band model),
and $\sigma_{v}=-1$ (VB) and $\sigma_{c}=1$ (CB). The model has electron-hole
symmetry, and the hopping constant enters only by scaling the band energies.
The square modulus of $s_{\mathbf{k}}$ is evaluated as

$\displaystyle|s_{\mathbf{k}}|^{2}$
$\displaystyle=3+2\left\\{\cos\left[\mathbf{k}\cdot(\bm{\delta}_{1}-\bm{\delta}_{2}+0\bm{\delta}_{3})\right]+\mathrm{c.p.}\right\\}$
(57)
$\displaystyle=3+2\cos(a_{0}k_{x})+4\cos(a_{0}k_{x}/2)\cos(\sqrt{3}a_{0}k_{y}/2)$
(58)

and vanishes at the two inequivalent $K$ points (e.g.,
$\mathbf{K}_{1,2}=\pm(1,0)^{\mathsf{T}}4\pi/3a_{0}$). Here, $\mathrm{c.p.}$
stands for cyclic permutation of the NN indices. For the connection elements
we find

$\displaystyle\bm{\xi}_{\eta\eta^{\prime}\mathbf{k}}$
$\displaystyle=\sigma_{\eta}\sigma_{\eta^{\prime}}\frac{i}{4}\frac{s_{\mathbf{k}}^{*}\partial_{\mathbf{k}}s_{\mathbf{k}}-s_{\mathbf{k}}\partial_{\mathbf{k}}s_{\mathbf{k}}^{*}}{|s_{\mathbf{k}}|^{2}}=-\sigma_{\eta}\sigma_{\eta^{\prime}}\frac{\left\\{\bm{\delta}_{3}\cos\left[\mathbf{k}\cdot(\bm{\delta}_{1}-\bm{\delta}_{2})\right]+\mathrm{c.p.}\right\\}}{2|s_{\mathbf{k}}|^{2}}$
(59)
$\displaystyle=\sigma_{\eta}\sigma_{\eta^{\prime}}\frac{1}{2|s_{\mathbf{k}}|^{2}}\begin{pmatrix}\sqrt{3}\sin(a_{0}k_{x}/2)\sin(\sqrt{3}a_{0}k_{y}/2)\\\
-\cos(a_{0}k_{x})+\cos(a_{0}k_{x}/2)\cos(\sqrt{3}a_{0}k_{y}/2)\end{pmatrix}.$
(60)

In this way, the intraband elements are just the negative of the (real)
interband elements,
$\xi_{vv\mathbf{k}}=\xi_{cc\mathbf{k}}=-\xi_{vc\mathbf{k}}=-\xi_{cv\mathbf{k}}$.
Around the $K$ points, the connection elements are a rotation field with a
singularity at $K$. The required $\mathbf{k}$-derivatives of
$\omega_{\eta\mathbf{k}}$ and $\bm{\xi}_{\eta\eta^{\prime}\mathbf{k}}$ are
straightforward to obtain.

Figure 3: (a) Energy contours for the monolayer graphene VB-CB transition,
$2\gamma_{0}|s_{\mathbf{k}}|$, cf. Eq. (56). (b) The CB spectrum of monolayer
(dashed) and bilayer (solid) graphene along the $k_{y}=0$ line, indicated by
the dashed line in (a). The inset shows a zoom close to the $K$ region. The
spectrum is electron-hole symmetric.

Using the analytically derived structure input, it is possible to directly
evaluate the $k_{\perp}$-integral (23), which is easier than the resonance
contour integral (24). We choose $\hat{\mathbf{k}}_{\parallel}=\hat{\bm{y}}$,
so we have (setting $f_{vc\mathbf{k}}=1$ for the ground state)

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$ $\displaystyle=\pi\tilde{K}_{n}\
\kappa_{2}\int\mathrm{d}k_{x}\sum_{k_{y,\mathrm{res}}}\frac{{\cal
F}^{(n)}_{vc\mathbf{k}}}{\left|\partial^{y}{\omega}_{vc\mathbf{k}}\right|}.$
(61)

For the following analysis it is convenient to introduce a dimensionless
frequency variable $\zeta\equiv\hbar\omega_{\Sigma}/2\gamma_{0}$. The
resonance condition $\omega_{cv\mathbf{k}}=\omega_{\Sigma}$ is equivalent to

$\displaystyle\cos(\sqrt{3}a_{0}k_{y}/2)$
$\displaystyle=\frac{\zeta^{2}-3-2\cos(a_{0}k_{x})}{4\cos(a_{0}k_{x}/2)},$
(62)

which, together with trigonometric identities, allows us to eliminate $k_{y}$
from Eq. (61). For $0<\zeta<1$ the resonance contours encircle the $K$ points,
while for $1<\zeta<3$ they encircle the $\Gamma$ point, see Fig. 3(a). Further
using multiple-angle formulae, the dependence on $k_{x}$ is effectively
through $u\equiv 2\cos(a_{0}k_{x}/2)$, which we substitute. This gives

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$ $\displaystyle=\pi g_{s}\tilde{K}_{n}\
\kappa_{2}{\cal I}^{(n)},$ (63) $\displaystyle{\cal I}^{(n)}$
$\displaystyle\equiv
4\int_{k_{x,\mathrm{min}}}^{k_{x,\mathrm{max}}}\mathrm{d}k_{x}\frac{{\cal
F}^{(n)}_{\lambda\mu\mathbf{k}}}{\left|\partial^{y}{\omega}_{\lambda\mu\mathbf{k}}\right|}=4\int_{u_{\mathrm{min}}}^{u_{\mathrm{max}}}\frac{\mathrm{d}u}{a_{0}\sqrt{1-u^{2}/4}}\frac{{\cal
F}^{(n)}_{\lambda\mu\mathbf{k}}}{\left|\partial^{y}{\omega}_{\lambda\mu\mathbf{k}}\right|},$
(64)

where the integration ranges have been chosen as

$\displaystyle[u_{\mathrm{min}},u_{\mathrm{max}}]=\left\\{\begin{array}[]{*{2}{l}}[-1-\zeta,-1+\zeta],&\quad
0<\zeta<1,\\\\[3.00003pt] [-1+\zeta,2],&\quad 1<\zeta<3.\end{array}\right.$
(67)

This corresponds to the displayed sections of contour lines below the straight
line connecting the $M$ points in Fig. 3(a). The factor 4 accounts for the
other sections / the other valley.

We start the evaluation with the linear absorption, $n=1$. The linear
conductivity for a hexagonal lattice is isotropic for normal incidence of
light onto the graphene plane, since
${\sigma}^{(1)}_{xx}={\sigma}^{(1)}_{yy}$. The direction $\hat{\mathbf{p}}$
can thus be chosen arbitrarily in the $xy$-plane, even though the integrand in
${\cal I}^{(n)}$ will differ. Here we choose $\hat{\mathbf{p}}=\hat{\bm{x}}$,
corresponding to the zigzag direction, and obtain

$\displaystyle{\cal I}^{(1)}$
$\displaystyle=\int_{u_{\mathrm{min}}}^{u_{\mathrm{max}}}\mathrm{d}u\frac{\sqrt{4-u^{2}}}{2\sqrt{3}\zeta^{2}u^{2}}\sqrt{4u^{2}-\left(1+u^{2}-\zeta^{2}\right)^{2}}.$
(68)

The integral can be transformed to a sum of well-known elliptical integrals,
which for shortness is not shown here explicitly. The low-frequent limit,
however, is given by a constant,

$\displaystyle\lim_{\zeta\to 0}{\cal I}^{(1)}$ $\displaystyle=\lim_{\zeta\to
0}\frac{1}{2\zeta^{2}}\int_{-1-\zeta}^{-1+\zeta}\mathrm{d}u\sqrt{4u^{2}-\left(1+u^{2}-\zeta^{2}\right)^{2}}$
(69) $\displaystyle=\lim_{\zeta\to
0}\int_{-1}^{1}\mathrm{d}x\sqrt{(x-\zeta^{-1})^{2}-\frac{1}{4}\left(\zeta^{-1}+\zeta(x-\zeta^{-1})^{2}-\zeta\right)^{2}}$
(70) $\displaystyle=\int_{-1}^{1}\mathrm{d}x\sqrt{1-x^{2}}=\frac{\pi}{2}.$
(71)

In the second equation we substituted $u=x\zeta-1$. The linear absorption is

$\displaystyle\alpha^{(1)}$
$\displaystyle=K_{1}\mathrm{Re}\left[{\sigma}^{(1)}_{xx}\right]=\pi
g_{s}K_{1}\frac{e_{0}^{2}}{\hbar}\kappa_{2}{\cal I}^{(1)}.$ (72)

If we now set $\kappa_{2}\equiv(2\pi)^{-{2}}$ and also substitute
$K_{1}=1/(n_{0{p}}\epsilon_{0}c)$, see Eq. (15), we obtain a 2D absorption
coefficient $\alpha^{(1)}_{\mathrm{2D}}$ and sheet conductance
$\sigma^{(1)}_{\mathrm{2D}}$ as

$\displaystyle\alpha^{(1)}_{\mathrm{2D}}$
$\displaystyle=\frac{1}{n_{0{p}}\epsilon_{0}c}\sigma^{(1)}_{\mathrm{2D}},\qquad\sigma^{(1)}_{\mathrm{2D}}=\frac{e_{0}^{2}}{4\hbar}\times\frac{{\cal
I}^{(1)}}{\pi/2},$ (73)

respectively. In the low-frequent (linear dispersion) regime we substitute Eq.
(71) and obtain the universal value
$\sigma^{(1)}_{\mathrm{2D}}=e_{0}^{2}/4\hbar$. [21] The 2D absorption
coefficient $\alpha^{(1)}_{\mathrm{2D}}$ is dimensionless and, since the thin-
film limit applies, corresponds to the relative intensity change
$\mathrm{d}I_{p}/I_{p}$. For a freely suspended graphene sheet,
$n_{0{p}}\approx 1$ is often used, and in the low-frequent limit results in
$\alpha^{(1)}_{\mathrm{2D}}\approx 0.0229$, in agreement with the literature.
[22, 21, 23] Obviously these low-frequent results could have been derived from
the same approach but using the graphene Hamiltonian valid near $K$ points,
${\cal H}_{\mathbf{k}}=\hbar v_{\mathrm{F}}\bm{\sigma}\cdot\mathbf{k}$, with
$v_{\mathrm{F}}=\sqrt{3}a_{0}\gamma_{0}/2\hbar$ the Fermi velocity and
$\bm{\sigma}$ the 2D vector of Pauli matrices in sublattice space.

The result (73), with the ${\cal I}^{(1)}$ integral carried out numerically,
is shown in Fig. 4. Starting from the universal value $e_{0}^{2}/4\hbar$ in
the linear dispersion regime, it increases towards a singularity at
$\hbar\omega_{p}=2\gamma_{0}$, corresponding to the van-Hove singularity
related with the saddle point at $M$. Above this peak, which becomes finite in
the presence of damping, the conductivity decays and reaches zero at the
maximal transition energy $\hbar\omega_{p}=6\gamma_{0}$.

The third-order conductivity tensor $\sigma^{(3)}$ of the hexagonal crystal
class (point group $6mm$) has 21 nonzero elements, of which 10 are
independent. [12] Only the ones in the $xy$-plane are relevant here:

$\displaystyle\sigma^{(3)}_{xxxx}$
$\displaystyle=\sigma^{(3)}_{yyyy}=\sigma^{(3)}_{xxyy}+\sigma^{(3)}_{xyxy}+\sigma^{(3)}_{xyyx}$
(74) $\displaystyle\sigma^{(3)}_{xxyy}$ $\displaystyle=\sigma^{(3)}_{yyxx},$
(75) $\displaystyle\sigma^{(3)}_{xyxy}$ $\displaystyle=\sigma^{(3)}_{yxyx},$
(76) $\displaystyle\sigma^{(3)}_{xyyx}$ $\displaystyle=\sigma^{(3)}_{yxxy}.$
(77)

We consider the degenerate 2PA only, $\omega_{p}=\omega_{e}\equiv\omega$, for
which intrinsic permutation symmetry further implies
$\sigma^{(3)}_{xyxy}=\sigma^{(3)}_{xxyy}$. The pump-probe absorption
coefficient for an arbitrary linearly polarized pump and probe
($\hat{\mathbf{p}}=\hat{\bm{x}}\cos(\theta_{p})+\hat{\bm{y}}\sin(\theta_{p})$)
is

$\displaystyle\beta=\frac{1}{2}\beta_{pe}$
$\displaystyle=\frac{1}{2\epsilon_{0}^{2}n_{0}^{2}c^{2}}\left(\sigma^{\prime(3)}_{xxxx}\cos^{2}(\theta_{p}-\theta_{e})+\sigma^{\prime(3)}_{xyxy}\sin^{2}(\theta_{p}-\theta_{e})\right),$
(78)

where single (double) prime denotes real (imaginary) part. This pump-probe
response is isotropic, i.e., unchanged by a rotation of the graphene sheet
with respect to the normal axis. An evaluation of the co- and cross-polarized
2PA is obviously sufficient to determine the complete third-order response,
which we do here. With one valence- and one conduction band in this model, the
third term in the bracket in Eq. (48) vanishes. Since both intraband elements
are identical ($\bm{\xi}_{vv\mathbf{k}}=\bm{\xi}_{cc\mathbf{k}}$), the
generalized derivative becomes the usual one. In the co-polarized case,
denoted by the subscript $\parallel$, we obtain

$\displaystyle{\cal I}^{(3)}_{\parallel}$
$\displaystyle=\int_{u_{\mathrm{min}}}^{u_{\mathrm{max}}}\mathrm{d}u\frac{a_{0}^{2}\hbar^{2}}{4\sqrt{3}\gamma_{0}^{2}u^{4}\zeta^{8}\sqrt{4-u^{2}}}\sqrt{4u^{2}-\left(1+u^{2}-\zeta^{2}\right)^{2}}$
$\displaystyle\qquad\qquad\times\left(8+2u^{4}-8\zeta^{2}+u^{2}\left(\zeta^{2}-10\right)\right)^{2}.$
(79)

Here we have chosen $\hat{\mathbf{p}}=\hat{\mathbf{e}}=\hat{\bm{x}}$; as
above, for other directions, the integrand differs but the integral remains
the same. As for 1PA, the integral can be reduced to a sum of elliptical
integrals, but we choose to evaluate it numerically. In the low-frequent
regime it scales as $\zeta^{-4}$ with a constant prefactor given by

$\displaystyle\lim_{\zeta\to 0}\zeta^{4}{\cal I}^{(3)}_{\parallel}$
$\displaystyle=\lim_{\zeta\to
0}\frac{a_{0}^{2}\hbar^{2}}{12\gamma_{0}^{2}\zeta^{4}}\int_{-1-\zeta}^{-1+\zeta}\mathrm{d}u\sqrt{4u^{2}-\left(1+u^{2}-\zeta^{2}\right)^{2}}$
$\displaystyle\qquad\times\left(8+2u^{4}-8\zeta^{2}+u^{2}\left(\zeta^{2}-10\right)\right)^{2}$
(80) $\displaystyle=\lim_{\zeta\to
0}\frac{a_{0}^{2}\hbar^{2}}{12\gamma_{0}^{2}}\int_{-1}^{1}\mathrm{d}x\sqrt{4(x-\zeta^{-1})^{2}-\left(\zeta^{-1}+\zeta((x-\zeta^{-1})^{2}-1)\right)^{2}}$
$\displaystyle\qquad\times\left(\frac{8}{\zeta}+2\zeta^{3}(x-\zeta^{-1})^{4}-8\zeta+\zeta(x-\zeta^{-1})^{2}\left(\zeta^{2}-10\right)\right)^{2}$
(81)
$\displaystyle=\frac{24a_{0}^{2}\hbar^{2}}{\gamma_{0}^{2}}\int_{-1}^{1}\mathrm{d}x\sqrt{1-x^{2}}x^{2}=\frac{3\pi
a_{0}^{2}\hbar^{2}}{\gamma_{0}^{2}},$ (82)

with the same substitution $u=x\zeta-1$ as for 1PA. The explicit expression
for the cross-polarized result ${\cal I}^{(3)}_{\perp}$ is skipped here. We
note, however, that its low-frequent limit is identical with the co-polarized
case.

From Eq. (63), the degenerate 2PA coefficient is

$\displaystyle\beta_{\mathrm{2D}}=\frac{1}{2}\alpha^{(3)}_{pe}$
$\displaystyle=g_{s}\frac{\pi}{2}\tilde{K}_{3}\ \kappa_{2}{\cal
I}^{(3)}=\frac{e_{0}^{4}}{8\pi
n_{0{p}}n_{0{e}}\hbar^{3}}\left(\frac{1}{\epsilon_{0}c}\right)^{2}{\cal
I}^{(3)}.$ (83)

In the low-frequent regime, substituting Eq. (82) we find the scaling
$\beta_{\mathrm{2D}}=\alpha^{(3)}_{pe}/2=\bar{K}/\omega^{4}$ with the roughly
constant prefactor

$\displaystyle\bar{K}$
$\displaystyle\equiv\frac{e_{0}^{4}}{n_{0{p}}n_{0{e}}\hbar^{5}}\left(\frac{1}{\epsilon_{0}c}\right)^{2}\frac{3a_{0}^{2}\gamma_{0}^{2}}{8}=\frac{1}{2n_{0{p}}n_{0{e}}\hbar^{3}}\left(\frac{v_{\mathrm{F}}e_{0}^{2}}{\epsilon_{0}c}\right)^{2}.$
(84)

Such scaling of the 2PA coefficient with the excitation frequency has been
found earlier [17, 19].

Multiplying the absorption coefficient $\beta_{\mathrm{2D}}$, with units
[area/power], by the pump intensity $I_{e}$ yields half the relative
absorption $\mathrm{d}I_{p}/I_{p}$ of the probe by the 2D material due to the
presence of the pump. Dividing $\beta_{\mathrm{2D}}$ by a thickness, i.e.,
setting $\kappa_{2}=(2\pi)^{-2}L_{\mathrm{c}}^{-1}$, yields a bulk 2PA
coefficient. For the thickness we substitute the layer spacing in graphite,
$L_{\mathrm{c}}=3.3$Å. The result for $\beta$ is shown in Fig. 5 (a), along
with $\beta_{\mathrm{2D}}\omega^{4}/\bar{K}$ in Fig. 5 (b). Here we used a
frequency-independent value $n_{0}=\sqrt{3}$. [19] The characteristic scaling
with $\omega^{-4}$ holds for photon energies smaller than $\gamma_{0}$, where
co- and cross-polarized 2PA is identical. In the vicinity of
$\hbar\omega=\gamma_{0}$ there is some enhancement. In this regime the co-
polarized 2PA becomes about twice as large as the cross-polarized one. At
higher energies the 2PA coefficients decrease and drop to zero at the band
edge, $\hbar\omega=3\gamma_{0}$.

Figure 4: Linear optical conductivity of mono- and bilayer graphene. The
monolayer result corresponds to Eq. (73).

Figure 5: Degenerate 2PA coefficients of mono- and bilayer graphene for co-
and cross-polarized excitation. These are isotropic, i.e., independent of the
light polarization with respect to the crystallographic axes. In (a), $\beta$
is shown as a function of photon wavelength, while in (b) the same data is
shown using a different scaling and as a function of energy. In (a), the
dotted line is $\beta=(\bar{K}/L_{\mathrm{c}})\omega^{-4}$, the dash-dotted
line is $\sim\omega^{-3}$. In (b), the 2D absorption coefficient (i.e., not
divided by $L_{\mathrm{c}}$), normalized by the long-wavelength result for the
monolayer, is shown. The horizontal lines correspond to the values 1 and 2.

####  IV.1.2 Graphene bilayer

A similar evaluation is now carried out for Bernal-stacked (AB-stacked)
bilayer graphene. The analogous model has four bands and serves as an example
of a multiband model. We consider only NN interlayer hopping characterized by
the constant $\gamma_{1}$. The bilayer Hamiltonian is

$\displaystyle{\cal H}_{\mathbf{k}}$
$\displaystyle=\left(\begin{array}[]{*{4}{c}}0&-\gamma_{0}s_{\mathbf{k}}&0&\gamma_{1}\\\
-\gamma_{0}s_{\mathbf{k}}^{*}&0&0&0\\\ 0&0&0&-\gamma_{0}s_{\mathbf{k}}\\\
\gamma_{1}&0&-\gamma_{0}s_{\mathbf{k}}^{*}&0\end{array}\right).$ (89)

For simplicity, we again set the overlap matrix to unity and keep the same
value for $\gamma_{0}$ as for monolayer in Sec. IV.1.1, while in reality it is
slightly larger. For $\gamma_{1}$ we use a value in the typical range, namely
0.4 eV. [24, 18] As shown in Fig. 3(b), the two VBs/CBs in this model follow
the monolayer bands with an offset of $\gamma_{1}/2$ in nearly the entire
Brillouin zone. In the vicinity of $K$ all four bands are parabolic instead of
linear.

An analytical solution (up to the final $k_{\parallel}$-integration) as for
monolayer graphene might be possible, but leads to lengthy expressions.
Instead we evaluate the absorption coefficients numerically, using the contour
integral method. The eigenvalues and -vectors of ${\cal H}_{\mathbf{k}}$ are
entered as analytical expressions so that no gauge fixing is needed. The
derivatives
($\xi_{\mathbf{k}},\partial\xi_{\mathbf{k}},\partial\omega_{\mathbf{k}}$) are
computed as balanced difference quotients (see Sec. III.3).

Even though the crystal structure of bilayer graphene belongs to a different
point group, its relevant tensor properties are identical with the monolayer
ones. Therefore the response in first and third order is isotropic, so that we
may again choose $\hat{\mathbf{p}}=\hat{\bm{x}}$ and
$\hat{\mathbf{e}}=\hat{\bm{x}}$ ($\hat{\mathbf{e}}=\hat{\bm{y}}$) for the co-
(cross-) polarized case. The bilayer 1PA and 2PA results are shown along with
the monolayer results in Figs. 4 and 5. To obtain a bulk absorption
coefficient with units [length/power], we use the same dielectric constant as
for monolayer, $n_{0}=\sqrt{3}$, and the doubled thickness,
$L_{\mathrm{c}}=6.6$Å.

For long wavelengths the linear conductivity, depicted in Fig. 4, is on
average twice as strong as that of a monolayer. However, a singularity along
with a jump exists at $\hbar\omega_{p}=\gamma_{1}$. There is no singularity at
$\hbar\omega_{p}=2\gamma_{0}$ as for the monolayer, but instead two
singularities occur at $2\gamma_{0}\pm\gamma_{1}$. Apart from that, the 1PA
spectra of mono- and bilayer look similar as expected from the band structure.

By contrast, the bilayer degenerate 2PA coefficient, plotted in Fig. 5,
differs from the monolayer result already in the low-frequent regime, where it
scales with $\omega^{-3}$ for $\hbar\omega\ll\gamma_{1}$. [17, 19] In the
range $\gamma_{1}/2<\hbar\omega<\gamma_{1}$ the co- and cross-polarized 2PA
strongly differ. At $\hbar\omega=\gamma_{1}$ there is a resonance peak, which
in our case of absent scattering is a singularity. These findings are in
agreement with the results of Refs. [17, 19]. For higher photon energies, the
bilayer 2PA (bulk) coefficients are very similar to the monolayer results.
This means that $\beta_{\mathrm{2D}}$ of the bilayer is about twice as large
as of the monolayer, which can be seen in Fig. 5(b). We thus cannot confirm
the conclusion of Ref. [19] that the bilayer 2PA is two orders of magnitudes
stronger than that of the monolayer in the VIS/NIR regime, however, this
finding seems to be also in disagreement with the results of Ref. [17].

###  IV.2 Zincblende semiconductors (${\bm{k\cdot}}{\bm{p}}$ theory)

In this section we apply our analytical approach to more complicated band-
structure models. In particular, we consider Zincblende semiconductors
described by well-established ${\bm{k\cdot}}{\bm{p}}$ models. First we analyze
the polarization dependence of $\beta$, which can be compared with the tensor
symmetry properties for the given crystal class, and later we compare $\beta$
with results from dynamical simulations.

####  IV.2.1 Polarization dependence of degenerate 2PA

10(4.1,-1.2) Fit parameters:

Figure 6:  2PA coefficient of Zincblende semiconductors as a function of the
probe polarization direction ($\theta_{p}$). The pump is co-polarized
($\theta_{e}\\!=\\!\theta_{p}$) or cross-polarized
($\theta_{e}\\!=\\!\theta_{p}\\!+\\!\pi/2$); the angle $\theta$ may be
interpreted here as either $\theta_{p}$, or, to match with Ref. [6],
$\theta_{e}$.

A set of 2PA measurements has been carried out and analyzed by Dvorak et al
[6] with bulk Zincblende semiconductors. From this set we consider the
polarization-dependent degenerate 2PA coefficient of GaAs(001) (at room
temperature), with $\hbar\omega_{p}=\hbar\omega_{e}=1.3$eV. The band-structure
input is obtained from the 30-band ${\bm{k\cdot}}{\bm{p}}$-model, [25] with
the band gap set to $\hbar\omega_{\mathrm{g}}=1.42$eV.

The third-order susceptibility tensor $\chi^{(3)}$ of Zincblende structures
(crystal group $\bar{4}3m$) has 21 nonzero elements, of which 4 are
independent. [12] Since we choose the propagation direction of both pulses in
$z$ (in crystal coordinates, i.e., [001]), the polarization unit vectors lie
in the $xy$-plane and we need to consider only elements therein, namely

$\displaystyle\sigma^{(3)}_{xxxx}$ $\displaystyle=\sigma^{(3)}_{yyyy},$ (90)
$\displaystyle\sigma^{(3)}_{xxyy}$ $\displaystyle=\sigma^{(3)}_{yyxx},$ (91)
$\displaystyle\sigma^{(3)}_{xyxy}$ $\displaystyle=\sigma^{(3)}_{yxyx},$ (92)
$\displaystyle\sigma^{(3)}_{xyyx}$ $\displaystyle=\sigma^{(3)}_{yxxy}.$ (93)

In the degenerate case studied here, we further have
$\sigma^{(3)}_{xyxy}=\sigma^{(3)}_{xxyy}$ due to intrinsic permutation
symmetry. The angle between the pump polarization and the $x$-axis ([100]) is
denoted by $\theta$. The probe is chosen as either co- or cross-polarized.
Following Ref. [6], the corresponding 2PA coefficients are found as

$\displaystyle\beta_{\parallel}$
$\displaystyle=\frac{1}{2}\beta_{pe}^{\parallel}=\frac{1}{2\epsilon_{0}^{2}n_{0}^{2}c^{2}}\sigma^{\prime(3)}_{xxxx}\left[1+\bar{\sigma}\left(\cos^{4}\theta+\sin^{4}\theta-1\right)\right],$
(94) $\displaystyle\beta_{\perp}$
$\displaystyle=\frac{1}{2}\beta_{pe}^{\perp}=\frac{1}{2\epsilon_{0}^{2}n_{0}^{2}c^{2}}\left(2\bar{\sigma}\sigma^{\prime(3)}_{xxxx}\cos^{2}\theta\sin^{2}\theta+\sigma^{\prime(3)}_{xyxy}\right),$
(95)

where the anisotropy parameter is given by

$\bar{\sigma}=1-\left(\sigma^{\prime(3)}_{xyyx}+2\sigma^{\prime(3)}_{xyxy}\right)/\sigma^{\prime(3)}_{xxxx}.$
(96)

Remember that the Cartesian indices are ordered as
${\sigma}^{(n)}_{\beta\alpha_{1}\dots\alpha_{n}}$, not
${\sigma}^{(n)}_{\beta\alpha_{n}\dots\alpha_{1}}$ as in Ref. [6].

Figure 6 shows the calculated 2PA coefficients as scatter plot (circles),
along with fits using Eqs. (94) and (95) by solid lines. The agreement of the
fit curves with the data is excellent, meaning that the model is in perfect
agreement with the crystal symmetry. The parameters extracted from the fit to
our numerical data are listed next to the plot in Fig. 6. The first and last
parameter have been obtained from $\beta_{\parallel}$, the second and third
from $\beta_{\perp}$. The dashed lines are the fits to experimental data from
Ref. [6]; the experimental data itself is not shown. The magnitudes of the
individual $\chi^{(3)}$-tensor elements agree with the experiment [6], at
least they are within the uncertainty ranges. The anisotropy parameter, even
though in rough agreement, is slightly too high compared with the experimental
value of $-0.76\pm 0.08$. In the plot in Fig. 6, our cross-polarized result
can be seen to be nearly identical with Ref. [6], while the co-polarized curve
is somewhat smaller than the experimental one but still in rather good
agreement.

####  IV.2.2 Dependence of 2PA on the non-degeneracy parameter

Here we carry out a comparison with dynamical simulations of non-degenerate
2PA in bulk GaAs and ZnSe. For ZnSe these simulations have been published in a
recent joint experimental and theoretical study [15], and we start with this
material. The dynamical simulations are just numerical solutions of the
perburbatively expanded SBE (1) for finite pulses. There, dephasing and
relaxation processes have been neglected as it is done here. Gaussian pulses
with FWHM of $\sim 30$ fs are shorter than the experimental ones but
sufficiently long so that the absorption spectrum does not change
significantly for longer pulses. The incidence of laser pulses is
(approximately) normal to the sample surface, so the label ‘(100)-sample’
refers to light propagation direction
$\hat{\bm{k}}_{p}(\approx\hat{\bm{k}}_{e})\parallel[100]$. The individual
polarization settings are as follows (illustrated in Ref. [15]):

  * •

(100)-sample: $\hat{\mathbf{p}}\parallel[011]$,
$\hat{\mathbf{e}}_{\perp}\parallel[0\bar{1}1]$

  * •

(110)-sample: $\hat{\mathbf{p}}\parallel[1\bar{1}0]$,
$\hat{\mathbf{e}}_{\perp}\parallel[001]$.

Here $\hat{\mathbf{e}}_{\perp}$ is defined as the direction perpendicular to
both $\hat{\mathbf{p}}$ and $\hat{\bm{k}}_{p}$. For co-polarized 2PA,
$\hat{\mathbf{e}}=\hat{\mathbf{p}}$, and for cross-polarized 2PA,
$\hat{\mathbf{e}}=\hat{\mathbf{e}}_{\perp}$. Both theoretical approaches are
based on the same eight-band Kane model, [26] where for convenience the
(100)-sample is replaced by the equivalent setting
$\hat{\bm{k}}_{p}\parallel[001]$, $\hat{\mathbf{p}}\parallel[110]$,
$\hat{\mathbf{e}}_{\perp}\parallel[\bar{1}10]$. Also, for the (110)-sample we
perform a rotation of the coordinate system within the
${\bm{k\cdot}}{\bm{p}}$-model, such that the $z$-axis corresponds to [110]. As
in Ref. [15] we model the frequency-dependent refractive index $n_{0}(\omega)$
for ZnSe by $n_{0}^{2}(\omega)=A+{B\lambda^{2}}/({\lambda^{2}-C^{2}})$, where
$\lambda$ is the wavelength in micrometers, and the parameters have the values
$A=4.00$, $B=1.90$, and $C^{2}=0.113$. [27]

Figure 7 shows the 2PA coefficient $\beta=\beta_{pe}/2$ in the pump-probe
setting with fixed frequency sum $\omega_{\Sigma}=3.098$eV$/\hbar$. The
results of the dynamical simulations can be seen to be nearly identical to the
analytical approach. The analytical approach describes a steady state response
to monochromatic plane waves, since we evaluate the 2PA result (49) at the
discrete carrier frequencies only. In Fig. 8 the corresponding ratio
$\beta_{\perp}/\beta_{\parallel}$ is plotted from the same data. Here one can
distinguish a slightly better fit of the analytical approach to the
measurement in the case of the (110)-sample.

Figure 7: Dependence of bulk ZnSe 2PA coefficients
$\beta_{\perp}$,$\beta_{\parallel}$ on the non-degeneracy parameter. The sum
frequency corresponds to $\hbar\omega_{\Sigma}=3.098$eV. The experimental data
labeled ‘data’ is the one shown in Ref. [15]; no error bars are available. The
polarization settings for each sample are also defined therein. The curves
labeled ‘cw’ are obtained from the analytical approach presented in this work,
using the same 8-band ${\bm{k\cdot}}{\bm{p}}$ model as in Ref. [15]. The
curves labeled ‘pulses’ are the results from dynamical simulations already
shown in Ref. [15], but corrected by a small post-processing mistake, a
constant factor $\pi/3$. All curves and data are for room temperature. The
curve labeled ‘2-band’ is a fit of $\beta_{\parallel}\sim
n_{0p}^{-1}n_{0e}^{-1}\omega_{p}^{-3}\omega_{e}^{-4}$ to the ‘cw’-curve.
Figure 8: Same data as in Fig. 7 but plotting the ratio of cross- and co-
polarized 2PA coefficients $\beta_{\perp}/\beta_{\parallel}$. See the caption
of Fig. 7 for details.

In Fig. 9 a similar comparison is carried out for bulk GaAs, also modeled by
the eight-band Kane model. [26] The four polarization settings are
(illustrated in Ref. [28]):

  * •

(100)-sample: $\hat{\mathbf{p}}\parallel[010]$,
$\hat{\mathbf{e}}_{\perp}\parallel[001]$

  * •

(110)-sample: $\hat{\mathbf{p}}\parallel[001]$,
$\hat{\mathbf{e}}_{\perp}\parallel[1\bar{1}0]$.

Again, the (100)-sample is replaced by an equivalent setting
$\hat{\bm{k}}_{p}\parallel[001]$, $\hat{\mathbf{p}}\parallel[100]$,
$\hat{\mathbf{e}}_{\perp}\parallel[010]$. Since the variation of refractive
index $n_{0}(\omega)$ in the considered light frequency range ($\sim
0.4-1.2$eV) is small, we have set it to a constant value of 3.4 for
simplicity. The sum frequency corresponds to $\hbar\omega_{\Sigma}=1.5694$eV.
Even though the Gaussian pulses have a duration of only $\sim 9$ fs (FWHM),
there is very good agreement between both approaches. Corresponding
experimental data can be found in Ref. [28] and qualitatively agrees with our
theory, both in order of magnitude and the relative dependence on
$\omega_{p}/\omega_{e}$. Surprisingly, the two sample orientations /
polarization settings give nearly identical theoretical curves for each
$\beta_{\parallel}$ and $\beta_{\perp}$. This is not quite the case in the
experiment [28], where for the (100)-sample the overall ratio
$\beta_{\perp}/\beta_{\parallel}$ remains at about $1/3$. On the other hand,
for the (110)-sample, both in experiment and theory this ratio starts at
around $1/2$ and slightly increases for larger $\omega_{p}/\omega_{e}$.

Figure 9: Comparison of the analytical approach with dynamical simulations as
in Fig. 7, for bulk GaAs with $\hbar\omega_{\Sigma}=1.5694$eV. The ‘pulses’
curves have some uncertainty and are spline fits to only three discrete data
points per curve, since the computational expense is large.

##  V Conclusions

We have investigated multi-photon absorption within a pump-probe setting by
semiconductors and graphene in the independent-particle approximation. The
basis of our numerical evaluations is a rather compact analytical expression
for the 2PA coefficient, which is expressed as a contour integral over the
resonance isosurfaces in $\mathbf{k}$-space. In principle, similar expressions
could be derived for $\ell>2$. Comparing Eqs. (28) and (48), we note that the
number of terms within the square modulus of ${\cal
F}^{(n)}_{\lambda\mu\mathbf{k}}$ just corresponds to $n$, so it might be that
3PA can be expressed by just five terms (some of which are sums over bands).

In order to verify the correctness of this approach, in the evaluation we have
focused on materials that have been analyzed previously by other means, both
experimentally and theoretically. The analytical result for the 1PA and 2PA
coefficients of monolayer (intrinsic) graphene in the nearest-neighbour tight-
binding approximation could be compared with the literature in the low-
frequent regime. The analysis for bulk Zincblende semiconductor has been based
on ${\bm{k\cdot}}{\bm{p}}$ theory. Here we have shown that the results are in
good agreement with dynamical simulations of the 2PA of finite-length pulses.
For sufficiently small spectral width of the pulses, the dynamical simulation
is directly comparable with the steady-state response at the carrier
frequencies as we have done here, instead of integrating over their spectra.
It should be possible to implement our algorithm in DFT programs, which would
open up its application to a wider range of materials.

## Acknowledgements

We acknowledge financial support from the Guangdong Province Science and
Technology Major Project (Future functional materials under extreme conditions
- 2021B0301030005).

## References

  * Aversa  _et al._ [1994] C. Aversa, J. E. Sipe, M. Sheik-Bahae, and E. W. van Stryland, Third-order optical nonlinearities in semiconductors: The two-band model, Phys. Rev. B 50, 18073 (1994).
  * Hannes and Meier [2019] W.-R. Hannes and T. Meier, Higher-order contributions and nonperturbative effects in the nondegenerate nonlinear optical absorption of semiconductors using a two-band model, Phys. Rev. B 99, 125301 (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).
  * Attaccalite  _et al._ [2018] C. Attaccalite, M. Grüning, H. Amara, S. Latil, and F. m. c. Ducastelle, Two-photon absorption in two-dimensional materials: The case of hexagonal boron nitride, Phys. Rev. B 98, 165126 (2018).
  * Hutchings and Wherrett [1994] D. C. Hutchings and B. S. Wherrett, Theory of anisotropy of two-photon absorption in zinc-blende semiconductors, Phys. Rev. B 49, 2418 (1994).
  * Dvorak  _et al._ [1994] M. D. Dvorak, W. A. Schroeder, D. R. Andersen, A. L. Smirl, and B. S. Wherrett, Measurement of the anisotropy of two-photon absorption coefficients in zincblende semiconductors, IEEE Journal of Quantum Electronics 30, 256 (1994).
  * Bolger  _et al._ [1993] J. A. Bolger, A. K. Kar, B. S. Wherrett, R. DeSalvo, D. C. Hutchings, and D. J. Hagan, Nondegenerate two-photon absorption spectra of ZnSe, ZnS and ZnO, Optics Communications 97, 203 (1993).
  * Hutchings and Stryland [1992] D. C. Hutchings and E. W. V. Stryland, Nondegenerate two-photon absorption in zinc blende semiconductors, J. Opt. Soc. Am. B 9, 2065 (1992).
  * Murayama and Nakayama [1994] M. Murayama and T. Nakayama, Two-photon-absorption spectra originating from higher-energy transitions, Phys. Rev. B 49, 5737 (1994).
  * Murayama and Nakayama [1995] M. Murayama and T. Nakayama, Ab initio calculations of two-photon absorption spectra in semiconductors, Phys. Rev. B 52, 4986 (1995).
  * Murayama and Nakayama [1997] M. Murayama and T. Nakayama, Symmetry-induced anisotropy of two-photon absorption spectra in zinc-blende semiconductors, Phys. Rev. B 55, 9628 (1997).
  * Boyd [2008] R. W. Boyd, in  _Nonlinear Optics (Third Edition)_, edited by R. W. Boyd (Academic Press, Burlington, 2008).
  * Aversa and Sipe [1995] C. Aversa and J. E. Sipe, Nonlinear optical susceptibilities of semiconductors: Results with a length-gauge analysis, Phys. Rev. B 52, 14636 (1995).
  * Passos  _et al._ [2021] D. J. Passos, G. B. Ventura, J. M. B. L. dos Santos, and J. M. V. P. Lopes, Nonlinear optical conductivity of a two-band crystal i, Journal of Physics: Condensed Matter 33, 465701 (2021).
  * Krauss-Kodytek  _et al._ [2021a] L. Krauss-Kodytek, W.-R. Hannes, T. Meier, C. Ruppert, and M. Betz, Nondegenerate two-photon absorption in znse: Experiment and theory, Phys. Rev. B 104, 085201 (2021a).
  * Vanderbilt [2018] D. Vanderbilt,  _Berry Phases in Electronic Structure Theory: Electric Polarization, Orbital Magnetization and Topological Insulators_ (Cambridge University Press, 2018).
  * Rioux  _et al._ [2011] J. Rioux, G. Burkard, and J. E. Sipe, Current injection by coherent one- and two-photon excitation in graphene and its bilayer, Phys. Rev. B 83, 195406 (2011).
  * McCann and Koshino [2013] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Reports on Progress in Physics 76, 056503 (2013).
  * Yang  _et al._ [2011] H. Yang, X. Feng, Q. Wang, H. Huang, W. Chen, A. T. S. Wee, and W. Ji, Giant two-photon absorption in bilayer graphene, Nano Letters 11, 2622 (2011).
  * 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, 235320 (2015).
  * Nair  _et al._ [2008] R. R. Nair, P. Blake, A. N. Grigorenko, K. S. Novoselov, T. J. Booth, T. Stauber, N. M. R. Peres, and A. K. Geim, Fine structure constant defines visual transparency of graphene, Science 320, 1308 (2008).
  * Kuzmenko  _et al._ [2008] A. B. Kuzmenko, E. van Heumen, F. Carbone, and D. van der Marel, Universal optical conductance of graphite, Phys. Rev. Lett. 100, 117401 (2008).
  * Holovský  _et al._ [2015] J. Holovský, S. Nicolay, S. De Wolf, and C. Ballif, Effect of the thin-film limit on the measurable optical properties of graphene, Scientific Reports 5, 15684 (2015).
  * Kuzmenko  _et al._ [2009] A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake, and K. S. Novoselov, Determination of the gate-tunable band gap and tight-binding parameters in bilayer graphene using infrared spectroscopy, Phys. Rev. B 80, 165406 (2009).
  * Richard  _et al._ [2004] S. Richard, F. Aniel, and G. Fishman, Energy-band structure of Ge, Si, and GaAs: A thirty-band $\mathbf{k}\mathbf{\cdot}\mathbf{p}$ method, Phys. Rev. B 70, 235204 (2004).
  * Winkler [2003] R. Winkler,  _Spin-orbit Coupling Effects in Two-Dimensional Electron and Hole Systems_ (Springer-Verlag Berlin Heidelberg, 2003).
  * Marple [1964] D. T. F. Marple, Refractive index of ZnSe, ZnTe, and CdTe, Journal of Applied Physics 35, 539 (1964).
  * Krauss-Kodytek  _et al._ [2021b] L. Krauss-Kodytek, C. Ruppert, and M. Betz, Near-infrared non-degenerate two-photon absorption coefficients of bulk gaas and si, Opt. Express 29, 34522 (2021b).
  * Furey  _et al._ [2021] B. J. Furey, R. M. Barba-Barba, R. Carriles, A. Bernal, B. S. Mendoza, and M. C. Downer, Im{$\chi^{(3)}$} spectra of 110-cut gaas, gap, and si near the two-photon absorption band edge, Journal of Applied Physics 129, 183109 (2021).
  * Note [1] The band indices $\mu\lambda$ will be skipped in this appendix.

##  Appendix A Relation between optical conductivity and absorption
coefficient

In this appendix we derive the relation between the optical conductivity of
particular order $(n)$ and the corresponding absorption coefficient, taking
care of obtaining the correct prefactors. For simplicity we do this for plane
wave (pump and probe) beams, which have a sharp carrier frequency, so we start
by setting

$\displaystyle\bm{E}_{\omega}$
$\displaystyle=\bm{E}_{\omega,p}+\bm{E}_{\omega,e}$ (97)
$\displaystyle\bm{E}_{\omega,p}$
$\displaystyle=2\pi\delta(\omega-\omega_{p})\boldsymbol{\cal
E}_{p}+2\pi\delta(\omega+\omega_{p})\boldsymbol{\cal E}_{p}^{*}$ (98)

and analogously for $\bm{E}_{\omega,e}$. The Fourier transform, Eq. (3), is
given by

$\displaystyle\bm{E}(\bm{r},t)$
$\displaystyle=\bm{E}_{p}(\bm{r},t)+\bm{E}_{e}(\bm{r},t)$ (99)
$\displaystyle\bm{E}_{p}(\bm{r},t)$ $\displaystyle=\boldsymbol{\cal
E}_{p}e^{i(\bm{k}_{p}\bm{r}-\omega_{p}t)}+\mathrm{c.c.}$ (100)

and analogously for $\bm{E}_{e}(\bm{r},t)$. Here, as for the perturbative
quantities $\rho_{\mathbf{k}}^{(n)}$ and $\bm{J}^{(n)}$ below, we explicitly
included the spatial dependence, which has been omitted in Sec. II. Note that
the real field amplitudes are $2|\boldsymbol{\cal E}_{p}|$ and
$2|\boldsymbol{\cal E}_{e}|$ for the probe and pump, respectively. The
intensities are thus

$\displaystyle I_{p}=2\epsilon_{0}n_{0p}c|\boldsymbol{\cal E}_{p}|^{2}$ (101)

and analogously for $I_{e}$. The following subsections correspond to two
alternative derivations leading to the same result.

###  A.1 $\langle\bm{J}\cdot\bm{E}\rangle$-method

$\ell$PA is determined by the $\ell$-th power of the intensity,

$\displaystyle(I_{p}+I_{e})^{\ell}=I_{p}^{\ell}+\ell
I_{p}^{\ell-1}I_{e}+\dots+\ell I_{p}I_{e}^{\ell-1}+I_{e}^{\ell}.$ (102)

For the weak-probe/strong-pump absorption considered here, only the second
last term is relevant. The rate of (pump-induced) energy absorption from the
probe beam is

$\displaystyle\langle\bm{J}^{(n)}\cdot\bm{E}_{p}\rangle$
$\displaystyle=\bm{J}^{(n)(+1|0)}\cdot\boldsymbol{\cal
E}_{p}^{*}e^{-i(\bm{k}_{p}\bm{r}-\omega_{p}t)}+$ $\displaystyle\
\quad\bm{J}^{(n)(-1|0)}\cdot\boldsymbol{\cal
E}_{p}e^{i(\bm{k}_{p}\bm{r}-\omega_{p}t)}$ (103)
$\displaystyle=2\,\mathrm{Re}\left(\bm{J}^{(n)(+1|0)}\cdot\boldsymbol{\cal
E}_{p}^{*}e^{-i(\bm{k}_{p}\bm{r}-\omega_{p}t)}\right)$ (104)
$\displaystyle=2\,\mathrm{Re}\left(J_{p}^{(n)(+1|0)}\cdot{\cal
E}_{p}^{*}e^{-i(\bm{k}_{p}\bm{r}-\omega_{p}t)}\right),$ (105)

with ${\cal E}_{p}=\boldsymbol{\cal E}_{p}\cdot\hat{\mathbf{p}}$ etc. Here,
after carrying out the spatial average $\langle\dots\rangle$ over a
macroscopic regime, where we employ the slowly varying amplitude
approximation, only the current components co- or counter-propagating with the
probe beam remain; these have the oscillatory dependence $exp[{\pm
i(\bm{k}_{p}\bm{r}-\omega_{p}t)}]$ denoted by the superscript $(\pm 1|0)$. We
obtain the absorption coefficient by normalizing by the intensities,

$\displaystyle\alpha^{(n)}_{p\\{e\\}}=\ell\alpha^{(n)}=\frac{2\,\mathrm{Re}\left(J_{p}^{(n)(+1|0)}\cdot{\cal
E}_{p}^{*}e^{-i(\bm{k}_{p}\bm{r}-\omega_{p}t)}\right)}{I_{p}I_{e}^{\ell-1}}.$
(106)

The factor $\ell$ is just the prefactor from the relevant intensity product
term above. At the same time, it is the number of relevant permutations of
field products leading to $\bm{J}^{(n)(+1|0)}$; for example, for 2PA ($n=3$),
the combinations are ${\cal E}_{p}{\cal E}_{e}{\cal E}_{e}^{*}$ and ${\cal
E}_{e}{\cal E}_{p}{\cal E}_{e}^{*}$. [6]

Now (see Eq. (11)),

$\displaystyle J_{p}^{(n)(+1|0)}(\bm{r},t)$
$\displaystyle=-\frac{e_{0}}{L^{3}}\sum_{\mathbf{k}}\mathrm{Tr}\left(v_{\mathbf{k}}^{p}\rho_{\mathbf{k}}^{(n)(+1|0)}(\bm{r},t)\right),$
(107)

and for the plane wave beams (see Eq. (4))

$\displaystyle\rho_{\mathbf{k}}^{(n)(+1|0)}(\bm{r},t)$
$\displaystyle=\left(\frac{e_{0}}{\hbar}\right)^{n}e^{i(\bm{k}_{p}\bm{r}-\omega_{p}t)}{\cal
E}_{p}|{\cal E}_{e}|^{n-1}\bar{\cal P}_{\mathbf{k}}^{(n)}$ (108)

with $\bar{\cal P}_{\mathbf{k}}^{(n)}$ as defined in Sec. II.3. Thus, Eq.
(106) becomes

$\displaystyle\alpha^{(n)}_{p\\{e\\}}=\ell\alpha^{(n)}=-2\frac{e_{0}}{L^{3}}\left(\frac{e_{0}}{\hbar}\right)^{n}\frac{\sum_{\mathbf{k}}\mathrm{Re}\left[\mathrm{Tr}\left(v_{\mathbf{k}}^{p}\bar{\cal
P}_{\mathbf{k}}^{(n)}\right)\right]}{(2\epsilon_{0}c)^{\ell}n_{0p}n_{0e}^{\ell-1}}$
(109)

which corresponds to Eq. (16). Finally using the symmetrized form of the
conductivity (13) we return to Eq. (15):

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$
$\displaystyle=\frac{2}{(2\epsilon_{0}c)^{\ell}}\frac{1}{n_{0p}n_{0e}^{\ell-1}}\mathrm{Re}\left[\sigma_{\mathrm{eff}}^{(n)}\right].$
(110)

###  A.2 Nonlinear wave equation method

This alternative derivation is based on Ref. [6] and the Supplementary
Material of Ref. [29]. The nonlinear wave equation is

$\displaystyle\left(\nabla^{2}-\frac{n_{0p}^{2}}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\right)\bm{E}_{p}(\bm{r},t)=\frac{1}{\epsilon_{0}c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\bm{P}^{(n)}(\bm{r},t).$ (111)

Using the plane-wave basis, and setting $\hat{\bm{k}}_{p}=\hat{\bm{z}}$,

$\displaystyle\left(\frac{\partial^{2}}{\partial
z^{2}}-\frac{n_{0p}^{2}}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\right){E}_{p}=\frac{1}{\epsilon_{0}c^{2}}\frac{\partial^{2}}{\partial
t^{2}}{P}_{p}^{(n)}=\frac{1}{\epsilon_{0}c^{2}}\frac{\partial}{\partial
t}{J}_{p}^{(n)}.$ (112)

Here ${J}_{p}^{(n)}={J}_{p}^{(n)(+1|0)}+{J}_{p}^{(n)(-1|0)}$ is the nonlinear
current density driving the optical response. The forward- and backward
propagating terms in the wave equation can be compared separately. Since they
are the complex conjugate of each other, we continue explicitly with the
forward-propagating part only, for which we find

$\displaystyle\left(\frac{\partial^{2}}{\partial
z^{2}}-\frac{n_{0p}^{2}}{c^{2}}\frac{\partial^{2}}{\partial t^{2}}\right){\cal
E}_{p}e^{i(k_{p}z-\omega_{p}t)}=\frac{1}{\epsilon_{0}c^{2}}\frac{\partial}{\partial
t}{J}_{p}^{(n)(+1|0)}.$ (113)

The left-hand side can thus be written

$\displaystyle\left(\frac{\partial^{2}{\cal E}_{p}}{\partial
z^{2}}+2ik_{p}\frac{\partial{\cal E}_{p}}{\partial z}-k_{p}^{2}{\cal
E}_{p}+\frac{n_{0p}^{2}\omega_{p}^{2}}{c^{2}}{\cal
E}_{p}\right)e^{i(k_{p}z-\omega_{p}t)}.$ (114)

In the slowly-varying envelope approximation the first term in the bracket is
neglected and the third and fourth terms cancel. This gives

$\displaystyle 2ik_{p}\frac{\partial{\cal E}_{p}}{\partial
z}e^{i(k_{p}z-\omega_{p}t)}=\frac{1}{\epsilon_{0}c^{2}}\frac{\partial}{\partial
t}{J}_{p}^{(n)(+1|0)},$ (115)

where we can cancel the oscillatory terms to obtain the reduced wave equation
(with ${J}_{p}^{(n)(+1|0)}={\cal J}_{p}e^{i(k_{p}z-\omega_{p}t)}$)

$\displaystyle 2ik_{p}\frac{\partial{\cal E}_{p}}{\partial
z}=-\frac{i\omega_{p}}{\epsilon_{0}c^{2}}{\cal J}_{p}$ (116)

or

$\displaystyle\frac{\partial{\cal E}_{p}}{\partial
z}=-\frac{1}{2n_{0p}\epsilon_{0}c}{\cal J}_{p}.$ (117)

Analogously to Eq. (108), the current (Eq. (12)), projected onto
$\hat{\mathbf{p}}$ and including relevant permutations, is

$\displaystyle J_{p}^{(n)(+1|0)}(z,t)$
$\displaystyle=e^{i(k_{p}z-\omega_{p}t)}{\cal E}_{p}|{\cal
E}_{e}|^{n-1}\sigma_{\mathrm{eff}}^{(n)}.$ (118)

We thus have

$\displaystyle\frac{\partial{\cal E}_{p}}{\partial
z}=-\frac{1}{2n_{0p}\epsilon_{0}c}{\cal E}_{p}|{\cal
E}_{e}|^{n-1}\sigma_{\mathrm{eff}}^{(n)}$ (119)

and correspondingly for the complex conjugate. This allows us to write the
derivative of the intensity (101) as

$\displaystyle\frac{\partial I_{p}}{\partial z}$
$\displaystyle=2n_{0p}\epsilon_{0}c\left({\cal E}_{p}\frac{\partial{\cal
E}_{p}^{*}}{\partial z}+\mathrm{c.c.}\right)$ (120) $\displaystyle=-2|{\cal
E}_{p}|^{2}|{\cal
E}_{e}|^{n-1}\mathrm{Re}\left[\sigma_{\mathrm{eff}}^{(n)}\right].$ (121)

On the other hand, the absorption coefficient is defined as

$\displaystyle\frac{\partial I_{p}}{\partial z}$
$\displaystyle=\alpha^{(n)}_{p\\{e\\}}I_{p}I_{e}^{\ell-1},$ (122)

so we get that

$\displaystyle\alpha^{(n)}_{p\\{e\\}}$
$\displaystyle=\frac{2}{(2\epsilon_{0}c)^{\ell}}\frac{1}{n_{0p}n_{0e}^{\ell-1}}\mathrm{Re}\left[\sigma_{\mathrm{eff}}^{(n)}\right],$
(123)

which corresponds to Eq. (15).

If instead of single plane waves, we have a spectrum,

$\displaystyle\bm{E}(z,t)=\int\frac{\mathrm{d}\omega}{2\pi}{\bm{E}}_{\omega}e^{i(k(\omega)z-\omega
t)}$ (124)

with $k(\omega)=n_{0}(\omega)\omega/c$, the derivation goes analogously, since
each frequency component can be treated independently, similar to the
splitting into forward/backward modes above.

##  Appendix B $\mathbf{k}_{\perp}$-integral

Here we show that111The band indices $\mu\lambda$ will be skipped in this
appendix.

$\int\mathrm{d}^{{D}-1}\mathbf{k}_{\perp}\sum_{k_{\parallel,\mathrm{res}}}\left|\frac{\partial{\omega}_{\mu\lambda\mathbf{k}}}{\partial
k_{\parallel}}\right|^{-1}\dots=\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}\left\|\nabla_{\mathbf{k}}{\omega}_{\mu\lambda\mathbf{k}}\right\|^{-1}\dots$
(125)

We first give a proof for ${D}=2$. Denoting the parallel and perpendicular
axes by $y$ and $x$, respectively, we can write for one specific resonance
path (there are two of them if the isoline is convex):

$\displaystyle\int\mathrm{d}k_{x}=\int\mathrm{d}s\frac{\mathrm{d}k_{x}}{\mathrm{d}s}=\int\mathrm{d}s\
\hat{t}_{\mathbf{k},x}=\int\mathrm{d}s\
\frac{\left|\frac{\partial{\omega}_{\mathbf{k}}}{\partial
k_{y}}\right|}{\left\|\nabla_{\mathbf{k}}^{xy}{\omega}_{\mathbf{k}}\right\|},$
(126)

where $s$ is the arc length along the path and $\hat{\bm{t}}_{\mathbf{k}}$ is
the tangential unit vector on the path. In the last equation,
$\hat{\bm{t}}_{\mathbf{k}}$ has been expressed by rotating the normal unit
vector $\hat{\bm{n}}_{\mathbf{k}}$ by $\pi/2$,

$\displaystyle\hat{\bm{t}}_{\mathbf{k}}=\pm\left(\begin{array}[]{*{2}{c}}0&1\\\
-1&0\end{array}\right)\hat{\bm{n}}_{\mathbf{k}},\quad\hat{\bm{n}}_{\mathbf{k}}=\frac{\nabla_{\mathbf{k}}^{xy}{\omega}_{\mathbf{k}}}{\left\|\nabla_{\mathbf{k}}^{xy}{\omega}_{\mathbf{k}}\right\|},\quad\nabla_{\mathbf{k}}^{xy}\equiv\left(\begin{array}[]{*{1}{c}}\partial/\partial{k_{x}}\\\
\partial/\partial{k_{y}}\end{array}\right)$ (131)

where $\pm=\mathrm{sgn}\left(\partial{\omega}_{\mathbf{k}}/\partial
k_{y}\right)$ in order to keep the direction of the integration
($\hat{t}_{\mathbf{k},x}>0$). In the last expression of (126) the direction
doesn’t matter any more, and we can also put together the two resonance paths
to form a closed contour. Then we directly arrive at Eq. (125), where for
${D}=2$ the gradient is understood as
$\nabla_{\mathbf{k}}=\nabla_{\mathbf{k}}^{xy}$.

Now we show that Eq. (125) also holds for ${D}=3$, where it reads (choosing
$\hat{\mathbf{k}}_{\parallel}=\hat{\bm{z}}$)

$\displaystyle\int\mathrm{d}k_{x}\mathrm{d}k_{y}\sum_{k_{z,\mathrm{res}}}\left|\frac{\partial{\omega}_{\mu\lambda\mathbf{k}}}{\partial
k_{\parallel}}\right|^{-1}\dots$
$\displaystyle=\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}\left\|\nabla_{\mathbf{k}}{\omega}_{\mu\lambda\mathbf{k}}\right\|^{-1}\dots$
(132)

We start from the right-hand side, where we first divide the closed surface
into ‘layers’ (in the $z$-direction) of individual surfaces
$S_{\mathrm{res}}$, so that each $S_{\mathrm{res}}$ is describable by a
function $k_{z,\mathrm{res}}(k_{x},k_{y})$ in a (simply connected) subspace of
the $k_{x},k_{y}$-plane. Then we express the surface integral as a plane 2D
integral using the appropriate surface element,

$\displaystyle\oint\mathrm{d}\mathbf{k}_{\mathrm{res}}$
$\displaystyle=\sum_{k_{z,\mathrm{res}}}\int\\!\\!\\!\\!\int\mathrm{d}S_{\mathrm{res}}$
(133)
$\displaystyle=\sum_{k_{z,\mathrm{res}}}\int\\!\\!\\!\\!\int\mathrm{d}k_{x}\mathrm{d}k_{y}\left\|\frac{\partial\mathbf{k}_{\mathrm{res}}}{\partial
k_{x}}\times\frac{\partial\mathbf{k}_{\mathrm{res}}}{\partial k_{y}}\right\|,$
(134)

with
$\mathbf{k}_{\mathrm{res}}^{\mathsf{T}}=(k_{x},k_{y},k_{z,\mathrm{res}}(k_{x},k_{y}))$.
Now we use that

$\displaystyle\frac{\partial k_{z,\mathrm{res}}}{\partial k_{x}}$
$\displaystyle=\frac{\hat{t}^{xz}_{z}}{\hat{t}^{xz}_{x}}=\frac{\partial_{k_{x}}{\omega}_{\mathbf{k}}}{\partial_{k_{z}}{\omega}_{\mathbf{k}}},\quad\frac{\partial
k_{z,\mathrm{res}}}{\partial
k_{y}}=\frac{\hat{t}^{yz}_{z}}{\hat{t}^{yz}_{y}}=\frac{\partial_{k_{y}}{\omega}_{\mathbf{k}}}{\partial_{k_{z}}{\omega}_{\mathbf{k}}},$
(135)

where, for instance, $\hat{\bm{t}}^{xz}$ is the tangent vector to the
resonance contour in the $xz$-plane, which is again expressed according to Eq.
(131) (the sign $\pm$ now doesn’t matter from the beginning). Therefore

$\displaystyle\left\|\frac{\partial\mathbf{k}_{\mathrm{res}}}{\partial
k_{x}}\times\frac{\partial\mathbf{k}_{\mathrm{res}}}{\partial
k_{y}}\right\|=\left\|\left(\begin{array}[]{*{1}{c}}1\\\ 0\\\
\hat{t}^{xz}_{z}/\hat{t}^{xz}_{x}\end{array}\right)\times\left(\begin{array}[]{*{1}{c}}0\\\
1\\\ \hat{t}^{yz}_{z}/\hat{t}^{yz}_{y}\end{array}\right)\right\|$ (142)
$\displaystyle=\left\|\left(\begin{array}[]{*{1}{c}}-\hat{t}^{xz}_{z}/\hat{t}^{xz}_{x}\\\
-\hat{t}^{yz}_{z}/\hat{t}^{yz}_{y}\\\
1\end{array}\right)\right\|=\frac{\left\|\nabla_{\mathbf{k}}{\omega}_{\mathbf{k}}\right\|}{|\partial_{k_{z}}{\omega}_{\mathbf{k}}|}$
(146)

Applying these transformations to the right-hand side of Eq. (132) completes
the proof for ${D}=3$.