# Band structures and contact points in phosphorene superlattice

Jilali Seffadi  Laboratory of Theoretical Physics, Faculty of Sciences,
Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Ilham
Redouani  Laboratory of Theoretical Physics, Faculty of Sciences, Chouaïb
Doukkali University, PO Box 20, 24000 El Jadida, Morocco  Youness Zahidi  MRI
Labortory, Polydisciplinary Faculty, Sultan Moulay Selimane University,
Khouribga, Morocco  Ahmed Jellal  Laboratory of Theoretical Physics, Faculty
of Sciences, Chouaïb Doukkali University, PO Box 20, 24000 El Jadida, Morocco
Canadian Quantum Research Center, 204-3002 32 Ave Vernon, BC V1T 2L7, Canada

###### Abstract

We study the band structures and the associated contact points for a
phosphorene superlattice made up of two periodic areas. We use the boundary
conditions to extract an equation describing the dispersion relation after
obtaining the eigen-wavefunctions. We show that energy transforms into linear
behavior near contact points, and fermions move at different speeds along
$x$\- and $y$\- directions. It was discovered that the periodic potential
caused additional Dirac points, which we located in $k$-space by establishing
their positions. We demonstrate that the barrier height and width can be used
to adjust the energy gap and modify the contact points. It might be that our
findings will be useful in the development of phosphorene-based electronic
devices.

###### pacs:

73.22.-f; 73.63.Bd; 72.10.Bg; 72.90.+y
Keywords: Phosphorene, superllatice potential, band stractures, Dirac points.

##  I INTRODUCTION

The electrical characteristics of 2D materials like graphene have received
enormous attention during the past few years. It has attracted great interest
in the scientific community since its realization in the laboratory [1, 2].
Such advantages result, especially, from the fact that a Hamiltonian of the
Dirac type governs the physics of the low-energy in graphene. Pure graphene
has no gaps and may be represented by the massless Dirac equation. The Klein
effect [3, 4, 5] prevents electrons in graphene from being constrained by
electrostatic potentials. Due to this characteristic, graphene cannot be used
in electronic devices. The most challenging task is to figure out how to use
electric fields to regulate the behavior of charge carriers in graphene. When
a gap is produced between the conduction and valence bands, Klein tunneling is
restrained, and the electronic properties of graphene change radically. To
contain the massless Dirac particle in the graphene sheet, a gap in the
spectrum’s energy can be opened [6, 7, 8]. Because of its straight bandgap,
high mobility, and anisotropic transport and mechanical properties,
phosphorene can compensate well for the lack of performance of other 2D
materials [9, 10].

In 2014, two research groups published the first report on the isolation of
phosphorene layers [11, 12], and today, there are two basic techniques for
making phosphorene: liquid exfoliation and scotch-tape cleavage [13, 14].
Unlike graphene, which has layers that are all plane, phosphorene is a weakly
van der Waals-bonded layered material with puckered surfaces on each layer
caused by the $sp^{3}$ hybridization of the $3s$ and $3p$ atomic orbitals
[15]. Phosphorene has an anisotropic puckered honeycomb band structure, which
distinguishes it from that of graphene. Because of this, graphene has a unit
cell with only two atoms, whereas there are four atoms in phosphorene [16].
Additionally, phosphorene’s strong anisotropy allows it to exhibit intriguing
direction-dependent transport and optical features. Phosphorene’s inherent
band gap makes it potentially useful for adjusting the transition between
conducting and insulating states. Due to its distinct electrical
characteristics and possible uses in optoelectronic, nanoelectronic, and other
applications, phosphorene has received a lot of attention [17, 18]. In fact, a
generation of anisotropic massless Dirac fermions and asymmetric Klein
tunneling in few-layer phosphorene superlattices was observed [19]. It is
shown that the highly anisotropic character of phosphorene in the presence of
bias allows for a field-induced semiconductor-metal transition [20].

Recent studies have used techniques including first-principles calculations,
$k\cdot p$ methods, and tight-binding models to determine the electronic
dispersion since it is crucial to have a complete understanding of the band
structure and charge carrier dynamics in phosphorene [21]. These studies
provided estimates of the energy gap for single-layer and multilayer
phosphorene, as well as evidence of significant anisotropy in effective mass.
In this respect, we study the band structures of phosphorene in the framework
of a continuum model, derived as the long-wavelength limit of a newly proposed
tight-binding mode. This study will be carried out by using a periodic
potential, resulting in a phosphorene superlattice. We describe the conditions
that allow this to occur and identify all Dirac points that appear when such
potential exists. In addition, we numerically analyze our results by
distinguishing between equal and unequal well and barrier widths under various
conditions of the physical parameters to underline the basic features of our
system. It turns out that the prediction that the periodic potential could be
used to manipulate the system band structure has increased interest in
phosphorene superlattices and opened up new avenues for the fabrication of
phosphorene-based electronic devices.

The present paper is organized as follows. In Sec. II, we provide the
theoretical model that underpins our approach and aims to derive the
dispersion relation for the entire cell explicitly. In Sec. III, we apply the
implicit function theorem to the dispersion relation to find the contact
points and identify their locations. As an illustration, we consider the first
Dirac point and derive the velocity components in the $x$\- and
$y$-directions, as well as show their behaviors against the barrier
parameters. In Sec. IV, for the equal and unequal well and barrier widths, we
present theoretical and numerical results based on various conditions. Our
final objective is to compare these two situations and demonstrate that they
are not equal. We conclude our results in Sec. V.

##  II Theoretical model

We analyze the band structure of charge carriers in phosphorene by taking into
account a periodic one-dimensional potential along the x-direction as depicted
in Fig. 1. We suggest that this structure can be created by giving the
phosphorene layer a local top gate voltage [22]. The one-dimensional
superlattice potential structure characterized by the heights $(V_{B},V_{W})$
and widths $(d_{B},d_{W})$. It is provided by the potential profile

$V_{j}(x)=\left\\{\begin{array}[]{llll}{V_{B}}&\mbox{if}&{id<x<id+d_{B}}\\\
{V_{W}}&\mbox{if}&{id+d_{B}<x<(i+1)d}\end{array}\right.$ (1)

where $i\in\mathbb{N}$ and $\left(j=B,W\right)$. The period is defined as a
series of wells and barriers that alternate in width $d$ $(d=d_{B}+d_{W})$.

Figure 1:  (Color online) Schematic representation of the periodic potentials
applied on the monolayer phosphorene.

The Hamiltonian for monolayer phosphorene can be calculated using the low-
energy five-hopping parameter tight-binding technique [23, 20, 16, 24]. This
process yields

$H(k)=\begin{pmatrix}u_{0}+\eta_{x}k_{x}^{2}+\eta_{y}k_{y}^{2}&\delta+\gamma_{x}k_{x}^{2}+\gamma_{y}k_{y}^{2}+i\chi
k_{x}\\\ \delta+\gamma_{x}k_{x}^{2}+\gamma_{y}k_{y}^{2}-i\chi
k_{x}&u_{0}+\eta_{x}k_{x}^{2}+\eta_{y}k_{y}^{2}\\\ \end{pmatrix}$ (2)

with the parameters quantities ${u}_{0}=-0.42$ eV, $\eta_{{x}}=0.58$ eV Å2,
$\eta_{{y}}=1.01$ eV Å2, $\delta=0.76$ eV, $\chi=5.25$ eV Å,
$\gamma_{{x}}=3.93$ eV Å2 and $\gamma_{{y}}=3.83$ eV Å. The phosphorene
superlattice is formed by n elementary cells positioned between the incidence
and transmission areas. The Hamiltonian described below can be used to
characterize our system in each $j$ region of the $n^{\text{th}}$ elementary
cell

$H_{j}(k)=H(k)+V_{j}(x)\mathbb{I}_{2}.$ (3)

By taking into account the conservation of the transverse wave vector $k_{y}$,
we can express the eigenspinors of these quasi-particles traveling along the
$\pm x$-directions in all areas as $\Phi_{j}(x,y)=\varphi_{j}(x)e^{ik_{y}y}$.
By resolving the eigenvalue equation $H_{j}\varphi_{j}=E_{j}\varphi_{j}$, the
associated eigenspinors can be found as

$\displaystyle\varphi_{j}=\Theta_{j}\cdot\Gamma_{j}$ (4)

by setting the quantities

$\displaystyle\Theta_{j}=\begin{pmatrix}e^{ik_{j}x}&e^{-ik_{j}x}\\\
z_{j}e^{ik_{j}x}&z_{j}^{-1}e^{-ik_{j}x}\end{pmatrix},\quad\Gamma_{j}=\begin{pmatrix}{a_{j}}\\\
{b_{j}}\\\ \end{pmatrix}$ (5)

as well as

$\displaystyle
z_{j}=s_{j}e^{-i\theta_{j}},\quad\theta_{j}=\arctan\left(\frac{\chi
k_{j}}{\delta+\gamma_{y}k_{y}^{2}}\right)$ (6)

where $s_{j}={\mbox{sign}}{\left(E_{j}-V_{j}-u_{0}-\eta_{y}k_{y}^{2}\right)}$.

We will then determine the probability that electrons will successfully pass
through the potential barriers. Specifically, for the transmission and
incidence regions, we have

$\displaystyle\Gamma_{\sf in}=\left(\begin{array}[]{c}{1}\\\
{r}\end{array}\right),\quad\Gamma_{\sf tr}=\left(\begin{array}[]{c}{t}\\\
{0}\end{array}\right)$ (11)

where $t$ and $r$ are the transmission and reflection coefficients, which can
be determined by imposing the continuity of the eigenspinors at interfaces.
The transfer matrix formalism is the most convenient way to explain this
process [25]. As a result, we obtain

$\displaystyle\begin{pmatrix}{1}\\\
{r}\end{pmatrix}=\Lambda\cdot\begin{pmatrix}{t}\\\ {0}\end{pmatrix}$ (12)

where the transfer matrix $\Lambda$ is given by

$\Lambda=\Theta_{\sf in}^{-1}[0]\cdot\lambda^{n}\cdot\Theta_{\sf
tr}[nd]=\begin{pmatrix}\Lambda_{11}&\Lambda_{12}\\\
\Lambda_{21}&\Lambda_{22}\end{pmatrix}$ (13)

and $\lambda$ takes the form

$\displaystyle\lambda=\Theta_{B}[0]\cdot\Theta_{B}^{-1}[d_{B}]\cdot\Theta_{W}[d_{B}]\cdot\Theta_{W}^{-1}[d].$
(14)

Now, using the Bloch theorem as well as the continuity at different interfaces
to get the result

$\displaystyle\varphi_{B}(0)\Gamma_{B}=\varphi_{W}(d)\Gamma_{W}e^{ik_{x}d}$
(15) $\displaystyle\varphi_{B}(d_{B})\Gamma_{B}=\varphi_{W}(d_{B})\Gamma_{W}$
(16)

where $d$ is the length of the unit cell and $k_{x}$ is the Bloch wave vector.
Thus, we can express $\Gamma_{W}$ in terms of $\Gamma_{B}$ as

$\displaystyle\Gamma_{W}$
$\displaystyle=\varphi_{W}^{-1}(d)\varphi_{B}(0)\Gamma_{B}e^{-ik_{x}d}$ (17)
$\displaystyle=\varphi_{W}^{-1}(d_{B})\varphi_{B}(d_{B})\Gamma_{B}.$ (18)

Using (5) along with (17) and (18) to end up with the relation

$\Xi\begin{pmatrix}{a_{B}}\\\ {b_{B}}\\\ \end{pmatrix}=0$ (19)

where the matrix $\Xi$ is provided by

$\Xi=\varphi_{W}^{-1}(d)\varphi_{B}(0)e^{-ik_{x}d}-\varphi_{W}^{-1}(d_{B})\varphi_{B}(d_{B}).$
(20)

Calculating $\det(\Xi)=0$ yields the dispersion relation

$\cos(k_{x}d)=\cos(d_{B}k_{B})\cos(d_{W}k_{W})+\Upsilon\sin(d_{B}k_{B})\sin(d_{W}k_{W})$
(21)

such that the wave vectors and the parameter $\Upsilon$ are

$\displaystyle
k_{B}=\frac{1}{\chi}\sqrt{(E-V_{B}-\mu_{0}-\eta_{y}k^{2}_{y})^{2}-(\delta+\gamma_{y}k^{2}_{y})^{2}}$
(22) $\displaystyle
k_{W}=\frac{1}{\chi}\sqrt{(E-V_{W}-\mu_{0}-\eta_{y}k^{2}_{y})^{2}-(\delta+\gamma_{y}k^{2}_{y})^{2}}$
(23)
$\displaystyle\Upsilon=\coth\theta_{B}\coth\theta_{W}-\frac{s_{B}s_{W}}{\sin{\theta_{B}}\sin{\theta_{W}}}.$
(24)

In the analysis that will follow, a theoretical study will be done to shed
light on the Dirac points. Also, we will see how the aforementioned findings
can be applied to understanding band structures. In fact, a numerical analysis
will be established to investigate the basic features of the present system.

##  III Dirac contact points

To establish the Dirac contact points, we apply the implicit function theorem
to the dispersion relation (21). Indeed, let us write (21) as an implicit
function of the form

$\begin{split}h(k_{x},k_{y},E)=0\end{split}$ (25)

which can be expanded near the contact point $(k_{xc},k_{yc},E_{c})$ as
follows

$\begin{split}h(k_{x},k_{y},E)=h(k_{xc},k_{yc},E_{c})+\Delta P\nabla
h(k_{xc},k_{yc},E_{c})+\frac{1}{2}\Delta
P^{t}\mathbb{H}h(k_{xc},k_{yc},E_{c})\Delta P\end{split}$ (26)

and we have

$\begin{split}h(k_{xc},k_{yc},E_{c})=0,\quad\nabla
h(k_{xc},k_{yc},E_{c})=0,\quad\Delta P=\begin{pmatrix}\Delta k_{x}\\\ \Delta
k_{y}\\\ \Delta E\end{pmatrix}\end{split}$ (27)

where $\mathbb{H}$ is the Hessian matrix of second-order partial derivatives
of $h$, acting as

$\begin{split}\mathbb{H}h(k_{xc},k_{yc},E_{c})=\begin{pmatrix}\rho&0&0\\\
0&\sigma&0\\\ 0&0&\tau\end{pmatrix}\end{split}$ (28)

and ($\rho$, $\sigma$, $\tau$) can be fixed in terms of the physical
parameters. Now injecting all into (26) to end up with the interesting result

$\begin{split}\rho\Delta k_{x}^{2}+\sigma\Delta k_{y}^{2}+\tau\Delta
E^{2}=0\end{split}$ (29)

or more explicitly

$E-E_{c}=\pm\sqrt{-\frac{\rho}{\tau}(k_{x}-k_{xc})^{2}-\frac{\sigma}{\tau}(k_{y}-k_{yc})^{2}}$
(30)

which is the energy of a generic Dirac Hamiltonian with anisotropic Fermi
velocity

$\begin{split}v_{y}=\frac{1}{\hbar}\sqrt{-\frac{\sigma}{\tau}},\quad
v_{x}=\frac{1}{\hbar}\sqrt{-\frac{\rho}{\tau}}.\end{split}$ (31)

where $\hbar$ is the Planck constant. This result is in agreement with that
obtained in [26] by studying the tunable electronic properties of multilayer
phosphorene and its nanoribbons.

As illustration, we consider the first contact point $(0,0,E_{0})$, and
therefore obtain

$\mathbb{H}h(0,0,E_{0})=\begin{pmatrix}d^{2}&0&0\\\ 0&\sigma_{1}&0\\\
0&0&\tau_{1}\end{pmatrix}$ (32)

such that $\sigma_{1}$ and $\tau_{1}$ are given by

$\displaystyle\sigma_{1}=4s_{B}s_{W}{\frac{\delta\sin\beta\left(\epsilon\cos\beta+\varsigma\sqrt{V_{B}^{2}-4{\delta}^{2}}\sin\beta\right)}{\chi\left(V_{B}^{2}-4{\delta}^{2}\right)^{5/2}}}$
(33)
$\displaystyle\tau_{1}=-s_{B}s_{W}{\frac{\zeta\cos\beta\sin\beta+\sqrt{V_{B}^{2}-4{\delta}^{2}}\left(\xi\cos^{2}\beta+\upsilon\right)}{2{\chi}^{2}\left(V_{B}^{2}-4\delta^{2}\right)^{7/2}}}$
(34)

and the quantities are

$\displaystyle\beta={\frac{d\sqrt{V_{B}^{2}-4{\delta}^{2}}}{4\chi}}$ (35)
$\displaystyle\epsilon=s_{B}s_{W}\left(V_{B}^{4}d\gamma_{y}-12V_{B}^{2}d\delta^{2}\gamma_{y}+32d{\delta}^{4}\gamma_{y}\right)+V_{B}^{4}d\gamma_{y}-4V_{B}^{2}d{\delta}^{2}\gamma_{y}$
(36)
$\displaystyle\varsigma=4s_{B}s_{W}\gamma_{y}V_{B}^{2}\chi-4\gamma_{y}V_{B}^{2}\chi$
(37) $\displaystyle\zeta=s_{B}s_{W}\left(-16V_{B}^{4}\chi
d{\delta}^{2}+192V_{B}^{2}\chi d{\delta}^{4}-512\chi
d{\delta}^{6}\right)-16V_{B}^{4}\chi d{\delta}^{2}+64V_{B}^{2}\chi
d{\delta}^{4}$ (38)
$\displaystyle\xi=s_{B}s_{W}\left(64V_{B}^{2}{\chi}^{2}{\delta}^{2}+256{\chi}^{2}{\delta}^{4}\right)-192V_{B}^{2}{\chi}^{2}{\delta}^{2}+256{\chi}^{2}{\delta}^{4}$
(39)
$\displaystyle\upsilon=\left(s_{B}s_{W}-1\right)\left(V_{B}^{6}{d}^{2}-4V_{B}^{4}{d}^{2}{\delta}^{2}\right)-s_{B}s_{W}\left(64V_{B}^{2}{\chi}^{2}{\delta}^{2}-256{\chi}^{2}{\delta}^{4}\right)+192V_{B}^{2}{\chi}^{2}{\delta}^{2}-256{\chi}^{2}{\delta}^{4}$
(40)

Thus, we get the two velocity components

$\begin{split}v_{y}=\frac{1}{\hbar}\sqrt{-\frac{\sigma_{1}}{\tau_{1}}},\quad
v_{x}=\frac{1}{\hbar}\sqrt{-\frac{d^{2}}{\tau_{1}}}.\end{split}$ (41)

(a) $V_{B}=2,2.5,3$ eV

(b) $d=2,4,6$ nm

Figure 2:  (color online) Velocity components ($v_{x}$, $v_{y}$) versus
barrier height $V_{B}$ near the first Dirac point for $k_{x}=k_{y}=V_{W}=0$,
$E_{0}=\frac{V_{B}}{2}+\mu_{{0}}$ and $d_{W}=d_{B}$. (a): $V_{B}=2$ eV (black,
purple), $V_{B}=2.5$ eV (red, yellow), $V_{B}=3$ eV (green, blue). (b): $d=2$
nm (black, purple), $d=4$ nm (red, yellow), $d=6$ nm (green, blue)

To establish a deeper connection between the barrier parameters ($V_{B}$, $d$)
and the group velocities along the $k_{x}$-direction $v_{x}$ (black, red,
green) and the $k_{y}$-direction $v_{y}$ (purple, yellow, blue), we plot Fig.
2 for the Dirac point $(0,0,E_{0})$. Fig. 2a depicts the behavior of $v_{x}$
and $v_{y}$ versus the barrier width $d$ for three different values of its
height $V_{B}$. It is clear that $v_{x}$ starts to decrease around $0.8\
10^{6}$ m/s, which is consistent with literature [27, 26]. $v_{x}$ decreases
as $d$ increases, eventually stabilizing at some fixed value dependent on
$V_{B}$. On the other hand, as $d$ increases, $v_{y}$ reaches a maximum around
$0.4\ 10^{6}$ m/s and then falls to zero. After a certain value of $d$,
$v_{y}$ exhibits a variety of small oscillations that are strongly influenced
by the barrier height $V_{B}$. The velocities are plotted against the barrier
height $V_{B}$ for three different values of its width $d$ in Fig. 2b. We can
see that $v_{x}$ increases as $V_{B}$ increases, eventually reaching a
constant value close to $0.8\ 10^{6}$ m/s that is highly dependent on $d$. In
contrast, we see that $v_{y}$ behaves similarly to that shown in Fig. 2a, but
with a slight variation in decreasing and oscillating as long as $V_{B}$ and
$d$ increase. This is not surprising given that $v_{x}$ is strongly dependent
on $d$ but $v_{y}$ is less so, according to (41). Since the Fermi velocities
control how rapidly electrons may flow between various components of nanoscale
devices like transistors or memory cells, they have an impact on
nanotechnology.

##  IV Analysis of band Structures

We will study the electrical band structures by distinguishing two cases.
First, consider equal well and barrier widths, and then consider asymmetrical
well and barrier widths. The valence and conductance mini-bands will be the
main topics of study in the sections that follow. We will assume that the
period of the phosphorene superlattice is constant $d=d_{W}+d_{B}$. We will
numerically investigate the dispersion relation for the phosphorene
superlattice in (21) in terms of the physical parameters $(V_{B},d)$
describing the applied potential in order to explain the symmetrical impact of
such a possibility.

###  IV.1 Equal well and barrier widths $d_{B}=d_{W}$

Before to numerically analyze the dispersion relation, let us investigate some
interesting cases. Indeed, It is worth noting that when $k_{W}=k_{B}$, we get

$\begin{split}E_{{0}}=\frac{V_{{B}}}{2}+\mu_{{0}}+\eta_{{y}}k_{{y}}^{2}\end{split}$
(42)

which is nothing but the quadratic energy part of the pure phosphorene [24]
that will be established later on. Now for $k_{x}=0$ and $d_{W}=d_{B}$, (21)
reduces to

$\begin{split}1=\cos^{2}(d_{W}k_{W})+\Upsilon\sin^{2}(d_{W}k_{W})\end{split}$
(43)

which is satisfied when $k_{W}=\frac{n\pi}{d_{W}}$, and $n\in\mathbb{Z}$,
resulting in quantized energy

$\begin{split}E_{n}=V_{W}+\mu_{0}+\eta_{y}k_{y}^{2}\pm\sqrt{\left(\frac{\pi\chi}{d_{W}}\right)^{2}n^{2}+\left(\delta+\gamma_{y}k_{y}^{2}\right)^{2}}\end{split}$
(44)

and the transverse momentum $k_{{y}}$

$\displaystyle
k_{{yn}}={\sqrt{\frac{d_{W}\delta\,\gamma_{{y}}+d_{W}E_{n}\eta_{{y}}-d_{W}\eta_{{y}}\mu_{{0}}+\sqrt{\left(4\eta_{\gamma}\pi^{2}{\chi}^{2}\right){n}^{2}+2\,d_{W}^{2}\delta\,\eta_{{y}}\gamma_{{y}}\left(E_{n}-\mu_{{0}}\right)+\kappa_{\mu}+d_{W}^{2}\kappa_{\eta}}}{{d_{W}\eta_{\gamma}}}}}$
(45)

where we have set $\eta_{\gamma}=\eta_{y}^{2}-\gamma_{y}^{2}$,
$\kappa_{\eta}=\eta_{y}^{2}{\delta}^{2}+\gamma_{{y}}^{2}{\mu_{{0}}}^{2}$, and
$\kappa_{\mu}=E_{n}d_{W}^{2}\gamma_{{y}}^{2}\left(E_{n}-2\,\mu_{{0}}\right)^{2}$.
One obtains from the zeros of (45)

$\begin{split}V_{B}=V_{Bn}=2\sqrt{\left(\frac{2\pi\chi}{d}\right)^{2}n^{2}+\delta^{2}}\end{split}$
(46)

which are the $V_{B}$ values that fill the minigap when $k_{y}=0$. It is worth
noting that $V_{B}$ has a critical value, which is $V_{Bc}=V_{B1}$. We make
the positions of the contact points obvious by noting that when $V_{B}$
exceeds $V_{Bc}$, two additional Dirac points appear at $k_{y}\neq 0$. The
energy minigap at $k_{y}=0$ closes at $V_{B}=V_{B2}$, creating three contact
points in the system. When $V_{B}>V_{B2}$, the minigap at $k_{y}=0$ opens, but
two additional Dirac points appear instead, leaving the system with four
contact points, and so on. Then, in the generic case, $E_{n}$ and $k_{yn}$
indicate the location of the contact point that is closest to $k_{y}=0$ for
$V_{B}=V_{Bn}$ or $V_{Bn}<V_{B}<V_{Bn+1}$. As a result, the contact points are
located as

$\displaystyle(E,k_{x},k_{y})=(E_{n},0,k_{yn}),(E_{n-1},0,k_{yn-1}),\cdots,(E_{1},0,k_{y1}).$
(47)

Now, by requiring $d_{B}=0$, we may rewrite (21) as follows:

$\begin{split}\cos(d_{W}k_{W})-\cos(k_{x}d_{w})=0.\end{split}$ (48)

which produces the energy of phosphorene as found in [24]

$\begin{split}E_{\text{phos}}=\mu_{0}+\eta_{y}k^{2}_{y}+V_{W}\pm\sqrt{k^{2}_{x}\chi^{2}+(\delta+\gamma_{y}k^{2}_{y})^{2}}.\end{split}$
(49)

For the Dirac points $\left(k_{x}=\frac{m\pi}{d_{W}},k_{y}=0\right)$, (48)
gives

$\begin{split}\cos\left(\frac{d_{W}}{\chi}\sqrt{(E-\mu_{0})^{2}-\delta^{2}}\right)=\cos(m\pi),\quad
m\in\mathbb{Z}\end{split}$ (50)

showing the quantized energy

$\begin{split}E_{m}=\mu_{0}\pm\sqrt{\left({\frac{\pi{\chi}}{{d_{W}}}}\right)^{2}{m}^{2}+{\delta}^{2}}.\end{split}$
(51)

The energies of holes and electrons are plotted as a function of the
transverse momentum $k_{y}$ with $k_{x}=0$ and $d_{B}=d_{W}=4.75$ nm for
$V_{B}=1.65$ eV (blue), 1.66 eV (red), and 1.68 eV (green) in Fig. 3a. We see
it is attainable to downgrade and control the gap until it is gapless by
elevating the height potential. When $V_{B}$ reaches a critical value
$V_{Bc}=1.68$ eV, we notice two Dirac points within the band structure, one on
the positive side and the other on the negative side. Same results are
obtained in studying the generation of anisotropic massless Dirac fermions and
asymmetric Klein tunneling in few-layer black phosphorus superlattices [19].
We can also see, that the Dirac points move up in energy as $V_{B}$ increases.
We show that the energy is symmetric for $k_{y}=0$. We can clearly see that
the critical value becomes $V_{Bc}=1.66$ eV and then $V_{Bc}=1.65$ eV in Figs.
3(b,c), where $d=10$ and 10.5 nm, respectively, instead of $d=9.5$ nm in Fig.
3a. The same results were attained for the electronic structure of a graphene
superlattice with massive Dirac fermions [28]. This is in contrast to the
electronic structure of a graphene superlattice with a modulated Fermi
velocity [29]. As we increase $V_{B}$ above $V_{Bc}$, a band gap appears at
$k_{y}=0$ and grows monotonically with increasing $V_{B}$, as seen in [30].

(a) $d=9.5$ nm

(b) $d=10$ nm

(c) $d=10.5$ nm

Figure 3: (color online) The dispersion relation (21) as a function of $k_{y}$
for $k_{x}=0$, $V_{B}=1.65$ eV (blue), $V_{B}=1.66$ eV (red), $V_{B}=1.68$ eV
(green) with (a): $d=9.5$ nm, (b): $d=10$ nm, (c): $d=10.5$ nm.

Fig. 4 demonstrates the absence of an original Dirac point and any vertical
Dirac points in our study, while we have a gap. However, it is evident that
the potential has an impact on how Dirac points occur. It is clearly seen that
the number of Dirac points increases as $V_{B}$ increases. Fig. 4a shows that
our system behaves like pure phosphorene, with a gap for a tiny barrier
$V_{B}$ ($V_{B}<V_{Bc}=1.66$ eV) and energy that is quadratic along the
$k_{y}$-direction. Fig. 4b shows the development of two Dirac points when
$V_{B}$ reaches the critical value $V_{Bc}$, and Fig. 4c supports our findings
by showing that when $V_{B}$ grows, the number of Dirac points similarly
increases. The band structure at Dirac points changed from seeming parabolic
to linear, as seen in Figs. 4(b,c).

(a) $V_{B}=1.642$ eV

(b) $V_{B}=1.68$ eV

(c) $V_{B}=3$ eV

Figure 4: (color online) The dispersion relation (21) as a function of $k_{y}$
for $k_{x}=V_{W}=0$, $d_{B}=d_{W}=\frac{d}{2}=5$ nm with (a): $V_{B}=1.642$
eV, (b): $V_{B}=1.68$ eV, (c): $V_{B}=3$ eV.

To emphasize the importance of barrier and well widths, Fig. 5 depicts the
effect of distance $q$ for $V_{B}=1.68$ eV and $d=10$ nm. We observe that by
increasing $q$, it is possible to downgrade and control the gap until the
system becomes gapless. Furthermore, the minigap at $k_{y}=0$ opens for
$q=\frac{1}{2}$ (green). However, at various values of $k_{y}$ and for the
parameters selected here, we notice there are extra Dirac points. As $q$
rises, the gap narrows until Dirac points are reached. We conclude that there
are two extreme cases to consider. Indeed, the system behaves like a pristine
phosphorene with a direct band gap for $q<\frac{2}{5}$, whereas
$q\geq\frac{2}{5}$ exhibits extra Dirac points.

(a) $V_{B}=1.68$ eV, $d=10$ nm Figure 5: The dispersion relation (21) as a
function of $k_{y}$ for $d_{W}=qd$, $d_{B}=(1-q)d$ $d=10$ nm, and $k_{x}=0$.

In Fig. 6a the plots of dispersion relation is shown as a function of $k_{x}$
for $k_{y}=0.054\ \AA$, $V_{B}=1.65$ eV (blue), $V_{B}=1.67$ eV (red),
$V_{B}=1.68$ eV (green) and $d=10$ nm. The gap can be adjusted until it
disappears by varying $V_{B}$, as seen when $V_{B}=1.68$ eV (green). The
electron and hole energies shift upward with increasing potential, as can be
seen. Nevertheless, The electron energy shift is not the same as the hole
energy shift, which means that for different values of $V_{B}$ there are
different gaps, some of which may be zero. Noteworthy is the fact that the
Dirac point is obtained at $k_{x}=0$. The electron and hole energies are
plotted as a function of $V_{B}$ in Fig. 6b for $k_{x}=V_{W}=0$, $k_{y}=0.054\
\AA^{-1}$, and $d_{B}=d_{W}=5$ nm. It is feasible to see that the energy
oscillates and exhibits gapless behavior at specific values of $V_{B}$, which
depend on different parameters. It is critical to remember that these findings
are similar to those reported in [31].

(a) $d=10$ nm

(b) $d=10$ nm

Figure 6: (color online) (a): The dispersion relation (21) as a function of
$k_{x}$ for $k_{y}=0.054\ \AA$, $V_{B}=1.65$ eV (blue), $V_{B}=1.67$ eV (red),
$V_{B}=1.68$ eV (green) and $d=10$ nm. (b): The dispersion relation (21) as a
function of $V_{B}$ for $d=10$ nm, $k_{x}=0$, and $k_{y}=0.054\ \AA^{-1}$.

The band structure around the Dirac points is depicted in Fig. 7a for $d=10$
nm and $V_{B}=1.68$ eV. There are two Dirac points, which they are isotropic
and symmetrically located with respect to the point $k_{y}=0$. The projection
of the low-energy band centered on the Dirac points with $V_{B}=1.68$ eV is
shown in Fig. 7b. We can have two Dirac points by selecting such a value for
$V_{B}$. The new Dirac cones have more clearly defined isotropic behavior.
Energy is symmetric for $k_{y}=0$ as shown more clearly by the valence and
conduction bands of the spectrum in Fig. 7c. Moreover, at normal incidence
there is no original Dirac point (ODP), which means that there is a gap in
contrast to the case for [31] where there is an ODP in pristine graphene. The
most affecting characteristic observed here is that instead of having an
original zero-energy Dirac point (DP) of the pristine graphene at $k_{y}=0$
there is the appearance of a pair of new contact points. This effect can also
be seen at contact points in the energy band structure of bilayer graphene
superlattices [32] and graphene superlattices with periodically modulated
Dirac gaps [33].

(a) $d=10$ nm, $V_{B}=1.68$ eV, $k_{x}=V_{W}=0$

(b) $d=10$ nm, $V_{B}=1.68$ eV

(c) $d=10$ nm, $V_{B}=1.68$ eV

Figure 7: (color online) (a): The dispersion relation (21) as a function of
$k_{y}$ for $k_{x}=V_{W}=0$ and $V_{B}=1.68$ eV. (b): Contour plot of the
energy band $E\in[0.38,0.45]$ eV with step $5.10^{-4}$, $d=10$ nm and
$V_{B}=1.68$ eV. (c): Contour plot 3D of the energy band as a function of
$k_{y}$ and $k_{x}$.

###  IV.2 Non equal well and barrier widths $d_{B}\neq d_{W}$

For the case $d_{B}\neq d_{W}$, we use also (43) to end up with

$\displaystyle
k_{W}=\frac{1}{\chi}\sqrt{(E-V_{W}-\mu_{0}-\eta_{y}k^{2}_{y})^{2}-(\delta+\gamma_{y}k^{2}_{y})^{2}}=\frac{m\pi}{d_{W}}$
(52) $\displaystyle
k_{B}=\frac{1}{\chi}\sqrt{(E-V_{B}-\mu_{0}-\eta_{y}k^{2}_{y})^{2}-(\delta+\gamma_{y}k^{2}_{y})^{2}}=\frac{m\pi}{d_{B}}$
(53)

where $m\in\mathbb{Z}$. We get by subtracting (53) from (52) the quantized
energy

$E_{m}={\frac{{\pi}^{2}{\chi}^{2}{m}^{2}\xi_{-}+2V_{{B}}d_{{B}}^{2}d_{{W}}^{2}\epsilon_{0}\mp\sqrt{{\pi}^{2}{\chi}^{2}{\eta_{{y}}}^{2}\left[{\pi}^{2}{\chi}^{2}\xi_{-}^{2}{m}^{2}-V_{{B}}^{2}d_{{B}}^{2}d_{{W}}^{2}\left(2\xi_{+}-V_{{B}}^{2}d_{{B}}^{2}d_{{W}}^{2}\eta_{{y}}^{2}\right)\right]}}{2\gamma_{{y}}d_{{W}}^{2}d_{{B}}^{2}V_{{B}}}}$
(54)

and we have set the quantities $\xi_{\pm}=d_{B}^{2}\pm d_{W}^{2}$ and
$\epsilon_{0}=-\delta\,\eta_{{y}}+\gamma_{{y}}\mu_{{0}}+\frac{V_{{B}}\gamma_{{y}}}{2}$.
Injecting (54) into (52) to find

$k_{{y}}=k_{ym}={\frac{\sqrt{d_{{W}}\eta_{\gamma}\left(d_{{W}}\left(\delta\,\gamma_{{y}}+E_{m}\eta_{{y}}-\eta_{{y}}\mu_{{0}}\right)+\sqrt{{\pi}^{2}{\chi}^{2}{m}^{2}\eta_{\gamma}+\kappa^{2}-\gamma_{{y}}d_{{W}}^{2}\mu_{{0}}\left(2\,\delta\,\eta_{{y}}+2\,E_{m}\gamma_{{y}}-\gamma_{{y}}\mu_{{0}}\right)}\right)}}{d_{{W}}\eta_{\gamma}}}$
(55)

where $\eta_{\gamma}=\ \eta_{y}^{2}-\gamma_{y}^{2}$ and
$\kappa=\delta\,d_{{W}}\eta_{{y}}+E_{m}d_{{W}}\gamma_{{y}}$. From the zeros of
(55), we get

$V_{{{\it B}}}=V_{{{\it
Bm}}}=\sqrt{{\left(\frac{\pi\chi}{d_{B}}\right)^{2}{m}^{2}}+{\delta}^{2}}-\sqrt{{\left(\frac{\pi\chi}{d_{W}}\right)^{2}{m}^{2}}+{\delta}^{2}}$
(56)

We will examine the impacts of a well and barrier with different widths, i.e.,
a well and barrier that aren’t exactly the same size. studying the electrical
band structure using the dispersion relation (21). In Fig. 8a, one observes
that it is feasible to reduce and control the gap until it disappears by
increasing the height barrier. This behavior is similar to that obtained for
the effect of one-dimensional superlattice potentials on the band gap of two-
dimensional materials [34]. When $V_{B}$ reaches the critical value of
$V_{Bc}=1.76$ eV, we can see two Dirac points in the band structure, one on
the positive side and one on the negative side, at $k_{y}\neq 0$. We show that
the energy is symmetric for $k_{y}=0$. In Fig. 8b for $d_{B}=2d_{W}=6$ nm, we
clearly see that when the potential rises, the electron and hole energies
shift up. Nevertheless, the critical value is different and increased in
comparison to the first case when we had $d_{B}=d_{W}$, but in this case,
$V_{Bc}=1.76$ eV is the same in both cases when $d_{B}=2d_{W}$ and
$d_{W}=2d_{B}$. We notice the same behavior as in the previous work with a bit
of difference in the critical potential. These behaviors are similar to those
obtained for the electronic structure of a graphene superlattice with massive
Dirac fermions [28].

(a) $d=9$ nm, $d_{W}=2d_{B}$

(b) $d=9$ nm, $d_{B}=2d_{W}$

Figure 8: (color online) The dispersion relation (21) as a function of $k_{y}$
for $k_{x}=V_{W}=0$, $V_{B}=1.72$ eV (blue), $V_{B}=1.74$ eV (red),
$V_{B}=1.76$ eV (green). (a): $d=9$ nm, $d_{W}=2d_{B}$ and (b): $d=9$ nm,
$d_{B}=2d_{W}$.

(a) $d_{B}=3$ nm, $d_{W}=2d_{B}$

(b) $d_{B}=3$ nm, $d_{W}=2d_{B}$

Figure 9: (color online) (a): The dispersion relation (21) as a function of
$k_{y}$ for $k_{x}=V_{W}=0$ and $V_{B}=1.76$ eV. (b): Contour plot of the
energy band $E\in[0.396,0.403]$ eV with step $1.5.10^{-4}$, $d=9$ nm and
$V_{B}=1.76$ eV .

The dispersion relation as a function of $k_{y}$ is plotted in Fig. 9a for
$k_{x}=0$, $d_{W}=2d_{B}=6$ nm and $V_{B}=1.76$ eV, in order to show the
impact of different well and barrier widths. We observe two Dirac points
emerge, one on the positive side and the other on the negative one. Fig. 9b
displays the low-energy band’s contour plot around the Dirac points. It is
clearly seen that the spectrum is symmetric about the normal incidence, there
is one Dirac point on each side of the normal incidence. It is obvious that
Dirac points appear when $d_{W}=2d_{B}$. When we increase the width of the
well, $d_{W}$, the contact points clearly shift down, as shown in [28].

(a) $d_{W}=3$ nm, $d_{B}=2d_{W}$

(b) $d_{W}=3$ nm, $d_{B}=2d_{W}$

Figure 10: (color online) (a): The dispersion relation (21) as a function of
$k_{y}$ for $k_{x}=V_{W}=0$ and $V_{B}=1.76$ eV. (b): Contour plot of the
energy band $E\in[0.52,0.53]$ eV with step $3.10^{-4}$, $d=9$ nm and
$V_{B}=1.76$ eV.

In Fig. 10a we show the dispersion relation as a function of $k_{y}$ for
$k_{x}=V_{W}=0$, $d_{B}=2d_{W}=6$ nm, and $V_{B}=1.76$ eV. We observe a pair
of Dirac points come out in the band structure, one on the positive side and
the other on the negative side. Fig. 10b demonstrates the low energy band’s
contour plot around the Dirac points. As can be seen, the spectrum is
symmetric about the normal incidence, there is one Dirac point on each side of
the normal incidence. When we increase the width barrier $d_{B}$, we can
clearly see that the contact points have moved upward, similar to the results
obtained in [28]. Because we obtained similar results to those in Fig. 8 with
a minor change in the placement of the Dirac points, it is clear that Dirac
points rise when $d_{B}>d_{W}$ but fall when the reverse is true.

(a) $d_{B}=2d_{W}$, $d_{W}=100/3$ nm

(b) $d_{W}=2d_{B}$, $d_{B}=100/3$ nm

Figure 11: (color online) The dispersion relation (21) as a function of
$k_{x}$ for $k_{y}=0.054$ Å-1, $V_{W}=0$, and $V_{B}=1.68$ eV (blue),
$V_{B}=1.7$ eV (red), $V_{B}=1.706$ eV (green). (a): $d=10$ nm,
$d_{B}=2d_{W}$, and (b): $d=10$ nm, $d_{W}=2d_{B}$.

Fig. 11 depicts the dispersion relation as a function of $k_{x}$ for
$k_{y}=0.054$ Å-1, $V_{W}=0$, $V_{B}=1.68$ eV (blue), $V_{B}=1.7$ eV (red),
and $V_{B}=1.706$ eV (green). In Fig. 11a, for $d_{B}=2d_{W}=200/3$ nm, we can
tune the gap until gapless phosphorene by varying $V_{B}$, as happens when
$V_{B}=1.706$ eV (green). This clearly demonstrates that the energy can be
modulated by the barrier height $V_{B}$. However, because the shift in
electron energy is not the same as the shift in hole energy, there are many
gaps that could all be zero for different values of $V_{B}$. It should be
noted that the Dirac points are obtained around $E=0.5$ eV. Fig. 11b
reproduces the same results as before for $d_{W}=2d_{B}=200/3$ nm, except that
the energy is reduced by showing different shifts depending on the values of
$V_{B}$. It is clear that when the potential increases, the electron and hole
energies alter. The Dirac points, which are located at $E=0.395$ eV, are
slightly different in their locations, but we still see the same behavior as
described in [35].

##  V Conclusion

Using a tight-binding effective low Hamiltonian and considering a constant
periodic potential, we have examined the electrical structure and the
corresponding contact points of a phosphorene superlattice. In the beginning,
we derived the dispersion relation governed by an equation resulting from the
boundary conditions. We were able to extract a linear behavior of the current
system near the contact points and find different velocity components
$(v_{x},v_{y})$ along the $x$\- and $y$-directions by inspecting this
relation. To highlight the fundamental feature of the system, we examined
$(v_{x},v_{y})$ for the original Dirac point $(0,0,E_{0})$. As a result, we
demonstrated that the barrier height and width can modulate $v_{y}$,
exhibiting different oscillations with small amplitudes up to some values of
these barrier parameters, whereas $v_{x}$ decreased toward a constant.

Two scenarios were studied: one with equal well and barrier widths while
maintaining a constant cell width, and the other with unequal well and barrier
widths. It was demonstrated that altering the well and barrier widths or
heights can be used to adjust the energy gap of the phosphorene superlattice.
When the barrier height is increased, the Dirac points are moved up in two
circumstances. It was discovered that as the barrier height or the well and
barrier widths increased, the energy gap shrank in both scenarios. The two
situations are not identical, though. When the barrier height $V_{B}$ exceeds
a critical value determined by the barrier width, Dirac points appear at
$k_{y}\neq 0$.Thus, the phosphorene superlattice becomes gapless under certain
conditions. We demonstrated that the spectrum is symmetric in both cases when
$k_{y}=0$. As a result, to control and modify the contact points in the
phosphorene superlattice, one may adjust either the well and barrier widths or
heights.

Dirac points for fermions in phosphorene under certain constraints are very
significant. In fact, the conduction and valence bands collide to generate
zero-energy points, whose positions are known as Dirac points in a material’s
band structure. This enables effective charge transfer because electrons can
readily flow between the two bands. As for the phosphorene superlattice, these
Dirac points are particularly intriguing for uses in electronics and
optoelectronics, as well as having potential uses in quantum computing and
communication technologies. They can also be employed as waveguides for
microwave or light communication systems because of their anisotropic
characteristic. Additionally, they can be utilized as energy harvesters in
thermoelectric devices due to their high density of states at the Fermi level.
Finally, Dirac points in phosphorene offer insight into basic physics
phenomena like topological insulators and Majorana fermions in addition to
these possible applications.

## References

  * [1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science 306, 666 (2004).
  * [2] F. Zhai and K. Chang, Phys. Rev. B 85, 155415 (2012).
  * [3] N. M. R. Peres, A. H. Castro Neto, and F. Guinea, Phys. Rev. B 73, 241403(R) (2006).
  * [4] M. I. Katsnelson, K. S. Novoselov, and A. K. Geim, Nat. Phys. 2, 620 (2006).
  * [5] C. W. J. Beenakker, Rev. Mod. Phys. 80, 1337 (2008).
  * [6] J. Reijniers, F. M. Peeters, and A. Matulis, Phys. Rev. B 64, 245314 (2001).
  * [7] S. Park and H. S. Sim, Phys. Rev. B 77, 075433 (2008).
  * [8] F. Zhai and K. Chang, Phys. Rev. B 77, 113409 (2008).
  * [9] M. Galluzzi, Y. Zhang, X. F. Yu, J. Appl. Phys. 128, 230903 (2020).
  * [10] X. J. Li, J. H. Yu, K. Luo, Z. H. Wu, and W. Yang, Nanotechnology 29, 174001 (2018).
  * [11] J. R. Brent, N. Savjani, E. A. Lewis, S. J. Haigh, D. J. Lewis, and P. O’Brien, Chemical Communications 50, 13338 (2014).
  * [12] Z. Guo, H. Zhang, S. Lu, Z. Wang, S. Tang, J. Shao, Z. Sun, H. Xie, H. Wang, X.-F. Yu, and K. Chu, Adv. Funct. Mater. 25, 6996 (2015).
  * [13] A. Carvalho, M. Wang, X. Zhu, A. S. Rodin, H. Su, and A. H. C. Neto, Nat. Rev. Mater. 1, 16061 (2016).
  * [14] R. Gusmao, Z. Sofer, and M. Pumera, Angewandte Chemie International Edition 56, 8052 (2017).
  * [15] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, Rev. Mod. Phys. 81, 109 (2009).
  * [16] M. Ezawa, New J. Phys. 16, 115004 (2014).
  * [17] S. Lin, Y. Li, J. Qian, and S. P. Lau, Materials Today Energy 12, 1 (2019).
  * [18] Y. Yi, Z. Sun, J. Li, P. K. Chu, and X. F. Yu, Small Methods 3, 1900165 (2019).
  * [19] Z. Li, T. Cao, M. Wu, and S. G. Louie, Nano letters 17, 2280 (2017).
  * [20] J. M. Pereira, Jr. and M. I. Katsnelson, Phys. Rev. B 92, 075437 (2015).
  * [21] A. N. Rudenko and M. I. Katsnelson, Phys. Rev. B 89, 201408(R) (2014).
  * [22] B. Huard, J. A. Sulpizio, N. Stander, K. Todd, B. Yang, and D. Goldhaber-Gordon, Phys. Rev. Lett. 98, 236803 (2007).
  * [23] Z. S. Popović, J. M. Kurdestany, and S. Satpathy, Phys. Rev. B 92, 035135 (2015).
  * [24] J. Seﬀadi, I. Redouani, Y. Zahidi, and A. Jellal, Solid State Commun. 351, 114777 (2022).
  * [25] J. Q. You, L. Zhang, and P. K. Ghosh, Phys. Rev. B 52, 17243 (1995).
  * [26] S. Soleimanikahnoj and I. Knezevic, J. Comput. Electron 16, 568 (2017).
  * [27] Zhesheng Chen, Jingwei Dong, Christine Giorgetti, Evangelos Papalazarou, Marino Marsi, Zailan Zhang, Bingbing Tian, Qingwei Ma, Yingchun Cheng, Jean-Pascal Rueff, Amina Taleb-Ibrahimi, and Luca Perfetti, 2D Mater. 7, 035027 (2020).
  * [28] J. R. Lima, J. Appl. Phys. 117, 084303 (2015).
  * [29] J. R. Lima, Phys. Lett. A 379, 1372 (2015).
  * [30] B. Ghosh, B. Singh, R. Prasad, and A. Agarwal, Phys. Rev. B 94, 205426 (2016).
  * [31] C. H. Pham, H. C. Nguyen, and V. L. Nguyen, J. Phys.: Condens. Matter 22, 425501 (2010).
  * [32] C. H. Pham and V. L. Nguyen, J. Phys.: Condens. Matter 26, 425502 (2014).
  * [33] G. M. Maksimova, E. S. Azarova, A. V. Telezhnikov, and V. A. Burdov, Phys. Rev. B 86, 205422 (2012).
  * [34] S. Ono, J. Appl. Phys. 121, 204301 (2017).
  * [35] Y. Lu, D. Zhou, G. Chang, S. Guan, W. Chen, Y. Jiang, and H. Lin, npj Computational Materials. 2, 16011 (2016).


#  What is the state of the art? Accounting for multiplicity in machine
learning benchmark performance

Kajsa Møllersen1 and Einar Holsbø2

(1Department of Community Medicine
2Department of Computer Science
UiT - The Arctic University of Norway
)

##  1 Introduction: Population vs sample state-of-the-art performance

Machine learning methods are commonly evaluated and compared by their
performance on data sets from public repositories, where methods can be
evaluated under identical conditions and across time. Theoretical results play
an important role in machine learning research, but for the last couple of
decades the success of increasingly complex methods has pushed the proof
towards the pudding: performance on benchmark data sets.

The numbers don’t lie: if method A performs better than method B on a specific
data set, that is what it did. But the truth that the numbers tell is limited.
We think of a data set, for example a collection of images, as a random sample
from an underlying population. In science we are interested in the performance
on the population, and the observed performance acts as an estimate. A good
estimator gives a reliable result, but overfitting to the data set,
originating from wrong use of cross-validation, data leakage and data reuse,
can give overly optimistic estimates. An independent test set, often referred
to as a hold-out set, is a sample that is not used for development of the
method, and gives an unbiased estimator. Competitions and challenges with
prize money typically withhold the test set so that it becomes impossible for
the participants to use it for development.

Although we know these things we have a tendency to extrapolate the
performance metric and interpret it as something about the world rather than
something about a data set. If we always appended “on MNIST” when we said so-
and-so method had accuracy so-and-so nobody would care. Who cares about MNIST?
We hold out data precisely because we want to estimate some sort of long-term
behavior on similar data.

An unbiased estimator solves only half of the problem; the variance tells just
as much as the bias regarding the quality of a method. For machine learning
methods, the sources of variance are plentiful. One is the random nature of
the test set, the training set and the validation set. In addition comes the
setting of hyperparameters, which can lead to substantial variation, and minor
variance sources such as random seeds and stochastic dropouts. The sources of
variance that are part of the training are often referred to as influencing
the method’s stability or robustness. With knowledge about the behaviour of
both the bias and the variance of an estimator, the error can be predicted,
often in terms of a confidence interval. It is certainly more useful to know
what the lower limit of the 95% confidence interval is, than the observed
performance on a specific test set.

The highest ranked performance on a problem is referred to as state-of-the-art
(SOTA) performance, and is used, among other things, as a reference point for
publication of new methods. As a specific method’s performance on a test set
is just an estimate of its true underlying performance, the reported SOTA
performances are estimates, and hence come with bias and variance.

State-of-the-art performance increasingly relies on benchmark data sets, often
in combination with a challenge or a competition. The practice of hold-out
test sets in competitions safe-guards against overoptimistic estimates for the
individual classifiers. Data is often released after the competition has
ended, and new methods can be compared to the state-of-the-art performance.

As we will show in the subsequent sections, using the highest-ranked
performance as an estimate for state-of-the-art is a biased estimator, giving
overoptimistic results. The mechanisms at play are those of multiplicity, a
topic that is well-studied in the context of multiple comparisons and multiple
testing, but has, as far as the authors are aware of, been nearly absent from
the discussion regarding state-of-the-art estimates.

Biased state-of-the-art has consequences both for practical implementations
and for research itself. A health care provider might want to use some machine
learning method for tumour detection, but only if the method correctly
classifies at least 80%. The health care provider closely follows updates on
state-of-the-art performance, and are ready to invest time and money in a
pilot project once the lower bound surpasses 80%. If the state-of-the-art
performance is an optimistic estimate, the health care provider has now lost
time and money, and perhaps also some faith in research results.

For the research community, the situation is perhaps even more dramatic. The
optimistic state-of-the-art estimate is used as a standard for evaluating new
methods. Not to say that a new method necessarily needs to beat the current
state-of-the-art to catch interest from reviewers and editors, but it
certainly helps if its performance is comparable. It has been pointed out (see
Section 2) that the incremental increase in performance can obstruct the
spread of new ideas that have substantial potential, but just haven’t reached
their peak yet.

In this article, we provide a probability distribution for the case of
multiple classifiers. We demonstrate the impact of multiplicity through
simulated examples. We show how classifier dependency impacts the variance,
but also that the impact is limited when the accuracy is high. Finally, we
discuss a real-world example; a Kaggle competition from 2020. All code is
available at https://github.com/3inar/ninety-nine.

##  2 Related work

Multiple comparisons is a well-studied subject within statistics, see, e.g.,
Tukey (1991) for an interesting discussion on some topics of multiple
comparisons, among them the confidence intervals. For a direct demonstration
of how multiple testing can lead to false discoveries, Bennett et al. (2009)
gives an example where brain activity is detected in a (dead) salmon “looking”
at photos of people. In the machine learning context, with multiple data sets
and classifiers, Demšar (2006) suggests adjusting for multiple testing,
according to well-known statistical theory. Our article focuses on the
estimate itself and its uncertainty, and not comparison of methods.

The re-use of benchmark datasets has been addressed in several publications.
Salzberg (1997) point to the problem that arises because these datasets are
static, and therefore run the risk of being exhausted by overfitting.
Exhausted public data sets still serve important functions, e.g., a new idea
can show to perform satisfactory, even when it is not superior. Thompson et
al. (2020) distinguishing between simultaneous and sequential multiple
correction procedures. In sequential testing, which is the reality of public
datasets, the number of tests is constantly updated. An interesting topic they
bring up is the conflict between sharing incentive, open access and stable
false positive rate. There are several demonstrations of the vulnerability of
re-use of benchmark datasets. Teney et al. (2020) performed a successful
attack on the VQA-CP benchmark dataset, that is, they are able to attain high
performance with a non-generalisable classifier. They point towards the risk
of ML methods capturing idiosyncrasies of a dataset. Torralba and Efros (2011)
introduced the game Name That Dataset! motivating their concern of methods
being tailored to a dataset rather than a problem. The lack of statistically
significant differences among top performing methods is demonstrated by, e.g.,
Everingham et al. (2015) and Fernández-Delgado et al. (2014).

The same mechanisms of overfitting are present in challenges through multiple
submissions to the validation set. Blum and Hardt (2015) introduced an
algorithm that releases the public leaderboard score only if there has been
significant improvement compared to the previous submission. This will prevent
overfitting by hindering adjustments to random fluctuations in the score, and
not true improvements of the method. Dwork et al. (2017) use the principles of
differential privacy to prevent overfitting to the validation set.

Several publications show that the overfitting is not as bad as expected, and
provide possible explanations. Recht et al. (2019) constructed new,
independent test sets of CIFAR-10 and ImageNet. They found that the accuracy
drop is dramatic, around 10%. However, the order of the best performing
methods is preserved, which can be explained by the fact that multiple classes
(Feldman et al. (2019)) and model similarity (Mania et al. (2019)) slows down
overfitting. Roelofs et al. (2019) compared the results from the private
leaderboard (one submission per team) and the public leaderboard (several
submissions per team) in Kaggle challenges. Several of their analyses indicate
overfitting even though their main conclusion is that there is little evidence
of overfitting.

Overfitting and its consequences are avoided in most challenges by withholding
the test set so that a method can only be evaluated once and not adapted to
the test set performance. Hold-out datasets also prevents wrong use of cross-
validation, see, e.g., Fernández-Delgado et al. (2014) and the seminal paper
of Dietterich (1998). Another contribution from public challenges is the
requirement of reproducibility to collect the prize money, which is a well-
known concern in science, see, e.g., Gundersen and Kjensmo (2018). Public
challenges typically rank the participants by a (single) performance measure,
and can fall victims of Goodhart’s law, paraphrased as “When a measure becomes
a target, it ceases to be a good measure.” by Teney et al. (2020). Ma et al.
(2021) offers a framework and platform for evaluation of NLP methods, where
other aspects such as memory use and robustness are evaluated.

The robustness/stability of a method, and the variability in performance are
crucial aspects that the hold-out test set of public challenges does not offer
an immediate solution to. Bouthillier et al. (2021) showed empirically that
hyperparameter choice and the random nature of the data are two large
contributors to variance, whereas data augmentation, dropouts, and weights
initialisations are minor contributors. Choice of hyperparameters are often
adjusted to validation set results, and the strategies of Blum and Hardt
(2015) and Dwork et al. (2017) can hinder that, if implemented by the
challenge hosts. Our contribution is a focus on better estimate for state-of-
the-art performance under the current conditions, where the approach of
Bouthillier et al. (2021) does not apply. Bousquet and Elisseeff (2002)
investigate how sampling randomness influences accuracy, both for regression
algorithms and classification algorithms, relying on the work of Talagrand
(1996) who stated that “A random variable that depends (in a “smooth” way) on
the influence of many independent variables (but not too much of any of them)
is essentially constant.”

There is a rich literature on various aspects regarding the use of public
challenges as standards for scientific publications, see, e.g., Varoquaux and
Cheplygina (2022) for a much broader overview than the one given here.
Concerns and solutions regarding multiplicity and better performance estimates
are active research fields. We contribute to this conversation by shining the
light on how multiplicity of classifiers’ performance on hold-out test sets
create optimistic estimates of state-of-the-art performance.

##  3 Multiple classifiers and biased state-of-the-art estimation

###  3.1 Notation

$\displaystyle n$ $\displaystyle:\text{number of trials/size of test set}$
$\displaystyle\theta$ $\displaystyle:\text{probability of success/probability
of correct prediction}$ $\displaystyle m$ $\displaystyle:\text{number of
experiments/number of classifiers}$ $\displaystyle\mathcal{B}ern(\theta)$
$\displaystyle:\text{Bernoulli distribution}$ $\displaystyle Y$
$\displaystyle:\text{random variable, indicates if prediction is correct}$
$\displaystyle\,\,\,\,Y\sim\mathcal{B}ern(\theta)$
$\displaystyle\mathcal{B}(n,\theta)$ $\displaystyle:\text{binomial
distribution}$ $\displaystyle X$ $\displaystyle:\text{random variable, number
of failures in one experiment}$
$\displaystyle\,\,\,\,X\sim\mathcal{B}(n,1\text{-}\theta)$
$\displaystyle\hat{\theta}(X)=\frac{n-X}{n}$ $\displaystyle:\text{estimator
for }\theta,\text{ unbiased}$ $\displaystyle P(X=x)$
$\displaystyle:\text{probability of $x$ failures in one experiment}$
$\displaystyle P_{x}=P(X\leq x)$ $\displaystyle:\text{probability of at most
$x$ failures in one experiment}$
$\displaystyle\mathbf{X}=\\{X_{1},X_{2},\ldots,X_{m}\\}$
$\displaystyle:\text{random sample, number of failures for classifiers
}j=1,\ldots,m$ $\displaystyle C_{x}$ $\displaystyle:\text{random variable,
number of experiments with at most $x$ failures}$
$\displaystyle\,\,\,\,C_{x}\sim\mathcal{B}(m,P_{x})$ $\displaystyle Z$
$\displaystyle:\text{random variable, the smallest number of failures among
all classifiers}$ $\displaystyle\theta_{SOTA}$ $\displaystyle:\text{highest
probability of correct prediction among all classifiers}$
$\displaystyle\hat{\theta}_{SOTA}(\mathbf{X})=\max_{j}\hat{\theta}_{j}$
$\displaystyle:\text{estimator for }\theta_{SOTA},\text{ biased}$
$\displaystyle\alpha$ $\displaystyle:\text{significance level}$

###  3.2 Two coin-flip examples

Multiplicity is, simply put, how the probability of an outcome changes when an
experiment is performed multiple times. Consider $n=20$ flips of a fair coin
where the outcome is the number of heads, referred to as the number of
successes in a binomial distribution, $\mathcal{B}(n,\theta)$, where the
probability of success is $\theta=0.5$. Let the random variable $X$ denote the
number of failures, $X\sim\mathcal{B}(n,1\text{-}\theta)$. The probability of
$15$ heads or more corresponds to observing $x=n-15$ tails or less, that is,
the cumulative distribution $P(X\leq x|n,1\text{-}\theta)$. This corresponds
to the cumulative distribution for the estimator
$\hat{\theta}(X)=\frac{n\text{-}X}{n}$, and we have that
$P(\hat{\theta}(X)\geq\frac{n\text{-}x}{n})=P(X\leq x)=0.02069$ for
$\hat{\theta}\geq 0.75$. In other words, if $\theta=0.5$, there is a
probability of $0.02069$ that the estimate is $0.75$ or above in a series of
$n=20$ trials. The probability of exactly $x=5$ can be calculated from the
probability mass function (pmf) $P(X=x|n,1\text{-}\theta)=0.01479$, but the
interest usually lies with the cumulative distribution. We focus on the case
where $X\leq x$, but it is straightforward to substitute $P(X\leq x)$ for
$P(X=x)$ in the following.

If we fix $x$ we can define the experiment described above as a Bernoulli
trial, where success is defined as $X\leq x$, with corresponding probability
$P(X\leq x)$, which we from now on write as $P_{x}$ for ease of notation. If
we perform the coin-flip experiment $m$ times, and success is defined as
$X\leq x$, we have a series of $m$ Bernoulli trials, and the number of
successes follow the binomial distribution $\mathcal{B}(m,P_{x})$. Let $C_{x}$
be the random variable that denotes the number of successful experiments. The
probability of, for example, $C_{x}=10$ experiments with at most $x=5$ tails,
is calculated from the pmf $P(C_{x}=10|m,P_{x})$. Of special interest is the
case where at least one experiment succeeds, $P(C_{x}>0)=1-P(C_{x}=0)$. In our
example, if the coin-flip experiment is performed $m=100$ times, the
probability of observing $15$ heads or more at least once is
$P(C_{x}>0|m,P_{x})=1-P(C_{x}=0|m,P_{x})$, with $P_{x}=0.02069$ from above. In
the estimation context, we have that although the probability of
$\hat{\theta}\geq 0.75$ is only $0.02069$ for the individual experiments, the
probability of at least one $\hat{\theta}\geq 0.75$ is $0.8765$ when the
experiment is performed $m=100$ times.

Consider an experiment where a coin is flipped $n=20$ times, and the number of
successes is the number of times a classifier correctly predicts the outcome.
An unbiased estimator for a classifier’s probability of correct prediction,
$\theta$, is ${\hat{\theta}=(n-X)/n}$, where $X$ is the number of failures.
The corresponding estimate is referred to as accuracy. The true state-of-the-
art for coin-flip classifiers is $\theta_{SOTA}=0.5$, since it is impossible
to train a classifier to predict coin flips. If $m=1,000$ classifiers
independently give their predictions for $n=20$ flips, we can consider this as
a binomial experiment with $1,000$ Bernoulli trials, where the probability of
success is the individual classifier’s probability of $n-X\leq n-x$ correct
predictions. The probability of at least one classifier achieving an accuracy
of at least $0.90$ is $1-P(C_{x}=0|m,P_{x})=0.1823$, where $P_{x}=P(X\leq
18|n,1\text{-}\theta)=0.00020$. This means that using the top-ranked method as
an estimate and ignoring multiplicity there is a substantial risk of
overestimating the state-of-the-art.

###  3.3 The probability distribution of $\hat{\theta}_{SOTA}$

We define the SOTA probability of correct prediction, $\theta_{SOTA}$, as the
highest probability of correct prediction among all classifiers applied to a
test set, in accordance with the common set-up for public competition. Of
interest is $P(\hat{\theta}_{SOTA}|m,n,\theta)$, the probability of
$\hat{\theta}_{SOTA}$, being the estimated state of the art from taking the
highest ranked classifier, given the number of classifiers, the size of the
test set, and the probabilities of correct prediction,
$\theta=\\{\theta_{1},\ldots,\theta_{m}\\}$.

We will in the following assume independent classifiers, and discuss
dependency in Section 4. We demonstrate the impact of multiplicity in the
simple case of identical probabilities of correct prediction,
$\theta_{j}=\theta$ for all classifiers, $j=1,\ldots,m$, and where
$\theta_{SOTA}=\theta$. We discuss non-identical classifiers in Section 4.2.

Let $n$ denote the size of the test set, $m$ be the number of classifiers, and
$\theta$ be the probability of correct prediction. Let $X$ be the random
variable that denotes the number of failures (wrong predictions) for a given
classifier. We have that $X\sim\mathcal{B}(n,1\text{-}\theta)$. We denote by
$P_{x}$ the cumulative distribution $P(X\leq x)$. The random sample,
$\mathbf{X}=\\{X_{1},X_{2},\ldots,X_{m}\\}$ consists of the number of wrong
predictions for each classifier $j=1,\ldots,m$ on a test set. Let $C_{x}$ be
the random variable that denotes the number of classifiers with at most $x$
failures. We have that $C_{x}\sim\mathcal{B}(m,P_{x})$, because each
classifier will independently have a probability $P_{x}$ of having at most $x$
wrong predictions.

We are interested in the probability of at least one classifier having at most
$x$ wrong predictions, because it corresponds to the probability function of
the $\hat{\theta}_{SOTA}$ estimator, $P(\hat{\theta}_{SOTA}>\frac{n-x}{n})$,
in other words, at least one $\hat{\theta}_{j}>\frac{n-x}{n}$. We can derive
this from $C_{x}$:

$\displaystyle P(C_{x}>0|m,n,\theta).$ (1)

Note that this corresponds to the cumulative distribution function (cdf)
$F(z)$ for the random variable $Z$, denoting the most wrong predictions among
all classifier. We have cdf $F(z)$, the probability that the most wrong
predictions among all classifiers is less than or equal to a number $z$, can
easily be derived from the binomial distribution

$\displaystyle F(z)$
$\displaystyle=P(C_{z}>0|m,n,\theta)=1-P(C_{z}=0|m,P_{z})$
$\displaystyle=1-(1-P_{z})^{m}.$ (2)

Let $f(z)$ denote the probability mass function of $Z$. We then have

$\displaystyle f(z)$
$\displaystyle=F(z)-F(z-1)=1-(1-P_{z})^{m}-\\{1-(1-P_{z-1})^{m}\\}$
$\displaystyle=(1-P_{z-1})^{m}-(1-P_{z})^{m}=$
$\displaystyle=(1-\sum_{\ell=0}^{z-1}P(Z=\ell))^{m}-(1-\sum_{\ell=0}^{z-1}P(Z=\ell)-P(Z=z))^{m}$
$\displaystyle=-\sum_{k=1}^{m}(-1)^{k}\binom{m}{k}(1-\sum_{\ell=0}^{z-1}P(Z=\ell))^{m-k}P(Z=z)^{k}.$
(3)

The last line is derived by binomial expansion. The practical relevance of the
binomial expansion is limited by the need to approximate binomial coefficient
for large $m$. From Eq. 3.3, the expected value and variance can easily be
calculated. An example of $F(z)$ can be seen in Fig. 2, and the corresponding
$f(z)$ in Fig. 3. The pmf $f(z)=P(Z=z)$ denotes the probability of at least
one classifier having exactly $z$ wrong predictions, and no classifiers having
fewer than $z$ wrong predictions.

###  3.4 A simulated public competition example

The multiplicity mechanisms described above apply to public competitions where
several teams submit their predictions on a test set. Winning teams are those
with highest accuracy, $\max_{j}\hat{\theta}_{j}$, where
$\hat{\theta}_{j}=(n-x_{j})/n$ is the accuracy of a classifier $j$ on the test
set. This is often referred to as the state of the art, and hence, we have
that $\hat{\theta}_{SOTA}(\mathbf{x})=\max_{j}\hat{\theta}_{j}$. As
exemplified in Section 3.2, multiplicity biases the estimation of
$\theta_{SOTA}$, and Section 3.3 gives us the tools to analyze its impact. In
the following, we will use a significance level of $\alpha=0.05$ and two-sided
confidence intervals.

Consider a classification problem with a test set of size $n=3,000$, and a
classifier with probability of correct prediction $\theta=0.90$. The pmf of a
single classifier is shown as a histogram in Fig. 1. The 95% confidence
interval for $\theta$ if we have observed $\hat{\theta}=0.90$ is
$(0.8887,0.9105)$, and we denote it by
$(\theta_{\alpha/2},\theta_{1\text{-}\alpha/2})$.

If there are $m=1,000$ classifiers, each with probability of correct
prediction $\theta=0.90$, the probability of at least one team achieving at
least $\theta_{1\text{-}\alpha/2}$ is close to 1. This comes as no surprise
given the definition of a confidence interval. It means that we can be almost
certain that the accuracy of the top ranked classifier, which is the default
SOTA, is above the upper limit of the confidence interval of the true
performance.

On the other hand, it is straightforward to show that if each individual
classifier has a probability of $P_{x}=0.00002532$ for at most $x$ failures,
then, with a probability of $\alpha/2$, at least one team will achieve at
least $\hat{\theta}=0.9213$. We denote this by
$\theta^{m}_{1\text{-}\alpha/2}$, the multiplicity-adjusted upper limit of the
$(1\text{-}\alpha)$ confidence interval.

Similarly, the expected value of $\hat{\theta}_{SOTA}$ can easily be
calculated from Eq. 3.3, and we have $E(\hat{\theta}_{SOTA})=0.9173$. Both are
optimistic estimates, knowing that $\theta_{SOTA}$ is $0.90$.

Within the competition context, this is not necessarily a problem, because
each team has the same probability of achieving the highest accuracy given
that they have the same $\theta$. In a scientific context, using
$\max_{j}\hat{\theta}_{j}$ as an estimate for $\theta_{SOTA}$ without
multiplicity corrections is problematic. Consider a new classifier with
probability of correct prediction $\theta_{m+1}=\theta_{1\text{-}\alpha/2}$
from above, displayed as the red curve in Fig. 1. By common scientific
standards, the expected performance, $E(\hat{\theta}_{m+1})=\theta_{m+1}$ is
considered significantly better than the expected performance of any of the
$1,000$ classifiers with $\theta=0.90$, and is of interest to the research
community. But the probability of this classifier beating the upper limit of
the multiplicity adjusted confidence interval is very low:
$P(\hat{\theta}_{m+1}\geq\theta^{m}_{1\text{-}\alpha/2})=0.0190$, and can be
seen as the red area in Fig. 1. Similarly, $P(\hat{\theta}_{m+1}\geq
E(\hat{\theta}_{SOTA}))=0.0996$, seen as the blue plus red area in Fig. 1.
This means that by using the maximum accuracy as SOTA, we can expect that 9
out of 10 classifiers that are in fact better will still be naively considered
worse.

The probability mass functions of single $\hat{\theta}$.

Figure 1: The grey histogram shows the pmf of a classifier with probability of
correct prediction $\theta$. The red curve is the pmf of a classifier with
probability of correct prediction $\theta_{1\text{-}\alpha/2}$, the upper
limit of the $95\%$ CI for $\hat{\theta}(\mathbf{x})=\theta$.
$E(\hat{\theta}_{SOTA})$ denotes the expected value of the highest ranked
among $m=1,000$ classifiers. $\theta^{m}_{1\text{-}\alpha/2}$ is the
multiplicity corrected upper limit of the $95\%$ CI. The shaded areas
corresponds to $P(\hat{\theta}_{m+1}\geq\theta^{m}_{1\text{-}\alpha/2})$ (red)
and $P\big{(}\hat{\theta}_{m+1}\geq E(\hat{\theta}_{SOTA})\big{)}$ (blue+red).

The cumulative distribution function, $F(z)$ from Eq. 3.3, is displayed in
Fig. 2, and describes the cumulative probability of at least one classifier
having at most $z$ failures. Or, in other words, the cumulative probability of
maximum accuracy. The number of failures are displayed on the top horizontal
axis, and the corresponding accuracy on the bottom horizontal axis. Note that
the horizontal axis is reversed compared to Fig. 1.

Cumulative distribution function, $F(z|m,n,\theta)$, for at least one
classifier having at most $z$ failures.

Figure 2: The figure shows the cumulative distribution function for at least
one classifier having at most $z$ failures. The horizontal axis denotes the
number of failures, $z$, on top and the accuracy at the bottom.

The corresponding pmf, $f(z)$ from Eq. 3.3, is displayed as a histogram in
Fig. 3.

Probability mass function, $f(z|m,n,\theta)$, for the smallest number of
failures, $z$, among all classifiers.

Figure 3: The figure shows the corresponding pmf for Fig. 2. The horizontal
axis denotes the number of failures, $z$, on top and the corresponding
accuracy at the bottom.

From Eq. 3.3 we get the expected value and variance. For $m=1,000$, $n=3,000$
and $\theta=0.9$ this corresponds to $E(\hat{\theta}_{SOTA})=0.9173$ with a
standard deviation of $0.001817$. We can see how different values for $m$, $n$
and $\theta$ affects the $\hat{\theta}_{SOTA}$ distribution in Table 1.

Expected values and standard deviations for $\hat{\theta}_{SOTA}$

m | n | $\theta$ | expected value | standard deviation
---|---|---|---|---
1,000 | 3,000 | 0.85 | 0.8707 | 0.002197
1,000 | 3,000 | 0.90 | 0.9173 | 0.001817
1,000 | 3,000 | 0.95 | 0.9624 | 0.001277
100 | 3,000 | 0.90 | 0.9135 | 0.002250
500 | 3,000 | 0.90 | 0.9163 | 0.001923
1,000 | 1,000 | 0.90 | 0.9294 | 0.003007
1,000 | 10,000 | 0.90 | 0.9096 | 0.001022

Table 1: The expectation and standard deviation for the SOTA estimator with
varying test set size, number of classifiers and probability of correct
prediction. The bold line is the example used in Fig. 3.

##  4 Non-i.i.d.

The probability distribution in Section 3.4 is based on independent
classifiers with identical $\theta$s. This assumption is rarely true in
practice. A classifier’s ability to perform better than random guesses is due
to correlation between features and labels in a training set. The teams will
pick up on many of the same features and correlations, regardless of which
method they choose. In addition, many teams will use different versions of
related methods. We investigate dependency and non-identical classifiers, and
the impact on the bias. Table 2 gives an overview, details are in the
following sections.

Upper bounds, expected values and standard deviations for
$\hat{\theta}_{SOTA}$

section | independent | identical | upper 95% CI | expected value | standard deviation
---|---|---|---|---|---
3.4 | YES | YES | 0.9213 | 0.9173 | 0.001817
4.1 | NO | YES | 0.9207 | 0.9140 | 0.003481
4.1 | NO∗ | YES | 0.9173 | 0.9140 | 0.001484
4.2 | YES | NO | 0.9177 | 0.9130 | 0.002129
4.3 | NO | NO | 0.9173 | 0.9101 | 0.003649

Table 2: The upper bound of the 95% CI, the expectation and standard deviation
for the SOTA estimator for different combinations of independency and
identicality.

###  4.1 Dependent, identical classifiers

The dependency between two classifiers $j,k$ can be described by their
correlation, $\rho_{jk}$. Let $Y_{j}$ be the Bernoulli variable with value $1$
if classifier $j$ predicts a data point correctly. We have
$P(Y_{j}=1)=\theta_{j}$, and $\theta_{j}=\theta$ for classifiers with equal
probability of success. Mania et al. (2019) calculated the similarity, defined
as the fraction of two models giving the same output, which is equivalent to

$\displaystyle
P\big{(}\\{Y_{j},Y_{k}\\}=\\{0,0\\}\big{)}+P(\\{Y_{j},Y_{k}\\}=\\{1,1\\})$ (4)
$\displaystyle=(1-\theta_{j})(1-\theta_{k})+\theta_{j}\theta_{k}$

for independent models. For dependent models with correlation $\rho_{jk}$,
this becomes

$\displaystyle(1-\theta_{j})(1-\theta_{k})+\theta_{j}\theta_{k}+2\rho_{jk}\sqrt{\theta_{j}\theta_{k}(1-\theta_{j})(1-\theta_{k})}.$
(5)

By setting $\theta_{j}=\theta$ and $\rho_{jk}=\rho$, $\rho$ can easily be
estimated from the estimations of $\theta$,
$P\big{(}\\{Y_{j},Y_{k}\\}=\\{0,0\\}\big{)}$ and
$P(\\{Y_{j},Y_{k}\\}=\\{1,1\\})$. In Mania et al. (2019), they used the
average accuracy for top-performing ImageNet models, $\hat{\theta}=0.756$. For
details, see Mania et al. (2019). Note that these are multi-class problems, so
a lower accuracy than for binary classification problems is to be expected.
From their calculated similarity of $0.85$, the estimated correlation is
$0.36$.

Boland et al. (1989) describes dependency in the majority systems context,
with $\rho_{0j}=\text{corr}(Y_{0},Y_{j}),\,j=1,\ldots m$, and the $Y_{j}$s are
independent given $Y_{0}$. This gives $\rho_{jk}=\rho_{0j}^{2}$. The
observations from $Y_{0}$ does not contribute directly to
$\hat{\theta}_{SOTA}$, but can be seen as, e.g., the influence that a common
training set has on the classifiers. For $\rho_{0j}=1$, we have $m$ identical
classifiers, and an unbiased E($\hat{\theta}_{SOTA})=\theta$. For
$\rho_{0j}=0$, we have $m$ independent classifiers, and the situation is
identical to the one described in Section 3.3.

The binomial distribution for dependent Bernoulli trials is not an easy
problem, see, e.g., Ladd (1975) and Hisakado et al. (2006). We here
demonstrate by simulation how dependency reduces the bias. We use the
estimated correlation from Mania et al. (2019), which corresponds to
$\rho_{0j}=0.6$ in the multiplicity set-up with an external $Y_{0}$.

We use the set-up from Section 3.4 with $m=1,000$ classifiers, test set size
$n=3,000$, and probability of correct prediction $\theta=0.90$. We simulate
two scenarios; one where $Y_{0}$ is random, and one where the $y_{0}$’s are
fixed with exactly $(n-x)/n=0.90$ correct predictions. The fixed scenario can
correspond to fixed training conditions, e.g., models pre-trained on ImageNet.
The random scenario picks up the randomness in, e.g., a training set. Although
all teams have the same training set, the training set is a random sample.

The upper bound of the $95\%$ confidence interval for the dependent
classifiers is at $0.9173$ for the fixed $y_{0}$ and $0.9207$ for the random
$Y_{0}$, but they have the same expected value.

Minimum number of failures in $m=1,000$ classifiers.

(a) $\rho_{0j}=0.6$, random $Y$

(b) $\rho_{0j}=0.6$, fixed $y$

Figure 4: The minimum number of failures in $100,000$ simulations, with
probability of success $\theta=0.9$.

Note that, with $\theta=0.90$, the similarity as described by Mania et al.
(2019) is $0.82$ for two independent classifiers. When the probability of
correct prediction is high, there is little room for variation between
classifiers, and the inherent similarity is high.

###  4.2 Non-identical, independent classifiers

The Poisson Binomial distribution, see, e.g., Wang1993, describes the sum of
$n$ non-identical independent Bernoulli trials with success probabilities
$p_{1},p_{2},\ldots,p_{n}$. The pmf for $K$ number of successes is

$\displaystyle P(K=k)=\sum\limits_{{A\in G_{k}}}\prod\limits_{{i\in
A}}p_{i}\prod\limits_{{j\in A^{c}}}(1-p_{j}),$ (6)

where $G_{k}$ is all subsets that can be selected from {1, 2, …, n}.

In the multiple classifier context, we are interested in generalising Eq. 3.3
to the non-identical case.

$\displaystyle
F(z)=P(C_{z}>0|m,n,\mathbf{\theta})=1-P(C_{z}=0|m,n,\mathbf{\theta}),$ (7)

where $P(C_{z}=0)$ corresponds to $0$ classifiers with non-identical
probabilities of success, $P(Z\leq z|\theta_{j})$. The cdf is

$\displaystyle F(z)$ $\displaystyle=1-\sum\limits_{{A\in
G_{0}}}\prod\limits_{{i\in A}}P(Z\leq z|\theta_{i})\prod\limits_{{j\in
A^{c}}}(1-P(Z\leq z|\theta_{j}))$
$\displaystyle=1-\prod\limits_{j=1}^{m}\big{(}1-P(Z\leq z|\theta_{j})\big{)}.$
(8)

The pmf is

$\displaystyle f(z)$ $\displaystyle=F(z)-F(z-1)$
$\displaystyle=1-\prod\limits_{j=1}^{m}\big{(}1-P(Z\leq
z|\theta_{j})\big{)}-\\{1-\prod\limits_{j=1}^{m}\big{(}1-P(Z\leq
z\text{-}1|\theta_{j})\big{)}\\}$
$\displaystyle=\prod\limits_{j=1}^{m}\big{(}1-P(Z\leq
z\text{-}1|\theta_{j})\big{)}-\prod\limits_{j=1}^{m}\big{(}1-P(Z\leq
z|\theta_{j})\big{)}$ (9)
$\displaystyle=\prod\limits_{j=1}^{m}\big{(}1-\sum\limits_{k=0}^{z\text{-1}}\binom{n}{k}(1-\theta_{j})^{k}\theta_{j}^{n\text{-}k}\big{)}-\prod\limits_{j=1}^{m}\big{(}1-\sum\limits_{k=0}^{z}\binom{n}{k}(1-\theta_{j})^{k}\theta_{j}^{n\text{-}k}\big{)}$

Analytically, it is more complex than the case of identical classifiers, but
it is easy to simulate.

We use a similar set-up as in Section 3.4 with $m=1,000$ classifiers, test set
size $n=3,000$, and probabilities of correct prediction
$\theta=\\{\theta_{1},\ldots,\theta_{m}\\}$, where
$\min\theta=0.875,\,\max\theta=0.90$, and they are equally spaced.

The simulated non-identical upper bound of the 95% confidence interval is
$0.9177$, with $100,000$ repetitions. The $E(\hat{\theta}_{SOTA})$ is
$0.9129$, with standard deviation of $0.0021$. The cdf and pmf are displayed
in Fig. 5 and Fig. 6.

Cumulative distribution function for non-identical classifiers.

Figure 5:

Probability mass function for non-identical classifiers.

Figure 6:

###  4.3 Non-identical, dependent classifiers

The most complex, but also the most realistic case, is where classifiers are
non-identical and dependent. An analytical solution is not known Hisakado et
al. (2006), but simulations are easy. The results of Boland et al. (1989) are
easily extended, and we have that

$\displaystyle P(Y_{j}|Y_{0}=1)$
$\displaystyle=\frac{\rho\sqrt{\theta_{j}(1-\theta_{j})\theta_{0}(1-\theta_{0})}+\theta_{j}\theta_{0}}{\theta_{0}}$
$\displaystyle P(Y_{j}|Y_{0}=0)$
$\displaystyle=\frac{-\rho\sqrt{\theta_{j}(1-\theta_{j})\theta_{0}(1-\theta_{0})}+\theta_{j}(1-\theta_{0})}{(1-\theta_{0})}$
(10)

We get the a the simulated dependent upper bound of the 0.95 confidence
interval of $0.91733$, with $100,000$ repetitions.

##  5 Discussion

We have addressed the mechanisms of multiplicity and how they affect the bias
of the SOTA estimator. When accuracy is used as a performance measure, the
exact probability distribution can be provided, and known analyses methods can
be engaged. The bias in SOTA is substantial in our example, and possibly
hinders new methods to attract interest from the scientific community. We have
also demonstrated that dependency might not have a big impact on the bias when
the probability of correct prediction is already high. The exception is when
the dependency is very strong.

The work presented here focuses on the simple case with an analytical
solution, and where the probabilities of correct prediction are known. In the
public competition or public data set context there are many factors that come
in to play, making the bias calculations everything but straightforward. We
discuss here an example, the SIIM-ISIC Melanoma Classification competition
from 2020 because it follows many of the common structures of public
challenges. The winning team chose many of the verified strategies from
Section 2.

The SIIM-ISIC Melanoma Classification competition consisted of more than
$40,000$ images of skin lesions, split into training ($33,126$), validation
($3,295$) and test set ($7,687$) from various sites. The data set is highly
imbalanced, with less than $2\%$ melanomas in the training set, which reflects
the reality in many cancer detection situations. A full description can be
found in Rotemberg et al. (2021). The size of the test set is smaller than the
recommendations of, e.g., Roelofs et al. (2019), but reflects the reality of
benchmark datasets and public challenges, as seen, e.g., in the overview of
Willemink et al. (2020). We want to emphasise that the shortcomings of this
competition are not unique, and that is precisely why it serves as a useful
example. Participants were allowed to use additional data for training if the
data was shared on the competition web site. The prize money was set to
$10,000$ USD for the best-performing method, measured as AUC on the hold-out
test set, published in the private leaderboard. Participants could get
intermediate results from the validation set, published on a public
leaderboard.

The team ‘All Data Are Ext’ consisting of three Kaggle Grandmasters won the
competition, achieving an AUC of $0.9490$ on the test set, and their method
and strategy is explained in Ha et al. (2020). Although we show that the AUC
is an optimistic estimate for the melanoma detection state-of-the-art, we do
not question their rank in the challenge. Quite the contrary; the research in
Section 2 supports their strategies. Their AUC decreased from the public to
the private leaderboard, but it was still within the confidence interval,
suggesting that they have avoided overfitting. They had a low rank on the
public leaderboard (880th out of 3300), and they have avoided the general
overfitting that can be seen in Fig. 7.

Figure 7: The line indicates public=private scores. There is a strong general
tendency of overfitting.

To avoid overfitting, they used additional data from previous years’
competitions, and they also did data augmentation. They chose ensemble
learning as a strategy to avoid overfitting, a strategy implicitly supported
by Talagrand (1996). They also trained on multiple classes, a strategy that
Feldman et al. (2019) showed will slow down the overfitting.

To estimate multiplicity effect on the AUC scores is not as straightforward as
with accuracy. The uncertainty of AUC can be calculated through bootstrapping,
but without access to the predictions for each data point in the test set, the
analyst must make some assumptions about the shape of the ROC curve. Fig. 8
shows different shapes of ROC curves and their confidence intervals. Fig. 9
shows corresponding plots for classifier predictions.

Figure 8: ROC curves with an AUC of .95 can look different. Above are three
ROC curves ordered from low variance AUC to high variance Figure 9: Scores
behind each of the above curves. The variability in the AUC curves has to do
with thhe way in which the distributions of scores for the two classes
overlap. A high-variance curve has mostly low overlap in scores but high
variance for the cases, meaning that a small handful of observations decide
the shape of the curve. (Recall that an ROC curve is built “backward:” the
threshold moves from left to right along the scores and the corresponding
operating characteristics are place in the plot from right to left)

However, with high AUC, the shape has less impact, and the analyst’s choices
are not suggestive of the overall conclusion. If we assume i.i.d. classifiers
and AUC = $0.9352$, the the upper bound of the $95\%$ confidence interval is
$0.9490$, same as the winning team. And since there are $427$ teams with AUC
$>0.9352$, we can expect observations outside the confidence interval. The
conclusion should be that although ‘All Data Are Ext’ did win the competition
according to the criteria, the AUC should not be used as the new SOTA. This
analysis show that with simple means a more sober SOTA estimate can be
communicated. A more elaborate analysis can be done through simulation set-
ups, and will decrease the SOTA even further. But the main conclusion remains
the same; highest-ranked performance is not the state-of-the-art, but rather
an optimistic estimate.

This work shows that a multiplicity-adjusted state-of-the-art estimate can be
calculated fairly easily through bootstrapping, and can serve as a better
reference point when introducing new methods to the scientific community.

## References

  * Tukey [1991] John W Tukey.  The philosophy of multiple comparisons.  _Statistical Science_ , 6:100–116, 1991.  URL https://www.jstor.org/stable/2245714#metadata_info_tab_contents.
  * Bennett et al. [2009] Craig M Bennett, Abigail A Baird, Michael B Miller, and George L Wolford.  Neural correlates of interspecies perspective taking in the post-mortem atlantic salmon: An argument for multiple comparisons correction.  _Neuroimage_ , 47:S125, 2009.
  * Demšar [2006] Janez Demšar.  Statistical comparisons of classifiers over multiple data sets.  _Journal of Machine Learning Research_ , 7:1–30, 2006.
  * Salzberg [1997] Steven L. Salzberg.  On comparing classifiers: Pitfalls to avoid and a recommended approach.  _Data Mining and Knowledge Discovery_ , 1:317–328, 1997\.  ISSN 13845810.  doi: 10.1023/A:1009752403260/METRICS.  URL https://link.springer.com/article/10.1023/A:1009752403260.
  * Thompson et al. [2020] William Hedley Thompson, Jessey Wright, Patrick G. Bissett, and Russell A. Poldrack.  Dataset decay and the problem of sequential analyses on open datasets.  _eLife_ , 9:1–17, 5 2020.  ISSN 2050084X.  doi: 10.7554/ELIFE.53498.
  * Teney et al. [2020] Damien Teney, Ehsan Abbasnejad, Kushal Kafle, Robik Shrestha, Christopher Kanan, and Anton van den Hengel.  On the value of out-of-distribution testing: An example of goodhart’s law.  _Advances in Neural Information Processing Systems_ , 33:407–417, 2020.
  * Torralba and Efros [2011] Antonio Torralba and Alexei A. Efros.  Unbiased look at dataset bias.  _Proceedings of the IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , pages 1521–1528, 2011.  ISSN 10636919.  doi: 10.1109/CVPR.2011.5995347.
  * Everingham et al. [2015] Mark Everingham, S. M.Ali Eslami, Luc Van Gool, Christopher K.I. Williams, John Winn, and Andrew Zisserman.  The pascal visual object classes challenge: A retrospective.  _International Journal of Computer Vision_ , 111:98–136, 1 2015.  ISSN 15731405.  doi: 10.1007/S11263-014-0733-5/FIGURES/27.  URL https://link.springer.com/article/10.1007/s11263-014-0733-5.
  * Fernández-Delgado et al. [2014] Manuel Fernández-Delgado, Eva Cernadas, Senén Barro, Dinani Amorim, and Amorim Fernández-Delgado.  Do we need hundreds of classifiers to solve real world classification problems?  _Journal of Machine Learning Research_ , 15:3133–3181, 2014\.  URL http://www.mathworks.es/products/neural-network.
  * Blum and Hardt [2015] Avrim Blum and Moritz Hardt.  The ladder: A reliable leaderboard for machine learning competitions.  pages 1006–1014. PMLR, 6 2015.  URL https://proceedings.mlr.press/v37/blum15.html.
  * Dwork et al. [2017] Cynthia Dwork, Vitaly Feldman, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Aaron Roth.  Guilt-free data reuse.  _Communications of the ACM_ , 60:86–93, 4 2017.  ISSN 15577317.  doi: 10.1145/3051088.
  * Recht et al. [2019] Benjamin Recht, Rebecca Roelofs, Ludwig Schmidt, and Vaishaal Shankar.  Do imagenet classifiers generalize to imagenet?  pages 5389–5400. PMLR, 5 2019.  URL https://proceedings.mlr.press/v97/recht19a.html.
  * Feldman et al. [2019] Vitaly Feldman, Roy Frostig, and Moritz Hardt.  The advantages of multiple classes for reducing overfitting from test set reuse.  pages 1892–1900. PMLR, 5 2019.  URL https://proceedings.mlr.press/v97/feldman19a.html.
  * Mania et al. [2019] Horia Mania, John Miller, Ludwig Schmidt, Moritz Hardt, and Benjamin Recht.  Model similarity mitigates test set overuse.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Roelofs et al. [2019] Rebecca Roelofs, Vaishaal Shankar, Benjamin Recht, Sara Fridovich-Keil, Moritz Hardt, John Miller, and Ludwig Schmidt.  A meta-analysis of overfitting in machine learning.  _Advances in Neural Information Processing Systems_ , 32, 2019.  URL https://www.kaggle.com/kaggle/meta-kaggle.
  * Dietterich [1998] Thomas G. Dietterich.  Approximate statistical tests for comparing supervised classification learning algorithms.  _Neural Computation_ , 10:1895–1923, 10 1998.  ISSN 0899-7667.  doi: 10.1162/089976698300017197.  URL https://direct.mit.edu/neco/article/10/7/1895/6224/Approximate-Statistical-Tests-for-Comparing.
  * Gundersen and Kjensmo [2018] Odd Erik Gundersen and Sigbjørn Kjensmo.  State of the art: Reproducibility in artificial intelligence.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 32:1644–1651, 4 2018.  ISSN 2374-3468.  doi: 10.1609/AAAI.V32I1.11503.  URL https://ojs.aaai.org/index.php/AAAI/article/view/11503.
  * Ma et al. [2021] Zhiyi Ma, Kawin, Ethayarajh, Tristan Thrush, Somya Jain, Ledell Wu, Robin Jia, Christopher Potts, Adina Williams, and Douwe Kiela.  Dynaboard: An evaluation-as-a-service platform for holistic next-generation benchmarking.  _Advances in Neural Information Processing Systems_ , 34:10351–10367, 12 2021.  URL https://eval.ai.
  * Bouthillier et al. [2021] Xavier Bouthillier, Pierre Delaunay, Mirko Bronzi, Assya Trofimov, Brennan Nichyporuk, Justin Szeto, Nazanin Mohammadi Sepahvand, Edward Raff, Kanika Madan, Vikram Voleti, Samira Ebrahimi Kahou, Vincent Michalski, Tal Arbel, Chris Pal, Gael Varoquaux, and Pascal Vincent.  Accounting for variance in machine learning benchmarks.  _Proceedings of Machine Learning and Systems_ , 3:747–769, 3 2021.
  * Bousquet and Elisseeff [2002] Olivier Bousquet and André Elisseeff.  Stability and generalization.  _Journal of Machine Learning Research_ , 2:499–526, 2002\.  URL http://sensitivity-analysis.jrc.cec.eu.int/.
  * Talagrand [1996] Michel Talagrand.  A new look at independence.  _The Annals of Probability_ , 24:1–34, 1996.  URL https://www.jstor.org/stable/2244830#metadata_info_tab_contents.
  * Varoquaux and Cheplygina [2022] Gaël Varoquaux and Veronika Cheplygina.  Machine learning for medical imaging: methodological failures and recommendations for the future.  _npj Digital Medicine 2022 5:1_ , 5:1–8, 4 2022.  ISSN 2398-6352.  doi: 10.1038/s41746-022-00592-y.  URL https://www.nature.com/articles/s41746-022-00592-y.
  * Boland et al. [1989] Philip J. Boland, Frank Proschan, and Y. L. Tong.  Modelling dependence in simple and indirect majority systems.  _Journal of Applied Probability_ , 26:81–88, 3 1989.  ISSN 0021-9002.  doi: 10.2307/3214318.  URL https://www.cambridge.org/core/journals/journal-of-applied-probability/article/abs/modelling-dependence-in-simple-and-indirect-majority-systems/070D6335BDDDDC7AF4D70BC9B21B0B7B.
  * Ladd [1975] Daniel W. Ladd.  An algorithm for the binomial distribution with dependent trials.  _Journal of the American Statistical Association_ , 70:333–340, 1975.  ISSN 1537274X.  doi: 10.1080/01621459.1975.10479867.
  * Hisakado et al. [2006] Masato Hisakado, Kenji Kitsukawa, and Shintaro Mori.  Correlated binomial models and correlation structures.  2006\.
  * Rotemberg et al. [2021] Veronica Rotemberg, Nicholas Kurtansky, Brigid Betz-Stablein, Liam Caffery, Emmanouil Chousakos, Noel Codella, Marc Combalia, Stephen Dusza, Pascale Guitera, David Gutman, Allan Halpern, Brian Helba, Harald Kittler, Kivanc Kose, Steve Langer, Konstantinos Lioprys, Josep Malvehy, Shenara Musthaq, Jabpani Nanda, Ofer Reiter, George Shih, Alexander Stratigos, Philipp Tschandl, Jochen Weber, and H. Peter Soyer.  A patient-centric dataset of images and metadata for identifying melanomas using clinical context.  _Scientific Data 2021 8:1_ , 8:1–8, 1 2021.  ISSN 2052-4463.  doi: 10.1038/s41597-021-00815-z.  URL https://www.nature.com/articles/s41597-021-00815-z.
  * Willemink et al. [2020] Martin J. Willemink, Wojciech A. Koszek, Cailin Hardell, Jie Wu, Dominik Fleischmann, Hugh Harvey, Les R. Folio, Ronald M. Summers, Daniel L. Rubin, and Matthew P. Lungren.  Preparing medical imaging data for machine learning.  _Radiology_ , 295:4–15, 2 2020.  ISSN 15271315.  doi: 10.1148/RADIOL.2020192224/ASSET/IMAGES/LARGE/RADIOL.2020192224.FIG5B.JPEG.  URL https://pubs.rsna.org/doi/10.1148/radiol.2020192224.
  * Ha et al. [2020] Qishen Ha, Bo Liu, and Fuxu Liu.  10 2020.  doi: 10.48550/arxiv.2010.05351.  URL https://arxiv.org/abs/2010.05351v1.