# Saturation of interband absorption in graphene

F.T. Vasko<EMAIL_ADDRESS>Institute of Semiconductor Physics, NAS of
Ukraine, Pr. Nauky 41, Kiev, 03028, Ukraine

###### Abstract

The transient response of an intrinsic graphene, which is caused by the
ultrafast interband transitions, is studied theoretically for the range of
pumping correspondent to the saturated absorption regime. Spectral and
temporal dependencies of the photoexcited concentration as well as the
transmission and relitive absotption coefficients are considered for mid-IR
and visible (or near-IR) spectral regions at different durations of pulse and
broadening energies. The characteristic intencities of saturation are
calculated and the results are compared with the experimental data measured
for the near-IR lasers with a saturable absorber. The negative absorption of a
probe radiation during cascade emission of optical phonons is obtained.

###### pacs:

78.47.jb, 78.67.Wj, 42.65.-k

##  I Introduction

The character of nonlinear response under the ultrafast interband excitation
of an intrinsic graphene is determined by a several physical processes which
are dependent on conditions of pumping. Under a low excitation level, when the
one-photon transitions take place, the energy relaxation and recombination of
photoexcited carriers were studied with the use of the time-resolved pump-
probe measurements, see experimental data and theoretical discussions in Refs.
1 and 2, respectively. Under an extremely high pumping, the multi-quantum
transitions, which cause the harmonics generation and the hybridization of
electron-hole states, take place. This regime is not investigated completely
for graphene, see general consideration for bulk materials or quantum wells in
Sects. 10 and 56 of Ref. 3. The nonlinear response is also possible within the
single-photon approach because the Rabi oscillations of coherent response
should take place in graphene if a pulse duration is less 100 fs. 4  Beside of
this, the nonlinear regime of energy relaxation and recombination due to the
Pauli blocking effect takes place if the photogeneration rate is comparable to
the energy relaxation or recombination rates. A saturation of transient
absorption, which have been investigated recently, 5 ; 6  is the most
important manifistation of such a regime. It is because this phenomenon was
exploited for realization of an ultrafast laser with a graphene saturable
absorber in the telecommunication spectral region. To the best of our
knowledge, a complete theoretical treatment of the saturation mechanism is not
performed yet and an investigation of this phenomena is timely now.

In this paper, we consider the temporal nonlinear evolution of carriers under
photoexcitation by the ultrafast pulses in mid-IR and visible (or near-IR)
spectral regions. For the mid-IR pump, one can neglect the relaxation and
recombination processes (the quasielastic relaxation due to acoustic phonons
remains uneffective up to nanoseconds) and the transient distribution of
nonequilibrium carriers is determined by the broadening of interband
transitions due to an elastic scattering and by the parameters of excitation.
For the excitation energies, $\hbar\Omega$, above the optical phonon energy,
$\hbar\omega_{\eta}$ ($\eta=\Gamma,~{}K$ labels the phonon modes correspondent
to the intra- and intervalley transitions), when a cascade emission of the
optical phonons dominates in relaxation, the transient distribution transforms
into a set of peaks. The effective electron-hole recombination takes place if
the lowest peak is placed around the half-energy of optical phonon,
$\hbar\omega_{\eta}/2$. So that, the character of response modifies
essentially if the frequency $\Omega$ varies over the $\omega_{\eta}$-range.

The saturation process is described within the framework of the temporally
local approach, when the decoherentization time (which determines the
broadening of the interband transitions) is shorter in comparision with the
duration of pumping. Spectral and temporal dependencies of the photoexcited
concentration and the response on the probe radiation (transmission and
absorption coefficients) are presented. The thresholds for saturation of
response are estimated to be about 0.2 MW/cm2, 60 MW/cm2, and 0.6 GW/cm2 for
$\hbar\Omega\sim$0.12, 0.8, and 1.5 eV, respectively (mid-IR, near-IR, and
visible spectral regions). These results are dependent on the
decoherentization and relaxation mechanisms and the electrodynamics
conditions. Their are discussed in comparison with the experimental data for
the near-IR pumping case 5 ; 6 . Conditions for the transient negative
absorption of a probe radiation during cascade emission of optical phonons are
also analyzed (this phenomenon under a steady-state pumping was considered
recently 7  in connection with a possibility of the THz lasing effect).

The consideration below is organized as follows. The temporally local approach
for description of the response of photoexcited carriers is developed in Sec.
II. Spectral and temporal dependencies of the response are described in Sects.
III and IV for the cases of excitation in mid-IR and visible (or near-IR)
spectral regions, respectively. A discussion of experimental data, the list of
assumptions used, and concluding remarks are given in the last section. In
Appendix we consider the mechanism of saturation caused by the collisionless
Rabi oscillations.

##  II Temporally local approach

Since the symmetry of the energy spectrum and scattering processes for
electrons and holes in an intrinsic graphene, we describe the phenomena under
consideration by the same distribution functions for the both types of
carriers, $f_{pt}$. According to Refs. 2b and 4, the kinetic equation for
$f_{pt}$ takes form:

$\frac{df_{pt}}{dt}=\nu_{pt}(1-2f_{pt})+J(f_{t}|p)$ (1)

and it should be solved with the initial condition $f_{pt\to\infty}=0$. Here
$J(f_{t}|p)$ is the collision integral, which is described the relaxation and
recombination processes, and $\nu_{pt}$ is the interband generation rate due
to the in-plane electric field ${\bf E}w_{t}\exp(-i\Omega t)+c.c.$, where
${\bf E}$ is the field strength, $\Omega$ is the pumping frequency, and
$w_{t}$ is the envelope form-factor of pulse with duration 2$\tau_{p}$
centered at $t=0$. Supposing that $\tau_{p}$ exceeds the dephasing time, we
have used in Eq. (1) the temporally local approach with the separated filling
factor, $(1-2f_{pt})$, and with the rate of photoexcitaion:

$\nu_{pt}=\nu_{R}w_{t}^{2}\Delta\left(\frac{2\upsilon
p-\hbar\Omega}{\gamma}\right),~{}~{}~{}~{}\nu_{R}=\frac{\pi(eE\upsilon/\Omega)^{2}}{\hbar\gamma}.$
(2)

Here $\upsilon=10^{8}$ cm/s is the velocity of neutrinolike quasiparticles,
$2\gamma$ is the broadening of the interband excitation described by the
phenomenological factor $\Delta(z)$. Below we consider the Lorentzian
lineshape of photoexcitation, when $\Delta(z)=\left[\pi(1+z^{2})\right]^{-1}$
and the Gaussian temporal envelope
$w_{\tau}=\sqrt[4]{2/\pi}\exp[-(t/\tau_{p})^{2}]$.

The solution of Eq. (1) determines both the photoinduced concentration, which
given by the standard formula

$n_{t}=\frac{2}{\pi\hbar^{2}}\int\limits_{0}^{\infty}{dpp}f_{pt},$ (3)

and the transient response on a probe radiation of frequency $\omega$
($\propto\exp(-i\omega t)$, which is described by the dynamic conductivity
$\sigma_{\omega t}$. For the collisionless case $\hbar\omega/\gamma\gg 1$,
when the parametric dependency on time takes place, 9  the real part of
$\sigma_{\omega t}$ is written as follows:

${\rm Re}\sigma_{\omega t}=\frac{e^{2}}{4\hbar}(1-2f_{p_{\omega},t}),$ (4)

where $p_{\omega}=\hbar\omega/\upsilon$. The imaginary part of $\sigma_{\omega
t}$ is determined through ${\rm Re}\sigma_{\omega t}$ with the use of the
dispersion relation and one can check that the carrier-induced contribution to
${\rm Im}\sigma_{\omega t}$ appears to be weak in comparison with (4) for the
peak-like distributions of carriers considered below. Thus, the only filling
factor in ${\rm Re}\sigma_{\omega t}$ is responsible for the nonlinear
behavior of the response under consideration.

We restrict ourselves by the the geometry of normal propagation of radiation.
The relative absorption of graphene sheet, $\xi_{\omega t}$, as well as the
reflection and transmission coefficients, $R_{\omega t}$ and $T_{\omega t}$,
are determined through $\sigma_{\omega t}$. Since the energy conservation
requirement, 8

$R_{\omega t}+T_{\omega t}+\xi_{\omega t}=1,$ (5)

we consider below only the absorption and transmission coefficients:

$\displaystyle\xi_{\omega t}\simeq\frac{16\pi}{\sqrt{\varepsilon}c}\frac{{\rm
Re}\sigma_{\omega t}}{|1+\sqrt{\epsilon}+4\pi\sigma_{\omega
t}/c|^{2}}\approx\xi_{m}(1-2f_{p_{\omega}t}),$ $\displaystyle T_{\omega
t}\simeq\frac{4\sqrt{\epsilon}}{\left|1+\sqrt{\epsilon}+4\pi\sigma_{\omega
t}/c\right|^{2}}\approx\frac{T_{m}}{(1-af_{p_{\omega}t})^{2}}.$ (6)

Here $\sqrt{\epsilon}$ is the refraction index of a thick substrate (for SiO2
substrate $\sqrt{\epsilon}\simeq$1.45 and dispersion of $\epsilon$ can be
neglected) and we approximately separated the carrier contributions using the
coefficients $\xi_{m}\approx
4\pi\alpha/\left[\sqrt{\epsilon}(1+\sqrt{\epsilon})\right]$, $T_{m}\approx
4\sqrt{\epsilon}/(1+\sqrt{\epsilon}+\pi\alpha)^{2}$, and $a\approx
2\pi\alpha/(1+\sqrt{\epsilon}+\pi\alpha)$ with $\alpha=e^{2}/\hbar c$. Notice,
that the negative absorption regime $\xi_{\omega t}<0$ takes place if
$f_{p_{\omega},t}>1/2$, under the population inversion condition (see
discussion in Sec. IV).

At $\omega=\Omega$ these relations describe the propagation of pumping pulse
with the time-dependent intensity $Sw_{t}^{2}$, where $S$ is the maximal
intensity. Performing the averaging of (6) over the pulse duration one obtains

$\left|\begin{array}[]{*{20}c}\xi_{S}\\\
T_{S}\end{array}\right|=\int\limits_{-\infty}^{\infty}\frac{dt}{\tau_{p}}w_{t}^{2}\left|\begin{array}[]{*{20}c}\xi_{\Omega
t}\\\ T_{\Omega t}\end{array}\right|,$ (7)

where we used $\int\limits_{-\infty}^{\infty}dtw_{t}^{2}/\tau_{p}=1$. Below we
solve Eq. (1) and analyze the responses (6) and (7) for different parameters
of pump and probe radiations.

##  III Mid-IR excitation

We consider here the mid-IR pumping case when the energy relaxation of
carriers is ineffective and $J(f_{t}|p)$ in Eq. (1) can be neglected. As a
result, the solution of the problem (1) takes form:

$f_{pt}=\int\limits_{-\infty}^{t}dt^{\prime}\nu_{pt^{\prime}}\exp\left(-2\int\limits_{t^{\prime}}^{t}d\tau\nu_{p\tau}\right).$
(8)

Evolution of such a distribution from zero value at $t\ll-\tau_{p}$ to the
saturated peak with the maximal value $f_{max}=1/2$ is shown in Fig. 1a versus
dimensionless time and energy at the pumping intensity $S=$1 MW/cm2. Temporal
dependencies of $f_{p_{\Omega}t}$ at different $S$ are shown in Fig. 1b. These
calculations were performed for $\hbar\Omega\simeq$120 meV (pumping by
CO2-laser), the pulse duration $2\tau_{p}\simeq$1 ps, and the broadening
energy $\gamma\simeq$6 meV which is in agreement with the mobility data for
the case of elastic scattering. 9  The temporally-dependent photoinduced
concentration $n_{t}$ is plotted Fig. 1c for the same parameters. The
saturated concentration versus intensity, which is attained at $t>\tau_{p}$,
is presented for $\gamma=$6 and 12 meV in Fig. 1d. These dependencies can be
fitted as

$n_{S}\approx\frac{bS}{1+S/S_{n}},$ (9)

where $b\simeq$6 or 12.2 MW-1 [$n_{S}$ is measured in 1011 cm-2] and
$S_{n}\simeq$1.76 or 10 MW/cm2 for and $\gamma=$6 or 12 meV, respectively.

Figure 1: (Color online) (a) Photoexcited distribution $f_{pt}$ versus energy
$\upsilon p$ and dimensionless time, $t/\tau_{p}$ at mid-IR pumping level
$S=$1 MW/cm2. (b) Temporal evolution of $f_{\Omega t}\equiv f_{p_{\Omega}t}$
at $S=$0.1, 0.3, 1, 3, and 6 MW/cm2 (marked as 1-5). (c) Potoinduced
concentration versus $t/\tau_{p}$ for the same conditions as in panel (b). (d)
Concentration $n_{S}$ at $t/\tau_{p}\to\infty$ versus $S$ for the different
broadening energies $\gamma$. Dotted curves are correspondent to the fit (9).

The relative absorption and transmission coefficients of a probe radiation of
frequency $\omega$ are determined through $f_{p_{\omega}t}$ according to Eqs.
(6). Spectral and temporal dependencies of $\xi_{\omega t}$ are shown in Fig.
2a for the conditions used in Fig. 1a. Since $af_{p_{\omega}t}\ll 1$, the peak
of relative transmission $T_{\omega t}/T_{m}$ resembles $f_{p_{\Omega}t}$
presented in Fig. 1a. Here $T_{m}=\simeq$0.95 is the transmission coefficient
without for non-doped graphene. The temporally-dependent relative absorption
and transmission at the pumping frequency $\Omega$ and at different $S$ are
presented in Figs. 2b and 2c, respectively. The saturated values of
$\xi_{S}/\xi_{m}$ and $T_{S}$ versus intensity are plotted in the upper and
lower panels of Fig. 2d for the parameters used in Fig. 1d ($\xi_{S}$ and
$T_{S}$ have only a weak dependency on $\gamma$). These curves can be fitted
as

$\xi_{S}\approx\frac{\xi_{m}}{1+S/\overline{S}},~{}~{}~{}~{}T_{S}\approx
T_{m}+\frac{hS}{1+S/\overline{S}},$ (10)

where $\overline{S}\approx$0.2 MW/cm2 and $h\approx$0.09 cm2/MW. The
saturation of $\xi_{S}$ and $T_{S}$ takes place at lower threshold intensities
in comparision to $n_{S}$, c. f. Figs. 1c and 2d, 2e. Thus, for the pumping
range $\geq$1 MW/cm2 the one-photon absorption is suppressed and a damage of
graphene by mid-IR radiation with $\tau_{p}\lesssim$1 ps is not possible.

Figure 2: (Color online) (a) Spectral and temporal dependencies of relative
absorption, $\xi_{\omega t}$ at $S=$1 MW/cm2. (b) Temporal evolution of
$\xi_{\Omega t}$ at $S$ used in Fig. 1b (marked). (c) The same as in panel (b)
for transmission, $T_{\Omega t}$. (d) Avaraged over pulse absorption and
transmission coefficients (upper and lower panels, respectively) versus $S$
[dashed curves are correspondent to Eq. (10)].

##  IV Cascade emission effect

In this section we consider the photoexcitation by visible and near-IR
radiation, when the cascade emission of optical phonons should be taken into
account in Eq. (1). For the temperatures below the optical phonon energies,
the spontaneous emission processes are only essential and the collision
integral is given by the finite-difference form (see evaluation in Refs. 2b
and 10)

$\displaystyle
J\left(f_{t}|p\right)=\sum_{\eta}\left[\nu_{p+p_{\eta}}\left(1-f_{pt}\right)f_{p+p_{\eta}t}\right.$
(11)
$\displaystyle\left.-\nu_{p-p_{\eta}}\left(1-f_{p-p_{\eta}t}\right)f_{pt}-\widetilde{\nu}_{p_{\eta}-p}f_{p_{\eta}-pt}f_{pt}\right].$

Here $\eta=\Gamma,~{}K$ is correspondent to the intra- and intervalley
transitions with the energy transfer, $\hbar\omega_{\eta}$, and the momentum
transfer, $p_{\eta}=\hbar\omega_{\eta}/\upsilon$. The last contribution of Eq.
(11) is responsible for the recombination process while the first and second
terms describe the interband cascade relaxation of carriers. The relaxation
rates $\nu_{p}$ and $\widetilde{\nu}_{p}$ are proportional to the density of
states,
$\nu_{p}\approx\widetilde{\nu}_{p}\approx\theta(p)\upsilon_{\eta}p/\hbar$,
where the characteristic velocities $\upsilon_{\Gamma,K}$ can be estimated
crudely as $\upsilon_{K}\approx 2\times 10^{6}$ cm/s and
$\upsilon_{\Gamma}\approx 10^{6}$ cm/s. 10  Thus, the $K$-mode emission gives
a dominant contribution to the relaxation process; moreover, the only
interband recombination is possible in the passive region, $0<\upsilon
p<\hbar\omega_{K}=$170 meV. Below we neglect other relaxation processes, so
that a peak-like transient distribution of carriers takes place due to the
negligible phonon dispersion and a narrow distribution of photoexcited
carriers, under the condition $\gamma\ll\hbar\omega_{K}$. For the sake of
simplicity, the cases of effective or suppressed recombination, when the lower
peak in the passive region is placed around or outside the energy
$\hbar\omega_{K}$ are considered. It is convenient to analyze calculations for
the near-IR and visible pumping cases separately.

Figure 3: (Color online) (a) Contour plots of photoexcited distributions
$f_{pt}$ versus energy and time for pulse duration $2\tau_{p}=$0.6 ps at
pumping level 200 MW/cm2 and $\hbar\Omega=$850 meV. (b) Transient evolution of
concentration $n_{t}$ and and populations of peaks around $\sim$43, 213, and
383 meV (marked as $n_{1}$, $n_{2}$, and $n_{3}$, respectively), for the same
conditions as in panel (a). (c) Evolution of $n_{t}$ for pumping levels
$S=$100, 200, and 300 MW/cm2 (marked as 1, 2, and 3). Solid and dashed curves
are plotted for $\hbar\Omega=$ 850 meV and 765 meV.

###  IV.1 Near-IR pumping

First, we consider the three-step cascade processes under the near-IR pumping
with wavelengths around $\sim 1.5~{}\mu$m and the pulse duration determined by
$\tau_{p}=$ 0.3 ps. We consider the regimes of the enhanced or suppressed
recombination supposing $\hbar\Omega=$ 850 meV or 765 meV. For this energy
region, the broadening of photoexcited peak is taken as $\gamma\simeq$18 meV,
so that $\hbar/\gamma\ll\tau_{p}$. The numerical solution of Eq. (1) with the
collision integral (11) is performed with the use of the temporal iterations
11  at different $S$. Figure 2a shows the contour plot of the three-peak
distribution function $f_{pt}$ for the case of efficient recombination
($\hbar\Omega=$ 850 meV) at $S=$200 MW/cm2. The carrier concentrations over
the peaks 1-3 and the total consentration $n_{t}$ given by Eq. (3) are shown
in Fig. 3b for the same parameters as in Fig. 1a. Since the relaxation rate in
Eq. (11) is proportional to the density of states, $\nu_{p}\propto p$, the
bottleneck effect takes place under the transition between the second and
third peaks and $n_{2}$ exceeds $n_{1,3}$. The transient evolution of
concentration for different $S$ is shown in Fig. 3c where the maximal
concentration exceeds 1012 cm-2 at $t\sim\tau_{p}$ and $S\geq$0.3 GW/cm2.
During the further evolution, $n_{t}$ decays due to the recombination process.
The case of the suppressed recombination ($\hbar\Omega=$765 meV) is different
because of, first, the peaks are shifted below (about 43 meV) and, second, the
decreasing of $f_{pt}$ and $n_{t}$ due to recombination is absent. The
saturated concentrations (dashed curves in Fig. 3c) exceed the peak
concentrations (solid curves in Fig. 3c) by factor $\sim$1.3. Note, that
$n_{t}$ decreases with increasing of $\gamma$ at fixed $S$ (not shown in Fig.
3).

Figure 4: (Color online) (a) Transient evolution of relative absorption for
$\hbar\Omega=$ 850 and 765 meV (solid and dotted curves, respectively) at
$S=$50, 100 and 200 MW/cm2 (marked as 1-3). (b) The same as in panel (a) for
transmission coefficient. (c) Avaraged over pulse absorption versus $S$
[dashed curve is correspondent to Eq. (9)]. (d) The same as in panel (c) for
transmission coefficient.

Transient evolutions of the absorption and transmission coefficients given by
Eq. (6) are shown in Figs. 4a and 4b at the different pumping frequencies
$\Omega$ (solid and dashed curves) and at different $S$. For
$t/\tau_{p}\leq$0.5, the temporal evolution of $\xi_{\Omega t}$ and $T_{\Omega
t}$ do not dependent on the character of recombination. For
$t/\tau_{p}\lesssim$1.5 this evolution is completely different: a quenching of
photoresponse or a steady-state contribution take place for the effective or
suppressed recombination cases. At $S\geq$300 MW/cm2 and $t/\tau_{p}\sim$0 one
obtains the saturated absorption around $\xi_{\Omega t}\sim$0.1. The negative
absorption takes place for a probe radiation with $\hbar\omega$ around the
first and second peaks. It is because $f_{p_{\omega}t}>1/2$, see Eq. (4) and
the contour plot in Fig. 3, where the regions of negative absorption are
separated by the thick (red) curves. Thus, the negative absorption (and a
possible stimulated emission of mid-IR radiation) is realized at $S\geq$100
MW/cm2 during time intervals $t\lesssim 5\tau_{p}$.

The absorption and transition coefficients averaged over pulse duration
according to Eq. (7) are shown in Figs. 4c and 4d. Since the transient
response at $|t|\lesssim\tau_{p}$ does not depend on the recombination
mechanism (see Figs. 3c, 4a, and 4b), the variation of $\xi_{S}$ and $T_{S}$
with $\hbar\Omega$ is less than 5%. These dependencies can be fitted by Eq.
(10) with the characteristic intensity $\overline{S}\approx$60 MW/cm2 and the
coefficient $h\approx$0.3 cm2/GW.

Figure 5: (Color online) (a) Transient evolution of concentration $n_{t}$ and
and populations of peaks around $\sim$0.09, 0.26, 0.43, 0.6, and 0.77 eV
(marked as $n_{1-5}$, respectively) for pulse duration $2\tau_{p}=$0.4 ps at
pumping level 0.4 GW/cm2 and $\hbar\Omega=$1.53 eV. (b) The same as in panel
(a) at $\hbar\Omega=$1.615 eV for peak’s positions $\sim$0.13, 0.3, 0.47,
0.64, and 0.81 eV marked as $n_{1-5}$. (c) Evolution of $n_{t}$ for pumping
levels $S=$0.2, 0.3, 0.4, 0.6 and 0.8 GW/cm2 (marked as 1-5, respectively) for
$\hbar\Omega=$1.53 eV. (d) The same as in panel (c) for $\hbar\Omega=$1.615
eV.

###  IV.2 Visible pumping

Next, we consider the visible light pumping, with wavelengths around $\sim
0.75~{}\mu$m, using the pulse duration $2\tau_{p}=$0.4 ps and the broadening
$\gamma\approx$34 meV (so that $\hbar/\gamma\ll\tau_{p}$). Supposing
$\hbar\Omega=$ 1.53 and 1.615 eV for the enhanced and suppressed recombination
regimes one arrive to the distribution function formed during the five-step
cascade process. Transient evolutions of the concentrations over the peaks 1-5
and of the total concentration $n_{t}$ at $S=$0.4 GW/cm2 are shown in Figs. 5a
and 5b for the cases of enhanced and suppressed recombination, respectively.
Similarly to the near-IR pumping case, the upper peak concentrations decrease
fast at $t>\tau_{p}$ and a maximal population of the second peak takes place
due to the bottleneck effect. Once again, at $t<\tau_{p}$ the shapes of
$n_{t}$ are the same for the both cases. At $t>\tau_{p}$ a quenching of
$n_{t}$ due to recombination takes place in Fig. 5a while there is no a
decreasing of $n_{t}$ in Fig. 5b. The temporal dependencies of concentration
for different $S$ are shown in Figs. 5c and 5d for the two recombination
regimes under consideration. The maximal concentration range up to 1013 cm-2
at $t\sim\tau_{p}$ and $S\approx$1 GW/cm2.

Figure 6: (Color online) (a) Transient evolution of relative absorption for
$\hbar\Omega=$1.53 and 1.615 eV, (solid and dotted curves respectively) at
$S=$0.2, 0.4, 0.6, 0.8 and 1.2 GW/cm2 (marked as 1-5). (b) The same as in
panel (a) for transmission coefficient. (c) Avaraged over pulse absorption
versus $S$ [dashed curve is correspondent to Eq. (9)]. (d) The same as in
panel (c) for transmission coefficient.

The temporal evolution of $\xi_{\Omega t}$ and $T_{\Omega t}$ at frequency
$\Omega$ (solid and dotted curves are correspondent to the two recombination
cases under consideration) are plotted in Figs. 5a and 5b. By analogy with
Sect. IVA, the negative absorption regime takes place at $S\geq$0.3 GW/cm2.
Beside of this, the conditions $\xi_{\omega t}<0$ take place around the peak
positions at $\omega<\Omega$ (not plotted, see a similar behavior in Fig. 3a);
for the first and second peaks the negative absorption regime is realized up
to $t\sim 5\tau_{p}$ at $S\geq$0.1 GW/cm2. In addition, at $t<0.5\tau_{p}$ the
response does not dependent on recombination and at $t>1.5\tau_{p}$ a damping
or time-independent response is realized for the effective or suppressed
recombination.

The averaged according to Eq. (7) absorption and transmission coefficients,
which do not depend on the recombination mechanism, are plotted in Figs. 6c
and 6d. Once again, $\xi_{S}$ and $T_{S}$ can be fitted by Eqs. (10) with the
characteristic intensity $\overline{S}\approx$0.56 GW/cm2 and the coefficient
$h\approx$0.03 cm2/GW. Since the departure rate from the photoexcited peak
increases if $\hbar\Omega$ grows, the characteristic intensity $\overline{S}$
is also increased in the visible spectral region in comparison with the near-
IR pumping case.

##  V Discussion and conclusions

To summarize, we have developed the nonlinear theory of transient response of
an intrinsic graphene under the ultrafast interband excitation. Within the
local time approach, the conditions of saturation of absorption were found in
the mid-IR, near-IR, and visible spectral regions. In addition, we have
demonstrated a possibility for the stimulated mid-IR radiation due to the
bottleneck effect during the cascade emission of optical phonons. Our
consideration is based on the set of assumptions about relaxation mechanisms.
First of all, the phenomenological model for the broadening with the
characteristic energy $\gamma$ is used for description of the intersubband
transitions. In Sects. III and IV we estimated $\gamma$ from the experimental
data for the departure relaxation rates. 1 ; 2 ; 9  Secondary, a simplified
description of energy relaxation is employed. We neglect the Coulomb
scattering which is not a dominant relaxation channel at $t\lesssim\tau_{p}$,
so that the results for $\xi_{S}$ and $T_{S}$ should not be modified
essentially. But a transient distribution at $t\gg\tau_{p}$ and a condition
for the negative absorption of a probe radiation in the mid-IR region can be
modified. Also, a possible contribution of the substrate vibration 12  is not
taken into account. These points require a special consideration but, anyway,
our calculation gives a lower bound of $S$. The other assumptions (parameters
for the electron-phonon coupling, conditions for the temporally-local
approach, and description of the interband response) are rather standard for
the calculations of the optical properties and the relaxation phenomena in
graphene. In addition, an inhomogenity of pumping, which causes the lateral
diffusion of carriers, 13  and a heating of phonons 14  may be essential;
these phenomena requre a special treatment, both experimental and theoretical.

We turn now to discussion of the experimental data available for the near-IR
spectral region. 5 ; 6  Numerical estimates for the saturation thresholds and
for the concentrations of the photoexcited carriers are in a qualitative
agreement with the consideration performed. But an accurate comparison with
the results presented is not possible for the two reasons. First, the graphene
structure was embedded into the laser cavity in 5 ; 6  so that the
electrodynamical conditions (for a propagated, reflected, and absorbed
radiation) were different from the simple geometry considered here. Second,
the multi-layer graphene or the graphene flakes were used, while a single-
layer graphene case was not under a detailed treatment. Thus, a special
measurements with the use of the simplest geometry of a well-characterized
sample placed over a semi-infinite substrate are necessary.

In closing, we have analyzed theoretically the conditions for realization of
an efficient graphene-based saturable absorber and have performed a comparison
with the experimental data. More extended treatment of this phenomena under
near-IR pumping, including an above-mentioned special measurements, in order
to improve an efficiency of the graphene based saturable absorber in the
lasers for telecommunications. An additional study in the mid-IR and visible
spectral regions should be useful for verification of different relaxation
mechanisms.

The author would like to thank E. I. Karp for insightful comments.

*

##  Appendix A Rabi oscillations regime

Below we describe the saturation of the averaged absorption and transmission
coefficients (7) under an ultrafast pumping for the case when the Rabi
oscillations conditions are satisfied. 4  The collisionless regime of response
is described by the $S$-dependent contribution to the distribution function

$1-2f_{pt}=\cos\left(\sqrt{\frac{S}{S_{R}}}\int\limits_{-\infty}^{t}\frac{dt^{\prime}}{\tau_{p}}w_{t^{\prime}}\right).$
(12)

Here the characteristic intensity is given by

$S_{R}=\frac{\sqrt{\epsilon}c}{4\pi}\left(\frac{\hbar\Omega}{e\tau_{p}\upsilon}\right)^{2}$
(13)

and $S_{R}\simeq$0.6 MW/cm2 for CO2 pumping with $\tau_{p}\simeq$0.1 ps. For
the near-IR or visible pumping with $\tau_{p}\simeq$30 fs, one obtains
$S_{R}\simeq$0.3 or 1 GW/cm2. Notice, that $S_{R}\propto(\Omega/\tau_{p})^{2}$
and (A.1) is not dependent on any other parameter if $\tau_{p}$ is shorter
than the dephasing relaxation time.

Figure 7: (Color online) Normalized absorption coefficient given by (A.3)
versus dimensionless intensity $S/S_{R}$. Dotted curve presents a monotonic
fit.

Substituting the distribution (A.1) into Eqs. (6) and (7) one obtains the
following analytical expression for the relative absorption

$\frac{\xi_{S}}{\xi_{m}}=\int\limits_{-\infty}^{\infty}\frac{dt}{\tau_{p}}w_{t}^{2}\cos\left(\sqrt{\frac{S}{S_{R}}}\int\limits_{-\infty}^{t}\frac{dt^{\prime}}{\tau_{p}}w_{t^{\prime}}\right)$
(14)

while the transmission coefficient is given by $T_{S}\approx
T_{m}(1-\widetilde{a}\xi_{S})$ with
$\widetilde{a}=\sqrt{\epsilon}(1+\sqrt{\epsilon})/2$. In Fig. 7 we plot the
function $\xi_{S}/\xi_{m}$ versus dimensionless intensity $S/S_{R}$ and the
oscillating character of response at $S/S_{R}\geq$4\. The oscillations appears
due to the dynamic inversion of transient population, see Ref. 4. The fit of
(A.3) at $S/S_{R}\leq 1$ is given by Eq. (10) with the characteristic
intensity $\overline{S}=2S_{R}$.

## References

  * (1) J. M. Dawlaty, S. Shivaraman, M. Chandrashekhar, F. Rana, and M. G. Spencer Appl. Phys. Lett. 92, 042116 (2008); D. Sun, Z.-K. Wu, C. Divin, X. Li, C. Berger, W. A. de Heer, P. N.First, and T. B. Norris, Phys. Rev. Lett. 101, 157402 (2008); R. W. Newson, J. Dean, B. Schmidt, and H. M. van Driel, Opt. Exp. 17, 2326 (2009).
  * (2) F. Rana, P. A. George, J. H. Strait, J. Dawlaty, S. Shivaraman, Mvs Chandrashekhar, and M. G. Spencer, Phys. Rev. B 79, 115447 (2009); P. N. Romanets and F.T. Vasko, Phys. Rev. B 81, 085421 (2010).
  * (3) F.T. Vasko and O.E. Raichev, Quantum Kinetic Theory and Applications (Springer, N.Y., 2005).
  * (4) P. N. Romanets and F.T. Vasko, Phys. Rev. B 81, 241411 (2010).
  * (5) Z. Sun, T. Hasan, F. Torrisi, D. Popa, G.Privitera, F. Wang, F. Bonaccorso, D. M. Basko, and A. C. Ferrari, ACS Nano 4, 803 (2010); Z. Sun, D. Popa, T. Hasan, F. Torrisi, F. Wang, E. J. R. Kelleher, J. C. Travers, and A. C. Ferrari, Nano Research 3, 653 (2010).
  * (6) H. Zhang, D.Y. Tang, L.M. Zhao, Q. Bao, K.P. Loh, Opt. Express, 17 17630 (2009); H. Zhang, D. Tang, R.J. Knize, L. Zhao, Q. Bao, K.P. Loh, Appl. Phys. Lett., 96, 111112 (2010).
  * (7) V. Ryzhii, M. Ryzhii, and T. Otsuji, J. Appl. Phys. 101, 083114 (2007); A. Satou, F. T. Vasko, and V. Ryzhii, Phys. Rev. B 78, 115431 (2008).
  * (8) L.A. Falkovsky, Phys. Usp. 51 887 (2008); M. V. Strikha and F.T. Vasko, Phys. Rev. B 81, 115413 (2010).
  * (9) N.M.R. Peres, Rev. Mod. Phys. 82, 2673 (2010); F. T. Vasko and V. Ryzhii, Phys. Rev. B 76, 233404 (2007); X. Hong, K. Zou, and J. Zhu, Phys. Rev. B 80, 241415 (2009).
  * (10) H. Suzuura and T. Ando, J. Phys. Soc. Japan, 77, 044703 (2008); S. Piscanec, M. Lazzeri, F. Mauri, A. C. Ferrari, and J. Robertson, Phys. Rev. Lett. 93, 185503 (2004).
  * (11) D. Potter, Computational Physics (J. Wiley, London, 1973).
  * (12) J. H. Chen, C. Jang, S. Xiao, M. Ishigami, and M. S. Fuhrer, Nature Nanotech. 3, 206 (2008); S. Fratini and F. Guinea, Phys. Rev. B 77, 195415 (2008).
  * (13) B. A. Ruzicka, S. Wang, L. K. Werake, B. Weintrub, K. P. Loh, and H. Zhao, arXiv:1005.3850.
  * (14) C. H. Lui, K. F. Mak, J. Shan, and T. F. Heinz, Phys. Rev. Lett. 105, 127404 (2010); H. Wang, J. H. Strait, P. A. George, S. Shivaraman, V. B. Shields, Mvs Chandrashekhar, J. Hwang, F. Rana, M. G. Spencer, C. S. Ruiz-Vargas, and J. Park, Appl. Phys. Lett. 96, 081917 (2010).


# Quantifying Distribution Shifts and Uncertainties for Enhanced Model
Robustness in Machine Learning Applications

Vegard Flovik<EMAIL_ADDRESS>Department of Electronic Systems, Norwegian
University of Science and Technology, Trondheim, Norway

###### Abstract

Distribution shifts, where statistical properties differ between training and
test datasets, present a significant challenge in real-world machine learning
applications where they directly impact model generalization and robustness.
In this study, we explore model adaptation and generalization by utilizing
synthetic data to systematically address distributional disparities. Our
investigation aims to identify the prerequisites for successful model
adaptation across diverse data distributions, while quantifying the associated
uncertainties. Specifically, we generate synthetic data using the van der
Waals equation for gases and employ quantitative measures such as Kullback-
Leibler divergence, Jensen-Shannon distance, and Mahalanobis distance to
assess data similarity. These metrics enable us to evaluate both model
accuracy and quantify the associated uncertainty in predictions arising from
data distribution shifts. Our findings suggest that utilizing statistical
measures, such as the Mahalanobis distance, to determine whether model
predictions fall within the low-error ”interpolation regime” or the high-error
”extrapolation regime” provides a complementary method for assessing
distribution shift and model uncertainty. These insights hold significant
value for enhancing model robustness and generalization, essential for the
successful deployment of machine learning applications in real-world
scenarios.

##  1 Introduction

In machine learning, ensuring model accuracy and reliability is paramount for
successful deployment in real-world applications. However, challenges arise
when the statistical properties of the data differ between training and test
datasets, a phenomenon commonly referred to as distribution shift [16]. This
presents significant challenges in dynamic systems where data distributions
evolve over time, as well as in transfer learning scenarios where models
encounter disparate data distributions [15]. Despite efforts to align source
and target domains, practical constraints often lead to distributional
disparities, impacting model accuracy and reliability in real-world
applications [3].

Numerous examples illustrate the implications of distribution shift for
machine learning models across various domains. Consider, for example, a
scenario where machine learning models monitor the performance of wind
turbines in a wind farm. Covariate shift occurs if environmental factors
affecting turbine performance, such as wind speed, wind direction, and
temperature, differ between historical training data and current operational
data due to changes in weather patterns or site conditions. Conversely, target
drift arises if the relationship between these environmental factors and
turbine performance changes over time due to factors like aging equipment or
modifications in operational protocols.

Similarly, in medical image analysis, models trained to detect tumors or
diseases from X-rays may fail when deployed in different hospitals with
varying equipment, leading to potential misdiagnoses and harmful outcomes for
patients [2]. Autonomous vehicles may encounter unforeseen scenarios not
adequately represented in training data, posing substantial risks to passenger
safety and public trust in autonomous driving technology [13]. Recommendation
systems deployed in online platforms may provide biased or inaccurate
suggestions when user preferences evolve over time or when operating in new
contexts, leading to sub-optimal user experiences and potential ethical
concerns [17]. These examples underscore the critical importance of addressing
distribution shift and its impact on model performance in real-world
applications spanning numerous industries and domains.

Previous studies by Shimodaira [18] and Liu et al. [10] have also shed light
on the detrimental effects of distributional disparities between training and
test datasets on model performance and generalization. These findings have
initiated the exploration of various domain adaptation methods aimed at
mitigating these disparities. For example, adversarial training [5] and
discrepancy-based approaches [11], have emerged as potential solutions to
align source and target domain, thereby reducing distributional disparities.
Furthermore, data augmentation techniques, including mixup [22] and generative
adversarial networks (GANs) [7], have also been explored to enhance model
robustness against distribution shift.

Additionally, Bayesian methods [6] and conformal prediction [1] offer
alternative approaches for uncertainty quantification in machine learning.
Bayesian methods provide a probabilistic framework for uncertainty estimation,
modeling uncertainty as a distribution over parameters or predictions.
Conformal prediction, on the other hand, offers a pragmatic approach by
providing prediction intervals with valid coverage guarantees. Integrating
these methods into machine learning pipelines can provide additional insights
into model uncertainty and enhance model robustness in real-world scenarios.

Previous work by Ovadia et al. [14] have also emphasized the importance of
quantifying distribution shift and its impact on model uncertainty. Proper
quantitative measures provide valuable insights into the extent of
distributional disparities, and these measures serve as critical tools for
assessing the robustness and generalization abilities of machine learning
models.

In this study, we aim to address a gap in understanding distribution shift and
its impact on model performance by investigating the role of synthetic data in
evaluating model robustness. Leveraging synthetic data generated using the Van
der Waals equation for gases, we systematically quantify distributional
disparities and assess their impact on model performance and uncertainty.

###  1.1 Research Questions

This study aims to answer the following fundamental questions surrounding the
challenges with data distribution shifts:

  * •

How can we quantify the degree of distribution shift and data similarity?

  * •

How well can machine learning models generalize across different data sets?

  * •

What uncertainties arise from differences in data distribution between train
and test data, and how can we quantify this?

By systematically investigating these questions through a comprehensive
analysis of data similarity measures and model uncertainty metrics, we aim to
advance our understanding of distribution shift and provide actionable
insights for enhancing the robustness and generalization of machine learning
models in real-world applications.

### Paper Structure

The paper is structured as follows:

  * •

Section 2 provides an overview of the main objectives and outlines the
methodology employed in the study.

  * •

Section 3 introduces the various measures and metrics used for quantifying
data similarity and their application within our study.

  * •

Section 4 introduces the architecture of our model and details the Monte Carlo
Dropout method for assessing model uncertainty.

  * •

Section 5 describes the process of generating our synthetic dataset using the
van der Waals equation.

  * •

Section 6 presents the technical details and results of our analysis.

  * •

Section 7 concludes the paper by summarizing the main outcomes and discussing
the conclusions drawn from our experiments.

##  2 Objectives and Methodology

To investigate these research questions in further detail, we perform two
different experiments, both based on generating synthetic data where we have
full control of the distributional properties across the different datasets.
Further details and methodology, including data generation and analysis, is
described in the following sections for experiments 1 and 2 respectively.

###  2.1 Experiment 1: Changes in Feature-Target Correlations

Objective: Experiment 1 aims to explore how variations in feature-target
correlations impact model accuracy and how data similarity can be quantified.
By systematically varying the feature-target correlations across different
datasets, we seek insights into the sensitivity of model performance to these
variations.

Significance: Understanding the effect of feature-target correlations and
distribution shift on model accuracy is crucial for developing robust models
and for assessing their ability to generalize to unseen data beyond the
training set.

Methodology: To investigate this, we employ the ”ideal gas” approximation to
generate synthetic training data. Subsequently, we generate datasets for
different gases, each with properties that deviate from ideal gas behavior,
approximated by the Van der Waals equation. In addition to visual comparisons
of the data distributions, we explore measures such as the Kullback-Leibler
Divergence and Jensen-Shannon Distance to quantify data similarity.
Subsequently, we train a machine learning model using the ideal gas data and
proceed to predict the properties of the other gases, assessing how
distribution shift impacts model accuracy.

###  2.2 Experiment 2: Feature Distribution Drift

Objective: Experiment 2 aims to investigate changes in the feature
distribution and their impact on model accuracy and uncertainty, with a focus
on quantifying these properties.

Significance: In real-world dynamical systems, variations in pressure,
temperature, or other properties over time can induce shifts in the feature
distribution, referred to as covariate shift. Detecting when the model
operates beyond the bounds of its training distribution and understanding the
implications for prediction accuracy and uncertainty are thus crucial for
maintaining model reliability in practical applications.

Methodology: We simulate changes in the feature distribution between the
training and test datasets using the ideal gas approximation. For each data
point in the test set, we measure the deviation of its feature values from the
training data distribution using the Mahalanobis distance. Subsequently, we
analyze model predictions and evaluate how distribution shift influences both
accuracy and uncertainty. During inference, we employ Monte Carlo Dropout as
an estimate of the model’s uncertainty, which is then correlated with the
degree of distribution shift quantified through the Mahalanobis distance.

##  3 Quantifying Data Similarity

A key aspect of our study involves quantifying the similarity between the
various datasets generated by the van der Waals equation. In this section, we
provide a brief theoretical background on the key metrics used to quantify
data similarity and distributional shifts, namely KL-Divergence, Jensen-
Shannon Distance, and Mahalanobis Distance. These metrics allow us to quantify
the extent to which one dataset diverges from another. We can thus evaluate
the overall difference or similarity of the various data distributions, and
how this affects the model’s ability to generalize beyond the training data.

###  3.1 Kullback-Leibler Divergence (KL-Divergence)

KL-Divergence serves as a measure of relative entropy between two probability
distributions [9]. Given two probability distributions $P(x)$ and $Q(x)$, KL-
Divergence is defined as:

$D_{\text{KL}}(P\|Q)=\sum_{x\in X}P(x)\log\frac{P(x)}{Q(x)}$ (1)

KL-Divergence quantifies the information lost when $Q(x)$ is used to
approximate $P(x)$. It is non-negative and equals zero if and only if the
distributions are identical.

In machine learning and statistics, it has been widely used in tasks such as
domain adaptation, transfer learning, and anomaly detection to quantify the
dissimilarity between data distributions [19]. In the context of comparing
data distribution shifts in this study, KL-Divergence thus provides a
quantitative measure of how much one dataset, $P(x)$, differs from another,
$Q(x)$.

###  3.2 Jensen-Shannon Distance

The Jensen-Shannon Distance is a symmetric measure of the similarity between
two probability distributions. It is derived from the KL-Divergence and is
defined as:

$D_{\text{JS}}(P,Q)=\frac{1}{2}D_{\text{KL}}(P\|M)+\frac{1}{2}D_{\text{KL}}(Q\|M)$
(2)

where $M$ is the average distribution of $P$ and $Q$.

Unlike KL-Divergence, which is a divergence measure, Jensen-Shannon Distance
is a proper metric. It combines elements of KL-Divergence to evaluate the
overall difference or similarity between datasets in a more balanced manner.
Furthermore, it is bounded between 0 and 1, where 0 indicates identical
distributions and 1 indicates completely different distributions. This bounded
nature of the Jensen-Shannon Distance provides some advantages when used as a
metric to quantify data distribution similarity. In the machine learning
literature, Jensen-Shannon Distance is often employed in tasks such as
clustering, generative modeling, and classification to assess the similarity
between probability distributions arising from different datasets or model
predictions [8].

###  3.3 Mahalanobis Distance

In statistical analysis and pattern recognition, the Mahalanobis distance is a
metric used to measure the distance between a point and a distribution [12].
Unlike Euclidean distance, which treats all dimensions equally, the
Mahalanobis distance adjusts for the covariance structure of the data,
providing a more accurate measure of similarity or dissimilarity.

Mathematically, the Mahalanobis distance $D_{M}(x)$ for a data point $x$ from
a distribution with mean vector $\mu$ and covariance matrix $S$ is defined as:

$D_{M}(x)=\sqrt{(x-\mu)^{T}S^{-1}(x-\mu)}$ (3)

Here:

  * •

$x$ is the vector representing the data point.

  * •

$\mu$ is the mean vector of the distribution.

  * •

$S$ is the covariance matrix of the distribution.

In this study, we utilize the Mahalanobis distance to quantify the similarity
of data points in the test set from the training distribution. By assessing
whether model predictions fall within the known data distribution
(interpolation) or extend beyond it (extrapolation), we obtain a deeper
understanding of the model’s ability to generalize and extrapolate beyond the
training data.

##  4 Model Architecture and Uncertainty Quantification

The selected model architecture in this study is deliberately kept simple,
allowing us to focus on the core research objectives related to data
distribution shifts and uncertainties rather than complex modeling approaches.
The model architecture includes the following components:

  1. 1.

Input Layer: This layer defines the shape of the input data.

  2. 2.

Dense Layers: The model includes three dense layers with 64, 64, and 32
neurons respectively. Each neuron in a dense layer is connected to every
neuron in the previous layer. The activation function used in these layers is
chosen as the Exponential Linear Unit (ELU).

  3. 3.

Custom Dropout Layer: This model incorporates a custom dropout layer with a
dropout rate of 0.1, randomly setting 10% of input units to zero. The custom
layer also enables dropout during inference to facilitate the Monte Carlo
Dropout technique, as discussed in Section 4.1.

  4. 4.

Output Layer: The final layer consists of a single neuron, representing the
output of the model. It does not apply any activation function, which is
standard for regression models where the goal is to predict continuous
numerical values.

For additional details about the technical implementation, readers are
directed to the code available on GitHub:
https://github.com/veflo/uncert_quant

###  4.1 Assessing Model Uncertainty with Monte Carlo Dropout

In this study, we employ Monte Carlo Dropout as a technique to estimate model
uncertainty during inference. This approach is particularly relevant for
capturing epistemic uncertainty arising from the model’s lack of knowledge or
understanding about certain regions of the input space [4].

Dropout is a regularization technique commonly used during model training. It
works by randomly ”dropping out” (i.e., setting to zero) a fraction of the
units in a layer during each training iteration. This prevents units from co-
adapting too much and ensures that the network learns a more robust and
generalizable representation of the data. Typically, dropout is turned off
during inference to obtain deterministic predictions. However, in Monte Carlo
Dropout, dropout is kept active during inference, and multiple forward passes
through the network are performed for the same input, resulting in different
predictions $y_{i}$ sampled from the model’s predictive distribution:

$y_{i}\sim p(y|x,\theta)$ (4)

where $\theta$ represents the model parameters.

$x_{1}$$x_{2}$$x_{3}$$h_{11}$$h_{12}$$h_{13}$$h_{14}$$h_{15}$$h_{21}$$h_{22}$$h_{23}$$h_{24}$$h_{25}$$\hat{y}$
Figure 1: Neural Network with Dropout, as indicated by red nodes

We can then aggregate this ensemble of predictions and calculate the mean,
$\mu$, and standard deviation, $\sigma$, as follows:

$\mu=\frac{1}{N}\sum_{i=1}^{N}y_{i}$ (5)

$\sigma=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\mu)^{2}}$ (6)

where $N$ is the number of forward passes.

By performing Monte Carlo Dropout during inference and investigating the
correlation between model uncertainty and the degree of distribution shift, as
quantified through the Mahalanobis distance, we thus aim to gain insights into
the model’s ability to generalize to unseen data with different statistical
properties.

##  5 Data Generation using the van der Waals Equation

The van der Waals equation is an equation of state for gases that provides a
more accurate description of the behavior of real gases, especially at high
pressures and low temperatures, where the assumptions of the ideal gas law
start to break down [21].

To generate synthetic datasets for our experiments, we utilize this equation
to simulate the behavior of gases under various conditions. Specifically, we
employ the equation to generate datasets for different gases, each exhibiting
deviations from ideal gas behavior, thus creating a range of feature-target
correlations.

The van der Waals equation is given by:

$\left(P+\frac{aN^{2}}{V^{2}}\right)(V-Nb)=NRT$ (7)

Where:

  * •

$P$ is the pressure of the gas (in atmospheres, atm)

  * •

$V$ is the volume of the gas (in liters, L)

  * •

$T$ is the temperature of the gas (in kelvin, K)

  * •

$N$ is the number of moles of the gas (in moles, mol)

  * •

$a$ and $b$ are van der Waals constants specific to each gas
($a=\text{[L}^{2}\,\text{atm/mol}^{2}$, $b=\text{[L/mol]}$)

  * •

$R$ is the gas constant
($0.0821\,\frac{\text{L}\cdot\text{atm}}{\text{mol}\cdot\text{K}}$)

The term $\frac{aN^{2}}{V^{2}}$ corrects for the attractive forces between gas
molecules, while the term $Nb$ corrects for the volume occupied by the gas
molecules themselves. It’s important to note that the van der Waals constants
$a$ and $b$ vary depending on the specific gas and are experimentally
determined.

In table 1 below, we have listed these constants for a selection of common
gases that we have used in this study.

Gas |  $a$ [L2 atm/mol2] |  $b$ [L/mol]
---|---|---
Ideal Gas | 0 | 0
Hydrogen (H2) | 0.244 | 0.0266
Helium (He) | 0.0346 | 0.0237
Neon (Ne) | 0.211 | 0.0174
Argon (Ar) | 1.355 | 0.032
Xenon (Xe) | 4.00 | 0.051
Nitrogen (N2) | 1.390 | 0.0391
Oxygen (O2) | 1.360 | 0.0318
Carbon Dioxide (CO2) | 3.610 | 0.0427
Methane (CH4) | 2.250 | 0.0428
Table 1: Van der Waals constants for gases used in this study. Adapted from
standard references on thermodynamics or physical chemistry.

For our data generation process, we first generate the variables for
Temperature ($T$), Volume ($V$), and Moles ($N$) from Gaussian distributions
with defined mean ($\mu$) and variance ($\sigma^{2}$). Specifically, we
generate $T$, $V$, and $N$ as follows:

$\displaystyle T$ $\displaystyle\sim\mathcal{N}(\mu_{T},\sigma_{T}^{2})$ (8)
$\displaystyle V$ $\displaystyle\sim\mathcal{N}(\mu_{V},\sigma_{V}^{2})$ (9)
$\displaystyle N$ $\displaystyle\sim\mathcal{N}(\mu_{N},\sigma_{N}^{2})$ (10)

We can then plug the gas constants and parameter distributions into the van
der Waals equation to calculate the corresponding Pressure ($P$), which serves
as our target variable for the synthetic dataset. We define the parameters for
generating synthetic datasets as indicated in the tables below for experiments
1 and 2, respectively:

Parameter | Mean | Std
---|---|---
Temperature ($T$) | $\mu_{T}=300\,\text{K}$ | $\sigma_{T}=25\,\text{K}\,$
Volume ($V$) | $\mu_{V}=50\,\text{L}$ | $\sigma_{V}=5\,\text{L}\,$
Moles ($N$) | $\mu_{N}=15\,\text{mol}$ | $\sigma_{N}=1\,\text{mol}\,$
Table 2: Parameters used for generating data for experiment 1 Dataset 1:
---
Temperature | $\mu_{T_{1}}=273\,\text{K}$ | $\sigma_{T_{1}}=50\,\text{K}$
Volume | $\mu_{V_{1}}=10\,\text{L}$ | $\sigma_{V_{1}}=1\,\text{L}$
Moles | $\mu_{N_{1}}=10\,\text{mol}$ | $\sigma_{N_{1}}=1\,\text{mol}$
Dataset 2:
Temperature | $\mu_{T_{2}}=300\,\text{K}$ | $\sigma_{T_{2}}=50\,\text{K}$
Volume | $\mu_{V_{2}}=9\,\text{L}$ | $\sigma_{V_{2}}=1.5\,\text{L}$
Moles | $\mu_{N_{2}}=11\,\text{mol}$ | $\sigma_{N_{2}}=1\,\text{mol}$
Table 3: Parameters used for generating data for experiment 2

  * •

For experiment 1, we use identical distributions for the input variables
(Temperature, Volume and Moles), but different van der Waals constant for the
various gases.

  * •

For experiment 2, we use the ideal gas approximation for both datasets (a = b
= 0), but rather change the distributions for the input variables, as
illustrated in table 3

##  6 Results and Discussion

This section presents a summary of the key findings from the study. For
comprehensive details and access to the code implementation of the analysis
methodology, readers are directed to the material available on GitHub:
https://github.com/veflo/uncert_quant

###  6.1 Experiment 1:

We begin by generating a synthetic dataset for the selection of gases, as
described in further detail in Section 5. The variables Temperature, Moles,
and Volume are uniformly generated for all gases, while the corresponding
pressure is calculated using the van der Waals equation. This equation
introduces corrections to the ideal gas approximation, resulting in varying
degrees of adjustments based on the properties of each gas. Figure 2(a)
illustrates the resulting pressure distribution, highlighting subtle
differences among the different gases. In addition to the generated datasets
for the other gases, we have also set aside a small subset of the ideal gas
data as an ”in distribution” test set.

(a) Pressure distribution, as calculated from van der Waals equation.

(b) t-SNE visualization

Figure 2: Pressure distribution and t-SNE plot for the various gases

To visualize the degree of similarity between the datasets, we employ t-SNE, a
popular technique for visualizing high-dimensional data [20]. The t-SNE plot,
as shown in Figure 2(b), reveals no clear clusters or outliers in this case.
However, given the subtle differences between the gases and that the ideal gas
approximation is a fairly accurate description in many scenarios, this result
is not unexpected. To perform a more quantitative investigation of data
similarity and distribution shift beyond visual comparison, we utilize KL-
divergence and Jensen Shannon Distance, as introduced in Section 3.

Next, we train our machine learning model using the ideal gas training data,
and then evaluate its generalization to the other datasets. Figure 3
illustrates the model’s predictions vs. actual pressure values for the
different gases.

Figure 3: Predicted vs Real Pressure, Experiment 1

In addition to the visual comparison of the predictions, we also compute the
Mean Absolute Percentage Error (MAPE). In Table 4, we summarize both the KL-
divergence, JS Distance, and MAPE between the ideal gas training dataset and
other gases.

Dataset | KL-Div | JS Distance | MAPE
---|---|---|---
Ideal Gas (subset) | 0.10 | 0.14 | 0.48
Nitrogen | 0.33 | 0.26 | 0.73
Oxygen | 0.57 | 0.33 | 0.86
Argon | 0.52 | 0.32 | 0.84
Carbon Dioxide | 5.49 | 0.73 | 3.30
Methane | 1.75 | 0.49 | 1.56
Hydrogen | 0.31 | 0.25 | 0.64
Helium | 0.55 | 0.33 | 0.76
Neon | 0.14 | 0.17 | 0.54
Xenon | 6.11 | 0.75 | 3.55
Table 4: KL-div, JS Distance, and prediction error between the training
dataset (Ideal Gas) and the datasets for the other gases

From these results, we note that smaller corrections to the van der Waals
equation typically yield lower KL-divergence and JS Distance, indicating
greater similarity to the training dataset. Consequently, datasets more
closely resembling the ideal gas approximation result in lower prediction
errors. Furthermore, it is worth noting that despite being drawn from the same
distribution, we observe a non-zero KL-divergence and JS distance between the
training data and the ”in-distribution” test set, referred to as ”Ideal Gas
(subset).” This observation underscores the inherent statistical variations
resulting from a finite number of generated data points, which in this case
were limited to 10,000 samples.

This phenomenon reflects a conventional practice in assessing machine learning
models, where a portion of the training data is reserved for testing purposes.
However, when the test data originates from the same distribution as the
training data, it often leads to overly optimistic estimates of model
accuracy. These discrepancies often become apparent upon deployment ”in the
wild,” where the model may encounter a different data distribution, resulting
in a corresponding degradation in accuracy.

(a) KL-Divergence vs. Prediction error

(b) JS-Distance vs. Prediction error

Figure 4: Comparison of KL-Divergence and JS-Distance vs. Prediction Error,
demonstrating a clear correlation between the degree of distribution shift and
the consequent decline in model accuracy.

As illustrated in Figure 4(a), we observe a linear trend between KL-divergence
and MAPE, while the relationship with JS-Distance in Figure 4(b) appears less
pronounced. This discrepancy may result from JS-Distance being a bounded
metric between [0,1], unlike the unbounded nature of the KL-divergence.
However, a deeper understanding of this scaling behavior would require more
detailed investigations, which are beyond the scope of this study.

In summary, the outcomes of Experiment 1 underscore the potential use of KL-
divergence and JS Distance as quantitative indicators of distribution shift,
serving as predictors of reduced model accuracy and increased uncertainty.
These metrics thus offer valuable insights for monitoring the performance and
robustness of machine learning models in real-world scenarios, where varying
degrees of distribution shift are frequently encountered.

###  6.2 Experiment 2

Experiment 2 extends the analysis of Experiment 1 by primarily focusing on
quantifying prediction uncertainty, thereby enhancing the assessment of model
reliability. While Experiment 1 simulated ”target drift”, where the
correlations between the input variables and the target differ across various
datasets, Experiment 2 focuses on ”covariate shift”. Here, the distribution of
the input variables themselves has changed, while the correlations to the
target variable remain constant

Specifically, we generate data using the ideal gas approximation for two
datasets, with the differences being determined by the parameters of the
distributions used for data generation. The resulting data distributions for
the generated data are illustrated in Figure 5. Further details on the data
generation process are also described in Table 3 in section 5.

(a) Volume

(b) Moles

(c) Temperature

(d) Pressure

Figure 5: Distribution of Volume, Moles, Temperature, and the resulting
calculated Pressure for the train (”in distribution”) vs. test (”out-of-
distribution”) datasets

We also here calculate the corresponding KL-divergence and JS-Distance between
the ”in distribution” training data vs. the ”out-of-distribution” test set, as
presented in Table 5. This is done by setting aside a subset of the training
data as an ”in distribution” test set, following a similar procedure as in
Experiment 1.

Dataset | KL-Div | JS Distance | MAPE
---|---|---|---
In distribution (subset) | 0.26 | 0.16 | 0.26
Out of distribution | 1.20 | 0.47 | 0.60
Table 5: KL-div, JS Distance, and prediction error for the ”in distribution”
vs ”out-of-distribution” datasets

We then proceed to train our machine learning model using the training data
and utilize it to make predictions for both the ”in-distribution” and ”out-of-
distribution” test sets. Figure 6 illustrates the predicted vs. real pressure
values. The largest deviations occur in regions of low and high pressure,
where the model encounters data outside the training distribution, as shown in
Figure 5(d). The associated prediction error (MAPE) is also provided in Table
5, for comparison between the ”in-distribution” and ”out-of-distribution” test
sets.

Figure 6: Predicted vs. Real Pressure, Experiment 2

####  6.2.1 Uncertainty Quantification

To evaluate model uncertainty, we employ the Monte Carlo Dropout technique, as
discussed in detail in section 4.1. This technique involves conducting N=100
forward passes through the network for each input, resulting in multiple
predictions $y_{i}$ sampled from the model’s predictive distribution. By
calculating the mean and standard deviation of these predictions for each
input data point, we thus obtain an estimate of the model’s epistemic
uncertainty.

We then compute the Mahalanobis distance, as introduced in section 3.3, which
serves as a metric to assess how far each datapoint deviates from the training
data distribution. Figure 7 illustrates the distribution of Mahalanobis
distances for the training data, including the limit indicating the 95th
percentile of this distribution. By computing the Mahalanobis distance for all
datapoints in the test set, we can then gain insights into the error and
uncertainty of model predictions relative to their distance from the training
distribution.

Figure 7: Mahalanobis distance for the ”In distribution” dataset

As illustrated in Figure 8(a), datapoints with a low Mahalanobis distance,
indicating proximity to the training distribution, generally exhibit lower
prediction errors. Conversely, an increase in the Mahalanobis distance,
particularly beyond the 95th percentile cutoff value, signifies a departure
from the training distribution. This correlates with a corresponding rise in
prediction errors and the spread of predictions, indicating increased
uncertainty.

This behavior is further explored in Figure 8(b), which illustrates the
relationship between prediction error and the standard deviation of
predictions, commonly used as an indicator of model uncertainty in the Monte
Carlo Dropout technique.

(a) Mahalanobis distance vs. prediction error

(b) Standard deviation vs. prediction error

Figure 8: Comparison of Mahalanobis distance and standard deviation as
indicators of prediction uncertainty.

While low standard deviation values often coincide with low prediction errors,
the correlation is less pronounced compared to the Mahalanobis distance.

These findings suggests that the Mahalanobis distance may offer superior, or
at least complementary, insights into prediction uncertainty and error
compared to relying solely on the Monte Carlo Dropout technique. Moreover, the
Mahalanobis approach offers significant computational advantages over Monte
Carlo methods. By requiring only a single forward pass per prediction instead
of multiple iterations, it greatly enhances the practicality and efficiency of
the approach.

Importantly, the Mahalanobis distance serves as a metric for assessing
prediction reliability on a per-data-point basis. Data points with Mahalanobis
distances well below the cutoff can be considered relatively accurate and
trustworthy, while those exceeding the cutoff may require further
investigation or precautionary measures.

Combining the Mahalanobis approach with Monte Carlo Dropout techniques, or
other methods for uncertainty quantification, holds promise for a more
comprehensive assessment of model uncertainty and determining when model
predictions can be trusted or not. This allows for making more informed
decisions based on the trustworthiness of the predictions, thereby enhancing
the robustness of machine learning systems deployed in dynamic real-world
environments.

##  7 Summary and Conclusions

This paper presents two experiments aimed at investigating the impact of
distribution shift on machine learning model performance and uncertainty
quantification.

In Experiment 1, synthetic datasets representing different gases were
generated, and a machine learning model was trained using the ideal gas
dataset. By quantifying the degree of distribution shift through the KL-
divergence and Jensen Shannon Distance between the various datasets, we
observed a clear correlation between distribution shift and prediction
accuracy. Higher errors were evident for datasets deviating further from the
ideal gas distribution, emphasizing the critical role of monitoring
distribution shift for maintaining model performance in real-world
applications. Furthermore, we illustrated how KL-divergence and JS Distance
serve as valuable indicators of distribution shift and predictors of model
prediction accuracy degradation.

Experiment 2 focused on uncertainty quantification, particularly in scenarios
of covariate shift. By employing the Monte Carlo Dropout technique for
estimating model uncertainty, and by calculating the Mahalanobis distance, we
correlated uncertainty with the distance of test data points from the training
distribution. Our findings revealed that as data points move away from the
training distribution, both prediction error and uncertainty tend to increase.
The Mahalanobis distance emerged as a promising metric for assessing
prediction reliability on a per-data-point basis, offering additional insights
into when model predictions can be trusted or not.

While uncertainty quantification methods like Bayesian methods [6] and
conformal prediction [1] offer alternative approaches, each presents its own
set of challenges and trade-offs. Bayesian methods provide a robust framework
for uncertainty estimation but face challenges such as computational
complexity and scalability limitations. Conversely, conformal prediction
offers a pragmatic alternative but can be susceptible to shifts in data
distribution over time, impacting prediction reliability.

The choice of methodology ultimately depends on the specific requirements of
the application and the inherent trade-offs between computational complexity,
interpretability, and performance. Rather than aiming to replace existing
methodologies, this study seeks to explore complementary approaches. The
findings presented here serve as a guide towards identifying appropriate
metrics and techniques for effective and robust monitoring and governance of
machine learning applications in real-world scenarios. By understanding the
strengths and limitations of different uncertainty quantification methods,
practitioners can make more informed decisions and tailor their approach to
suit the demands of their specific application domains.

## References

  * [1] Larry Brown and Alice White.  Conformal prediction: A pragmatic approach to prediction reliability.  Machine Learning, 28(4):543–567, 2016.
  * [2] Randall J Ellis, Ryan M Sander, and Alfonso Limon.  Twelve key challenges in medical machine learning and solutions, 2022.
  * [3] Abolfazl Farahani, Sahar Voghoei, Khaled Rasheed, and Hamid R Arabnia.  A brief review of domain adaptation.  Advances in data science and information engineering: proceedings from ICDATA 2020 and IKE 2020, pages 877–894, 2021.
  * [4] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  arXiv preprint arXiv:1506.02142, 2016.
  * [5] Yaroslav Ganin, Evgeniya Ustinova, Hana Ajakan, Pascal Germain, Hugo Larochelle, François Laviolette, Mario March, and Victor Lempitsky.  Domain-adversarial training of neural networks.  Journal of machine learning research, 17(59):1–35, 2016.
  * [6] Jakob Gawlikowski, Cedrique Rovile Njieutcheu Tassi, Mohsin Ali, Jongseok Lee, Matthias Humt, Jianxiang Feng, Anna Kruspe, Rudolph Triebel, Peter Jung, Ribana Roscher, et al.  A survey of uncertainty in deep neural networks.  Artificial Intelligence Review, 56(Suppl 1):1513–1589, 2023.
  * [7] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 27. Curran Associates, Inc., 2014.
  * [8] Sharu Theresa Jose and Osvaldo Simeone.  Information-theoretic bounds on transfer generalization gap based on jensen-shannon divergence.  In 2021 29th European Signal Processing Conference (EUSIPCO), pages 1461–1465. IEEE, 2021.
  * [9] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  The annals of mathematical statistics, 22(1):79–86, 1951.
  * [10] Jiashuo Liu, Zheyan Shen, Yue He, Xingxuan Zhang, Renzhe Xu, Han Yu, and Peng Cui.  Towards out-of-distribution generalization: A survey.  arXiv preprint arXiv:2108.13624, 2021.
  * [11] Mingsheng Long, Yue Cao, Jianmin Wang, and Michael Jordan.  Learning transferable features with deep adaptation networks.  In International conference on machine learning, pages 97–105. PMLR, 2015.
  * [12] Prasanta Chandra Mahalanobis.  On the generalized distance in statistics.  Sankhyā: The Indian Journal of Statistics, Series A (2008-), 80:S1–S7, 2018.
  * [13] Khan Muhammad, Amin Ullah, Jaime Lloret, Javier Del Ser, and Victor Hugo C de Albuquerque.  Deep learning for safe autonomous driving: Current challenges and future directions.  IEEE Transactions on Intelligent Transportation Systems, 22(7):4316–4336, 2020.
  * [14] Yaniv Ovadia, Emily Fertig, Jie Ren, Zachary Nado, D. Sculley, Sebastian Nowozin, Joshua Dillon, Balaji Lakshminarayanan, and Jasper Snoek.  Can you trust your model's uncertainty? evaluating predictive uncertainty under dataset shift.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [15] Sinno Jialin Pan and Qiang Yang.  A survey on transfer learning.  IEEE Transactions on Knowledge and Data Engineering, 22(10):1345–1359, 2010.
  * [16] Joaquin Quiñonero-Candela, Masashi Sugiyama, Anton Schwaighofer, and Neil D Lawrence.  Dataset shift in machine learning.  Mit Press, 2022.
  * [17] Idris Rabiu, Naomie Salim, Aminu Da’u, and Akram Osman.  Recommender system based on temporal models: a systematic review.  Applied Sciences, 10(7):2204, 2020.
  * [18] Hidetoshi Shimodaira.  Improving predictive inference under covariate shift by weighting the log-likelihood function.  Journal of Statistical Planning and Inference, 90(2):227–244, 2000.
  * [19] Masashi Sugiyama, Song Liu, Marthinus Christoffel du Plessis, Masao Yamanaka, Makoto Yamada, Taiji Suzuki, and Takafumi Kanamori.  Direct divergence approximation between probability distributions and its applications in machine learning.  Journal of Computing Science and Engineering, 7(2):99–111, 2013.
  * [20] Laurens Van der Maaten and Geoffrey Hinton.  Visualizing data using t-sne.  Journal of machine learning research, 9(11), 2008.
  * [21] Johannes Diderik Van Der Waals and John Shipley Rowlinson.  On the continuity of the gaseous and liquid states.  Courier Corporation, 2004.
  * [22] Hongyi Zhang, Moustapha Cisse, Yann N Dauphin, and David Lopez-Paz.  mixup: Beyond empirical risk minimization.  arXiv preprint arXiv:1710.09412, 2017.