# Optical nonlinearities in the excited carrier density of atomically thin
transition metal dichalcogenides

D. Erben  Institute for Theoretical Physics, University of Bremen, P.O. Box
330 440, 28334 Bremen, Germany<EMAIL_ADDRESS>A. Steinhoff
Institute for Theoretical Physics, University of Bremen, P.O. Box 330 440,
28334 Bremen, Germany  M. Lorke  Institute for Theoretical Physics, University
of Bremen, P.O. Box 330 440, 28334 Bremen, Germany  F. Jahnke  Institute for
Theoretical Physics, University of Bremen, P.O. Box 330 440, 28334 Bremen,
Germany

###### Abstract

In atomically thin semiconductors based on transition metal dichalcogenides,
photoexcitation can be used to generate high densities of electron-hole pairs.
Due to optical nonlinearities, which originate from Pauli blocking and many-
body effects of the excited carriers, the generated carrier density will
deviate from a linear increase in pump fluence. In this paper, we use a
theoretical approach that combines results from ab-initio electronic-state
calculations with a many-body treatment of optical excitation to describe
nonlinear absorption properties and the resulting excited carrier dynamics. We
determine the validity range of a linear approximation for the excited carrier
density vs. pump power and identify the role and magnitude of optical
nonlinearities at elevated excitation carrier densities for MoS2, MoSe2, WS2,
and WSe2 considering various excitation conditions. We find that for above-
band-gap photoexcitation, the use of a linear absorption coefficient of the
unexcited system can strongly underestimate the achievable carrier density for
a wide range of pump fluences due to many-body renormalizations of the two-
particle density-of-states.

##  1 Introduction

Monolayers of transition metal dichalcogenide (TMDC) semiconductors exhibit
strong Coulomb interaction of their charge carriers that occupy a rich valley
structure in reciprocal space. As a new degree of freedom these valleys are
selectively optically addressable 1, 2, giving rise to a variety of Coulomb-
bound, bright and dark excitons, trions and biexcitons.3, 4, 5, 6, 7, 8, 9, 10
Furthermore, giant band-gap renormalization 11, 12, efficient carrier
scattering13, 14, 15, and the Mott transition of excitons at elevated carriers
densities 16 have been discussed. For a better understanding of the underlying
physics and the prospects of TMDC semiconductors as active materials in future
optoelectronic devices like light-emitting diodes 17, 18, 19, 20, solar cells
17, 18 or lasers 21, 22, 23, 24, 25, experimental techniques such as
photoluminescence and pump-probe spectroscopy are used. 26, 27, 28, 29, 30
These experiments involve photoexcitation of electron-hole pairs, frequently
above the quasi-particle band gap, and benefit from reliable estimates of the
present excited carrier density. The simplest approach to the excitation
density is based on the frequency-dependent linear absorption coefficient. 31,
32, 33 Normalizing the pump fluence to the laser spot size and taking the
absorption coefficient at the excitation energy together with the photon
energy yields a reasonable approximation for low intensities. However, as soon
as the exciting laser pulse is intense enough, effects such as phase space
filling are expected to reduce the absorption in a nonlinear way. Furthermore,
excited carriers will induce screening of the Coulomb interaction, thereby
reducing excitonic effects and modifying interband Coulomb enhancement in the
interacting density of states entering the optical absorption. Therefore, it
is desirable to quantify the involved nonlinearities microscopically and
calculate the resulting excitation densities for typical experimental
situations. In this paper we provide a numerical analysis of the photo-excited
charge carrier dynamics as well as optical nonlinearities in absorption and
excited carrier density. Our approach combines material-realistic electronic
state and interaction matrix element calculations with semiconductor Bloch
equations including many-body effects.

##  2 Photoexcited charge carrier densities in TMDCs

Figure 1: Optically excited charge carrier density vs. pump fluence for
excitation of MoS2 with a 150 fs laser pulse tuned to the band-gap energy (a),
above the C exciton (b), and well above band gap corresponding to the Ti-
sapphire laser emission wavelength (c). Solid lines and symbols represent
calculated carrier densities including optical nonlinearities. Dashed lines
correspond to carrier densities obtained from the calculated linear absorption
coefficients of the unexcited system at the respective energies. Absorption
spectra of MoS2 for a given excited carrier density in thermal equilibrium at
300K are shown in (d) together with the energetic position and spectral width
of the pump pulses corresponding to (a)-(c). The carrier density increases
from zero (black line) to $3\times 10^{14}$ cm-2 (gray line) corresponding to
the data points in (e). The nonlinear behaviour of the extracted absorption
coefficients vs. carrier density at the three pump energies is depicted in
(e).

A microscopic theory for the photoexcited charge carrier density in TMDC
semiconductors requires addressing several subtopics and suitable interfacing
of the methods used in these parts. For electronic state calculations, we
start from density functional theory (DFT) with many-body perturbation theory
corrections in GW approximation and perform a basis set reduction to the
relevant part of the band structure by constructing a six-band lattice
Hamiltonian.11 Four conduction and two valence bands are used to describe the
dominant contributions to the photoexcitation process and the subsequent
excited carrier dynamics. Light-matter interaction and the resulting excited
carrier dynamics are treated using the semiconductor Bloch equations (SBE)11,
34 which include excitonic effects, Pauli-blocking due to excited carriers,
energy renormalizations (shifts of the band structure) due to the Coulomb
interaction of excited carriers, as well as carrier scattering (equilibration
and thermalization) processes. The SBE can be used for calculating
photoexcited carrier densities and optical absorption spectra under the
influence of the above discussed effects. The theory includes a consistent
calculation of Coulomb and dipole interaction matrix elements, the description
of screening effects (background lattice contributions including encapsulation
layers and excited carrier contributions), as well as dephasing and light
propagation effects in the system. Details are provided in the Methods
section.
The density of photoexcited charge carriers at a specific fluence is obtained
from a time-domain numerical solution of the SBE into steady state without the
presence of recombination processes. Assuming that the excited-carrier
dynamics is faster than typical recombination processes, our results reflect
the achievable photoexcitation densities in a time window between pulsed
excitation and recombination.
Our findings for the dependence of excited-carrier density on detuning and
pump fluence are provided in Figs. 1 and 2. For better interpretation of the
results, we also show calculated linear absorption spectra for different
constant densities of photoexcited carriers in thermal equilibrium at room
temperature. Only for these calculations of the spectra, the SBE are solved
for fixed quasi-equilibrium occupation probabilities $F_{\mathbf{k}}^{e/h}$
and a weak electric probe field is used to extract the linear optical
susceptibility of the TMDC including many-body effects. The susceptibility is
translated into absorption using the transfer matrix method for a thin layer
in an hBN environment.35 We note, that for optical pumping considered in Fig.
1 (a) - (c) and Fig. 2, photoexcited carrier density and nonequilibirum
population functions $f^{e/h}_{\mathbf{k}}(t)$ vary continuously during the
pulse, which is fully considered in the dynamical calculations.
Results for MoS2 are collected in Fig. 1. For pumping at the single-particle
band gap, a nonlinear relation between density and fluence is found, which can
be directly related to the Pauli blocking effect at the band edge.
Nevertheless, this effect is partially compensated by quasi-particle energy
renormalizations that spectrally shift the band gap away from the pump pulse
to lower energies, thereby opening up phase space at the pump energy. This
leads to a nearly linear increase of the excited carrier density for fluences
up to $20$ $\mu$J cm-2 and corresponding carrier densities up to $10^{13}$
cm-2. As shown in the inset of Fig. 1 (a), band-structure renormalizations
even overcompensate phase space filling in this fluence range, so that a
slighty superlinear increase of the density is achieved.

Figure 2: Optically excited charge carrier density vs. pump fluence for
excitation of MoSe2 (left), WS2 (center), and WSe2 (right) with a 150 fs laser
pulse tuned to the band-gap energy (red), above the C exciton (blue), and well
above band gap corresponding to the TiSa emission wavelength (green). Solid
lines and symbols represent calculated carrier densities including optical
nonlinearities. Dashed lines correspond to carrier densities obtained from the
calculated linear absorption coefficients of the unexcited system at the
respective energies.

For pumping above the band gap, the linear regime extends only to $5$ $\mu$J
cm-2 and excited-carrier densities are typically larger than for pumping at
the band gap. This is a result of the reduced Pauli blocking as the carrier
scattering rapidly distributes excited carriers away from the pump laser
resonance to lower energies. In comparison to the carrier densities obtained
from extrapolating the zero-density linear absorption (dotted lines), larger
values are obtained. We attribute this to increasing absorption at the pump
laser resonance with increasing excited carrier density, as can be seen in
Fig. 1 (d) and (e). Energy renormalization shifts regions of larger interband
density-of-states from higher energies $h\nu$ onto the pump laser resonance.
Furthermore, the increasing dephasing also contributes to the effect.
A similar behaviour for pumping at the single-particle band gap is found for
other TMDC materials, as shown in Fig. 2. While the superlinear increase of
carrier density due to overcompensation of the phase-space-filling reduction
by energy renormalization is quantitatively comparable in the Mo-based
materials, the initial overshooting extends to larger pump fluences for the
W-based materials. This effect originates from stronger spin-orbit
interaction: the splitting between A and B excitons is increased so that the B
exciton is located above the single-particle gap at zero density. Under
photoexcitation, band-gap renormalization, screening, and Pauli blocking cause
a red shift and bleaching of A and B excitons. The remnant of the B exciton is
shifted into resonance with the pump laser for excitation at the single-
particle band gap, thereby providing an additional boost of the two-particle
density-of-states, see the density-dependent absorption spectra of Fig. S1 in
the Supporting Information. For above band-gap pumping, the dependence of
photoexcited carrier density on fluence is superlinear already at fluences of
several $\mu$J cm-2 for the sulfides, while it is closer to or even below the
linear estimate for the selenides. This can be traced back to the single-
particle band gap of the selenides being smaller by about $0.3$ eV, see Tab.
S1. As a consequence, pumping above the C exciton or at the Ti-sapphire laser
frequency drives band-to-band transitions that are correspondingly further
above the gap. These spectral regions profit less from a renormalization of
the two-particle density-of-states than in the sulfides, see Fig. S1 and Fig.
S2 in the Supporting Information.
It should be noted, that the calculated excited carrier densities in Figs. 1
and 2 follow from a self-consistent solution of the SBE. This solution
includes coupled polarization and population dynamics and, as such, does not
use the absorption spectrum of the system. For the used parameters and
excitation conditions, the system is close to the adiabatic regime, where the
carrier density evolution follows approximately the pulse area. During
different stages of the pump pulse, the population increase is approximately
determined by the absorption of the system for the momentary carrier
distribution. Correspondingly, the absorption spectrum of the system in Fig. 1
(d) and Fig. S1 can be used to interprete the results. In the realized
incoherent regime with ultra-fast equilibration of carriers, absorption
changes during the pump pulse reflect the excited-carrier nonlinearities.

Figure 3: Time evolution of the optically excited carrier density in MoS2 for
spin-up electrons (a) and holes (c) resolved for the different valleys in the
Brillouin zone when pumping at the band gap. The corresponding band-structure
renormalizations $\Delta
E_{\nu}^{\lambda}=|\varepsilon_{\nu}^{\lambda}(0)-\varepsilon_{\nu}^{\lambda}(t)|$
for these valleys are provided in (b) and (d). Curves for K${}^{\prime}_{e}$
and $\Sigma^{\prime}_{e}$ have been omitted for better visibility as they show
the same overall behaviour as Ke and $\Sigma_{e}$ except that
$\Sigma^{\prime}_{e}$ is renormalized weaker than K${}^{\prime}_{e}$. The
excitation pulse is represented by the gray shaded area.

We also note, that radiative coupling of the induced interband polarization
$\mathbf{P}_{\textrm{TMDC}}(t)$ within the TMDC layer reduces the effective
electric field $\mathbf{E}(t)$ according to Eq. 2, see Methods section. Thus,
less coherence is converted into photoexcited carrier density. The net effect
is approximately a 10% excitation-density reduction.
Additional insight into the build-up of the excited-carrier densities in the
different valleys for pumping at the single-particle band gap is provided in
Fig. 3. Panels (a) and (c) show the rise of carrier density in the most
important band-structure valleys (K, K’, $\Sigma$, $\Sigma^{\prime}$). The
large densities at $\Sigma/\Sigma^{\prime}$ reflect a drain of electrons from
the K-valley, where they are excited by the pump pulse. This is due to the
valley energy renormalizations induced by the excited carriers, see Fig. 3 (b)
and (d). A transition to an indirect band gap takes place, as the conductions-
band $\Sigma$-valley experiences a stronger energy lowering than the K-valley
and eventually the $\Sigma$-valley energy is below that of the $K$-valley36.
Since the band-structure renormalizations almost directly follow the pump
pulse and the relaxation of carriers, carrier drain takes place faster than
typical radiative recombination times, which are on the order of $1-10$ ps37
depending on substrate and temperature38. Loss of carriers into the
$\Sigma$-valley has been associated with a sublinear increase of
photoluminescence with pump fluence27.

##  3 Conclusion

This paper provides a critical analysis, to which extent linear absorption
coefficients can be used to determine the excited carrier density in TMDC
semiconductors. A combination of material-realistic electronic-state
calculations with a theory for optical properties under the influence of
excited-carrier many-body effects has been used to model optical pulse
excitation and the subsequent carrier dynamics. Results are different for
pumping at or above the quasi-particle band gap. In the former case, a linear
regime is found for surprisingly large excited carrier densities up to
$1\times 10^{13}$ cm-2 for MoSe2 and up to $2\times 10^{12}$ cm-2 for MoS2,
WS2, and WSe2 due to a compensation of Pauli-blocking nonlinearities in the
absorption and band-structure renormalizations. For above-gap excitation,
excited carrier densities can strongly exceed the linear extrapolation values
since absorption reduction due to Pauli blocking remains weak while stronger
absorption regions at higher energies can be shifted onto the pump resonance
due to giant band-gap shrinkage.

##  4 Methods

###  4.1 Microscopic theory of photoexcited TMDC monolayers

A microscopic description of TMDC semiconductors photoexcited with intense
short laser pulses needs to adress various aspects: (i) The material specific
properties of the band structure as well as dipole and Coulomb interaction
matrix elements. The latter include dielectric screening effects due to the
TMDC lattice as well as possible encapsulation layers and exhibit a specific
momentum dependence owing to the atomically thin material 4. (ii) The dynamics
of the light-matter interaction following a short-pulse excitation and
subsequent carrier relaxation. The effects are reflected in the excited
carrier population dynamics. (iii) Many-body interaction effects of the
excited carriers with a variety of contributions leading to energy
renormalizations of the electronic single-particle states, excitonic effects,
carrier scattering processes, as well as dephasing influencing optical
properties and providing dissipation in the coherent excitation dynamics. (iv)
Light-propagation effects, where the applied electromagnetic laser field
induces a coherent material polarization in the atomically thin TMDC layer
which acts back on the electromagnetic field thus leading to an effective
field within the material.35 While the comprehensive solution of the full
problem is a challenging task, the goal of this paper is to use suitable
approximation schemes for the individual problems and construct interfaces
between treatments of the subproblems.
Our treatment is based on well-established ab-initio electronic-state
calculations of TMDC monolayers using the DFT-GW scheme. The excited carrier
dynamics will be described with semiconductor Bloch equations (SBE). 34, 11 To
interface both approaches, we use a Wannier function lattice Hamiltonian with
a reduced basis set involving the relevant bands for direct coupling to the
optical field. This interface provides material-realistic dipole and Coulomb
interaction matrix elements. In the coherent excitation dynamics, excitonic
effects are included by solving the SBE with these matrix elements entering
the electron-hole interaction terms.
For the treatment of many-body effects we focus on above-bandgap excitation,
as used in many experiments. A typical scenario involves two-pulse excitation,
where the intense pump pulse generates a large density of excited carriers in
the sample and a delayed weak test pulse probes the reflectivity or
transmittivity in a broad spectral range, that can include also the excitonic
resonances. The pump pulse will induce a coherent polarization in the
material, which will undergo a rapid transition from coherent to incoherent
regime as determined by dephasing processes. Owing to the above-bandgap
excitation, the pump pulse photoexcites unbound charge carriers in the form of
an electron-hole-plasma. Initially one expects a nonequilibrium hot-carrier
distribution, which equilibrates due to carrier-carrier Coulomb scattering and
cools down to the lattice temperature by means of carrier-phonon interaction.
In the course of these processes, excited carriers will occupy the valleys of
the band structure. 14, 39, 40, 41 Typically a quasi-equilibrium situation is
established as the relaxation processes occur before recombination sets in38,
37.
The total density of excited electrons and holes $(e,h)$ can be calculated as

$\begin{split}n^{e/h}(t)=\frac{1}{\mathcal{A}}\sum\limits_{\mathbf{k}}f_{\mathbf{k}}^{e/h}(t)\,,\end{split}$
(1)

where
$f_{\mathbf{k}}^{e/h}(t)=\big{\langle}a^{\dagger,e/h}_{\mathbf{k}}(t)\,a^{\vphantom{\dagger}e/h}_{\mathbf{k}}(t)\big{\rangle}$
are the occupation probabilities of electron and hole Bloch states and
$\mathcal{A}$ is the crystal area. The carrier dynamics is governed by the
SBE, which are coupled equations of motion (EOM) for $f_{\mathbf{k}}^{e/h}(t)$
and the microscopic inter-band polarization
$\psi_{\mathbf{k}}^{he}(t)=\big{\langle}a^{\vphantom{\dagger
h}}_{\mathbf{k}}(t)\,a^{\vphantom{\dagger}e}_{\mathbf{k}}(t)\big{\rangle}$
driven by an electric field $\mathbf{E}(t)$ within the atomically thin layer.
The field is given by the applied electromagnetic laser field
$\mathbf{E}_{0}(t)$ together with the reflected field reradiated by the
induced coherent TMDC polarization, as determined by a solution of Maxwell’s
equations 35:

$\displaystyle\mathbf{E}(t)$ $\displaystyle=$
$\displaystyle\mathbf{E}_{0}(t)-\frac{\mu_{0}c_{0}}{2}\frac{\partial}{\partial
t}\mathbf{P}_{\textrm{TMDC}}(t)$ (2) $\displaystyle\approx$
$\displaystyle\mathbf{E}_{0}(t)+\frac{i\omega_{0}}{2\varepsilon_{0}c_{0}n}\mathbf{P}_{\textrm{TMDC}}(t)\,,$

where
$\mathbf{P}_{\textrm{TMDC}}(t)=\frac{1}{\mathcal{A}}\sum_{\mathbf{k},h,e}\mathbf{d}_{\mathbf{k}}^{he}\psi_{\mathbf{k}}^{he}(t)$
is the macroscopic polarization of the TMDC layer, $\omega_{0}$ is the pump
laser frequency, $n$ is the refractive index of the dielectric environment and
$\mathbf{d}_{\mathbf{k}}^{he}$ denotes the dipole matrix element. The approach
takes into account the embedding of the TMDC into encapsulation layers.
Specifically, $\mathbf{E}_{0}(t)$ is the externally applied optical field
within the encapsulation layer. The SBE are given by:

$\begin{split}\frac{\mathrm{d}}{\mathrm{d}t}f_{\mathbf{k}}^{e/h}(t)=\frac{2}{\hbar}\mathrm{Im}\big{[}(\Omega_{\mathbf{k}}^{he})^{*}\psi_{\mathbf{k}}^{he}\big{]}-\frac{f_{\mathbf{k}}^{e/h}(t)-F_{\mathbf{k}}^{e/h}}{\tau_{\textrm{relax}}}\,,\end{split}$
(3) $\displaystyle
i\hbar\frac{\mathrm{d}}{\mathrm{d}t}\psi_{\mathbf{k}}^{he}(t)$
$\displaystyle=$
$\displaystyle\left(\tilde{\varepsilon}_{\mathbf{k}}^{h}+\tilde{\varepsilon}_{\mathbf{k}}^{e}-i\Gamma(n)\right)\psi_{\mathbf{k}}^{he}(t)$
(4)
$\displaystyle-\left(1-f_{\mathbf{k}}^{e}(t)-f_{\mathbf{k}}^{h}(t)\right)\Omega_{\mathbf{k}}^{he}\,.$

Here,
$\Omega_{\mathbf{k}}^{he}=\mathbf{d}_{\mathbf{k}}^{eh}\cdot\mathbf{E}(t)+\frac{1}{\mathcal{A}}\sum_{h^{\prime},e^{\prime},\mathbf{k}^{\prime}}W_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}^{eh^{\prime}he^{\prime}}\psi_{\mathbf{k}^{\prime}}^{h^{\prime}e^{\prime}}(t)$
is the Rabi energy renormalized by the screened Coulomb interaction between
electrons and holes. The inter-band Coulomb interaction give rise to excitonic
resonances below the quasi-particle band gap as well as a redistribution of
oscillator strength between band-to-band transitions. Screening of the Coulomb
interaction consists of two fundamentally different contributions.
“Background” screening due to virtual excitation of valence electrons into
empty conduction-band states as well as the polarizability of surrounding
dielectric material in a heterostructure is included via the dielectric
function $\varepsilon^{-1}_{\textrm{b}}$ contained in Coulomb matrix elements
$V_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}^{eh^{\prime}he^{\prime}}$.
In addition to this, screening induced by photoexcited carriers is captured by
the dielectric function $\varepsilon^{-1}_{\textrm{exc}}$, such that
$W_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}^{eh^{\prime}he^{\prime}}=\varepsilon^{-1}_{\textrm{exc}}(|\mathbf{k}-\mathbf{k}^{\prime}|)V_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}^{eh^{\prime}he^{\prime}}$
is the fully screened Coulomb matrix element. The renormalized single-particle
energies $\tilde{\varepsilon}_{\mathbf{k}}^{\lambda}$ are discussed in the
following section. The Pauli exclusion principle enters via the blocking term
$1-f_{\mathbf{k}}^{e}(t)-f_{\mathbf{k}}^{h}(t)$. Excitation-power-dependent
dephasing processes are described by $\Gamma(n)$. To model the equilibration
and cooling of excited carriers, we use the relaxation time approximation
given by the last term of Eq. (3). The non-equilibrium carrier population
develops into a Fermi distribution $F_{\mathbf{k}}^{e/h}$ for the respective
species of charge carriers on a characteristic timescale
$\tau_{\textrm{relax}}$. As the carrier relaxation does not change the density
of excited carriers, we use for $F_{\mathbf{k}}^{e/h}$ a Fermi-Dirac function
with the same carrier density as the momentary population
$f_{\mathbf{k}}^{e/h}(t)$. For the temperature of $F_{\mathbf{k}}^{e/h}$ the
lattice temperature is taken to model efficient carrier-phonon scattering. For
the relaxation time due to carrier-carrier14 and carrier-phonon40, 42
scattering we use $\tau_{\textrm{relax}}=100$ fs.

###  4.2 Coulomb interaction and band-structure renormalizations

We combine the above theory of photoexcitation with band structures, dipole
matrix elements and bare as well as screened Coulomb matrix elements on a G0W0
level as input for the SBE (3)-(4). Our hybrid approach is designed to yield a
numerically tractable and yet material-realistic model, using a minimal basis
to represent band structures and interaction matrix elements on the same
footing. A six-band model (one valence band and two conduction bands for each
spin degree of freedom) is considered to describe the photoexcitation physics
around the fundamental band gap. We use a large set of momentum states
spanning the full Brillouin zone, thereby including all band-structure valleys
that can be occupied by excited electrons and holes, which directly influence
the optical nonlinearities of the photoexcited carrier density. A second key
ingredient is the realistic description of the long-range Coulomb interaction
of excited carriers to accurately quantify exciton binding energies as well as
band-gap renormalizations. We utilize a lattice Hamiltonian
$H_{\mathbf{k}}^{\alpha\beta}$ formulated in a three-dimensional localized
basis $\big{|}\alpha\big{>}$ consisting of Wannier functions with dominant
d-orbital character from the transition metal (Mo or W). The valence- and
conduction-band splitting caused by spin-orbit interaction is included in the
lattice Hamiltonian along the lines of 43, 11. Diagonalization of the
Hamiltonian yields the band structure $\varepsilon_{\mathbf{k}}^{\lambda}$ and
the Bloch states
$\big{|}\psi_{\mathbf{k}}^{\lambda}\big{>}=\sum_{\alpha}c^{\lambda}_{\alpha,\mathbf{k}}\big{|}\mathbf{k},\alpha\big{>}$,
where the coefficients $c^{\lambda}_{\alpha,\mathbf{k}}$ describe the
momentum-dependent contribution of the orbital $\alpha$ to the Bloch band
$\lambda$. Further technical details on the approach are given in 11. Using
the so-called Peierls approach44, 11, dipole matrix elements can be calculated
directly from the lattice Hamiltonian:

$\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k}}=\frac{e}{i}\frac{1}{\varepsilon^{\lambda}_{\mathbf{k}}-\varepsilon^{\lambda^{\prime}}_{\mathbf{k}}}\sum_{\alpha\beta}(c^{\lambda}_{\alpha,\mathbf{k}})^{*}c^{\lambda^{\prime}}_{\beta,\mathbf{k}}\nabla_{\mathbf{k}}H^{\alpha\beta}_{\mathbf{k}}.$
(5)

Coulomb matrix elements for TMDC monolayers embedded in a dielectric
environment are obtained from bare Coulomb matrix elements and the RPA
dielectric function of a TMDC monolayer calculated in the localized basis
$\big{|}\alpha\big{>}$.16 The bare matrix elements and the dielectric function
$\varepsilon_{\textrm{b}}$ including environmental screening effects are
parametrized as a function of $|\mathbf{q}|$ using the Wannier function
continuum electrostatics approach 45. This combines a macroscopic
electrostatic model for the screening by the dielectric environment in a
heterostructure with a localized description of Coulomb interaction and yields
background-screened matrix elements $V^{\alpha\beta}_{|\mathbf{q}|}$. The
actual parametrization is provided in Ref. 16. We obtain Coulomb matrix
elements in the Bloch-state representation by a unitary transformation using
the coefficients $c^{\lambda}_{\alpha,\mathbf{k}}$:

$\begin{split}V^{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}_{\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}\mathbf{k}_{4}}&=\\\
&\\!\\!\\!\\!\\!\\!\\!\\!\sum_{\alpha,\beta}\Big{(}c_{\alpha,\mathbf{k}_{1}}^{\lambda_{1}}\Big{)}^{*}\Big{(}c_{\beta,\mathbf{k}_{2}}^{\lambda_{2}}\Big{)}^{*}c_{\beta,\mathbf{k}_{3}}^{\lambda_{3}}c_{\alpha,\mathbf{k}_{4}}^{\lambda_{4}}V^{\alpha\beta}_{|\mathbf{k}_{1}-\mathbf{k}_{4}|}\,.\end{split}$
(6)

The encapsulation of TMDC monolayers in hexagonal boron nitride (hBN) is
frequently used to reduce inhomogeneous contributions to the linewidth,
stemming from surface wrinkling or doping. 41 For such a situation we assume a
dielectric environment for all four investigated TMDCs (MoS2, MoSe2, WS2,
WSe2) with a dielectric constant of $\varepsilon_{r}=4.5$.46 In addition a
narrow gap of $0.3$ nm between the monolayer and the surrounding hBN layers
has been taken into account.8 Due to the dielectric screening induced by hBN
encapsulation, the single-particle band gap of the unexcited monolayer
exhibits a shrinkage, that can be described in the $G\Delta W$-formalism as
shown in Refs. 8, 47, 48. In the limit of static screening, the
renormalization of the electron and hole bands $\lambda$ can be expressed as:

$\displaystyle\Sigma^{\lambda,\textrm{G}\Delta\textrm{W}}_{\mathbf{k}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\frac{1}{\mathcal{A}}\sum\limits_{\mathbf{k}^{\prime}}\Delta
V^{\textrm{G}\Delta\textrm{W},\lambda\lambda\lambda\lambda}_{\mathbf{k}\mathbf{k}^{\prime}\mathbf{k}\mathbf{k}^{\prime}}\,,$
(7)

where $\Delta V=V^{\textrm{HS}}-V^{\textrm{free}}$ is the difference of the
Coulomb interaction macroscopically screened by the heterostructure and the
freestanding monolayer.8 The band structure of a TMDC monolayer in a
dielectric environment is therefore given by

$\displaystyle\varepsilon_{\mathbf{k}}^{\lambda,\textrm{G}\Delta\textrm{W}}=\varepsilon_{\mathbf{k}}^{\lambda}+\Sigma^{\lambda,\textrm{G}\Delta\textrm{W}}_{\mathbf{k}}\,.$
(8)

As an example, we find a band-gap shrinkage of $360$ meV for an hBN/MoS2/hBN
heterostructure.
Our description focusses on photoexcitation of unbound electrons and holes on
short time scales before a significant amount of excitons is formed, which
involves slower acoustic-phonon-assisted relaxation processes. On these
grounds, screening due to excited electrons and holes on top of the background
screening is calculated using the plasma dielectric function
$\varepsilon^{-1}_{\textrm{exc}}(|\mathbf{q}|)$. In static approximation and
in the long-wavelength limit,36 the dielectric function assumes the well-known
form
$\varepsilon^{-1}_{\textrm{exc}}(|\mathbf{q}|)=1+\kappa_{|\mathbf{q}|}/|\mathbf{q}|$.
Under photoexcitation, the band structure experiences momentum-dependent
renormalizations $\Sigma_{\mathbf{k}}^{\lambda,\textrm{exc}}$ due to Coulomb
interaction between photoexcited carriers:

$\displaystyle\tilde{\varepsilon}_{\mathbf{k}}^{\lambda}=\varepsilon_{\mathbf{k}}^{\lambda,\textrm{G}\Delta\textrm{W}}+\Sigma_{\mathbf{k}}^{\lambda,\textrm{exc}}\,.$
(9)

With increasing excitation density, a pronounced shrinkage of the band gap is
obtained in connection with the strong Coulomb interaction in TMDCs.49, 50,
51, 52, 29, 51, 53, 54, 55 Recently, it has been predicted that different
shifts for the various band-structure valleys can lead to a transition from
direct to indirect band gaps.36 In general, band-structure renormalizations
are composed of Hartree-Fock contributions and correlation terms. Owing to the
complexity of the used description in the full Brillouin zone, we rely on the
screened exchange Coulomb hole (SXCH) approximation that uses a quasi-static
approximation to correlation terms on a GW-level 36. The corresponding self-
energy of excited carriers is given by

$\displaystyle\Sigma_{\mathbf{k}}^{\lambda,\textrm{exc}}$ $\displaystyle=$
$\displaystyle\Sigma_{\mathbf{k}}^{\lambda,\textrm{H}}+\Sigma_{\mathbf{k}}^{\lambda,\textrm{U}}+\Sigma_{\mathbf{k}}^{\lambda,\textrm{SX}}+\Sigma_{\mathbf{k}}^{\lambda,\textrm{CH}}\,,$
(10)

where $\Sigma_{\mathbf{k}}^{\lambda,\textrm{H}}$ is the renormalization due to
the Hartree interaction, $\Sigma_{\mathbf{k}}^{\lambda,\textrm{U}}$ is the
unscreened electron-hole exchange and
$\Sigma_{\mathbf{k}}^{\lambda,\textrm{SX}}+\Sigma_{\mathbf{k}}^{\lambda,\textrm{CH}}$
constitutes the SXCH contribution. As discussed in Ref. 36,
$\Sigma_{\mathbf{k}}^{\lambda,\textrm{U}}$ and
$\Sigma_{\mathbf{k}}^{\lambda,\textrm{SX}}+\Sigma_{\mathbf{k}}^{\lambda,\textrm{CH}}$
induce the direct-indirect transition of the band gap. Since we take into
account two conduction bands per spin, the self-energy (10) includes Coulomb
interaction between electrons in different conduction bands $e$ and
$e^{\prime}$. Special care has to be taken when going beyond the Hartree-Fock
level for these band combinations. In the SXCH scheme, the Hartree-Fock self-
energy is augmented by including correlations on the GW level in static
approximation. It is discussed in 36 that inter-band correlation terms on the
GW level are sensitive to dynamical screening at roughly the energy difference
between the involved bands. Since photoexcited electrons reside at the K and
$\Sigma$ points, where the energy difference between the conduction bands is
large, relevant GW correlation terms would be sensitive to plasma screening at
high frequencies on the order of $1$ eV. We assume that plasma screening,
which is mainly caused by intra-band polarization, is weak at such
frequencies, so that the inter-band correlations can be neglected. Hence self-
energy terms connecting different conduction bands are treated on the Hartree-
Fock level.

###  4.3 Influence of excitation-induced dephasing

The optical absorption spectrum is directly influenced by the spectral width
of excitonic and band-to-band transitions. The spectral HWHM is determined by
the dephasing $\Gamma$ of inter-band polarizations entering the SBE (4). Due
to the intrinsic connection between dephasing and carrier-carrier scattering,
the dephasing rate is in general power-dependent and can therefore contribute
to a nonlinear dependence of excited carrier density on pump fluence.

Figure 4: Calculated carrier-density dependence of the excitation-induced
dephasing for MoS2, MoSe2, WS2, and WSe2 obtained from a model that is
described in detail in the text.

For the low-excitation regime (small pump fluence), the dominant contribution
to dephasing stems from the coupling of charge carriers to phonons. This
effect strongly depends on the lattice temperature and on the TMDC material.
It has been quantified for the A exciton in Ref. 41 and can be used to fix the
value of the excitation-density dependent dephasing $\Gamma(n)$ at $n=0$. To
model the power-dependent contribution, we refer to Ref.31, where the increase
of the A-exciton linewidth in WS2 has been extracted from a pump-probe
experiment for carrier densities up to several $10^{12}$ cm-2. We assume that
the Mott density, where excitons are fully dissociated into unbound carriers,
marks a turning point in the density dependence of excitation-induced
dephasing. Moreover, it can be expected that the dephasing rate saturates for
large densities due to phase space filling in the carrier-carrier scattering
rates. These effects are phenomenologically captured by the expression

$\displaystyle\Gamma(n)$ $\displaystyle=$
$\displaystyle\left[\arctan\left(\frac{n-n_{\textrm{Mott}}}{\eta}-1\right)+\frac{\pi}{2}\right]$
(11) $\displaystyle\times\frac{\gamma_{0}(f-1)}{\pi}+\gamma_{0}\,,$

where $\gamma_{0}$ is taken from 41 to be $22.5$ meV for MoS2, $17.5$ meV for
MoSe2, $12$ meV for WS2 and $16.5$ meV for WSe2. The Mott density is
approximated as $n_{\textrm{Mott}}=1\times 10^{13}$cm-216. The coefficients
$\eta=1\times 10^{13}$cm-2 and $f=2.5$ can be estimated using the data from
Ref. 31. We assume $n_{\textrm{Mott}}$, $\eta$ and $f$ to be the same for all
TMDC materials. For MoS2, this yields a maximum dephasing rate of $56$ meV.
The results for all considered materials are shown in Fig. 4.

###  4.4 Pump pulse characteristics

The fluence of the pump pulse propagating in a dielectric material with
refractive index $n$ corresponds to the transmitted electromagnetic energy per
area given by

$\begin{split}F=\int dt\left|\mathbf{S}(t)\right|=\varepsilon_{0}c_{0}n\int
dt\left|\mathbf{E}_{0}(t)\right|^{2}\,,\end{split}$ (12)

with S denoting the Poynting vector. Here, $\mathbf{E}_{0}(t)$ is the pump
laser field inside the hBN encapsulation layer with refractive index
$n=\sqrt{\varepsilon_{r}}$. In specific experimental situations, the electric
field inside the sample can be estimated from the field outside the sample,
using a transfer matrix approach. We describe the laser field
$\mathbf{E}_{0}(t)$ by a Gaussian pulse envelope with $150$ fs FWHM duration,
circular light polarization, and incidence perpendicular to the TMD layer. We
assume that the field intensity is given by the field intensity outside the
sample, reduced by the reflectance $R=\left|\frac{1-n}{1+n}\right|^{2}$ of the
outer boundary of the hBN capping layer. Throughout the paper, we refer to the
fluence outside the sample that can be directly compared to experimental
parameters. This work has been supported by the Deutsche
Forschungsgemeinschaft through the graduate school Quantum Mechanic Materials
Modelling (QM3, RTG 2247) and project JA 619/18-1. We acknowledge a grant for
CPU time at the HLRN (Göttingen/Berlin).

## References

  * Xu et al. 2014 Xu, X.; Yao, W.; Xiao, D.; Heinz, T. F. Spin and pseudospins in layered transition metal dichalcogenides. _Nature Physics_ 2014, _10_ , 343–350
  * Ugeda et al. 2014 Ugeda, M. M.; Bradley, A. J.; Shi, S.-F.; Jornada, F. H. d.; Zhang, Y.; Qiu, D. Y.; Ruan, W.; Mo, S.-K.; Hussain, Z.; Shen, Z.-X.; Wang, F.; Louie, S. G.; Crommie, M. F. Giant bandgap renormalization and excitonic effects in a monolayer transition metal dichalcogenide semiconductor. _Nature Materials_ 2014, _13_ , 1091
  * Qiu et al. 2013 Qiu, D. Y.; da Jornada, F. H.; Louie, S. G. Optical Spectrum of ${\mathrm{MoS}}_{2}$: Many-Body Effects and Diversity of Exciton States. _Physical Review Letters_ 2013, _111_ , 216805
  * Chernikov et al. 2014 Chernikov, A.; Berkelbach, T. C.; Hill, H. M.; Rigosi, A.; Li, Y.; Aslan, O. B.; Reichman, D. R.; Hybertsen, M. S.; Heinz, T. F. Exciton Binding Energy and Nonhydrogenic Rydberg Series in Monolayer ${\mathrm{WS}}_{2}$. _Physical Review Letters_ 2014, _113_ , 076802
  * Ye et al. 2014 Ye, Z.; Cao, T.; O’Brien, K.; Zhu, H.; Yin, X.; Wang, Y.; Louie, S. G.; Zhang, X. Probing excitonic dark states in single-layer tungsten disulphide. _Nature_ 2014, _513_ , 214–218
  * Tang et al. 2019 Tang, Y.; Mak, K. F.; Shan, J. Long valley lifetime of dark excitons in single-layer WSe 2. _Nature Communications_ 2019, _10_ , 1–7
  * Plechinger et al. 2016 Plechinger, G.; Nagler, P.; Arora, A.; Schmidt, R.; Chernikov, A.; Águila, A. G. d.; Christianen, P. C. M.; Bratschitsch, R.; Schüller, C.; Korn, T. Trion fine structure and coupled spin–valley dynamics in monolayer tungsten disulfide. _Nature Communications_ 2016, _7_ , 1–9
  * Florian et al. 2018 Florian, M.; Hartmann, M.; Steinhoff, A.; Klein, J.; Holleitner, A. W.; Finley, J. J.; Wehling, T. O.; Kaniber, M.; Gies, C. The Dielectric Impact of Layer Distances on Exciton and Trion Binding Energies in van der Waals Heterostructures. _Nano Letters_ 2018, _18_ , 2725–2732
  * Hao et al. 2017 Hao, K.; Specht, J. F.; Nagler, P.; Xu, L.; Tran, K.; Singh, A.; Dass, C. K.; Schüller, C.; Korn, T.; Richter, M.; Knorr, A.; Li, X.; Moody, G. Neutral and charged inter-valley biexcitons in monolayer MoSe 2. _Nature Communications_ 2017, _8_ , 1–7
  * Steinhoff et al. 2018 Steinhoff, A.; Florian, M.; Singh, A.; Tran, K.; Kolarczik, M.; Helmrich, S.; Achtstein, A. W.; Woggon, U.; Owschimikow, N.; Jahnke, F.; Li, X. Biexciton fine structure in monolayer transition metal dichalcogenides. _Nature Physics_ 2018, _14_ , 1199–1204
  * Steinhoff et al. 2014 Steinhoff, A.; Rösner, M.; Jahnke, F.; Wehling, T. O.; Gies, C. Influence of Excited Carriers on the Optical and Electronic Properties of MoS2. _Nano Letters_ 2014, _14_ , 3743–3748
  * Meckbach et al. 2018 Meckbach, L.; Stroucken, T.; Koch, S. W. Giant excitation induced bandgap renormalization in TMDC monolayers. _Applied Physics Letters_ 2018, _112_ , 061104
  * Kaasbjerg et al. 2014 Kaasbjerg, K.; Bhargavi, K. S.; Kubakaddi, S. S. Hot-electron cooling by acoustic and optical phonons in monolayers of ${\mathrm{MoS}}_{2}$ and other transition-metal dichalcogenides. _Phys. Rev. B_ 2014, _90_ , 165436
  * Steinhoff et al. 2016 Steinhoff, A.; Florian, M.; Rösner, M.; Lorke, M.; Wehling, T. O.; Gies, C.; Jahnke, F. Nonequilibrium carrier dynamics in transition metal dichalcogenide semiconductors. _2D Materials_ 2016, _3_ , 031006
  * Danovich et al. 2017 Danovich, M.; Aleiner, I. L.; Drummond, N. D.; Fal’ko, V. I. Fast Relaxation of Photo-Excited Carriers in 2-D Transition Metal Dichalcogenides. _IEEE Journal of Selected Topics in Quantum Electronics_ 2017, _23_ , 168–172
  * Steinhoff et al. 2017 Steinhoff, A.; Florian, M.; Rösner, M.; Schönhoff, G.; Wehling, T. O.; Jahnke, F. Exciton fission in monolayer transition metal dichalcogenide semiconductors. _Nature Communications_ 2017, _8_ , 1166
  * Pospischil et al. 2014 Pospischil, A.; Furchi, M. M.; Mueller, T. Solar-energy conversion and light emission in an atomic monolayer p–n diode. _Nature Nanotechnology_ 2014, _9_ , 257–261
  * Baugher et al. 2014 Baugher, B. W. H.; Churchill, H. O. H.; Yang, Y.; Jarillo-Herrero, P. Optoelectronic devices based on electrically tunable p–n diodes in a monolayer dichalcogenide. _Nature Nanotechnology_ 2014, _9_ , 262–267
  * Ross et al. 2014 Ross, J. S.; Klement, P.; Jones, A. M.; Ghimire, N. J.; Yan, J.; Mandrus, D. G.; Taniguchi, T.; Watanabe, K.; Kitamura, K.; Yao, W.; Cobden, D. H.; Xu, X. Electrically tunable excitonic light-emitting diodes based on monolayer WSe2 p–n junctions. _Nature Nanotechnology_ 2014, _9_ , 268–272
  * Withers et al. 2015 Withers, F.; Pozo-Zamudio, O. D.; Mishchenko, A.; Rooney, A. P.; Gholinia, A.; Watanabe, K.; Taniguchi, T.; Haigh, S. J.; Geim, A. K.; Tartakovskii, A. I.; Novoselov, K. S. Light-emitting diodes by band-structure engineering in van der Waals heterostructures. _Nature Materials_ 2015, _14_ , 301–306
  * Wu et al. 2015 Wu, S.; Buckley, S.; Schaibley, J. R.; Feng, L.; Yan, J.; Mandrus, D. G.; Hatami, F.; Yao, W.; Vučković, J.; Majumdar, A.; Xu, X. Monolayer semiconductor nanocavity lasers with ultralow thresholds. _Nature_ 2015, _520_ , 69–72
  * Ye et al. 2015 Ye, Y.; Wong, Z. J.; Lu, X.; Ni, X.; Zhu, H.; Chen, X.; Wang, Y.; Zhang, X. Monolayer excitonic laser. _Nature Photonics_ 2015, _9_ , 733–737
  * Salehzadeh et al. 2015 Salehzadeh, O.; Djavid, M.; Tran, N. H.; Shih, I.; Mi, Z. Optically Pumped Two-Dimensional MoS2 Lasers Operating at Room-Temperature. _Nano Letters_ 2015, _15_ , 5302–5306
  * Li et al. 2017 Li, Y.; Zhang, J.; Huang, D.; Sun, H.; Fan, F.; Feng, J.; Wang, Z.; Ning, C. Z. Room-temperature continuous-wave lasing from monolayer molybdenum ditelluride integrated with a silicon nanobeam cavity. _Nature Nanotechnology_ 2017, _12_ , 987–992
  * Lohof et al. 2019 Lohof, F.; Steinhoff, A.; Florian, M.; Lorke, M.; Erben, D.; Jahnke, F.; Gies, C. Prospects and Limitations of Transition Metal Dichalcogenide Laser Gain Materials. _Nano Letters_ 2019, _19_ , 210–217
  * He et al. 2014 He, K.; Kumar, N.; Zhao, L.; Wang, Z.; Mak, K. F.; Zhao, H.; Shan, J. Tightly Bound Excitons in Monolayer ${\mathrm{WSe}}_{2}$. _Physical Review Letters_ 2014, _113_ , 026803
  * Steinhoff et al. 2015 Steinhoff, A.; Kim, J.-H.; Jahnke, F.; Rösner, M.; Kim, D.-S.; Lee, C.; Han, G. H.; Jeong, M. S.; Wehling, T. O.; Gies, C. Efficient Excitonic Photoluminescence in Direct and Indirect Band Gap Monolayer MoS2. _Nano Letters_ 2015, _15_ , 6841–6847
  * Klots et al. 2014 Klots, A. R.; Newaz, A. K. M.; Wang, B.; Prasai, D.; Krzyzanowska, H.; Lin, J.; Caudel, D.; Ghimire, N. J.; Yan, J.; Ivanov, B. L.; Velizhanin, K. A.; Burger, A.; Mandrus, D. G.; Tolk, N. H.; Pantelides, S. T.; Bolotin, K. I. Probing excitonic states in suspended two-dimensional semiconductors by photocurrent spectroscopy. _Scientific Reports_ 2014, _4_ , 6608
  * Chernikov et al. 2015 Chernikov, A.; Ruppert, C.; Hill, H. M.; Rigosi, A. F.; Heinz, T. F. Population inversion and giant bandgap renormalization in atomically thin WS2 layers. _Nature Photonics_ 2015, _9_ , 466
  * Wang et al. 2019 Wang, J.; Ardelean, J.; Bai, Y.; Steinhoff, A.; Florian, M.; Jahnke, F.; Xu, X.; Kira, M.; Hone, J.; Zhu, X.-Y. Optical generation of high carrier densities in 2D semiconductor heterobilayers. _Science Advances_ 2019, _5_ , eaax0145
  * Sie et al. 2017 Sie, E. J.; Steinhoff, A.; Gies, C.; Lui, C. H.; Ma, Q.; Rösner, M.; Schönhoff, G.; Jahnke, F.; Wehling, T. O.; Lee, Y.-H.; Kong, J.; Jarillo-Herrero, P.; Gedik, N. Observation of Exciton Redshift–Blueshift Crossover in Monolayer WS2. _Nano Letters_ 2017, _17_ , 4210–4216
  * Bataller et al. 2019 Bataller, A. W.; Younts, R. A.; Rustagi, A.; Yu, Y.; Ardekani, H.; Kemper, A.; Cao, L.; Gundogdu, K. Dense Electron–Hole Plasma Formation and Ultralong Charge Lifetime in Monolayer MoS2 via Material Tuning. _Nano Letters_ 2019, _19_ , 1104–1111
  * Li et al. 2019 Li, Y.; Shi, J.; Chen, H.; Mi, Y.; Du, W.; Sui, X.; Jiang, C.; Liu, W.; Xu, H.; Liu, X. Slow Cooling of High-Energy C Excitons Is Limited by Intervalley-Transfer in Monolayer MoS2. _Laser & Photonics Reviews_ 2019, _13_ , 1800270
  * Haug and Koch 2004 Haug, H.; Koch, S. _Quantum theory of the optical and electronic properties of semiconductors_ ; World Scientific, 2004
  * Jahnke et al. 1997 Jahnke, F.; Kira, M.; Koch, S. Linear and nonlinear optical properties of excitons in semiconductor quantum wells and microcavities. _Z. Phys. B_ 1997, _104_ , 559–572
  * Erben et al. 2018 Erben, D.; Steinhoff, A.; Gies, C.; Schönhoff, G.; Wehling, T. O.; Jahnke, F. Excitation-induced transition to indirect band gaps in atomically thin transition-metal dichalcogenide semiconductors. _Physical Review B_ 2018, _98_ , 035434
  * Fang et al. 2019 Fang, H. H.; Han, B.; Robert, C.; Semina, M. A.; Lagarde, D.; Courtade, E.; Taniguchi, T.; Watanabe, K.; Amand, T.; Urbaszek, B.; Glazov, M. M.; Marie, X. Control of the Exciton Radiative Lifetime in van der Waals Heterostructures. 2019,
  * Palummo et al. 2015 Palummo, M.; Bernardi, M.; Grossman, J. C. Exciton Radiative Lifetimes in Two-Dimensional Transition Metal Dichalcogenides. _Nano Letters_ 2015, _15_ , 2794–2800
  * Bertoni et al. 2016 Bertoni, R.; Nicholson, C.; Waldecker, L.; Hübener, H.; Monney, C.; De Giovannini, U.; Puppin, M.; Hoesch, M.; Springate, E.; Chapman, R.; Cacho, C.; Wolf, M.; Rubio, A.; Ernstorfer, R. Generation and Evolution of Spin-, Valley-, and Layer-Polarized Excited Carriers in Inversion-Symmetric ${\mathrm{WSe}}_{2}$. _Physical Review Letters_ 2016, _117_ , 277201
  * Selig et al. 2016 Selig, M.; Berghäuser, G.; Raja, A.; Nagler, P.; Schüller, C.; Heinz, T. F.; Korn, T.; Chernikov, A.; Malic, E.; Knorr, A. Excitonic linewidth and coherence lifetime in monolayer transition metal dichalcogenides. _Nature Communications_ 2016, _7_ , 1–6
  * Cadiz et al. 2017 Cadiz, F.; Courtade, E.; Robert, C.; Wang, G.; Shen, Y.; Cai, H.; Taniguchi, T.; Watanabe, K.; Carrere, H.; Lagarde, D.; Manca, M.; Amand, T.; Renucci, P.; Tongay, S.; Marie, X.; Urbaszek, B. Excitonic Linewidth Approaching the Homogeneous Limit in ${\mathrm{MoS}}_{2}$-Based van der Waals Heterostructures. _Physical Review X_ 2017, _7_ , 021026
  * Molina-Sánchez et al. 2016 Molina-Sánchez, A.; Palummo, M.; Marini, A.; Wirtz, L. Temperature-dependent excitonic effects in the optical properties of single-layer ${\mathrm{MoS}}_{2}$. _Physical Review B_ 2016, _93_ , 155435
  * Liu et al. 2013 Liu, G.-B.; Shan, W.-Y.; Yao, Y.; Yao, W.; Xiao, D. Three-band tight-binding model for monolayers of group-VIB transition metal dichalcogenides. _Phys. Rev. B_ 2013, _88_ , 085433
  * Tomczak and Biermann 2009 Tomczak, J. M.; Biermann, S. Optical properties of correlated materials: Generalized Peierls approach and its application to ${\text{VO}}_{2}$. _Phys. Rev. B_ 2009, _80_ , 085117
  * Rösner et al. 2015 Rösner, M.; Şaşıoğlu, E.; Friedrich, C.; Blügel, S.; Wehling, T. O. Wannier function approach to realistic Coulomb interactions in layered materials and heterostructures. _Physical Review B_ 2015, _92_ , 085102
  * Geick et al. 1966 Geick, R.; Perry, C. H.; Rupprecht, G. Normal Modes in Hexagonal Boron Nitride. _Phys. Rev._ 1966, _146_ , 543–547
  * Rohlfing 2010 Rohlfing, M. Electronic excitations from a perturbative $\text{LDA}+GdW$ approach. _Physical Review B_ 2010, _82_ , 205127
  * Winther and Thygesen 2017 Winther, K. T.; Thygesen, K. S. Band structure engineering in van der Waals heterostructures via dielectric screening: the G$\Delta$W method. _2D Materials_ 2017, _4_ , 025059
  * Ulstrup et al. 2016 Ulstrup, S.; Čabo, A. G.; Miwa, J. A.; Riley, J. M.; Grønborg, S. S.; Johannsen, J. C.; Cacho, C.; Alexander, O.; Chapman, R. T.; Springate, E.; Bianchi, M.; Dendzik, M.; Lauritsen, J. V.; King, P. D. C.; Hofmann, P. Ultrafast Band Structure Control of a Two-Dimensional Heterostructure. _ACS Nano_ 2016, _10_ , 6315–6322
  * Gao and Yang 2017 Gao, S.; Yang, L. Renormalization of the quasiparticle band gap in doped two-dimensional materials from many-body calculations. _Physical Review B_ 2017, _96_ , 155410
  * Liu et al. 2019 Liu, F.; Ziffer, M. E.; Hansen, K. R.; Wang, J.; Zhu, X. Direct Determination of Band-Gap Renormalization in the Photoexcited Monolayer ${\mathrm{MoS}}_{2}$. _Physical Review Letters_ 2019, _122_ , 246803
  * Cunningham et al. 2017 Cunningham, P. D.; Hanbicki, A. T.; McCreary, K. M.; Jonker, B. T. Photoinduced Bandgap Renormalization and Exciton Binding Energy Reduction in WS2. _ACS Nano_ 2017, _11_ , 12601–12608
  * Fan et al. 2017 Fan, X.; Zheng, W.; Liu, H.; Zhuang, X.; Fan, P.; Gong, Y.; Li, H.; Wu, X.; Jiang, Y.; Zhu, X.; Zhang, Q.; Zhou, H.; Hu, W.; Wang, X.; Duan, X.; Pan, A. Nonlinear photoluminescence in monolayer WS 2 : parabolic emission and excitation fluence-dependent recombination dynamics. _Nanoscale_ 2017, _9_ , 7235–7241
  * Xie et al. 2019 Xie, Y.; Zhang, S.; Li, Y.; Dong, N.; Zhang, X.; Wang, L.; Liu, W.; Kislyakov, I. M.; Nunzi, J.-M.; Qi, H.; Zhang, L.; Wang, J. Layer-modulated two-photon absorption in MoS2: probing the shift of the excitonic dark state and band-edge. _Photonics Research_ 2019, _7_ , 762–770
  * Kang et al. 2017 Kang, M.; Kim, B.; Ryu, S. H.; Jung, S. W.; Kim, J.; Moreschini, L.; Jozwiak, C.; Rotenberg, E.; Bostwick, A.; Kim, K. S. Universal Mechanism of Band-Gap Engineering in Transition-Metal Dichalcogenides. _Nano Letters_ 2017, _17_ , 1610–1615


Srijan<EMAIL_ADDRESS>Swapnaneel
<EMAIL_ADDRESS>Indian Statistical Institute
Kolkata, India Indian Statistical Institute
Kolkata, India

A Machine Learning approach of Ecological Modeling:
A New method to find Similarity Index

Srijan Chattopadhyay
<EMAIL_ADDRESS>
Swapnaneel Bhattacharyya
<EMAIL_ADDRESS>
Indian Statistical Institute, Kolkata, India

###### Abstract

In many scientific research, it’s often imperative to determine whether pairs
of entities have similarities in themselves or not. There are standard
approaches to this problem such as Jaccard${}^{\href
https://doi.org/10.1111/j.1469-8137.1912.tb05611.x}$, Sørensen-Dice[2], and
Simpson${}^{\href https://www.nature.com/articles/163688a0}$. Recently, a
better index for analysis of Co-Occurrence and similarity $(\hat{\alpha})$ was
developed[4] and it reversed all the results obtained by Standard Indices and
supported theoretical predictions. In this paper, we propose a new measure of
similarity using Maximum Likelihood Estimation${}^{\href
https://projecteuclid.org/journals/statistical-science/volume-12/issue-3/RA-
Fisher-and-the-making-of-maximum-
likelihood-1912-1922/10.1214/ss/1030037906.full}$, Principal Component
Analysis${}^{\href https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4792409/\href
https://stats.oarc.ucla.edu/sas/output/principal-components-analysis/}$,
Linear Discriminant Analysis${}^{\href
https://www.sciencedirect.com/science/article/pii/S0167865504001825?via3Dihub}$,
and clustering. Our index depends strongly on the data before introducing
randomness in prevalence. Then we proposed a new method of randomization,
which changed the whole pattern of the results. Before randomization, it was
strongly dependent on the prevalance and hence was following the pattern of
the Jaccard index. So, we introduced the new randomization technique and hence
the whole results reversed and followed that of $\hat{\alpha}$. Also, we will
show some limitations of $\hat{\alpha}$ which we try to resolve through
different pathways.

##  1 Introduction

Similarity or Co-occurrence analysis has been an interesting area of research
since the 18th century${}^{\href https://www.worldcat.org/title/christiani-
hugenii-libellus-de-ratiociniis-in-ludo-ale-or-the-value-of-all-chances-in-
games-of-fortune-cards-dice-wagers-lotteries-c-mathematically-
demonstrated/oclc/1006229641?referer=di&ht=edition}$ as it plays an important
role in many fields like Ecology[10], Medical Science${}^{\href
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4135978/}$,
Biodiversity${}^{\href
https://www.science.org/doi/10.1126/science.1248484\href
https://www.hindawi.com/journals/ijfr/2020/4865845/}$, Genetics${}^{\href
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3569618/}$, Bio-
geography${}^{\href https://www.nature.com/articles/35016066\href
https://pubmed.ncbi.nlm.nih.gov/31624208/}$, Epidemiology${}^{\href
https://projecteuclid.org/Proceedings/berkeley-symposium-on-mathematical-
statistics-and-probability/proceedings-of-the-third-berkeley-symposium-on-
mathematical-statistics-and/toc/bsmsp/1200501640}$, Neuroscience${}^{\href
https://pubmed.ncbi.nlm.nih.gov/23798414/}$,
Evolution${}^{\href https://www.nature.com/articles/nature13827}$,
Economics${}^{\href https://www.jstor.org/stable/2232173}$, Geology${}^{\href
https://www.scirp.org/journal/paperinformation.aspx?paperid=71033}$, etc.
Scientists use similarity indices and Co-occurrence Analyses to determine
distributional similarities between two populations in the same time period or
the same population in different time periods. Also in medical sciences or
genetics, it is important to observe whether two drugs or two genes behave
similarly or not.

In the recent past three of the most used similarity indices are Jaccard,
Sorensen-Dice, and Simpson indices.

  * •

Jaccard Index: This index is due to Paul Jaccard which measures the similarity
between finite sample sets and is defined as the size of the intersection
divided by the size of the union of the sample sets, $J(A,B)=\frac{|A\cap
B|}{|A\cup B|}$. So by construction $0\leq J\leq 1.$

  * •

Sørensen-Dice Index: This index is due to Thorvald Sørensen and Lee Raymond
Dice defined as $\frac{2|X\cap Y|}{|X|+|Y|}$ for two given finite sets $X,Y$.

  * •

Simpson Index: Simpson’s Diversity Index is a statistic to measure the
diversity of species in a community. In a population with $N$ many entities if
there are $R$ communities and $n_{i}$ denotes the number of entities belonging
to the ith community then the Simpson Index is defined as
$\frac{\sum_{i=1}^{R}{n_{i}(n_{i}-1)}}{N(N-1)}$ which can be interpreted as if
two entities are randomly selected from the population then the probability of
two entities coming from the same community gives the Simpson Index.

All these indices interpret smaller values (close to 0) as a highly negative
association and those near one to be highly positive.
Recently a new measure of similarity named $(\hat{\alpha})$ was discovered by
Mainali, Slud, Singer, and Fagan. They claimed that the three above-mentioned
metrics and numerous related measures suffer from major Statistical flaws that
not only bias the interpretation but completely invalidate them. There they
provided an affinity model depending on some unknown parameters
(probabilities) $p_{1}$ and $p_{2}$ based on which the proposed theoretical
measure $\alpha=log(\frac{p_{1}}{1-p_{1}}/\frac{p_{2}}{1-p_{2}})$ and the MLE
estimation of $\alpha$ is their proposed similarity index. There the following
two propositions were derived.

  * •

Proposition 1: (null hypothesis distribution). If the numbers $m_{A}$,
$m_{B}$, and $N$ are fixed and the $m_{A}$ balls colored A are placed in boxes
(with at most one ball to a box) in any way whatever, and, given the placement
of the A-colored balls, the $m_{B}$ balls colored B are placed equiprobably in
the N boxes, then the probability distribution of the count X of boxes
containing balls of both colors is hypergeometric $(N,m_{A},m_{B})$

  * •

Proposition 2: (affinity model distribution). Suppose that the numbers
$m_{A},m_{B}$, and $N$ are fixed, and the $m_{A}$ balls colored A are placed
in boxes (with at most one ball in each box) in any way whatever, and given
the placement of the A-colored balls, balls colored B are placed independently
in boxes with probability $p_{1}$ for boxes containing A-colored balls and
with probability $p_{2}$ for boxes not containing A-colored balls. Then,
conditionally given that a total of $m_{B}$ balls are placed in boxes, either
conditionally given the configuration of the A-colored balls or
unconditionally, the distribution of the count X of co-occurrences is extended
hypergeometric with parameters $N,m_{A},m_{B}$, and
$\alpha=log(\frac{p_{1}}{1-p_{1}}/\frac{p_{2}}{1-p_{2}})$

Hence

$P(X=k)=\frac{{m_{A}\choose k}{N-m_{A}\choose m_{B}-k}e^{\alpha
k}}{\sum_{j=s}^{t}{m_{A}\choose j}{N-m_{A}\choose m_{B}-j}e^{\alpha j}}$ (1)

for $s=max(m_{A}+m_{B}-N,0)\leq k\leq min(m_{A},m_{B})=t$. Here $X$ is the
count of co-occurrences. This is Extended Hypergeometric or Non-Central
Hypergeometric Distribution. So, for observed $X=k$ with fixed
$(N,m_{A},m_{B})$ they estimated $\alpha\in(-\infty,\infty)$ by its ML
Estimation $\hat{\alpha}$ and concluded it to be the unique solution of the
equation when a finite solution exists

$X=\frac{\sum_{j=s}^{t}j{m_{A}\choose j}{N-m_{A}\choose m_{B}-j}e^{\alpha
j}}{\sum_{j=s}^{t}{m_{A}\choose j}{N-m_{A}\choose m_{B}-j}e^{\alpha j}}$ (2)

But when $X=t$ they defined $\hat{\alpha}$ to be $-\infty$ and to be $+\infty$
when $X=s$.
The most interesting feature of $\alpha$ was the distribution of it hardly
depends on the prevalence i.e. on the $(N,m_{A},m_{B})$, whereas the most used
three do. Our measure also doesn’t directly depend on the prevalence, but it
depends on the whole data, i.e. for the same prevalence set different data may
result in different values of the index and we will later discuss why it is
like that. And that is why we call it a Joint to Marginal Approach, because,
we first use the whole information at a time and then cut whatever is needed.
The whole process will be explained later in detail. We will also give some
methods to apply when there are missing data in the data set. And,
interestingly, the patterns of our results exactly match with $\hat{\alpha}$
when we apply it to the data. We used $Probit$ and $Logit$ links to make the
data continuous if it is not, but that is totally data specific.

##  2 An Issue of $\hat{\alpha}$

We point out an issue of $\hat{\alpha}$. Firstly the approach to compute
$\hat{\alpha}$ is classical i.e. no prior distribution of $p_{1},p_{2}$ is
assumed. But practically $p_{1},p_{2}$ can’t be assumed to be degenerate.
That’s why one can assume some empirical prior distribution for $p_{1},p_{2}$
and try to obtain the Bayes estimator for $\hat{\alpha}$. For example, we know
for given $p_{1},p_{2}$,

$P(X=k|p_{1},p_{2})=\frac{{m_{A}\choose k}{N-m_{A}\choose
m_{B}-k}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{k}}{\sum_{j=s}^{t}{m_{A}\choose
j}{N-m_{A}\choose m_{B}-j}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{j}}$ (3)

for $s=max(m_{A}+m_{B}-N,0)\leq k\leq min(m_{A},m_{B})=t$. Now, $p_{1}$ is the
probability of species 2 occupying a site already occupied by species 1, and
$p_{2}$ is the probability of species 2 occupying a site already not occupied
by species 1. Since $p_{1},p_{2}$ will depend on time and space, it’s
necessary to assume a prior distribution for them.

##  3 A Possible Bayesian Approach to calculate $\hat{\alpha}$

  * •

Assume the prior distribution of $p_{1}$ to be $TN(\mu_{1},1,0,1)^{\href
https://en.wikipedia.org/wiki/Truncated_normal_distribution}$ and of $p_{2}$
to be $TN(\mu_{2},1,0,1)$. So, by the law of total probability

$P(X=k)=\int_{0}^{1}\int_{0}^{1}{\frac{{m_{A}\choose k}{N-m_{A}\choose
m_{B}-k}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{k}}{\sum_{j=s}^{t}{m_{A}\choose
j}{N-m_{A}\choose
m_{B}-j}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{j}}\frac{e^{\frac{-(p_{1}-\mu_{1})^{2}}{2}}}{\int_{0}^{1}{e^{\frac{-(p_{1}-\mu_{1})^{2}}{2}}}dp_{1}}\frac{e^{\frac{-(p_{2}-\mu_{2})^{2}}{2}}}{\int_{0}^{1}{e^{\frac{-(p_{2}-\mu_{2})^{2}}{2}}}dp_{2}}dp_{1}dp_{2}}$

So, we have to maximize this equation $g(\mu_{1},\mu_{2})$ with respect to
$\mu_{1}$ and $\mu_{2}$, for finding the MLE of $\mu_{1}$ and $\mu_{2}$, which
is barely impossible to solve analytically, We can try to implement an EM
algorithm over here, but that also will be very complicated. We tried to do it
using R using a grid search method, but with high prevalence values, space
complexity becomes drastically high along with a high time complexity. We did
check using GMP package, but due to inner algorithmic issues, the results were
unexpectedly wrong. So, it becomes tremendously challenging to solve using
Truncated Normal prior.

  * •

So, as uniform density is the simplest one, we try to obtain Bayes Estimate of
$\alpha$ assuming Uniform Priors i.e. both $p_{1}$ and $p_{2}$ follow
$Unif(0,\alpha)$ and we will then try to maximize again the likelihood. So,
here,

$P(X=k)=\int_{0}^{\alpha}\int_{0}^{\alpha}{\frac{{m_{A}\choose
k}{N-m_{A}\choose
m_{B}-k}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{k}}{\sum_{j=s}^{t}{m_{A}\choose
j}{N-m_{A}\choose
m_{B}-j}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{j}}dp_{1}dp_{2}}$ (4)

This is a function of $\alpha$,

$f(\alpha)=\frac{1}{\alpha^{2}}\int_{0}^{\alpha}\int_{0}^{\alpha}{g(p_{1},p_{2})dp_{1}dp_{2}}\\\
$ (5)

where, $g(p_{1},p_{2})=\frac{{m_{A}\choose k}{N-m_{A}\choose
m_{B}-k}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{k}}{\sum_{j=s}^{t}{m_{A}\choose
j}{N-m_{A}\choose m_{B}-j}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{j}}$

$f^{\prime}(\alpha)=\frac{-2}{\alpha^{3}}\int_{0}^{\alpha}\int_{0}^{\alpha}{g(p_{1},p_{2})dp_{1}dp_{2}}+\frac{1}{\alpha^{2}}\frac{d}{d\alpha}{\int_{0}^{\alpha}\int_{0}^{\alpha}{g(p_{1},p_{2})dp_{1}dp_{2}}}$
(6)

Now,

$\frac{d}{d\alpha}{\int_{0}^{\alpha}\int_{0}^{\alpha}{g(p_{1},p_{2})dp_{1}dp_{2}}}=\int_{0}^{\alpha}{g(\alpha,p_{2})dp_{2}}+\int_{0}^{\alpha}{g(p_{1},\alpha)dp_{1}}$
(7)

Hence, $f^{\prime}(\alpha)$ becomes:

$\displaystyle
f^{\prime}(\alpha)=\frac{1}{\alpha^{2}}[\int_{0}^{\alpha}{\frac{{m_{A}\choose
k}{N-m_{A}\choose
m_{B}-k}(\frac{\alpha(1-p_{2})}{p_{2}(1-\alpha)})^{k}}{\sum_{j=s}^{t}{m_{A}\choose
j}{N-m_{A}\choose
m_{B}-j}(\frac{\alpha(1-p_{2})}{p_{2}(1-\alpha)})^{j}}}+\frac{{m_{A}\choose
k}{N-m_{A}\choose
m_{B}-k}(\frac{p_{1}(1-\alpha)}{\alpha(1-p_{1})})^{k}}{\sum_{j=s}^{t}{m_{A}\choose
j}{N-m_{A}\choose m_{B}-j}(\frac{p_{1}(1-\alpha)}{\alpha(1-p_{1})})^{j}}]$
$\displaystyle-\frac{2}{\alpha^{3}}\int_{0}^{\alpha}\int_{0}^{\alpha}{\frac{{m_{A}\choose
k}{N-m_{A}\choose
m_{B}-k}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{k}}{\sum_{j=s}^{t}{m_{A}\choose
j}{N-m_{A}\choose
m_{B}-j}(\frac{p_{1}(1-p_{2})}{p_{2}(1-p_{1})})^{j}}dp_{1}dp_{2}}$

Now $\hat{\alpha}$ is the root of the equation $f^{\prime}(\alpha)=0$. From
the above equations, it’s clear that it’s almost computationally impossible to
find a numerical value of $\hat{\alpha}$ even when we use both prior
distributions to be the same with $Unif(0,\alpha)$. So, considering the need
and difficulty for the assumption of the prior distribution we introduce a new
method to deal with the Similarity or Co-Occurrence Index.

##  4 Our Method

Assumption and Model Description: Let’s say we have data of the species
$S_{1},...,S_{k}$, and the data was collected from the islands
$I_{1},...,I_{n}$ over the time periods $t_{1},...,t_{l}$. Now our assumption
is For a fixed time period the species $S_{i}$ can occur in any of the $n$
islands independently with the same probability i.e.

  1. 1.

$P(S_{i}$ is found in $I_{j}$ in time period $t_{r}$) = $p_{i,r}$

So each $S_{i}$ in $r$ th time period follows $Bernoulli(p_{i,r})$ So each
entry of the data is either 1 (if present ) or 0. So first we have to obtain
the estimate $\hat{p}_{i,r}$.
Necessary Transformtions: For that, we first do some transformations. For each
time period, if, $S_{i,j}<n-1$, we add 1 to it, and for the other cases, we
subtract 1 from it. Now the most convenient choice is the MLE based on this
transformed data we have. And the MLE of the Bernoulli parameter is the sample
mean and by SLLN it’s a consistent Estimator of the population mean as well
which is $p_{i,r}$. So we have got the estimated distributions. Now since the
matrix is discrete in nature we make necessary monotone transformations to get
some information from the matrix which are continuous but contain most of the
information of the actual data. The two possible transformations are

$T^{r}_{i,j}=\begin{cases}log(\frac{\hat{p}_{i,r}}{1-\hat{p}_{i,r}}),&\text{if
$S^{r}_{i,j}=1,\hat{p}_{i,r}\geq\frac{1}{2}$}\\\
log(\frac{1-\hat{p}_{i,r}}{\hat{p}_{i,r}}),&\text{if
$S^{r}_{i,j}=0,\hat{p}_{i,r}\geq\frac{1}{2}$}\\\
log(\frac{\hat{p}_{i,r}}{1-\hat{p}_{i,r}}),&\text{if
$S^{r}_{i,j}=0,\hat{p}_{i,r}<\frac{1}{2}$}\\\
log(\frac{1-\hat{p}_{i,r}}{\hat{p}_{i,r}}),&\text{if
$S^{r}_{i,j}=1,\hat{p}_{i,r}<\frac{1}{2}$}\\\ \end{cases}$

Or,

$T^{r}_{i,j}=\begin{cases}\Phi^{-1}(\hat{p}_{i,r}),&\text{if
$S^{r}_{i,j}=1,\hat{p}_{i,r}\geq\frac{1}{2}$}\\\
\Phi^{-1}(1-\hat{p}_{i,r}),&\text{if
$S^{r}_{i,j}=0,\hat{p}_{i,r}\geq\frac{1}{2}$}\\\
\Phi^{-1}(\hat{p}_{i,r}),&\text{if
$S^{r}_{i,j}=0,\hat{p}_{i,r}<\frac{1}{2}$}\\\
\Phi^{-1}(1-\hat{p}_{i,r}),&\text{if
$S^{r}_{i,j}=1,\hat{p}_{i,r}<\frac{1}{2}$}\\\ \end{cases}$

where $\Phi$ denotes the CDF of $N(0,1)$.
Applying PCA to transformed data and introducing "score" matrix: So
corresponding to each time period we get a matrix of Continuous data entries.
So we have the matrices $T^{1},...,T^{l}$ where $T^{i}$ corresponds to the
time period $t_{i}$.
Now we apply Principal Component Analysis${}^{\href
https://stats.oarc.ucla.edu/sas/output/principal-components-analysis/}$ to the
transformed data. Now we define a new matrix X as

$\displaystyle X$ $\displaystyle=\begin{bmatrix}{T^{1}}^{T}\\\ {T^{2}}^{T}\\\
\vdots\\\ {T^{l}}^{T}\end{bmatrix}$ (8)

So X has $k$ columns and $nl$ many rows. Now we define another matrix $Y$ as
follows

$\displaystyle\begin{bmatrix}Y_{1,i}\\\ Y_{2,i}\\\ \vdots\\\
Y_{nl,i}\end{bmatrix}$ $\displaystyle=\begin{bmatrix}X_{1,i}-\bar{X}_{i}\\\
X_{2,i}-\bar{X}_{i}\\\ \vdots\\\ X_{nl,i}-\bar{X}_{i}\end{bmatrix}$ (9)

where

$\bar{X}_{i}=\frac{X_{1,i}+...+X_{nl,i}}{nl}$

Now we want to find the eigenvectors of $Y^{T}Y$ which of dimension $k\times
k$. Now usually what happens is $nl$ becomes much smaller than $k$. Now we
know that non-zero eigenvalues of $AB$ and $BA$ are the same. Now

$YY^{T}\vec{v}=\lambda\vec{v}$

$Y^{T}Y(Y^{T}\vec{v})=\lambda(Y^{T}\vec{v})$

So if we can find the eigenvectors of $YY^{T}$ then right multiplying by
$Y^{T}$ we can get those of $Y^{T}Y$. Now let the eigenvalues of $YY^{T}$ be

$\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{n}\geq 0$

Now taking only the "significant" eigenvalues ( say $\geq 10^{-5})$ we get $r$
many eigenvalues. Then taking only those eigenvectors corresponding to those
eigenvalues we get a bunch of eigenvectors and adjoining them column-wise into
a matrix we get $P^{\prime}$. Then define $P=Y^{T}P^{\prime}$. Clearly, $P$ is
the set of eigenvectors of $Y^{T}Y$ corresponding to the "significant"
eigenvalues. Then we normalize each column and call the new matrix $Q$. So
clearly we got an Orthonormal Basis for the row space of $Y$. Now we express
the rows of $Y$ in terms of this orthonormal basis, hence we get a "Score"
vector for each of the rows of $Y$ i.e. each of the columns of the actual
data. So the score will be

$\textit{score}=YQ$ (10)

Hence we are now in a space, the dimension of which is less than
$\min\\{k,nl\\}$. so we are done with our PCA.
Applying LDA to score matrix: Now we apply Linear Discriminant Analysis on the
score  matrix. If necessary we can cut the last few columns of score matrix to
reduce again the number of Eigenvectors corresponding to significant
Eigenvalues. say we take the first $s$ many columns of scores and call that
$R$. Then define the matrices $U_{1},...,U_{l}$ such that
$U_{1}=R[1:n,],U_{2}=R[n+1:2n,],..,U_{i}=R[(i-1)n+1:in,]$ etc. Then define
$Y_{1},..,Y_{l}$ in the same fashion we defined $Y$ from $X$, i.e. by
centering the data. Then define $m_{1},.,m_{l}$ such that $m_{i}$ is the
vector consisting of the mean of the elements in each column of $U_{i}$, i.e.
the $jth$ element of $m_{i}$ will be the $\frac{j\text{ th column sum of
$U_{i}$}}{\text{number of rows of $U_{i}$ }}$. Then construct $Z$ such that

$\displaystyle Z$ $\displaystyle=\begin{bmatrix}{U_{1}}\\\ {U_{2}}\\\
\vdots\\\ {U_{l}}\end{bmatrix}$ (11)

Then by C.R.Rao’s formula, if there are $C$ classes and each of the classes,
has to mean $\mu_{i}$ and the same covariance $\Sigma$ then the scatter
between class variability may be defined by the sample covariance of the class
means

$\Sigma_{b}=\frac{1}{C}\sum_{i=1}^{C}(\mu_{i}-\mu)(\mu_{i}-\mu)^{T}$ (12)

where $\mu$ is the mean of the class means. The class separation in a
direction $\vec{x}$ is given by

$S=\frac{\vec{x}^{T}\Sigma_{b}\vec{x}}{\vec{x}^{T}\Sigma\vec{x}}$ (13)

So we have to find the direction which maximizes $S$ because in that direction
the variation between the same class is minimized and variation between
different classes is maximized. So in our case $W=cZ^{T}Z$ where $c$ is the
normalization constant, which represents $\Sigma_{b}$ and $B=d\hat{\Sigma}$
where $\hat{\Sigma}$ is Co-variance matrix of $m_{1},...,m_{l}$ and $d$ is
Normalising Constant. So we have to find the

$\underset{\vec{x}\in\mathbb{R}^{s}\setminus{\vec{0}}}{\text{arg max :
}}\frac{\vec{x}^{T}W\vec{x}}{\vec{x}^{T}B\vec{x}}$

Now we use the following theorem. The proof is given in the Appendix.
Theorem: Let $W$ be Positive Definite.Then

$\underset{\vec{x}\not=\vec{0}}{\text{max}}\frac{\vec{x}^{T}B\vec{x}}{\vec{x}^{T}W\vec{x}}=\lambda_{1}(W^{-1}B)$

where $\lambda_{1}(W^{-1}B)$ is the largest eigenvalue of $W^{-1}B$ and the
maximum is attained at $\vec{x}_{0}$ iff $\vec{x}_{0}$ is an eigenvector of
$W^{-1}B$ corresponding to $\lambda_{1}(W^{-1}B)$.

Now in the practical situation, two things can happen:

  * •

From the proof, the eigenvalues of $B^{-1}W$ must be real. But the algorithmic
approaches done by different languages (e.g. R) don’t always give real
eigenvalues as an output. So in that case following the proof, we do the
Cholesky Decomposition of $B$ and ask to return the corresponding eigenvalue.
So Using the Cholesky Decomposition

$W=M^{T}M$

$(M^{-1})^{T}BM^{-1}\vec{y}_{0}=\lambda_{1}\vec{y}_{0}$

$M^{-1}(M^{-1})^{T}B(M^{-1}\vec{y}_{0})=\lambda_{1}(M^{-1}\vec{y}_{0})$

$W^{-1}B\vec{x}_{0}=\lambda_{1}\vec{x}_{0}$

where $\vec{x}_{0}=M^{-1}\vec{y}_{0}$. Instead, we consider
$\vec{z}_{0}=M\vec{y}_{0}=M^{2}\vec{x}_{0}$ and look at the scores with
respect to $\vec{z}_{0}$.

  * •

Also in real situations $w$ may not always be invertible. Then we follow the
following technique We take some small $\lambda$ greater than 0 and add
$\lambda I$ to $W$ to make it invertible. Note that such a $\lambda$ always
exists since there exists finitely many $\lambda$ for which $(A+\lambda I)$ is
non-invertible as $det(A+\lambda I)$ is a polynomial of $\lambda$.

The Similarity Index: After all these steps, we got a matrix consisting of
eigenvectors of $W^{-1}B$ as columns. Then we want to express $R$ in terms of
$M^{2}S$. So we introduce a new matrix $scores_{new}$ as:

$scores_{new}=R(S^{T}_{1})^{-1}$ (14)

where $S_{1}=M^{2}S$. Hence the first column of $scores_{new}$ will be the
scores with respect to $\vec{z}_{0}$. Call $N=scores_{new}[,1]$. Then $N$ is a
vector of dimension $nl$.
So our similarity index for $i^{\text{ th}}$ island will be

1 - standard deviation of $N_{i},N_{i+n},..,N_{i+n(l-1)}$.

So the possible range of this index is $(-\infty,1]$ and the more it is close
to 1 indicates more association and the more negative indicates less
association or similarity.

##  5 Results and Discussion

Below, we show our results when we apply our non-random algorithm to some
data. Firstly, we have applied to bio-geographical data. Previously, the
Jaccard index and $\hat{\alpha}$ have been applied to it and the results are
also there. It was shown that Jaccard and $\hat{\alpha}$ are behaving
oppositely, and we will show our non-random measure shows the same pattern
with Jaccard, and it happens because we did the whole thing based on whatever
data is given to us up to some parameter estimation. That implies that this
index is highly dependent on the prevalence because of its non-stochastic
behavior and hence behaves like Jaccard.

###  5.1 Biogeography: Time Similarity between different islands

Beta Diversity is a measure that indicates the diversity of ecological
communities across space or compositional change in a single community over
different time periods. So, low beta diversity signifies more similarities or
less significant changes over space or time. On the other hand, more beta
diversity means high values of species turnover (across space) or high levels
of local change (over time). The data${}^{\href
https://www.nature.com/articles/s41598-017-05114-5#:~:text=The%20flora%20of%20the%20Tuscan,indicate%20very%20high%20turnover%20rates.}$
consists of data from 16 islands of the Tuscan Archipelago, Italy in 1830-1950
and 1951-2015. All the latitudes of the islands are in the range 42.25280 and
43.43093 and the longitudes are in the range 9.81988 and 11.31667, so the
islands are in very comparable environments and hence we made an assumption
that, in any two islands, the probability of getting a species is same in the
same time period. Now, we already have the results using Jaccard and
$\hat{\alpha}$ with us. The dataset has a total of 1842 species, each of which
has their data in 16 islands in 2 different aforementioned time periods.
Firstly, we assume that the probability of getting a species on any island is
the same for a fixed time period. So, say if $X_{i,j}$ is the random variable
specifying the existence of the ith species in the jth time period, then
$X_{i,j}$ basically follows $Bernoulli(p_{ij})$, where $p_{ij}$ is the
probability of finding the ith species in the jth time period on any of the
islands. Now, firstly we estimate $p_{ij}$ by its MLE (mean here), calling it
$\hat{p}_{i,j}$. We first change some boundary values and add(or, subtract) 1
to(or, from) each. The whole process is described later. Then on that data, we
did the following. We do the following transformation to make the dataset
continuous from discrete. We did it in two different ways and the results will
be similar for both, hence can use any one of them and the results won’t
differ much. If the ith species is present in the kth island in the jth time
period if $p_{i,j}\geq\frac{1}{2}$, we replace the corresponding entry of the
actual data by $F(\frac{p_{i,j}}{1-p_{i,j}})$ and if $p_{i,j}<\frac{1}{2}$, we
replace by $F(\frac{1-p_{i,j}}{p_{i,j}})$ and it is finite because of the
previous boundary value considerations and the corresponding manipulations.
And if it is not present, we reverse the previous. Here $F$ is any of the
following two $\log$ or $\Phi^{-1}$, i.e. the logarithm(Logit transform) or
the inverse standard normal CDF(Probit transform). So, we have now converted
the Bernoulli categorical variables into a continuous data set. Then we apply
PCA and LDA to the data and we obtain some scores for each island for two
different time periods. For kth island, we define the similarity index between
two time periods to be the standard deviation of $S_{k,2}$ and $S_{k,1}$,
where $S_{k,j}$ is the score obtained in the jth period for the kth island
after the whole process. We call that our proposed index, $\hat{\beta_{L}}$ if
Logit transform was applied and $\hat{\beta_{N}}$ if Probit transform was
applied.
$\hat{\alpha}$ shows that large islands have more stable communities than
small ones, but the Jaccard index shows the opposite result, and so does our
index. Jaccard is positively biased for species pairs with higher prevalence,
because high prevalence allows for higher co-occurrence (and higher J values)
even in null scenarios, same thing happens for our non-random index also, so
they yield the same result. Similarly, we see that as the distance to the
nearest bigger island increases, $\hat{\alpha}$ decreases, Jaccard and
$\hat{\beta}_{N}$, $\hat{\beta}_{L}$ increases. Similarly increase in species
richness yields a decrease in $\hat{\alpha}$, an increase in Jaccard,
$\hat{\beta}_{N}$ and $\hat{\beta}_{L}$

Figure 1: Relation of beta(L) and beta(N) with island area

Figure 2: Relation of beta(L) and beta(N) with distance to nearest bigger
island

Figure 3: Relation of beta(L) and beta(N) with species richness

Also, There is evidence that similar to species richness, beta-diversity also
decreases with increasing elevation${}^{\href
https://www.hindawi.com/journals/ijfr/2020/4865845/}$Kraft et al.${}^{\href
https://www.science.org/doi/10.1126/science.1208584}$ concluded that beta-
diversity decreases with increasing elevation and latitude just as a
consequence of the decreasing size of the species pool. Here we got the
evidence in favor of the results using $\hat{\beta_{L}}$ and
$\hat{\beta_{N}}$.

Figure 4: Relation of beta(L) and beta(N) with elevation

###  5.2 Biogeography: Time Similarity between different islands and Space
Similarity between different birds

Now we will first consider two beta diversities, first over time, then between
different species. We have data from 88 bird species of which each has
data${}^{\href https://www.nature.com/articles/s41559-017-0457-3}$ from 96
islands in three different years 2011, 2012, and 2013. We want first want to
see the time similarities between various species of an island. Then we want
to see the similarities between different bird species time-wise and as a
whole. Then we apply all the processes as before on the data. And, for this
case, we define time similarity as a whole of 3 years as
$sd(S_{k,1},S_{k,2},S_{k,3})$ and those between any two time periods to be the
same definition as before. We found $\hat{\beta_{N}}$ and $\hat{\beta_{L}}$
both and found amazing similarities (with some shift) between their density
based on the kernel density estimation method${}^{\href
https://en.wikipedia.org/wiki/Kernel_density_estimation}$. We also plotted a
scatter plot between these, and as expected, we got a diagonal line. We also
plotted the QQ plot between $\hat{\beta}_{L}$ and $\hat{\beta}_{N}$ and
observed a diagonal trend which indicates matching of their sample quantiles,
hence of the distributions. Now, from the Score plot of three different
periods, we can see that the scores can be distinguished from the
$\hat{\beta}_{L}$, whether it is not from $\hat{\beta}_{N}$. We can see three
clusters clearly in the case of $\hat{\beta}_{L}$. When we considered the
joint similarity over the whole time, the plots of $\hat{\beta}_{L}$ and
$\hat{\beta}_{N}$ look very similar. That implies that the diversity was
almost the same in the years 2011 and 2012.

Figure 5: Density Plot for $\hat{\beta_{L}}$ and $\hat{\beta_{N}}$

Figure 6: Scatter and QQ plot of $\hat{\beta_{L}}$ and $\hat{\beta_{N}}$

Figure 7: Plot of the of $\hat{\beta_{L}}$ and $\hat{\beta_{N}}$ the islands

Figure 8: Scores for L and N for different Time periods and islands

Now let us move on to the analysis of similarity analyses between different
bird species. We projected the whole similarity matrix into several
heatmaps${}^{\href https://en.wikipedia.org/wiki/Heat_map}$ time-wise and as a
whole to visualize the whole thing properly. In each of the heatmaps, we see
there is a white off-diagonal line which indicates a value close to 1 because
it is self-similarity. The white portions indicate relatively strong
similarity to the blue portions. As the score is bounded by 1, hence there is
no red box in the figures according to the color scale.

  * •

In 2013, Golden Whistler, Brown Quail, Varied Sittella, and Southern Whiteface
these species behaved most similarly as we can see from the central white part
of the heat map of 2013.

  * •

In 2012, Little Corella, Golden Whistler, Apostlebird, Brown Quail, and
Sittella these species behaved very similarly as we can see from the lower
white part of the heat map of 2012.

  * •

In 2011, Little Corella, Apostlebird, Little Raven, etc behaved very similarly
as in the figure of the heat map of 2011.

  * •

In each period, more groups of birds behaved very similarly as we can see from
the heatmaps, there are various white portions in each of them. There are some
deep blue boxes in each heat map also which indicates strong dissimilarity.

Figure 9: Heat map of the similarity of the species for all three years Figure
10: Heat map of the similarity of the species for 2011 Figure 11: Heat map of
the similarity of the species for 2012 Figure 12: Heat map of the similarity
of the species for 2013

##  6 Randomization

In different time periods, we are keeping the marginals of each species fixed.
Now our index highly depends on the prevalence since there is nothing
stochastic. So in order to get rid of the issue we introduce the method of
randomization: The initial distribution of $X_{0}$ which denotes the presence
of each species in each time period is given by
$\mu=\begin{bmatrix}p&1-p\end{bmatrix}$ where $p$ = 1- $\hat{p}_{i,r}$[If it
is done for $i$ th species and for $r$th time period]. We will change the
entries 0 or 1 randomly with the probabilities

$P(X_{1}=0|x_{0}=0)=a$

$P(X_{1}=1|x_{0}=0)=1-a$

$P(X_{1}=0|x_{0}=1)=b$

$P(X_{1}=1|x_{0}=1)=1-b$

where $X_{1}$ denotes the entries after the randomisation. So the marginal
distribution of $X_{1}$ is given by

$X_{1}\sim\begin{bmatrix}p&1-p\end{bmatrix}\begin{bmatrix}a&1-a\\\
b&1-b\end{bmatrix}$

Now our transformation should have the following properties:

  * •

The marginal distribution of $X_{0}$ and $X_{1}$ should be same.

  * •

The Kullback Leibler Distance${}^{\href
https://projecteuclid.org/journals/annals-of-mathematical-
statistics/volume-22/issue-1/On-Information-and-
Sufficiency/10.1214/aoms/1177729694.full}$ between $X_{0}$ and $X_{1}$ should
be minimum.

To satisfy the first condition we need to select $a,b$ so that the following
equation holds:

$\begin{bmatrix}p&1-p\end{bmatrix}=\begin{bmatrix}p&1-p\end{bmatrix}\begin{bmatrix}a&1-a\\\
b&1-b\end{bmatrix}$

So we have the equation

$ap+b(1-p)=p$ (15)

which gives the following relations between $a,b$

$b=\frac{p}{1-p}(1-a),0<p<1,0<a<1$ (16)

Now we choose a so that the Kullback Leiblar Divergence${}^{\href
https://projecteuclid.org/journals/annals-of-mathematical-
statistics/volume-22/issue-1/On-Information-and-
Sufficiency/10.1214/aoms/1177729694.full}$ between $X_{0}$ and $X_{1}$ is
minimal. So for that we at first calculate the Kullback Distance between
$X_{0}$ and $X_{1}$. The joint distribution of $(X_{0},X_{1})$ is given by

$p(0,1)=b(1-p)$

$p(1,0)=(1-a)p$

$p(0,0)=ap$

$p(1,1)=(1-b)(1-p)$

So the Kullback Leibler Divergence $D$ is given by

$D=p(0,1)\textbf{log}(\frac{p(0,1)}{p(0)p(1)})+p(1,0)\textbf{log}(\frac{p(1,0)}{p(1)p(0)})+p(0,0)\textbf{log}(\frac{p(0,0)}{p(0)p(0)})+p(1,1)\textbf{log}(\frac{p(1,1)}{p(1)p(1)})$

Hence,

$\displaystyle
D=b(1-p)log(\frac{b(1-p)}{p(1-p)})+(1-a)plog(\frac{(1-a)p}{p(1-p)})+$
$\displaystyle
aplog(\frac{ap}{p^{2}})+(1-b)(1-p)log(\frac{(1-b)(1-p)}{(1-p)^{2}})$
$\displaystyle=b(1-p)log(\frac{b}{p})+(1-a)plog(\frac{(1-a)}{(1-p)})+aplog(\frac{a}{p})+$
$\displaystyle(1-b)(1-p)log(\frac{(1-b)}{(1-p)})$

Since from (16) it can be seen that optimum $b$ is a function of $a$, hence
this expression is a function of $a$. So

$D(a)=2(1-a)p\textbf{log}(\frac{1-a}{1-p})+ap\textbf{log}(\frac{a}{p})+(1-2p+ap)\textbf{log}(\frac{1-2p+ap}{(1-p)^{2}})$
(17)
$D^{\prime}(a)=-2p(log(\frac{a-1}{p-1})+1)+p(log(\frac{a}{p})+1)+p(log(\frac{p(a-2)+1}{(p-1)^{2}})+1)$
(18)

$D^{\prime}(a)=0$ yields the solution $a=p$, since $D(p)=0$ and by Gibbs’
Inequality $D(a)\geq 0$ $\forall a\in(0,1)$ so

$p=\underset{a\in(0,1)}{\text{arg min}}D(a)$

We apply this randomization $N$ times where $N$ is a large natural number.
Then for each step, we get a score vector and store them in a matrix row-wise.
Then we normalize each row (i.e. scale down the norm of every row to 1) and
then take the mean column-wise to get an average score vector of all the
randomized datasets. We call this vector $N$, a vector of dimension $nl$.
So our similarity index (after randomisation) for $i^{\text{ th}}$ island will
be

1 - standard deviation of $N_{i},N_{i+n},..,N_{i+n(l-1)}$.

So the possible range of this index is $(-\infty,1]$ and the more it is close
to 1 indicates more association and the more negative indicates less
association or similarity.

##  7 Results after Randomization

The main goal of randomization was to vary through all the prevalence sets
keeping the marginal of each species in each time period fixed. So, we
performed 4 sets of randomizations with $N=1000,1000,2000$, and $5000$
respectively. We noticed a strong consistency in the results. That is, in each
such randomization set, we got very similar results. And surprisingly, We
reanalyzed the island dataset using the randomization method with our index.
And that describes the importance of randomization. Randomization means we are
going through all possible prevalences but keeping the distribution of each
species in different periods fixed. That’s why now the whole trend goes in
favor of $\hat{\alpha}.$ Here we only described the results using
$\hat{\beta}_{L}$, as we have seen before $\hat{\beta}_{L}$ and
$\hat{\beta}_{N}$ yields very similar results. We have seen the pattern of
$\hat{\alpha}$ now matches with $\hat{\beta}_{L}$ after randomization.

Figure 13: Results after first set of (1k) randomizations

Figure 14: Results after second set of (1k) randomizations

Figure 15: Results after third set of (2k) randomizations

Figure 16: Results after fourth set of (5k) randomizations

##  8 Handling Missing Data Case

In many practical situations, it’s often very challenging to collect
ecological data which leads to the case of missing data. In case we use the
usual Bootstrap way of re-sampling to deal with the missing data case:
We first estimate based on whatever we have for each species in each time
period and then randomly generate from that estimated distribution for the
missing parts. And then apply the whole process as before. We here present the
result after making the data missing and the actual data without
randomization. We worked on the island data here also. We have artificially
randomized the process of data removal and in each time period, we removed 16
data points, i.e. a total of 32 data points from the whole data, and then
applied bootstrap to account for those. And then we plotted those in the
figure and clearly, we can see that in the same scale, it goes along the line
$X=Y$. Hence this method of handling missing data performs well here. Hence as
long as there is sufficient data relative to missing, the method should
perform well.

Figure 17: Plot of score based on missing and full data

##  9 References

  1. 1.

The Distribution Of The Flora In The Alpine Zone: Paul Jaccard.

  2. 2.

A method of establishing groups of equal amplitude in plant sociology based on
similarity of species and its application to analyses of the vegetation on
Danish commons: Sørensen.

  3. 3.

Measurement of Diversity: E. H. Simpson

  4. 4.

A better index for analysis of co-occurrence and similarity: Kumar P. Mainali,
Eric Slud, Michael C. Singer, And William F. Fagan

  5. 5.

The value of all chances in games of fortune, cards, dice, wagers, lotteries,
c.: mathematically demonstrated, Christiaan Huygens, William Browne, Augustus
De Morgan

  6. 6.

R.A. Fisher and the making of maximum likelihood 1912-1922: John Aldrich

  7. 7.

Principal component analysis: a review and recent developments: Ian T.
Jolliffe and Jorge Cadima

  8. 8.

Principal Components Analysis | SAS Annotated Output

  9. 9.

Online local learning algorithms for linear discriminant analysis: G.K. Demir,
K. Ozmehmet

  10. 10.

Species Co-Occurrence: A Meta-Analysis Of J. M. Diamond’s Assembly Rules
Model: Nicholas J. Gotelli, Declan J. McCabe.

  11. 11.

Antagonism between Bacteriostatic and Bactericidal Antibiotics Is Prevalent:
Paolo S. Ocampo, Viktória Lázár, Balázs Papp, Markus Arnoldini, Pia Abel zur
Wiesch, Róbert Busa-Fekete, Gergely Fekete, Csaba Pál, Martin Ackermann, and
Sebastian Bonhoeffera

  12. 12.

Assemblage Time Series Reveal Biodiversity Change but Not Systematic Loss:
Maria Dornelas, Nicholas J. Gotelli, Brian Mcgill, Hideyasu Shimadzu, Faye
Moyes, Caya Sievers, And Anne E. Magurran.

  13. 13.

On measures of association among genetic variables: Daniel Gianola, Eduardo
Manfredi, and Henner Simianer.

  14. 14.

Diversity and endemism of the benthic seamount fauna in the southwest Pacific:
Bertrand Richer de Forges, J. Anthony Koslow and G. C. B. Poore.

  15. 15.

The geography of biodiversity change in marine and terrestrial assemblages:
Shane A Blowes, Sarah R Supp, Laura H Antão, Amanda Bates, Helge Bruelheide,
Jonathan M Chase, Faye Moyes, Anne Magurran, Brian McGill, Isla H Myers-Smith,
Marten Winter, Anne D Bjorkman, Diana E Bowler, Jarrett E K Byrnes, Andrew
Gonzalez, Jes Hines, Forest Isbell, Holly P Jones, Laetitia M Navarro, Patrick
L Thompson, Mark Vellend, Conor Waldock, Maria Dornelas

  16. 16.

J. Cornfield, in Proceedings of the Third Berkeley Symposium on Mathematical
Statistics and Probability (University of California Press, 1956), vol. 4, pp.
135–148

  17. 17.

Cognitive relevance of the community structure of the human brain functional
coactivation network: Nicolas A Crossley, Andrea Mechelli, Petra E Vértes,
Toby T Winton-Brown, Ameera X Patel, Cedric E Ginestet, Philip McGuire, Edward
T Bullmore

  18. 18.

Long-term phenotypic evolution of bacteria: Germán Plata, Christopher S. Henry
and Dennis Vitkup

  19. 19.

The Export Similarity Index: Some Structural Tests: Mitchell Kellman and Tim
Schroder

  20. 20.

Employing Diversity and Similarity Indices to Evaluate Geopedological Soil
Mapping in Miyaneh, East Azerbaijan Province, Iran: Naser Nazari, Shahla
Mahmoodi, Mohammad Hassan Masihabadi

  21. 21.

Comparative Analyses of Diversity and Similarity Indices of West Bank Forest
and Block A Forest of the International Institute of Tropical Agriculture
(IITA) Ibadan, Oyo State, Nigeria: Ariyo Oluyinka Christopher

  22. 22.

Of Beta-Diversity And Its Decrease With Increasing Elevation

  23. 23.

Disentangling the Drivers of $\beta$ Diversity Along Latitudinal and
Elevational Gradients

  24. 24.

Interpretation of Heatmaps, Wikipedia

  25. 25.

On Information and Sufficiency : S. Kullback, R. A. Leibler

  26. 26.

Truncated Normal Distribution, Wikipedia

  27. 27.

Bootstrap Re-sampling Technique, Wikipedia

  28. 28.

Kernel Density Estimation, Wikipedia

  29. 29.

Species co-occurrence analysis predicts management outcomes for multiple
threats, Ayesha I. T. Tulloch, Iadine Chadès and David B. Lindenmayer

  30. 30.

Plant recording across two centuries reveals dramatic changes in species
diversity of a Mediterranean archipelago, Alessandro Chiarucci, Simone
Fattorini, Bruno Foggi, Sara Landi, Lorenzo Lazzaro, János Podani and Daniel
Simberloff

##  10 Appendix

  1. 1.

Theorem: Let $W$ be Positive Definite.Then

$\underset{\vec{x}\not=\vec{0}}{\text{max}}\frac{\vec{x}^{T}B\vec{x}}{\vec{x}^{T}W\vec{x}}=\lambda_{1}(W^{-1}B)$

where $\lambda_{1}(W^{-1}B)$ is the largest eigenvalue of $W^{-1}B$ and the
maximum is attained at $\vec{x}_{0}$ iff $\vec{x}_{0}$ is an eigenvector of
$W^{-1}B$ corresponding to $\lambda_{1}(W^{-1}B)$.

###### Proof.

Let $\mu=\lambda_{1}(W^{-1}B)$ and $W=C^{T}C$ where $C$ is non-singular.
Writing $C\vec{x}=\vec{y}$ we get

$\underset{\vec{x}\not=\vec{0}}{\text{max}}\frac{\vec{x}^{T}B\vec{x}}{\vec{x}^{T}W\vec{x}}=\underset{\vec{y}\not=\vec{0}}{\text{max}}\frac{\vec{y}^{T}(C^{-1})^{T}BC^{-1}\vec{y}}{\vec{y}^{T}\vec{y}}=\lambda_{1}((C^{-1})^{T}BC^{-1})$

Now the characteristic roots of $(C^{-1})^{T}BC^{-1}$ are all real and, by the
next Theorem are the same as the characteristic roots of
$C^{-1}(C^{-1})^{T}B=W^{-1}B$ Hence the theorem follows. Also,
$\frac{\vec{x}^{T}B\vec{x}}{\vec{x}^{T}W\vec{x}}=\mu$ iff $C\vec{x}$ is an
eigenvector of $(C^{-1})^{T}BC^{-1}$ corresponding to $\mu$ which is
equivalent to: $\vec{x}$ is an eigenvector of $W^{-1}B$ corresponding to
$\mu$. ∎

  2. 2.

Theorem : Let $A,B$ are $n\times n$ matrices and $A$ is non-singular
(invertible). Then the characteristic polynomial of $AB$ and $BA$ are the
same.

###### Proof.

Observe since $A$ is invertible, we can write

$A^{-1}(AB)A=BA$

So $AB$ and $BA$ are similar. Since similar matrices have the same
characteristic polynomial( proved as the next theorem) we have characteristic
polynomials of $AB$ and $BA$ are the same. ∎

  3. 3.

Theorem : Characteristic polynomials of similar matrices are the same.

###### Proof.

Say, $A$ and $B$ are similar matrices. So, there exists nonsingular $P$ such
that $A=PBP^{-1}$. So,

$det(A-\lambda I)=det(PBP^{-1}-\lambda I)=det(P(B-\lambda
I)P^{-1})=det(P)det(B-\lambda I)det(P^{-1})=det(B-\lambda I)$

∎