# Early kinetic decoupling effect on the forbidden dark matter annihilations
into standard model particles

Yu Liu<EMAIL_ADDRESS>Xuewen Liu<EMAIL_ADDRESS>Bin Zhu
<EMAIL_ADDRESS>Department of Physics, Yantai University, Yantai
264005, China

###### Abstract

The early kinetic decoupling (eKD) effect is an inevitable ingredient in
calculating the relic density of dark matter (DM) for various well-motivated
scenarios. It appears naturally in forbidden dark matter annihilation, the
main focus of this work, which contains fermionic DM and a light singlet
scalar that connects the DM and standard model (SM) leptons. The strong
suppression of the scattering between DM and SM particles happens quite early
in the DM depletion history, where the DM temperature drops away from the
thermal equilibrium, $T_{\chi}<T_{\rm SM}$, leading to the decreased kinetic
energy of DM. The forbidden annihilation thus becomes inefficient since small
kinetic energy cannot help exceed the annihilation threshold, naturally
leading to a larger abundance. To show the eKD discrepancy, we numerically
solve the coupled Boltzmann equations that govern the evolution of DM number
density and temperature. It is found that eKD significantly affects the DM
abundance, resulting in almost an order of magnitude higher than that by the
traditional calculation. We also discuss the constraints from experimental
searches on the model parameters, where the viable parameter space shrinks
when considering the eKD effect.

##  I Introduction

Relic density is an essential topic for dark matter (DM) physics. The
classical scenario that explains the observed abundance in the present
Universe is thermal particle production in the early Universe, which is the
so-called freeze-out mechanism for the weakly interacting massive particle
(WIMP) Lee and Weinberg (1977); Hut (1977); Sato and Kobayashi (1977); Dicus
_et al._ (1977); Wolfram (1979). The DM particles were initially in thermal
equilibrium with the heat bath via the intense interactions among them. The
number density dilutes along with the expansion of the Universe and finally
freezes out of the heat bath once the annihilation rate falls behind the
cosmic expansion rate, resulting in a comoving constant.

Gondolo et al. Gondolo and Gelmini (1991); Edsjo and Gondolo (1997) have
developed the renowned treatment of calculating the DM relic density by
solving the Boltzmann equation of the number density with high accuracy,
called the “standard” method. One hypothesis entering this treatment is that
DM keeps in local kinetic equilibrium with the thermal plasma during or even
after the freeze-out process. The scatterings with the standard model (SM)
particles have been at a much more intense level Bringmann (2009).

However, this is not always true for many well-motivated mechanisms, where
kinetic equilibrium might decouple earlier than assumed, leading to the early
kinetic decoupling (eKD) around the freeze-out period. The eKD effect has been
extensively studied in the literature van den Aarssen  _et al._ (2012); Binder
_et al._ (2017); Brümmer (2020); Ala-Mattinen and Kainulainen (2020); Abe
(2020); Binder  _et al._ (2021); Zhu and Liu (2022); Hryczuk and Laletin
(2021); Abe (2021); Du  _et al._ (2022); Ala-Mattinen  _et al._ (2022);
Hryczuk and Laletin (2022). The most influenced scenarios include 1) resonant
annihilation of dark matter Binder  _et al._ (2017); Abe (2020); Binder  _et
al._ (2021), 2) Sommerfeld-enhanced annihilation Binder  _et al._ (2021), and
3) subthreshold annihilation (also known as the forbidden annihilation) Binder
_et al._ (2021). In these regimes, the eKD occurs because the elastic
scattering processes are suppressed; DM and the SM particles experience
different temperatures. All these cases demonstrated the actual DM abundance
can be affected by up to 1 order of magnitude compared to the traditional
method, at least in some parts of the parameter space.

Here, we study the eKD in the forbidden DM scenario D’Agnolo and Ruderman
(2015). In such a scenario, DM dominantly annihilates into heavier final
states, which can proceed at finite temperatures in the early Universe,
relying on the thermal tail with the high velocity of DM. Many works were
devoted to studying the forbidden annihilations with different theoretical
models and a variety of phenomenological topics Delgado  _et al._ (2017);
D’Agnolo  _et al._ (2021); Wojcik and Rizzo (2022); Yang (2022); Cheng  _et
al._ (2022); Herms  _et al._ (2022a, b). In this work, we employ the model
where forbidden DM annihilations into the SM leptons are mediated by a singlet
scalar. Such channels were studied in Ref.D’Agnolo  _et al._ (2021), which is
experimentally viable and predicts a very narrow mass range for DM, that can
be tested at future beam dump experiments. However, the eKD effect has not
been studied in the context of this model.

Actually, the scatterings of DM against SM particles are strongly suppressed
in the early times because of the mass splitting between the DM and SM
leptons. The DM temperature drops away from the thermal bath temperature,
$T_{\chi}<T_{\rm SM}$, which leads to the decreased kinetic energy of DM. The
forbidden annihilations thus become inefficient since small kinetic energy
cannot help exceed the annihilation threshold. So this naturally leads to a
larger abundance and the eKD effects in such a model should not be neglected.
On the contrary, eKD will cause significant impacts on DM abundance. We
investigate the relic density beyond the standard treatment used in Ref.
D’Agnolo  _et al._ (2021) by considering the coupled Boltzmann differential
equations, where the temperature evolution of the dark sector could be taken
into account. We use the public code DRAKE Binder  _et al._ (2021) to perform
the numerical calculations. We find a DM relic density that differs by up to
an order of magnitude from the standard treatment and leaves a reduced
feasible parameter space under the various experimental constraints.

The rest of this paper is structured as follows. In Sec. II, we start with a
general description of the coupled Boltzmann equations that govern the
evolution of DM number density and temperature. In addition, we discuss the DM
model for forbidden annihilations and analyze the occurrence of early kinetic
decoupling. Section III is devoted to a thorough study of DM relic density for
the forbidden channels, and makes a detailed comparison between the new
treatment and the traditional one. We further discuss various constraints from
collider searches and astrophysical observations on the parameter space in
Sec. IV. Finally, we conclude in Sec. V.

##  II Early kinetic decoupling effects on forbidden annihilations

###  II.1 Basic formulas

Keeping kinetic equilibrium during and even after the freeze-out epoch is one
underlying assumption for traditional relic density calculations. However,
this is not always the case for various scenarios, where kinetic decoupling
happens earlier than the chemical decoupling process. To study the DM relic
density by taking into account the early kinetic decoupling effect, we should
consider the following Boltzmann equation for DM phase-space distribution
Binder  _et al._ (2017, 2021)

$\displaystyle E\left(\frac{\partial}{\partial
t}-H\vec{p}\cdot\frac{\partial}{\partial\vec{p}}\right)f_{\chi}(t,\vec{p})=C_{\rm
ann.}[f_{\chi}]+C_{\rm el.}[f_{\chi}],$ (1)

where $E$ is the energy of the DM, $H$ is the Hubble constant, $\vec{p}$ is
the momentum of DM, and $f_{\chi}$ is the DM phase-space density. The
collision term $C_{\rm ann.}$ represents the annihilation of DM particles into
thermal bath particles, and $C_{\rm el.}$ is for elastic scattering processes
between DM and SM scattering partners. For two-body processes,

$\displaystyle C_{\rm ann.}$ $\displaystyle=$
$\displaystyle\frac{1}{2g_{\chi}}\sum\int\frac{d^{3}p^{\prime}}{(2\pi)^{3}2E_{p^{\prime}}}\int\frac{d^{3}k}{(2\pi)^{3}2E_{k}}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}2E_{k^{\prime}}}$
$\displaystyle\times$
$\displaystyle(2\pi)^{4}\delta^{4}(p+p^{\prime}-k-k^{\prime})$
$\displaystyle\times$ $\displaystyle\Bigl{(}-|{\cal M}_{\chi\chi\to{\cal
B}{\cal B}^{\prime}}|^{2}f_{\chi}(\vec{p})f_{\chi}(\vec{p^{\prime}})(1\pm
f^{eq}_{\cal B}(\vec{k}))(1\pm f^{eq}_{\cal B^{\prime}}(\vec{k^{\prime}}))$
$\displaystyle+$ $\displaystyle|{\cal M}_{{\cal B}{\cal
B^{\prime}}\to\chi\chi}|^{2}f^{eq}_{\cal B}(\vec{k})f^{eq}_{{\cal
B}^{\prime}}(\vec{k^{\prime}})(1\pm f_{\chi}(\vec{p}))(1\pm
f_{\chi}(\vec{p^{\prime}}))\Bigr{)},$ $\displaystyle C_{\rm el.}$
$\displaystyle=$
$\displaystyle\frac{1}{2g_{\chi}}\sum\int\frac{d^{3}p^{\prime}}{(2\pi)^{3}2E_{p^{\prime}}}\int\frac{d^{3}k}{(2\pi)^{3}2E_{k}}\int\frac{d^{3}k^{\prime}}{(2\pi)^{3}2E_{k^{\prime}}}$
$\displaystyle\times$
$\displaystyle(2\pi)^{4}\delta^{4}(p+p^{\prime}-k-k^{\prime})$
$\displaystyle\times$ $\displaystyle\Bigl{(}-|{\cal M}_{\chi{\cal
B}\to\chi{\cal B}}|^{2}f_{\chi}(\vec{p})f^{eq}_{\cal B}(\vec{k})(1\pm
f_{\chi}(\vec{p^{\prime}}))(1\pm f^{eq.}_{\cal B}(\vec{k^{\prime}}))$
$\displaystyle+$ $\displaystyle|{\cal M}_{\chi{\cal B}\to\chi{\cal
B}}|^{2}f_{\chi}(\vec{p^{\prime}})f^{eq.}_{\cal B}(\vec{k^{\prime}})(1\pm
f_{\chi}(\vec{p}))(1\pm f^{eq.}_{\cal B}(\vec{k}))\Bigr{)},$

where ${\cal B}$ and ${\cal B}^{\prime}$ stand for particles in the thermal
bath such as SM leptons, $g_{\chi}$ is the number of internal degrees of
freedom of DM, and $f_{\cal B}^{eq}$ is given by the Fermi-Dirac or Bose-
Einstein distribution depending on the spin of ${\cal B}$. The summation
should be taken for all the internal degrees of freedom for all the particles.
For the nonrelativistic DM, $C_{\rm el.}$ can be simplified as the Fokker-
Planck operator Binder  _et al._ (2016); Bertschinger (2006); Bringmann and
Hofmann (2007); Bringmann (2009) 111As pointed out by Ref.Binder  _et al._
(2021), for forbidden DM, the Fokker-Planck approximation is not that
accurate, but the dominant eKD effect on the relic density can in many cases
still be fairly well captured by the Fokker-Planck approximation. For more
precise treatment, we leave it to future work.:

$\displaystyle C_{\mathrm{el}}$
$\displaystyle\simeq\frac{E}{2}\gamma(T){\Bigg{[}}TE\partial_{p}^{2}\\!+\left(2T\frac{E}{p}\\!+\\!p\\!+\\!T\frac{p}{E}\right)\partial_{p}+3{\Bigg{]}}f_{\chi}\,.$
(4)

In the above, the momentum transfer rate $\gamma(T)$ is given by (see also
Ref. Gondolo  _et al._ (2012))

$\displaystyle\gamma=\frac{1}{3g_{\chi}m_{\chi}T}\\!\int\\!\frac{\text{d}^{3}k}{(2\pi)^{3}}f_{\cal
B}^{\pm}(E_{k})\left[1\\!\mp\\!f_{\cal
B}^{\pm}(E_{k})\right]\\!\\!\\!\int\limits^{0}_{-4k_{\mathrm{cm}}^{2}}\\!\\!\\!\text{d}t(-t)\frac{\text{d}\sigma}{\text{d}t}v\,,$
(5)

where the differential cross section can be expressed as
$({\text{d}\sigma}/{\text{d}t})v$ $\equiv$ $|\mathcal{M}|^{2}_{\chi
f\leftrightarrow\chi f}$ /$(64\pi{k}E_{k}m_{\chi}^{2})$, and $k_{\rm cm}^{2}$
is given by

$\displaystyle k_{\rm cm}^{2}=\frac{m_{\chi}^{2}(E_{k}^{2}-m_{\cal
B}^{2})}{m_{\chi}^{2}+m_{\cal B}^{2}+2m_{\chi}E_{k}}.$ (6)

Here $E_{k}$ is the energy of heat bath particle ${\cal B}$. Note that $k_{\rm
cm}^{2}\neq E_{k}^{2}-m_{\cal B}^{2}=|\vec{k}|^{2}$.

During the chemical decoupling, the scattering processes may not be frequent
enough to maintain the kinetic equilibrium, which means that DM particles own
a different temperature $T_{\chi}$ from the thermal plasma in their following
evolution. A common definition of the DM temperature is

$\displaystyle T_{\chi}$ $\displaystyle=$
$\displaystyle\frac{g_{\chi}}{3n_{\chi}}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{\vec{p}^{2}}{E}f_{\chi}(\vec{p})\equiv\frac{s^{2/3}}{m_{\chi}}y,$
(7)

which is also a function of the thermal bath temperature $T$. In this
definition $n_{\chi}$ is the number density of the DM, and $s$ is the entropy
density. Here $y$ is a dimensionless version in analogy to the DM yield
$Y(=n_{\chi}/s)$.

To reach a suitable description of the DM temperature evolution and then
explore the eKD effect on the chemical decoupling process, we should consider
the second moment of $f_{\chi}$ as a dynamical degree of freedom. By
integrating Eq. (1) with $g_{\chi}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{E}$
and
$g_{\chi}\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{E}\frac{\vec{p}^{2}}{E^{2}}$,
one obtains the zeroth and second moments of the Boltzmann equation,
respectively. This leads to a relatively simple coupled system of Boltzmann
differential equations ( denoted as the cBE method hereafter),

$\displaystyle\frac{Y^{\prime}}{Y}$ $\displaystyle=$
$\displaystyle\frac{sY}{x\tilde{H}}\left[\frac{Y_{\rm
eq}^{2}}{Y^{2}}\left\langle\sigma v\right\rangle_{T}-\left\langle\sigma
v\right\rangle_{T_{\chi}}\right]\,,$ (8) $\displaystyle\frac{y^{\prime}}{y}$
$\displaystyle=$ $\displaystyle\frac{1}{x\tilde{H}}\langle
C_{\text{el}}\rangle_{2}+\frac{sY}{x\tilde{H}}\left[\left\langle\sigma
v\right\rangle_{T_{\chi}}-\left\langle\sigma
v\right\rangle_{2,T_{\chi}}\right]$
$\displaystyle+\frac{sY}{x\tilde{H}}\frac{Y_{\rm
eq}^{2}}{Y^{2}}\left[\frac{y_{{\rm eq}}}{y}\left\langle\sigma
v\right\rangle_{2,T}\\!-\\!\left\langle\sigma
v\right\rangle_{T}\right]+2(1-w)\frac{H}{x\tilde{H}}\,,$

where $x$ is defined as usual $x=m_{\chi}/T$ and $Y_{\rm eq}(x)\equiv n_{\rm
eq}(T)/s$. $\tilde{H}\equiv H/\left[1+(1/3)d(\log g^{s}_{\rm eff})/d(\log
T)\right]$, with $g^{s}_{\rm eff}$ being the entropy degrees of freedom of the
background plasma. $w(T_{\chi})\equiv 1-{\langle
p^{4}/E^{3}\rangle_{T_{\chi}}}/({6T_{\chi}})$, with $\langle
p^{4}/E^{3}\rangle=\frac{g_{\chi}}{n_{\chi}^{eq}(T_{\chi})}\int\frac{d^{3}p}{(2\pi)^{3}}{\left(\vec{p}\cdot\vec{p}\right)^{2}}/{E^{3}}e^{-\frac{E}{T_{\chi}}}$.
Note that the elastic scattering term given in Eq. (4) does not contribute to
the zeroth moment term. This is a natural consequence because the elastic
scattering processes do not change the number density of DM.

The above compact form of the differential equations contains the following
thermally averaged cross sections,

$\displaystyle\langle
C_{\text{el}}\rangle_{2}\equiv\frac{g_{\chi}}{3nT_{\chi}}\int\frac{\text{d}^{3}p}{(2\pi)^{3}}\frac{p^{2}}{E^{2}}C_{\text{el}}\;,$
(10) $\displaystyle\left\langle\sigma v\right\rangle_{T_{\chi}}$
$\displaystyle\equiv$
$\displaystyle\frac{g_{\chi}^{2}}{(n_{\chi}^{eq})^{2}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}q}{(2\pi)^{3}}\left(\sigma
v\right)_{\chi\chi\to{\cal B}{\cal B}^{\prime}}$ (11) $\displaystyle\times$
$\displaystyle
f_{\chi}^{eq}(\vec{p},T_{\chi})f_{\chi}^{eq}(\vec{q},T_{\chi}),$

The thermal average $\langle\sigma v\rangle_{2,T}$ is a variant of the
commonly used thermal average $\langle\sigma v\rangle_{T}$, and is explicitly
stated in Ref. Binder  _et al._ (2017) and introduced as

$\displaystyle\left\langle\sigma v\right\rangle_{2,T_{\chi}}$ $\displaystyle=$
$\displaystyle\frac{g_{\chi}^{2}}{(n_{\chi}^{eq})^{2}T_{\chi}}\int\frac{d^{3}p}{(2\pi)^{3}}\int\frac{d^{3}q}{(2\pi)^{3}}\frac{\vec{p}\cdot\vec{p}}{3E}\left(\sigma
v\right)_{\chi\chi\to{\cal B}{\cal B}^{\prime}}$ (12) $\displaystyle\times$
$\displaystyle
f_{\chi}^{eq}(\vec{p},T_{\chi})f_{\chi}^{eq}(\vec{q},T_{\chi}).$

For $\left\langle\sigma v\right\rangle_{T}$ and $\left\langle\sigma
v\right\rangle_{2,T}$, replace $T_{\chi}$ by $T$ in $\left\langle\sigma
v\right\rangle_{T_{\chi}}$ and $\left\langle\sigma
v\right\rangle_{2,T_{\chi}}$, respectively. And $n_{\chi}^{eq}(T_{\chi})$ is
given by

$\displaystyle
n_{\chi}^{eq}(T_{\chi})=g_{\chi}\int\frac{d^{3}p}{(2\pi)^{3}}f_{\chi}^{eq}(\vec{p},T_{\chi})=g_{\chi}\int\frac{d^{3}p}{(2\pi)^{3}}e^{-\frac{E_{p}}{T_{\chi}}}.$
(13)

In this work, we use the numerical routine DRAKE to solve the coupled
Boltzmann equations. The measured value of $\Omega h^{2}$ by the Planck
Collaboration is $\Omega h^{2}=0.120\pm 0.001$ Aghanim  _et al._ (2020). The
viable parameter space is determined by matching this value.

###  II.2 Model and discussion on forbidden channels

We have adopted a simple model that only takes into account DM annihilations
into SM leptons. The DM is a Dirac fermion coupled to the SM sector via the
scalar portal $\phi$. After the electroweak symmetry breaking, the effective
Lagrangian can be written as

$-\mathcal{L}\supset
g_{ij}\phi\bar{l}_{i}l_{j}+g^{A}_{ij}\phi\bar{l}_{i}\gamma_{5}l_{j}+g_{\chi}\phi\bar{\chi}\chi+g^{A}_{\chi}\phi\bar{\chi}\gamma_{5}\chi,$
(14)

where the indices on the couplings $i,j=e,\mu,\tau$. This model has been
studied thoroughly for the forbidden mechanism (for detail, refer to Ref.
D’Agnolo  _et al._ (2021)). The merits include the following: 1) DM mass is
limited in quite a small window close to the masses of the SM leptons, which
is a strong prediction that can be tested soon by colliders or beam-dump
experiments. 2) Kinematically forbidden DM naturally evade the stringent
constraints from the energy injections into the cosmic microwave background
(CMB) D’Agnolo and Ruderman (2015). In the forbidden DM scenario, the energy
injection processes suffer the Boltzmann suppression at $T\lesssim$ eV, so
that sub-GeV thermal relics are consistent with the experiment, making
annihilations to SM leptons with DM masses $m_{\chi}\ll 10$ GeV still viable.

In this scenario, the DM relic density should be carefully scrutinized, as the
eKD effect appears generic. The actual relic density receives a significant
correction compared with the conventional method, as shown in the following
sections.

Following Refs. D’Agnolo and Ruderman (2015); D’Agnolo  _et al._ (2021), we
consider a pair of DM particles dominantly annihilating into two SM particles
$2\ell$ with mass $m_{\chi}<m_{\ell}$. Cosmological constraints make forbidden
annihilations into electrons unfeasible, including big bang nucleosynthesis
(BBN) and CMB Sabti  _et al._ (2020). So, for simplicity, we only consider the
$\mu^{+}\mu^{-}$ and $\tau^{+}\tau^{-}$ channels, with abbreviated couplings
$g^{(A)}_{\mu},~{}g^{(A)}_{\tau},$ and $g_{\chi}^{(A)}$, which allows us to
explore all the relevant DM phenomenology systematically.

Applying the detailed balance condition for the DM number-changing process,
the cross section of the forbidden channels is exponentially suppressed,

$\langle\sigma_{\chi}v\rangle=\langle\sigma_{\ell}v\rangle\frac{(n_{\ell}^{\rm
eq})^{2}}{(n_{\chi}^{\rm eq})^{2}}\simeq\langle\sigma_{\ell}v\rangle
e^{-2\Delta x}\,,$ (15)

where $\Delta\equiv(m_{\ell}-m_{\chi})/m_{\chi}$.
$\sigma_{\chi}\equiv\sigma(\chi\chi\to\ell\ell_{)}$, while $\sigma_{\ell}$ is
the cross section for the inverse process. When the annihilation rate becomes
slower than the Hubble expansion, DM is no longer in equilibrium with the SM
thermal bath, resulting in chemical decoupling.

What about the scattering between DM and SM particles during this period? From
Eq. (5), the momentum transfer rate $\gamma$ is proportional to an exponential
factor:

$\gamma(x)\propto e^{-(\Delta+1)x}.$ (16)

The full expression of the rates is listed in Appendix A. In the forbidden
scenario, $\Delta+1>1$ implies the scattering frequency experiences a strong
suppression at a much earlier period.

It is known that DM kinetically decouples out of the SM thermal bath as long
as the momentum transfer rate $\gamma$ is smaller than the Hubble expansion,
$\gamma<H$. For illustration, we show the momentum transfer rate $\gamma$ for
the $\chi\mu^{\pm}\to\chi\mu^{\pm}$ scattering process, in Fig. 1. In
comparison, we plot the Hubble parameter $H(x)$ as a function of $x$. The red
lines stand for the evolution of $\gamma(x)$ and $H(x)$ in the forbidden DM
case where we take $m_{\chi}=0.1{\rm GeV}$. To demonstrate the distinctiveness
of the forbidden DM, we also provide the results of a nonforbidden case where
DM mass is larger than that of the annihilation products ($m_{\chi}$=1 GeV).
The most remarkable finding is that $\gamma$ of the forbidden case becoming
comparable with $H(x)$ happens much earlier than that of the nonforbidden
case. The kinetic decoupling starts at around $x=20$, which is usually the
same time as DM chemically decoupled from the thermal bath. The reason for
this very early kinetic decoupling is straightforward to understand as the
result of an exponential suppressed momentum transfer rate, as derived in Eq.
(16). One can obtain similar results for the
$\chi\bar{\chi}\to\tau^{+}\tau^{-}$ forbidden channel. We conclude that eKD
exists in the forbidden DM scenario.

Figure 1: Evolution of momentum transfer rate $\gamma(x)$ for
$\chi\mu^{\pm}\to\chi\mu^{\pm}$ scattering, and comparing with Hubble constant
H(x). The red lines stand for the forbidden scenario, where $m_{\chi}=0.1{\rm
GeV}$ and $m_{\phi}=0.26{\rm GeV}$. The gray lines correspond to the
nonforbidden case, that we take as $m_{\chi}=1{\rm GeV}$.

The next step is to find the effect of the eKD more concretely. We
systematically study DM relic density in both cBE and traditional methods
(denoted as nBE as in Ref. Binder  _et al._ (2021)) and then discuss the
phenomenological possibilities.

##  III Relic density: comparison between cBE and nBE approaches

In this section, we compute the relic density in both the standard method
(nBE) and the cBE approach. We restrict our study to the same range of
$m_{\chi}$ that $0.9m_{\mu}\lesssim m_{\chi}\lesssim m_{\mu}$ as derived in
Ref.D’Agnolo  _et al._ (2021), in which forbidden annihilation into
$\mu^{+}\mu^{-}$ is experimentally viable. For the
$\chi\bar{\chi}\to\tau^{+}\tau^{-}$ case, the corresponding DM mass is
$0.8m_{\tau}\lesssim m_{\chi}\lesssim m_{\tau}$. The detailed annihilation
cross sections and the momentum transfer rates for scatterings are presented
in Appendix A.

Table 1: Benchmark points, selected for $\chi\bar{\chi}\to\mu^{+}\mu^{-}$ and $\chi\bar{\chi}\to\tau^{+}\tau^{-}$ channels. We set $g_{\chi}=0$ and $g^{A}_{\chi}=1.121$ for all cases as in Ref.D’Agnolo  _et al._ (2021) which can make a direct comparison. Benchmark ($\mu^{+}\mu^{-}$) |  | $m_{\chi}$ |  | $m_{\phi}$ |  | $g_{\mu}$ |  | $g_{\mu}^{A}$ |  | $\Omega_{\rm cBE}h^{2}/\Omega_{\rm nBE}h^{2}$
---|---|---|---|---|---|---|---|---|---|---
BP$\mu$1 |  | 0.1 GeV |  | 0.26 GeV |  | 0.00343326 |  | 0 |  | 0.12/0.018
BP$\mu$2 |  | 0.1 GeV |  | 0.26 GeV |  | 0.00091757 |  | $g_{\mu}$ |  | 0.12/0.0186
BP$\mu$3 |  | 0.1055 GeV |  | 0.3 GeV |  | 0.00145082 |  | 0 |  | 0.12/0.036
BP$\mu$4 |  | 0.1055 GeV |  | 0.3 GeV |  | 0.00026751 |  | $g_{\mu}$ |  | 0.12/0.064
Benchmark ($\tau^{+}\tau^{-}$) |  | $m_{\chi}$ |  | $m_{\phi}$ |  | $g_{\tau}$ |  | $g_{\tau}^{A}$ |  | $\Omega_{\rm cBE}h^{2}/\Omega_{\rm nBE}h^{2}$
BP$\tau$1 |  | 1.6 GeV |  | 5 GeV |  | 1.09261 |  | 0 |  | 0.12/0.0088
BP$\tau$2 |  | 1.6 GeV |  | 5 GeV |  | 0.268806 |  | $g_{\tau}$ |  | 0.12/0.0078
BP$\tau$3 |  | 1.77 GeV |  | 6 GeV |  | 0.0565632 |  | 0 |  | 0.12/0.032
BP$\tau$4 |  | 1.77 GeV |  | 6 GeV |  | 0.0096236 |  | $g_{\tau}$ |  | 0.12/0.058

To reveal the effects of the eKD and the differences between the cBE and nBE
approaches, we first find several benchmark points for the two forbidden
channels, shown in Table 1. We also set
$g_{\chi}=0,~{}(g^{A}_{\chi})^{2}/4\pi=0.1$ as in Ref.D’Agnolo  _et al._
(2021), for making a rough but straight comparison. The values of the rest
model parameters are fixed to obtain the observed DM relic density for
$\Omega_{\rm cBE}h^{2}$. With these inputs, we can compute the relic density
in the nBE method. It can be seen that the ratio of $\Omega_{\rm
cBE}/\Omega_{\rm nBE}$ is sizable, even reaching an order of magnitude.

The significant difference between the cBE and nBE results exactly comes from
the eKD effects. In Fig. 2, we show the temperature and abundance evolution
for selected benchmarks in Tab. 1. From the left panel, the green lines are
the evolution curves of $y_{\chi}$, namely the temperature of DM, which depart
from the thermal bath temperature (the gray line) at around $x=20$.
Qualitatively, DM needs higher momenta to overcome the annihilation threshold,
leading to a self-cooling phase as soon as it is no longer kinetically coupled
to the muons. This is why there is a drop and the temperatures evolve
separately for the dark sector and the SM sector since then. DM annihilation
becomes less efficient much earlier just because of this cooling, which
results in a higher DM abundance than in the nBE approach, as shown in the
right panel of Fig. 2. Note that the same cooling phenomena also have been
found in Ref. Binder  _et al._ (2021).

Figure 2: Evolution of DM temperature $y$ and DM abundance $Y$. Left: the
green lines depict the DM temperature evolution in the cBE approach, with
benchmarks BP$\mu$1 (solid), BP$\mu$2 (dashed), BP$\mu$3 (dot-dashed), and
BP$\mu$4 (dotted). The gray line stands for the temperature of the thermal
bath. Right: the red lines are the DM yields in the cBE approach, for
different parameter settings. The gray lines represent the yields in the
traditional nBE approach. The convention of the line styles is the same as in
the left diagram.

In Fig. 3, we show a global picture of the eKD effect for the forbidden cases
of $\chi\bar{\chi}\to\mu^{+}\mu^{-}$ and $\chi\bar{\chi}\to\tau^{+}\tau^{-}$,
where we define the deviations of the relic density in cBE and nBE approaches
by ${\rm Deviation}\equiv(\Omega_{\rm cBE}-\Omega_{\rm nBE})/\Omega_{\rm
cBE}$. We display the results in the $(m_{\phi},~{}g_{\mu/\tau})$ plane by the
density plotting method, with setting $m_{\chi}=0.1{\rm GeV}$ and
$m_{\chi}=1.48{\rm GeV}$ for illustration. The sizeable deviations appear in
most parameter spaces from 20% up to almost 100%. The maximum deviations seem
to emerge in the resonance region. However, it should be noted that eKD
effects already exist in the resonant annihilations of DM Binder  _et al._
(2017); Abe (2020); Binder  _et al._ (2021). So in this region, one should
study the eKD for resonance and forbidden channels together. It is beyond the
scope of our work, as we mainly focus on the forbidden annihilations.

Figure 3: Density distribution of the deviation between cBE and nBE
approaches, which is defined as $(\Omega_{\rm cBE}-\Omega_{\rm
nBE})/\Omega_{\rm cBE}$. The left panel is for the
$\chi\bar{\chi}\to\mu^{+}\mu^{-}$ channel and the right is for the
$\chi\bar{\chi}\to\tau^{+}\tau^{-}$ channel; we take $m_{\chi}=0.1$ GeV and
$m_{\chi}=1.48$ GeV for the left and right panels respectively.

Figure 4: Relic density comparison between $\Omega_{\rm cBE}h^{2}$ and
$\Omega_{\rm nBE}h^{2}$ for the $\chi\bar{\chi}\to\mu^{+}\mu^{-}$ and
$\chi\bar{\chi}\to\tau^{+}\tau^{-}$ channels. Here, we choose the parameters
to match the observed abundance in the nBE approach. The left diagram
corresponds to the $\mu^{+}\mu^{-}$ mode, and we take $m_{\phi}=0.26$ GeV, and
the right diagram presents the $\tau^{+}\tau^{-}$ mode, for which we take
$m_{\phi}=5$ GeV.

To emphasize the importance of the improved treatment of the decoupling
history near the mass threshold, we plot in Fig. 4 the ratio of the resulting
relic density to that of the standard nBE approach. Here the parameters
satisfy the requirement for $\Omega_{\rm nBE}h^{2}=0.12$. The different
choices of coupling $g^{A}_{l}$ correspond to different curves, as labeled in
the plots. The mass ratio $m_{\ell}/m_{\chi}$ can be divided into three
regions. At lower mass ratios, dark matter evolves as ordinary WIMPs where
kinetic equilibrium is maintained. The relic density derived by the cBE and
nBE methods agrees with each other. The gray-shaded region is known as the
resonance region, where $2m_{\chi}\simeq m_{\phi}$. Here, the eKD effect
arises due to the distinct cooling and heating effects of dark matter.For a
more detailed discussion of the origin of these features, please refer to
Ref.Binder  _et al._ (2017). In the forbidden DM mass region (brown-shaded),
we can see that the cBE results are larger than nBE several times with the
same parameters. For $\mu^{+}\mu^{-}$ case, $3\lesssim\frac{\Omega_{\rm
cBE}}{\Omega_{\rm nBE}}\lesssim 10$, and $2\lesssim\frac{\Omega_{\rm
cBE}}{\Omega_{\rm nBE}}\lesssim 15$ for the $\tau^{+}\tau^{-}$ forbidden case.

As the mass ratio $m_{\mu/\tau}/m_{\chi}$ increases, the coupling
($g_{\mu/\tau}$), which is needed to obtain the correct relic density, rises
rapidly (see Fig. 5). We find the upper bounds of the mass ratio, beyond which
the couplings become non-perturbative (i.e., $g_{\mu/\tau}<\sqrt{4\pi}$).
Since a larger coupling is required in the cBE approach, a smaller mass ratios
are allowed compared to the standard nBE treatment.

Figure 5: Correct relic density as function of mass ratio and coupling
$g_{\mu/\tau}$ in both nBE and cBE approaches. We take $m_{\phi}=0.26$ GeV for
the $\mu^{+}\mu^{-}$ mode, and $m_{\phi}=5$ GeV for the $\tau^{+}\tau^{-}$
mode.

Speaking overall, DM relic density in the cBE method is larger than that by
using the standard nBE method. When the elastic scattering is strongly
suppressed, the temperature of DM particles drops below that of the thermal
bath ($T_{\chi}<T_{\rm SM}$), meaning DM particles do not have enough kinetic
energy. The forbidden annihilation thus becomes ineffective due to the
inability to overcome the annihilation threshold, leading to a larger
abundance in cBE treatment.

##  IV Parameter space and experimental constraints

In this section, we will find the feasible parameter space for the cBE
approach in the $(m_{\phi},~{}g_{\mu,\tau})$ plane, by requiring the correct
DM relic density. The viable parameter space further displays more accurate
results compared with nBE. Of course, there are numerous constraints on the
model parameters that are imposed by the collider searches and the
astrophysical observations across a wide range.

###  IV.1 $\chi\bar{\chi}\to\mu^{+}\mu^{-}$

We first study the $\mu^{+}\mu^{-}$ forbidden channel. The numerical results
are shown in Fig. 6. We set $g_{\mu}^{A}=0$ for the left panel, while
$g_{\mu}^{A}=g_{\mu}$ for the right panel. The different choices of the
couplings to muons do not strongly affect the phenomenology.

To elaborate on our findings, we start by searching the parameter boundary of
the forbidden annihilation, in which $\Delta\to 0$ and the right relic density
are fulfilled. The green line in the plots shows the boundary of the cBE
approach; the green shaded region, denoted by the expression $\Delta(\rm{\bf
cBE})\leqslant 0$, represents the unforbidden space. In the remaining
parameter space, where the forbidden annihilations dominate the DM depletion,
$m_{\chi}$ is chosen at each point to match the correct relic density. The
required coupling becomes larger when $m_{\mu}/m_{\chi}>1$, leading to the
corresponding curves lying inside the funnel area with a similar shape as the
boundary. For comparison, we also repeated the parameter boundary of the
standard nBE approach in the same figure, as shown by the gray lines. The gray
regions stand for the non-forbidden regions.

The allowed parameter space of the two approaches is noticeably distinctive.
The eKD effect reduces the parameter space, compared to that of nBE. Or we can
say that to satisfy the relic density requirement, the larger coupling
$g_{\mu}$ is required for the cBE scenario. The reason is straightforward as
already pointed out in the last section. Around the freeze-out stage, the
temperature of DM decreases resulting in reduced kinetic energy and then
inefficient forbidden annihilations. To maintain the DM dilution process, a
larger coupling is required for $T_{\chi}<T$.

Figure 6: Results and constraints on forbidden channel
$\chi\bar{\chi}\to\mu^{+}\mu^{-}$. The green line represents the parameter
boundary of the forbidden DM region with $m_{\mu}/m_{\chi}=1$; and the green
shaded region stands for the non-forbidden space that $m_{\chi}>m_{\mu}$. The
rest space is allowed for the forbidden annihilations to get the correct DM
relic density. The gray region labeled $\Delta({\rm nBE})\leqslant 0$ is
obtained in the nBE approach. Also, the non-forbidden space that is dominated
by the $\chi\bar{\chi}\to\phi\phi$ process is shown. The brown region displays
the Planck bounds on the energy injection process Slatyer (2016). The orange
region is excluded by the E137 electron beam-dump experiment Bjorken  _et al._
(1988). The purple region is excluded by energy loss process of SN1987A Hirata
_et al._ (1987); Jaeckel  _et al._ (2018); Dolan  _et al._ (2017). The
projected sensitivities for future beam-dump experiments from BDX Bondí (2017)
(light blue), M3 Berlin  _et al._ (2019); Kahn  _et al._ (2018) (dashed blue
line), $\mathrm{NA62}$ Cortina Gil  _et al._ (2017) (red), and NA64-$\mu$ Chen
_et al._ (2018); Gninenko  _et al._ (2015); Gninenko and Krasnikov (2001)
(dashed magenta line) are also shown here.

In the small mediator mass region, the annihilation into pairs of mediators
dominates the relic density. We depict it in gray with the label
$\chi\bar{\chi}\to\phi\phi$ as the non-forbidden DM region. For experimental
constraints, the most crucial parameter is the mediator mass $m_{\phi}$. In
the following, we show the constraints one by one:

  * •

Planck (brown region).

DM annihilations into SM electromagnetically interacting particles could
modify the anisotropies of the CMB Adams  _et al._ (1998); Chen and
Kamionkowski (2004); Padmanabhan and Finkbeiner (2005). The measurements of
the CMB by the Planck satelliteAghanim  _et al._ (2020) can thus put robust
constraints on such energy injection processes Slatyer (2016). In this model,
the photon pairs can be produced in DM annihilations via a muon loop, leading
to the injection into the CMB. We recast the corresponding constraints from
Ref. D’Agnolo  _et al._ (2021) which are shown in brown.

As pointed out in Ref. D’Agnolo  _et al._ (2021), when $m_{\phi}$ near
$2m_{\chi}$ which is also close to but smaller than $2m_{\mu}$, the
annihilation cross section of $\chi\bar{\chi}\to\gamma\gamma$ is enhanced due
to $\sigma(\chi\bar{\chi}\to\gamma\gamma)\sim
1/[(m_{\phi}^{2}-4m_{\chi}^{2})^{2}+m_{\phi}^{2}\Gamma_{\phi}^{2}]$. So the
CMB constraints exclude the most parameter space for $m_{\phi}<2m_{\mu}$ and
only leave a small part of the allowed room when $m_{\phi}>2m_{\mu}$.

From Fig. 6, we find that the eKD effects narrow down the allowed parameter
space that obtains the right amount of DM relic abundance. Especially the CMB
constraints exclude the most parameter space with cBE in the pure scalar
interaction scenario (i.e. $g_{\mu}^{A}=0$).

  * •

E137 (for orange region) and BDX (blue region).

Secondary muons are produced from the electron beam-dump experiments, such as
SLAC E137 Bjorken  _et al._ (1988) and Jefferson Lab BDX experimentsBondí
(2017), which can be used to explore the signals of light scalar emission and
the muon-scalar coupling via muon-nucleon scattering process
$\mu+N\rightarrow\mu+N+\phi$. The E137 experiment’s null result established
exclusion limits on the parameter space, and the upcoming BDX experiment can
likewise yield a predicted limit.

In Fig. 6, we display, as shaded areas, constraints from the E137 electron
beam-dump experiment Marsicano  _et al._ (2018), and projections from BDX
Bondí (2017). The constraints exclude part of the available parameter space in
the ($m_{\phi}-g_{\mu}~{}(g_{\mu}^{A})$) plane.

  * •

$\mathrm{NA62}$ (red shaded region).

NA62 is a fixed-target experiment at the CERN Super Proton Synchrotron (SPS)
that is dedicated to measurements of kaon rare decays, including projected
searches on $K\to\mu\nu\phi$ Cortina Gil  _et al._ (2017). Such a decay
channel is an excellent probe of new light scalars that couples preferentially
to muons. Ref. Krnjaic  _et al._ (2020) has derived the probe sensitivity for
this process, which can be used to test our parameter space, as shown in red.

  * •

NA64-$\mu$ and M3 (magenta and blue dashed lines).

Similar to the above beam-dump experiments, NA64-$\mu$ Chen  _et al._ (2018);
Gninenko  _et al._ (2015); Gninenko and Krasnikov (2001) and ${\rm M}^{3}$
Berlin  _et al._ (2019); Kahn  _et al._ (2018) are designed to search light
scalars in the muon-nucleon scattering process $\mu N\rightarrow\mu N\phi$,
using muon beams. It’s worth noting that the allowed parameter space can be
tested in the coming future.

  * •

SN 1987A (purple region).

At last, for the parameter space involving new light scalars, we should also
consider the constraints from the observation of supernovae cooling. The most
famous constraints arise from the energy loss process via the Primakoff effect
$\gamma p\to p\phi$ in SN1987AHirata  _et al._ (1987); Jaeckel  _et al._
(2018); Dolan  _et al._ (2017). The excluded region is displayed in purple.

###  IV.2 $\chi\bar{\chi}\to\tau^{+}\tau^{-}$

Figure 7: Results and constraints on forbidden channel
$\chi\bar{\chi}\to\tau^{+}\tau^{-}$. The green and gray shaded regions have
the same meaning as in Fig. 6, and we use mass ratio $m_{\tau}/m_{\chi}=1$. We
show constraints from searches at BaBar for $e^{+}e^{-}\rightarrow\phi\gamma$
Dolan  _et al._ (2017); Lees  _et al._ (2017); Chen  _et al._ (2018) (orange
region), searches at LEP for $Z\rightarrow\bar{\tau}\tau+$MET Chen  _et al._
(2018); Tanabashi  _et al._ (2018) (brown region), and Planck constraints on
DM annihilations (blue). We include projections for Belle II Dolan  _et al._
(2017); Abe  _et al._ (2010) (red region) and future $Z$-factories Liu  _et
al._ (2018); d’Enterria (2016, 2017); Ahmad  _et al._ (2015) (purple region).

The forbidden annihilation in the di-tau case is discussed here. Annihilation
to $\tau^{+}\tau^{-}$ shares several qualitative features with annihilations
to muons. The same computations are performed, including the searches for the
forbidden annihilation parameter region, the comparison of cBE and nBE
treatments, and the variety of limitations from experimental searches. The
results are displayed in Fig. 7 with different parameter settings. Note that
the consideration of the experimental constraints is rather different from the
$\chi\bar{\chi}\to\mu^{+}\mu^{-}$ case due to the different mediator mass.

  * •

BaBar (orange region) and Belle II (red region)

Searching for $e^{+}e^{-}\to\gamma+{\rm invisible}$ at the $e^{+}e^{-}$
colliders, such as BaBar and the future Belle II experiments, sets existed and
projected constraints on the parameter space. As depicted in Fig. 7, the
orange line and region show the constraints arising from BaBar Dolan  _et al._
(2017); Lees  _et al._ (2017); Chen  _et al._ (2018), and the red region shows
the exquisite sensitivity from Belle II using 50 ab-1 integrated luminosity
Dolan  _et al._ (2017); Abe  _et al._ (2010), which indicates that a large
portion of the viable parameter space will be tested.

  * •

LEP (brown region) and Tera Z (purple region)

Precision measurements of $Z$ boson decay width can place constraints on the
model parameters, namely on the exotic $Z$ decays. In this scenario,
$Z\to\tau^{+}\tau^{-}\phi~{}(\phi\to{\rm invisible})$ contributed to the
measured $Z\to\tau^{+}\tau^{-}$ width. Large Electron Positron Collider (LEP)
has set corresponding limits on such channels Chen  _et al._ (2018); Tanabashi
_et al._ (2018) which are shown in brown.

Additionally, there have been several proposals for future $Z$-factories to
search the same processes Liu  _et al._ (2018); d’Enterria (2016, 2017); Ahmad
_et al._ (2015), based on Circular Electron Positron Collider (CEPC) and the
Future Circular Collider $e^{+}e^{-}$ (FCC-ee) for instance. The projections
of Tera-Z options (with accumulated $10^{12}$ Z’s ) provide leading
sensitivities in the tens GeV range, shown in purple.

  * •

Planck (blue region)

The CMB constraints are similar with the $\chi\bar{\chi}\to\mu^{+}\mu^{-}$
forbidden channel. The parameter space at $m_{\phi}<2m_{\tau}$ is much
constrained as always, as shown in Fig. 7, due to the enhancement of
$\chi\bar{\chi}\to\gamma\gamma$ when $m_{\phi}\simeq 2m_{\chi}$.

##  V Conclusion

In the forbidden DM scenario, consideration of the early kinetic decoupling is
not only a correction to the DM relic density but an indispensable ingredient.
In this work, we investigate the early kinetic decoupling effect in forbidden
channels, where DM annihilation to SM leptons is kinetically forbidden.
Specifically we focus on the $\chi\bar{\chi}\to\mu^{+}\mu^{-}$ and
$\chi\bar{\chi}\to\tau^{+}\tau^{-}$ modes. By analyzing the scattering
momentum transfer rate, we found the kinetic equilibrium breaks at about
$x\simeq 20$, which is the same stage of chemical decoupling. So eKD should be
taken seriously into account. With different benchmark points, we found that
there is a cooling phase during the evolution that causes the DM temperature
to deviate from the thermal bath and evolve solely. The decreased kinetic
energy of DM particles suppresses the forbidden annihilation rate, which gives
rise to larger abundances. The difference between the cBE and traditional nBE
methods is significant, within part of the parameter space showing a deviation
up to an order of magnitude larger, for both $\mu^{+}\mu^{-}$ and
$\tau^{+}\tau^{-}$ forbidden channels. We also considered the experimental
constraints from beam-dump experiments, collider searches, and astrophysical
observations. The viable parameter space in the forbidden DM model has been
reduced when using the cBE treatment. Most of the parameter space will be
tested by the forthcoming experimental searches.

## ACKNOWLEDGEMENTS

We thank Murat Abdughani for the helpful discussions. This work was supported
by the National Natural Science Foundation of China under Grants No. 12005180
and No. 12275232 and by the Natural Science Foundation of Shandong Province
under Grant No. ZR2020QA083.

##  Appendix A Cross sections for annihilation and scattering processes

The annihilation cross section for $\chi\bar{\chi}\rightarrow l^{+}l^{-}$ is

$\sigma=\frac{\Big{(}\sqrt{s-4m_{l}^{2}}\big{(}g^{A2}_{l}s+g_{l}^{2}(s-4m_{l}^{2})\big{)}\big{(}g^{A2}_{\chi}s+g_{\chi}^{2}(s-4m_{\chi}^{2})\big{)}\Big{)}}{\Big{(}16\pi
s\sqrt{s-4m_{\chi}^{2}}\big{(}(s-m_{\phi}^{2})^{2}+m_{\phi}^{2}\Gamma_{\phi}^{2}\big{)}\Big{)}},$
(17)

where the total decay rate of scalar $\phi$ is calculated as

$\begin{split}\Gamma_{\phi}={}&\frac{1}{16\pi
m_{\phi}}\sqrt{1-\frac{4m_{l}^{2}}{m_{\phi}^{2}}}\times\\\
&\big{(}(g_{l}-g^{A}_{l})^{2}\big{(}m_{\phi}^{2}-2m_{l}^{2}\big{)}+(g^{A}_{l}+g_{l})^{2}\big{(}m_{\phi}^{2}-2m_{l}^{2}\big{)}-4m_{l}^{2}(g^{A}_{l}+g_{l})(g_{l}-g^{A}_{l})\big{)}\\\
&+\frac{1}{16\pi m_{\phi}}\sqrt{1-\frac{4m_{\chi}^{2}}{m_{\phi}^{2}}}\times\\\
&\big{(}(-g^{A}_{\chi}-g_{\chi})^{2}\big{(}m_{\phi}^{2}-2m_{\chi}^{2}\big{)}+(g^{A}_{\chi}-g_{\chi})^{2}\big{(}m_{\phi}^{2}-2m_{\chi}^{2}\big{)}+4m_{\chi}^{2}(g^{A}_{\chi}-g_{\chi})(g^{A}_{\chi}+g_{\chi})\big{)}\end{split}$
(18)

For the elastic scattering process $\chi\mu^{\pm}\rightarrow\chi\mu^{\pm}$,
the amplitude is

$|\mathcal{M}|^{2}=\frac{(g^{A2}_{l}t+g_{l}^{2}(t-4m_{l}^{2}))(g^{A2}_{\chi}t+g_{\chi}^{2}(t-4m_{\chi}^{2}))}{(m_{\phi}^{2}-t)^{2}}.$
(19)

And the momentum transfer rate can be written as,

$\displaystyle\gamma=\frac{1}{3g_{\chi}m_{\chi}T}\\!\int\\!\frac{\text{d}^{3}k}{(2\pi)^{3}}f_{\cal
B}^{\pm}(E_{k})\left[1\\!\mp\\!f_{\cal
B}^{\pm}(E_{k})\right]\\!\\!\\!\int\limits^{0}_{-4k_{\mathrm{cm}}^{2}}\\!\\!\\!\text{d}t(-t)\frac{\text{d}\sigma}{\text{d}t}v\,,$
(20)

with

$\int_{-4k_{cm}^{2}}^{0}dt(-t)\frac{d\sigma}{dt}v=\frac{1}{64\pi
kE_{k}m_{\chi}^{2}}\int_{-4k_{cm}^{2}}^{0}dt(-t){|\mathcal{M}|^{2}}=\frac{1}{64\pi
kE_{k}m_{\chi}^{2}}\times t{\rm Amp},$ (21)

and

$\displaystyle t{\rm Amp}$
$\displaystyle=4k^{2}m_{\chi}^{2}\Bigg{(}\frac{2(g^{A2}_{l}+g_{l}^{2})(g^{A2}_{\chi}+g_{\chi}^{2})k^{2}m_{\chi}^{2}}{(m_{\chi}^{2}+m_{\mu}^{2}+2m_{\chi}E_{k})^{2}}+{}$
(22)
$\displaystyle\frac{4(g^{A2}_{l}+g_{l}^{2})g_{\chi}^{2}m_{\chi}^{2}+4g_{l}^{2}(g^{A2}_{\chi}+g_{\chi}^{2})m_{\mu}^{2}-2(g^{A2}_{l}+g_{l}^{2})(g^{A2}_{\chi}+g_{\chi}^{2})m_{\phi}^{2}}{m_{\chi}^{2}+m_{l}^{2}+2m_{\chi}E_{k}}+{}$
$\displaystyle\frac{\big{(}4g_{l}^{2}m_{l}^{2}-(g^{A2}_{l}+g_{l}^{2})m_{\phi}^{2}\big{)}\big{(}-4g_{\chi}^{2}m_{\chi}^{2}+(g^{A2}_{\chi}+g_{\chi}^{2})m_{\phi}^{2}\big{)}}{4k^{2}m_{\chi}^{2}+m_{\phi}^{2}(m_{\chi}^{2}+m_{l}^{2}+2m_{\chi}E_{k})}\Bigg{)}-{}$
$\displaystyle\Big{(}16g_{l}^{2}g_{\chi}^{2}m_{\chi}^{2}m_{\mu}^{2}-8\big{(}(g^{A2}_{l}+g_{l}^{2})g_{\chi}^{2}m_{\chi}^{2}+g_{l}^{2}(g^{A2}\chi+g_{\chi}^{2})m_{l}^{2}\big{)}m_{\phi}^{2}+{}$
$\displaystyle
3(g^{A2}_{l}+g_{l}^{2})(g^{A2}_{\chi}+g_{\chi}^{2})m_{\phi}^{4}\Big{)}\Bigg{(}\log
m_{\phi}^{2}-\log(m_{\phi}^{2}+\frac{4k^{2}m_{\chi}^{2}}{m_{\chi}^{2}+m_{l}^{2}+2m_{\chi}E_{k}})\Bigg{)}.$

## References

  * Lee and Weinberg (1977) B. W. Lee and S. Weinberg, Phys. Rev. Lett. 39, 165 (1977).
  * Hut (1977) P. Hut, Phys. Lett. B 69, 85 (1977).
  * Sato and Kobayashi (1977) K. Sato and M. Kobayashi, Prog. Theor. Phys. 58, 1775 (1977).
  * Dicus  _et al._ (1977) D. A. Dicus, E. W. Kolb, and V. L. Teplitz, Phys. Rev. Lett. 39, 168 (1977), [Erratum: Phys.Rev.Lett. 39, 973 (1977)].
  * Wolfram (1979) S. Wolfram, Phys. Lett. B 82, 65 (1979).
  * Gondolo and Gelmini (1991) P. Gondolo and G. Gelmini, Nucl. Phys. B 360, 145 (1991).
  * Edsjo and Gondolo (1997) J. Edsjo and P. Gondolo, Phys. Rev. D 56, 1879 (1997), arXiv:hep-ph/9704361 .
  * Bringmann (2009) T. Bringmann, New J. Phys. 11, 105027 (2009), arXiv:0903.0189 [astro-ph.CO] .
  * van den Aarssen  _et al._ (2012) L. G. van den Aarssen, T. Bringmann, and Y. C. Goedecke, Phys. Rev. D 85, 123512 (2012), arXiv:1202.5456 [hep-ph] .
  * Binder  _et al._ (2017) T. Binder, T. Bringmann, M. Gustafsson, and A. Hryczuk, Phys. Rev. D 96, 115010 (2017), [Erratum: Phys.Rev.D 101, 099901 (2020)], arXiv:1706.07433 [astro-ph.CO] .
  * Brümmer (2020) F. Brümmer, JHEP 01, 113 (2020), arXiv:1910.01549 [hep-ph] .
  * Ala-Mattinen and Kainulainen (2020) K. Ala-Mattinen and K. Kainulainen, JCAP 09, 040 (2020), arXiv:1912.02870 [hep-ph] .
  * Abe (2020) T. Abe, Phys. Rev. D 102, 035018 (2020), arXiv:2004.10041 [hep-ph] .
  * Binder  _et al._ (2021) T. Binder, T. Bringmann, M. Gustafsson, and A. Hryczuk, Eur. Phys. J. C 81, 577 (2021), arXiv:2103.01944 [hep-ph] .
  * Zhu and Liu (2022) B. Zhu and X. Liu, Sci. China Phys. Mech. Astron. 65, 231011 (2022), arXiv:2104.03238 [hep-ph] .
  * Hryczuk and Laletin (2021) A. Hryczuk and M. Laletin, JHEP 06, 026 (2021), arXiv:2104.05684 [hep-ph] .
  * Abe (2021) T. Abe, Phys. Rev. D 104, 035025 (2021), arXiv:2106.01956 [hep-ph] .
  * Du  _et al._ (2022) Y. Du, F. Huang, H.-L. Li, Y.-Z. Li, and J.-H. Yu, JCAP 04, 012 (2022), arXiv:2111.01267 [hep-ph] .
  * Ala-Mattinen  _et al._ (2022) K. Ala-Mattinen, M. Heikinheimo, K. Kainulainen, and K. Tuominen, Phys. Rev. D 105, 123005 (2022), arXiv:2201.06456 [hep-ph] .
  * Hryczuk and Laletin (2022) A. Hryczuk and M. Laletin, Phys. Rev. D 106, 023007 (2022), arXiv:2204.07078 [hep-ph] .
  * D’Agnolo and Ruderman (2015) R. T. D’Agnolo and J. T. Ruderman, Phys. Rev. Lett. 115, 061301 (2015), arXiv:1505.07107 [hep-ph] .
  * Delgado  _et al._ (2017) A. Delgado, A. Martin, and N. Raj, Phys. Rev. D 95, 035002 (2017), arXiv:1608.05345 [hep-ph] .
  * D’Agnolo  _et al._ (2021) R. T. D’Agnolo, D. Liu, J. T. Ruderman, and P.-J. Wang, JHEP 06, 103 (2021), arXiv:2012.11766 [hep-ph] .
  * Wojcik and Rizzo (2022) G. N. Wojcik and T. G. Rizzo, JHEP 04, 033 (2022), arXiv:2109.07369 [hep-ph] .
  * Yang (2022) K.-C. Yang, JHEP 11, 083 (2022), arXiv:2209.10827 [hep-ph] .
  * Cheng  _et al._ (2022) Y. Cheng, S.-F. Ge, X.-G. He, and J. Sheng, (2022), arXiv:2211.05643 [hep-ph] .
  * Herms  _et al._ (2022a) J. Herms, S. Jana, V. P. K., and S. Saad, Phys. Rev. Lett. 129, 091803 (2022a), arXiv:2203.05579 [hep-ph] .
  * Herms  _et al._ (2022b) J. Herms, S. Jana, V. P. K., and S. Saad, in  _14th International Workshop on the Identification of Dark Matter 2022_ (2022) arXiv:2209.14882 [hep-ph] .
  * Binder  _et al._ (2016) T. Binder, L. Covi, A. Kamada, H. Murayama, T. Takahashi, and N. Yoshida, JCAP 11, 043 (2016), arXiv:1602.07624 [hep-ph] .
  * Bertschinger (2006) E. Bertschinger, Phys. Rev. D 74, 063509 (2006), arXiv:astro-ph/0607319 .
  * Bringmann and Hofmann (2007) T. Bringmann and S. Hofmann, JCAP 04, 016 (2007), [Erratum: JCAP 03, E02 (2016)], arXiv:hep-ph/0612238 .
  * Gondolo  _et al._ (2012) P. Gondolo, J. Hisano, and K. Kadota, Phys. Rev. D 86, 083523 (2012), arXiv:1205.1914 [hep-ph] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Sabti  _et al._ (2020) N. Sabti, J. Alvey, M. Escudero, M. Fairbairn, and D. Blas, JCAP 01, 004 (2020), arXiv:1910.01649 [hep-ph] .
  * Slatyer (2016) T. R. Slatyer, Phys. Rev. D 93, 023527 (2016), arXiv:1506.03811 [hep-ph] .
  * Bjorken  _et al._ (1988) J. Bjorken, S. Ecklund, W. Nelson, A. Abashian, C. Church, B. Lu, L. Mo, T. Nunamaker, and P. Rassmann, Phys. Rev. D 38, 3375 (1988).
  * Hirata  _et al._ (1987) K. Hirata _et al._ (Kamiokande-II), Phys. Rev. Lett. 58, 1490 (1987).
  * Jaeckel  _et al._ (2018) J. Jaeckel, P. C. Malta, and J. Redondo, Phys. Rev. D 98, 055032 (2018), arXiv:1702.02964 [hep-ph] .
  * Dolan  _et al._ (2017) M. J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer, and K. Schmidt-Hoberg, JHEP 12, 094 (2017), [Erratum: JHEP 03, 190 (2021)], arXiv:1709.00009 [hep-ph] .
  * Bondí (2017) M. Bondí (BDX), EPJ Web Conf. 142, 01005 (2017).
  * Berlin  _et al._ (2019) A. Berlin, N. Blinov, G. Krnjaic, P. Schuster, and N. Toro, Phys. Rev. D 99, 075001 (2019), arXiv:1807.01730 [hep-ph] .
  * Kahn  _et al._ (2018) Y. Kahn, G. Krnjaic, N. Tran, and A. Whitbeck, JHEP 09, 153 (2018), arXiv:1804.03144 [hep-ph] .
  * Cortina Gil  _et al._ (2017) E. Cortina Gil _et al._ (NA62), JINST 12, P05025 (2017), arXiv:1703.08501 [physics.ins-det] .
  * Chen  _et al._ (2018) C.-Y. Chen, J. Kozaczuk, and Y.-M. Zhong, JHEP 10, 154 (2018), arXiv:1807.03790 [hep-ph] .
  * Gninenko  _et al._ (2015) S. N. Gninenko, N. V. Krasnikov, and V. A. Matveev, Phys. Rev. D 91, 095015 (2015), arXiv:1412.1400 [hep-ph] .
  * Gninenko and Krasnikov (2001) S. N. Gninenko and N. V. Krasnikov, Phys. Lett. B 513, 119 (2001), arXiv:hep-ph/0102222 .
  * Adams  _et al._ (1998) J. A. Adams, S. Sarkar, and D. W. Sciama, Mon. Not. Roy. Astron. Soc. 301, 210 (1998), arXiv:astro-ph/9805108 .
  * Chen and Kamionkowski (2004) X.-L. Chen and M. Kamionkowski, Phys. Rev. D 70, 043502 (2004), arXiv:astro-ph/0310473 .
  * Padmanabhan and Finkbeiner (2005) N. Padmanabhan and D. P. Finkbeiner, Phys. Rev. D 72, 023508 (2005), arXiv:astro-ph/0503486 .
  * Marsicano  _et al._ (2018) L. Marsicano, M. Battaglieri, A. Celentano, R. De Vita, and Y.-M. Zhong, Phys. Rev. D 98, 115022 (2018), arXiv:1812.03829 [hep-ex] .
  * Krnjaic  _et al._ (2020) G. Krnjaic, G. Marques-Tavares, D. Redigolo, and K. Tobioka, Phys. Rev. Lett. 124, 041802 (2020), arXiv:1902.07715 [hep-ph] .
  * Lees  _et al._ (2017) J. P. Lees _et al._ (BaBar), Phys. Rev. Lett. 119, 131804 (2017), arXiv:1702.03327 [hep-ex] .
  * Tanabashi  _et al._ (2018) M. Tanabashi _et al._ (Particle Data Group), Phys. Rev. D 98, 030001 (2018).
  * Abe  _et al._ (2010) T. Abe _et al._ (Belle-II), (2010), arXiv:1011.0352 [physics.ins-det] .
  * Liu  _et al._ (2018) J. Liu, L.-T. Wang, X.-P. Wang, and W. Xue, Phys. Rev. D 97, 095044 (2018), arXiv:1712.07237 [hep-ph] .
  * d’Enterria (2016) D. d’Enterria, Frascati Phys. Ser. 61, 17 (2016), arXiv:1601.06640 [hep-ex] .
  * d’Enterria (2017) D. d’Enterria, in  _17th Lomonosov Conference on Elementary Particle Physics_ (2017) pp. 182–191, arXiv:1602.05043 [hep-ex] .
  * Ahmad  _et al._ (2015) M. Ahmad _et al._ , in  _CEPC-SPPC Preliminary Conceptual Design Report. 1. Physics and Detector_ (2015).


# Quantum-confinement and Structural Anisotropy result in Electrically-Tunable
Dirac Cone in Few-layer Black Phosphorous

Kapildeb Dolui and Su Ying Quek  Department of Physics, Centre for Advanced 2D
Materials and Graphene Research Centre, National University of Singapore, 6
Science Drive 2, Singapore 117546, Singapore.

###### Abstract

Abstract

2D materials are well-known to exhibit interesting phenomena due to quantum
confinement. Here, we show that quantum confinement, together with structural
anisotropy, result in an electric-field-tunable Dirac cone in 2D black
phosphorus. Using density functional theory calculations, we find that an
electric field, $E_{\rm ext}$, applied normal to a 2D black phosphorus thin
film, can reduce the direct band gap of few-layer black phosphorus, resulting
in an insulator-to-metal transition at a critical field, $E_{c}$. Increasing
$E_{\rm ext}$ beyond $E_{c}$ can induce a Dirac cone in the system, provided
the black phosphorus film is sufficiently thin. The electric field strength
can tune the position of the Dirac cone and the Dirac-Fermi velocities, the
latter being similar in magnitude to that in graphene. We show that the Dirac
cone arises from an anisotropic interaction term between the frontier orbitals
that are spatially separated due to the applied field, on different halves of
the 2D slab. When this interaction term becomes vanishingly small for thicker
films, the Dirac cone can no longer be induced. Spin-orbit coupling can gap
out the Dirac cone at certain electric fields; however, a further increase in
field strength reduces the spin-orbit-induced gap, eventually resulting in a
topological-insulator-to-Dirac-semi-metal transition.

2D materials, Black phosphorus, Electric field effects, Band gap tuning, Dirac
semimetal.

## Introduction

Two-dimensional (2D) layered materials, where interlayer interactions are
dominated by weak van der Waals (vdW) forces, has attracted tremendous
attention in nanoelectronics Wang . Of particular interest are quantum
confinement effects in which the electronic structure of the 2D layered
material changes abruptly as a function of the number of layers. A well-known
example is the indirect to direct band gap transition when 2D MoS2 is reduced
to a single layer Splendiani . Here, we show that quantum confinement results
in the emergence of an electrically tunable Dirac cone in a recently
discovered 2D material - 2D black phosphorus.

Black phosphorus, a layered semiconductor material, is the thermodynamically
stable form of phosphorus. Different experimental techniques, such as
mechanical Li  and liquid Coleman  exfoliation, have been employed to thin
down the bulk to a monolayer. Experimentalists have succeeded also in
measuring a thickness-dependent transport gap, ranging from $\sim$0.30 eV
(bulk) to 1.0 eV (monolayer) Saptarshi , and in fabricating field effect
transistors based on ultrathin black phosphorus Li . Such transistors exhibit
high mobility of $\sim 1000~{}{\rm cm^{2}/V\cdot s}$ and appreciably high
on/off ratios, up to 104, opening up a possible potential application in
nanoelectronics Li . In contrast to previously studied 2D materials such as
graphene and MoS2, 2D black phosphorus layers are highly anisotropic ( see
Figure 1), with zigzag chains in one direction (y in the figure) and armchair
chains in the other direction (x). This anisotropy adds another dimension of
interest to 2D materials, resulting in direction-dependent optical and
transport properties Xia .

Figure 1: (Color online) The optimized structure of a bilayer (2L) black
phosphorus. (a) and (b) show the top and side view of 2L black phosphorus. An
external vertical electric field, $E_{\mathrm{ext}}$ is applied along $z$-axis
and $2U$ denotes the corresponding potential difference across the layers.
$\Delta$ is the potential difference between the two vertically inequivalent
sites in each layer. For the sake of simplicity we use $|\Delta|=U/2$ for our
tight-binding model of the bilayer system. In (a), hopping parameters of the
tight-binding model are indicated. $t$ and $t^{\prime}$ denote intralayer and
interlayer hopping parameters respectively. In (b), dark and light yellow
correspond to two different layers of black phosphorus. Figure 2: (Color
online) Electronic structure of 4L black phosphorus as a function of an
external electric field, $E_{\mathrm{ext}}$. In (a) we present: (upper panel)
the band structure and (lower panel) the charge density, respectively at the
CBM (top) and VBM (bottom) for $E_{\mathrm{ext}}=0$. In (b) and (c), the same
quantities are shown for $E_{\mathrm{ext}}$ = 0.31 and 0.72 V/Å, respectively.
Note that the critical field, $E_{c}$, is $\sim$ 0.48 V/Å for the 4L black
phosphorus. In (c), where $E_{\rm ext}>E_{c}$, a Dirac cone appears along the
$k_{y}$ direction, and a band gap of $\sim$39 meV opens along the $k_{x}$
direction. Inset in (c): zoomed-in bandstructure around the $\Gamma$ point
near the Fermi level, $E_{F}$. Band inversion is evident in (c).

In this work, we explore from first principles calculations the effect of
external electric fields applied perpendicular to thin films of black
phosphorus. Similar to the Stark effect previously predicted for boron nitride
nanotubes, two-dimensional transition metal dichalcogenides, and other 1D and
2D materials Ashwin ; Khoo ; Kapildeb ; Park ; Son ; Gava  we find that the
electric field localizes the valence band maximum (VBM) and conduction band
minimum (CBM) states at opposite surfaces of the slab. The resulting potential
difference between the VBM and CBM states results in a decrease in the band
gap, eventually leading to an insulator-to-metal transition at a critical
applied field, $E_{c}$. However, in contrast to other materials, we also find
that for black phosphorus films below a certain critical thickness, a further
increase in applied field strength results in a highly anisotropic opening in
the band gap that leads to a formation of a Dirac cone. This result is not
observed for black phosphorus films that are too thick. We show that the
emergent Dirac cone physics arises from anisotropic interaction terms between
the electric-field-induced quantum-confined VBM and CBM states in few layer
black phosphorus.

## RESULTS AND DISCUSSION

Electronic structure at zero $\bf E_{\rm ext}$ We begin our study by
calculating the electronic structure of pristine black phosphorus at zero
external electric field. Electronic structure calculations are performed by
using density functional theory (DFT) within the Perdew-Burke-Ernzerhof
parametrization (PBE) PBE  of the generalized gradient approximation (GGA) of
the exchange and correlation potential, as implemented in the Quantum Espresso
QE  package. To treat van der Waals (vdW) interactions, we employ the vdW+DF
approach for the exchange-correlation functional vdW (See Methods for further
details).

Bulk black phosphorus has an orthorhombic crystal structure; the optimized
lattice parameters of bulk black phosphorus are $a$=3.35 Å, $b$=4.68 Å, and
$c$=11.42 Å, in good agreement with the corresponding experimental values
($a$=3.31 Å, $b$=4.37 Å, and $c$=10.47 Å) expt-lattice . Bulk black phosphorus
is a semiconductor and we predict a direct band gap of 0.43 eV at the $\Gamma$
point. Although this band gap value is slightly larger than the experimental
value of 0.35 eV Warschauer , it is in good agreement with other vdW-corrected
density functional theory calculations Qiao . For thin film black phosphorus,
we find that the structural parameters do not change significantly. On the
other hand, the band gap increases monotonically with decreasing thickness,
reaching 0.99 eV for the monolayer. Interestingly, the band gap remains direct
at the $\Gamma$ point when the thickness reduces from bulk to monolayer limit.

As previously mentioned, black phosphorus is a layered compound consisting of
puckered phosphorene layers with highly anisotropic bonding within each layer
(Figure 1). We denote multi-layered black phosphorus by $n$L-BP, where $n$ is
the number of layers. 2(a) shows the electronic structure of 4L-BP, which
exhibits a direct band gap of 0.60 eV at the $\Gamma$ point. The valence band
maximum (VBM) and conduction band minimum (CBM) charge density are both
delocalized over the layers. Our analysis of the projected density of states
indicates that the charge density at the VBM and the CBM originate from
phosphorus 3$s$, 3$p_{x}$ and 3$p_{z}$ orbitals (see Supplementary Figure S1).
The $p_{z}$ orbitals have the largest contribution ($\textgreater$ 70%) to the
VBM and CBM among all the orbitals.

Figure 3: (Color online) Variation of the (PBE+vdW) band gap with the applied
electric field, $E_{\mathrm{ext}}$, for 2L (black circles), 3L (red squares),
4L (green triangles), 5L (blue stars), 6L (cyan triangles), 7L (magenta
diamonds), and 8L (orange triangles) -BP. Inset: variation of the Stark effect
coefficient, $S_{nL}$, with number of layers. Note that we consider the linear
region of the band gap versus $E_{\rm ext}$ curve to compute $S_{nL}$, and the
fitted linear part is indicated by the corresponding dashed line.

Electronic structure for $\bf E_{\rm ext}<E_{c}$ We now discuss the response
of the electronic structure of $n$L-BP to a static external electric field
with magnitude $E_{\rm ext}$, applied normal to the layers. As $E_{\rm ext}$
increases, the band gap decreases monotonically while remaining direct at
$\Gamma$, and eventually closes at a critical field, $E_{c}$ (see Figure 2(b)
and Figure 3). This phenomenon essentially arises from the Stark effect: the
applied electric field lifts the degeneracy among distinct layers in the thin
film. After the degeneracy is lifted, states that are localized in regions
with higher potential will have a higher energy due to electrostatic
interaction with the external field, while those localized in regions with
lower potential will have a lower energy. Consequently, the VBM becomes
localized on the surface layer with higher potential, while the CBM becomes
localized on the opposite surface layer with lower potential, resulting in a
decrease in band gap. We note that, strictly speaking, the VBM and CBM states
also have tails into the interior layers (Figure 2(b-c)). These tails are
important for the VBM and CBM states to interact when they become energy
degenerate at the critical field (see later).

Figure 3 shows the modulation of band gap as a function of $E_{\rm ext}$ for
different $n$L-BP, where $n$ runs from 2 to 8. We find that the bandgap
reduces very rapidly with $E_{\rm ext}$ as the thickness gets larger. Such a
thickness dependence can be easily understood by considering the potential
drop across the layers. If one simply assumes that the potential drops
linearly across the layers, a larger thickness requires a smaller electric
field to sustain the same potential difference between the two edges.
Furthermore, the band gap at zero field also decreases with increasing
thickness. As such, the magnitude of the critical field, $E_{c}$, decreases
with thickness. Indeed we observe that the value of $E_{c}$ drops to 0.29 eV/Å
for 8L-BP. In practice, the potential drop is not strictly linear in space due
to screening within the layers; however, an analysis of the screening effect
shows that the potential drop in the interlayer region is much larger than
that in the intralayer region, and the average potential in each layer drops
linearly across the layers (see Supplementary Figure S2).

We quantify the above analysis as follows. Under the assumption that the
potential drops linearly across the layers, the external field, $E_{\rm ext}$,
induces a potential difference, $\Delta V$, which can be expressed as $\Delta
V=(eE_{\rm ext}/\kappa)\left<Z\right>$, where $\kappa$ is the dielectric
constant along the direction of the applied field and $\left<Z\right>$ is
defined as $\left<Z\right>=\left<Z\right>_{\rm V}-\left<Z\right>_{\rm C}$.
Here $\left<Z\right>_{\rm V}$ and $\left<Z\right>_{\rm C}$ denote the center
positions (along the $z$-axis) of the localised VBM and CBM states,
respectively Fawei . Since the VBM and CBM localize at opposite surfaces for
large $E_{\rm ext}$, $\left<Z\right>=n$c$-2\delta$, where $c$ is the layer-
layer distance and $\delta$ is the position of $\left<Z\right>_{\rm V,C}$ from
the surfaces Fawei . We note that we are assuming here that $\delta$ is
independent of the applied field and film thickness. Although this is not
likely to be valid, deviations should be small compared to $\left<Z\right>$.
Under this model, we can write the change in band gap, $\Delta E_{\rm g}$, as
follows:

$\displaystyle\Delta E_{g}=e\Delta V=\frac{eE_{\rm
ext}}{\kappa}\left<Z\right>.$ (1)

In general for the critical field, $E_{c}$, we obtain from equation 1 that

$\displaystyle-E_{\rm
g}=\frac{eE_{c}}{\kappa}\left<Z\right>=\frac{eE_{c}}{\kappa}(nc-2\delta).$ (2)

Indeed, we find that $E_{g}/E_{c}$ versus $n$ is approximately linear (see
Supplementary Figure S3), with deviations resulting from the assumption that
$\delta$ is constant.

In the initial response to $E_{\rm ext}$, the band gap varies non-linearly
(likely quadratic) with $E_{\rm ext}$, and the slope of the curve vanishes at
$E_{\rm ext}=0$ (Figure 3). This is consistent with the fact that the
variation in band gap should be independent of the polarity of the applied
field, so that in this low field regime, the above linear analysis breaks
down. However, for larger electric fields, for instance $(E_{\rm
ext}$/$E_{c}$)> 0.1 for a specific thickness, $\Delta E_{\rm g}$ is an
approximately linear function of $E_{\rm ext}$ (Figure 3) and this linear
slope can be described as

$\displaystyle\frac{dE_{\rm g}}{dE_{\rm ext}}=-eS_{n\rm L},$ (3)

where $e$ is the electron charge and $S_{n\rm L}$ is defined as the linear
giant stark effect (GSE) coefficient. Combining Equation 1 with Equation 3,
the GSE coefficient can be written as

$\displaystyle S_{n\rm
L}=\frac{\left<Z\right>}{\kappa}=\frac{(nc-2\delta)}{\kappa}.$ (4)

Indeed, the inset of Figure 3 shows that the Stark coefficient increases
approximately linearly as a function of the number of layers.

Figure 4: (Color online) (a) The PBE+vdW bandstructure of 4L-BP is plotted in
the 2D brillouin zone for $E_{\rm ext}$ = 0.72 V/Å, above the critical field.
The Dirac cone appears at $\pm\Lambda$. (b) The PBE+vdW bandstructure along
Y-$\Gamma$-X, for the three different $E_{\mathrm{ext}}$ above the critical
value, namely $E_{\mathrm{ext}}$ = 0.51 (black), 0.72 (red), and 0.93 (green)
V/Å. (c) The bandgap that opens along $\Gamma$-X direction, $E_{g}^{\Gamma-X}$
, is plotted against ($E_{\rm ext}-E_{c}$) for 4L-BP (blue line) and 8L-BP
(red line). Note that the values of $E_{c}$ for 4L-BP and 8L-BP are 0.48 and
0.28 V/Å, respectively. (d) The Fermi velocity along $k_{x}$, $v_{F}^{x}$ (red
dashed line), and along $k_{y}$, $v_{F}^{y}$ (red solid line) are plotted as a
function of ($E_{\rm ext}-E_{c}$) for 4L-BP. $v_{\rm F}$ is calculated as
$v_{\rm F}=\delta E/\delta k$ at the Dirac point.

At the critical field, the system becomes metallic when the VBM and CBM touch
at $\Gamma$. Unlike the electrically driven Stoner magnetism previously
reported for MoS2 nanoribbons Kapildeb , there is no spin polarization in this
system (Supplementary Figure S4). This is because the density of states at the
Fermi level is much smaller in the 2D system than in the 1D nanoribbons.

Electronic structure beyond critical field An interesting and unique effect is
observed in very thin-layer BPs (for example in 4L-BP), when an electric field
larger than the critical field is applied. Specifically, a Dirac cone appears
at a $k$-point along $\Gamma$-Y, which we shall call $\Lambda$, as shown in
Figure 2(c) and Figure 4(a-b). The valence band manifold and conduction band
manifold touch at exactly two points, $\pm\Lambda$, in the 2D Brillouin Zone
[Figure 4(a)]. By plotting the charge density of the highest occupied and
lowest unoccupied states at k-points close to the Dirac point(see
Supplementary Figure S5), we see that the bands are inverted at the Dirac
point, indicating that we have a new topological phase. As black phosphorus
films exhibit a direct band gap even in the presence of the external electric
fields, the band inversion in itself may not be too surprising. However, the
fact that the valence and conduction band manifolds touch at only one point,
the Dirac point, is intriguing and deserves further analysis.

As $E_{\rm ext}$ increases beyond the critical value, the Dirac point,
$\Lambda$, shifts progressively from $\Gamma$ toward Y in the 2D Brillouin
zone [Figure 4(b)]. This shifts occurs simultaneously with an enhancement of
the corresponding Dirac Fermi velocity, $v_{\rm F}$ [Figure 4(c)]. The
existence of the Dirac point is associated with an asymmetry between the
$\Gamma$-Y and $\Gamma$-X directions - so that the valence and conduction
bands touch only at the Dirac points and not in a ring in the 2D Brillouin
Zone; neither does a band gap open up throughout the 2D Brillouin Zone. We
quantify this asymmetry and therefore also the robustness of the Dirac
electronic structure by the band gap along $\Gamma$-X, $E_{\rm g}^{\Gamma-X}$.
As shown in Figure 4(c), $E_{\rm g}^{\Gamma-X}$ increases initially as $E_{\rm
ext}$ increases beyond the critical field; however, as $E_{\rm ext}$ continues
to increase, $E_{\rm g}^{\Gamma-X}$ decreases. For 4L-BP, $E_{\rm
g}^{\Gamma-X}$ reaches a maximum value of $\sim$39 meV at $E_{\rm ext}$=0.72
V/Å , where the corresponding $x~{}{\rm and}~{}y$-components of $v_{F}$ are
0.84$\times~{}10^{6}$ and 0.56$\times 10^{6}$ m/s, respectively. These Fermi
velocities are as large as the predicted value of $\sim 0.86\times 10^{6}$ m/s
in graphene Daniel . Interestingly, we find that the maximum value of $E_{\rm
g}^{\Gamma-X}$ reduces to 8 meV when the thickness of black phosphorus
increases to 8L [Figure 4(c)], and vanishes for the 12L case (see
Supplementary Figure S6). Salient features that arise from the above analysis
are that $E_{\rm g}^{\Gamma-X}$ cannot be arbitrarily increased with electric
field strength, and that $E_{\rm g}^{\Gamma-X}$ is smaller for thicker black
phosphorus films.

Figure 5: (Color online) The bandstructure of bilayer black phosphorus
calculated by using a simple tight-binding model. In (a) the bandstructure is
plotted for $E_{\mathrm{ext}}$ = 0 (black line) and $E_{\mathrm{ext}}$
<$E_{c}$ (red line). In (b) the bandstructures are plotted for three different
$E_{\mathrm{ext}}$ above the critical value, where $U_{1}$<$U_{2}$ <$U_{3}$
<$U_{4}$. For the applied bias $U_{2}$ (green lines) and $U_{3}$ (magenta
lines), the interlayer hopping, $t^{\prime}$, is finite and kept fixed,
whereas $t^{\prime}=0$ for the applied bias $U_{4}$ (cyan dashed lines). The
horizontal dotted line corresponds to the Fermi level. (c) The $k_{x}-k_{y}$
plot of the magnitude of $<\phi_{\rm V}|H(t^{\prime})|\phi_{\rm C}>$, where
$\phi_{\rm V}$ and $\phi_{C}$ are the eigenvectors of the Hamiltonian $H$ at
$t^{\prime}=0$ for the highest occupied and the lowest unoccupied energy
eigenvalues, respectively. We use $t^{\prime}$=0.2 eV for the 2D plot.

To understand the above intriguing phenomena, we develop a simple tight-
binding model that captures the essential physics of black phosphorus thin
films in the presence of an external electric field.

Tight-binding model For simplicity we consider a tight-binding model for BP
films that are two-layer thick (bilayer BP). Since we are concerned with the
valence band (VB) and conduction band (CB) which stem from $p_{z}$ orbitals,
we include only $p_{z}$ orbitals in our model. In the presence of large
external electric fields, the VB and CB are mainly localized at the opposite
surface layers of the thin film; therefore, the bilayer model can effectively
also represent the VB and CB electronic structure in thicker films.

We begin with a standard tight-binding Hamiltonian:

$\displaystyle
H=\sum_{i}(\varepsilon_{i}+U_{i})~{}c_{i}c_{i}^{\dagger}+\sum_{<i,j>}t_{ij}~{}c_{i}c_{j}^{\dagger}+\sum_{<i,j>}t^{\prime}_{ij}~{}c_{i}c_{j}^{\dagger}$
(5)

where $\varepsilon_{i}$ is the energy of $i$-th site, $U_{i}$ is the potential
shift at the $i$-th site upon the applied vertical field (see figure 1),
$c_{i}^{\dagger}$($c_{j}$) is the creation (annihilation) operator of
electrons at site $i(j)$, and $t_{ij}$ and $t^{\prime}_{ij}$ are respectively
the intralayer and interlayer inter-site hopping integrals as depicted in
figure 1. In the first instance, we include hopping terms up to the fourth and
fifth nearest neighbor for the interlayer and intralayer interactions
respectively, using the parameters reported previously in Ref. Rudenko . This
provides an accurate description of the VB and CB, as given by the $GW$
approximation. Figure 5(a) shows that our model nicely reproduces the band
dispersion of the VB and CB in comparison with that of the $GW$ bands Rudenko
.

When an external field is applied, the on-site energy of each site,
$\varepsilon_{i}$, is rigidly shifted by the corresponding potential
difference, $U_{i}$ (see Figure 1). For small $U_{i}$ [ Figure 5(a)] we see
that the band gap closes as described above. When $U_{i}$ increases above a
critical value, $U_{c}$, the Dirac cone emerges along the $\Gamma$-Y
direction, just as was predicted by DFT, with the position of the Dirac cone
shifting with the size of $U_{c}$ (Figure 5(b), green and magenta curves).

However, interestingly, when the interlayer interaction terms $t^{\prime}$ are
set to zero, $U>U_{c}$ causes the band gap to vanish all over the two-
dimensional Brillouin Zone (Figure 5(b), cyan curve). This finding indicates
that the interlayer interaction is critical to the emergence of a Dirac cone.
Intuitively, we expect the VB and CB states to interact via these interlayer
terms, since the states are localized at opposite surfaces of the thin film.
The opening of a band gap everywhere in the Brillouin Zone except along
$\Gamma$-Y thus results from interactions between the VB and CB states,
mediated by the interlayer hopping integrals. This also explains why the
maximum value of $E_{\rm g}^{\Gamma-X}$ decreases as the thickness of the BP
film increases and the interactions between the VB and CB become smaller.
Thus, while the critical field $E_{c}$ may become smaller with thicker films
due to a decrease in the band gap with increasing thickness, the Dirac cone
phenomenon will be absent beyond a critical thickness even with gap closure.
We find that within DFT, when the film thickness increases to 12L, $E_{\rm
g}^{\Gamma-X}$ becomes negligible (equal to the precision of our
calculations)(see Supplementary Figure S6). Quantum confinement is therefore
important for the Dirac cone phenomenon.

Figure 6: (Color online) Effect of spin-orbit coupling on the band structure
of 4L black phosphorus, zoomed-in at the $\Lambda$ point. Black and red lines
denote the bandstructure calculated without and with spin-orbit coupling,
respectively, at $E_{\rm ext}$ = (a) 0.56, (b) 0.61 and (c) 0.66 V/Å.

The unique Dirac cone band structure implies also that for $U>U_{c}$, there is
no interaction between the VB and CB states along the $\Gamma$-Y direction,
but non-zero interactions along all other directions. This is a direct
consequence of the anisotropic interlayer interactions in black phosphorus.
Indeed, as shown in Figure 5(c), the $k$-dependent Hamiltonian overlap matrix
elements between the VB and CB states are highly anisotropic, going to zero
along $k_{x}=0$, i.e. along the $\Gamma$-Y direction, but non-zero elsewhere
in the Brillouin Zone.

From Figure 5(c), we can also see that for a fixed electric field strength,
the interaction between the VB and CB states increases as we move along
$\Gamma$-X away from the $\Gamma$ point. As the electric field strength is
increased, the VB and CB cross further away from $\Gamma$. This explains why
$E_{\rm g}^{\Gamma-X}$ initially increases as $E_{\rm ext}$ is increased [
Figure 4(c)]. However, a larger electric field strength also implies greater
localization of the VB and CB states, and therefore a smaller overall
interaction term. These two effects compete with each other, resulting in an
initial increase and then a subsequent decrease of $E_{\rm g}^{\Gamma-X}$ as
$E_{\rm ext}$ increases [ Figure 4(c)].

Thus, we have shown that the Dirac cone phenomenon arising for $E_{\rm
ext}>E_{c}$ stems from a combination of quantum confinement and anisotropy
effects in black phosphorus.

Figure 7: (Color online) Total energy per unit cell plotted against $d-d_{0}$
for different values of $E_{\mathrm{ext}}$ in V/Å. $d$ is the inter-layer
distance and $d_{0}$=3.59 Åis the equilibrium inter-layer distance. Inset
shows the variation of binding-energy (BE), calculated within the PBE+vdW
functional, as a function of the applied electric field, $E_{\mathrm{ext}}$,
for 2L black phosphorus.

Bandstructure with spin-orbit correction As already mentioned before, the band
gap remains direct at the $\Gamma$ point until the the gap closes completely.
As a consequence, the band gap does not change significantly by the inclusion
of spin-orbit coupling (SOC). However, the role of SOC becomes preeminent when
$E_{\rm ext}>E_{c}$. Interestingly, Figure 6(a) shows that the Dirac cone is
gapped out in the band structure by including the SOC, which demonstrates that
the system is likely to be a quantum spin Hall insulator. The prediction can
be explicitly confirmed by calculating the $Z_{2}$ topological invariant of
the system TI-zahid , as shown in Ref Zunger . For 4L-BP, the maximum SOC-
induced gap-opening at the Dirac point is as large as 3 meV at $E_{\rm ext}$=
0.56 V/Å. When $E_{\rm ext}$ is further increased, the SOC-induced gap
decreases (to 2meV at $E_{\rm ext}$= 0.61 V/Å) and subsequently, the Dirac
semimetal phase appears at $E_{\rm ext}$=0.66 V/Å along the $\Gamma$-Y
direction (Figure 6). Thus, the external electric field can be used to tune
the opening and closure of the gap along $\Gamma$-Y when SOC is included. It
is interesting to mention that the band inversion is controlled by
electrically driven band overlap between the VB and CB, unlike other
topological insulators, for example Bi2Se3, where band inversion is caused by
spin-orbit coupling.

Robustness of results against functional Our results above are robust against
the exact choice of functional. Specifically, with the PBE functional, the
band gap also decreases monotonically with thickness (see supplementary figure
S7), but the rate of reduction of the PBE-bandgap is different from that of
the vdW-corrected one. However, the qualitative physics, i.e. both the
thickness-dependent bandgap, and the external field induced bandgap-closure
and Dirac cone, consistently appear regardless of the functional. Notably, the
strength of critical field required for the insulator-to-metal transition of
bilayer BP is almost the same for both PBE and PBE+vdW functionals. Moreover,
although many-electron $GW$ calculations are in principle required for
accurate calculation of the band gap, the PBE and vdW corrected functionals
underestimate and overestimate the bulk band gap respectively, and in both
cases, the same physics emerges.

Figure 8: (Color online) Tuning of optical properties by thickness and applied
$E_{\mathrm{ext}}$. (a) The optical conductivity of 2L (black), 4L (red), 6L
(green), and 8L (blue) BP, plotted as a function of frequency. Solid and
dashed lines denote two different components of optical conductivity,
$\sigma_{xx}$ and $\sigma_{yy}$, respectively. (b) $\sigma_{xx}$ (solid lines)
and $\sigma_{yy}$ (dashed lines) of 4L-BP are plotted against the frequency
for two different applied fields, $E_{\mathrm{ext}}$ = 0.0 (black) and 0.21
(red) V/Å.

Binding energy We also calculate the total energy as a function of interlayer
distance in bilayer BP for different $E_{\rm ext}$, as shown in Figure 7. We
find that the binding energy reduces by increasing the strength of $E_{\rm
ext}$ and eventually vanishes for a field of $E_{\rm ext}\sim$ 1.6 V/Å, at
which point the layers becomes unbounded. This field is much larger than the
critical fields at which the Dirac cone phenomenon was predicted. Notably the
equilibrium distance between the layers is almost independent of $E_{\rm
ext}$. This result suggests that an external field can trigger the exfoliation
process of black phosphorus as an alternative to mechanical Li  and chemical
exfoliation Coleman . Note that few-layer graphene has been successfully
exfoliated in the presence of $E_{\rm ext}$ Liang .

Optical Properties Given the predicted changes in band gap as a function of
thickness and electric field strength, we also consider a related quantity -
the optical conductivity of the system. Figure 8(a) shows both the $xx$ and
$yy$-components of the optical conductivity, $\sigma_{xx}$ and $\sigma_{yy}$,
as a function of frequency for 2L, 4L, 6L and 8L structures. We find that the
onset of $\sigma_{xx}$ and $\sigma_{yy}$ depend on the dipole transitions at a
point near R, and at the $\Gamma$ point, respectively. As the bandgap reduces
with increasing thickness, the onset of $\sigma_{yy}$ also decreases, while
that of $\sigma_{xx}$ decreases slightly. By applying an external field to
4L-BP, we find that the electric-field induced reduction in direct bandgap
decreases the onset of $\sigma_{yy}$ whereas that of $\sigma_{xx}$ remains
almost unchanged [ Figure 8(b)].

## CONCLUSIONS

In summary, we have shown using first principles calculations that the
application of an external electric field to few-layer BP results in a
reduction in the direct band gap of the material, eventually resulting in an
insulator-to-metal transition at a critical field. Above this critical field,
a Dirac cone emerges at $\Lambda$ along the $\Gamma-Y$ direction, provided
that the BP film is sufficiently thin. We quantify the robustness of the Dirac
cone by the energy gap in the $\Gamma-X$ direction, $E_{\rm g}^{\Gamma-X}$.
Using a tight-binding model, we show that quantum confinement and structural
anisotropy are critical for the emergent Dirac cone phenomenon. These
conclusions further explain the subtle dependence of $E_{\rm g}^{\Gamma-X}$ on
the electric field strength and film thickness, providing us with a handle to
tune the resulting electronic structure. The electric field strength can be
used to tune the position of $\Lambda$ and the Dirac-Fermi velocities, the
latter being similar in magnitude to those in graphene. We also show that by
including spin-orbit coupling, the Dirac point can be gapped out at certain
electric fields, but increasing the electric field strength further can change
the material from a topological insulator to a Dirac semi-metal. We have thus
shown that the electronic band structure of few-layer black phosphorus can be
tuned in very interesting ways by an applied electric field. These results
will motivate the search for other direct band gap materials where quantum
confinement and structural anisotropy may lead to tunable emergent Dirac cone
physics.

## METHODS

Electronic structure calculations are performed by using density functional
theory within the generalized gradient approximation (GGA) of the exchange and
correlation potential Perdew-Burke-Ernzerhof parametrization (PBE) PBE  as
implemented in the Quantum Espresso QE  package. To treat van der Waals (vdW)
interactions, we employ the vdW+DF approach for the exchange-correlation
functional vdW . The wavefunction and the charge density are expanded using
energy cutoffs of 50 and 300 Ry, respectively. Brillouin zone sampling is done
by using a (10$\times$14$\times$1) Monkhorst-Pack $k$-grid. Periodic boundary
conditions have been included and a vacuum layer of at least 20 Å is included
to suppress the interaction between the periodic images. The conjugate
gradient algorithm is used to obtain optimized geometries, where all the atoms
in the unit cell are allowed to relax until the forces on each atom are less
than 0.002 eV/Å, and the coordinates are fixed in the calculations using the
finite external electric field. We have checked that for 4L-BP, the results
are qualitatively similar when the atoms are relaxed in the presence of the
applied electric field (see Supplementary Figure S8).

Simulations with an external electric field are performed using a periodic
sawtooth-type potential perpendicular to the layers swatooth . In order to
remove the effect of electrostatic interactions between the periodic
supercells, a dipole correction scheme suggested by L. Bengtsson dipole , is
included.

The real part of the optical conductivity is computed, using the following
formula: $\Re(\sigma)=\omega\times\Im(\epsilon)$. Here the imaginary part of
the dielectric tensor, $\Im(\epsilon_{\alpha\beta}(\omega))$, can be expressed
as

$\displaystyle\Im(\epsilon_{\alpha\beta}(\omega))=\frac{4\pi e^{2}}{\Omega
N_{k}m^{2}}\sum_{v,c}\int_{BZ}\frac{\hat{\textbf{M}}_{\alpha\beta}\times
f(E_{\textbf{k},v})}{\nabla_{k}(E_{v}-E_{c})}...$
$\displaystyle...\bigg{[}\delta(E_{\textbf{k},c}-E_{\textbf{k},v}+\hbar\omega)+\delta(E_{\textbf{k},c}-E_{\textbf{k},v}-\hbar\omega)\bigg{]},$
(6)

where $e~{}{\rm and}~{}m$ are the electronic charge and mass, respectively,
$\Omega$ is the lattice volume, $N_{k}$ is the total number of k-points, the
indices $c$ and $v$ denote the conduction and the valance band, respectively,
and $f$ is the Fermi-Dirac distribution. $\hat{\textbf{M}}_{\alpha\beta}$ is
defined as follows: $\hat{\textbf{M}}_{\alpha\beta}=\langle
u_{\textbf{k},c}|\hat{\textbf{p}}_{\alpha}|u_{\textbf{k},v}\rangle\langle
u_{\textbf{k},v}|\hat{\textbf{p}}_{\beta}^{\dagger}|u_{\textbf{k},c}\rangle$,
where $|u_{\bf k}\rangle$ and $\hat{\bf P}$ are the single particle Bloch
function and the dipole moment operator, respectively. An uniform k mesh of
36$\times$26$\times$1 is used for self consistent steps in order to calculate
the optical conductivity. Note that our calculation does not include many-body
effects which are required to obtain an accurate description of the optical
conductivity.

## ACKNOWLEDGMENTS

We gratefully acknowledge support from the National Research Foundation,
Singapore, for funding under the NRF Fellowship (NRF-NRFF2013-07) and under
the NRF medium-sized centre programme. We thank the Graphene Research Center
for the computational resources provided.

## Author contribution

KD performed the calculations. KD and SYQ planned the study. SYQ provided
supervision. KD and SYQ wrote the manuscript.

## Additional Information

The authors declare no competing financial interests.

## References

  * (1) Wang, Q.H., Kalantar-Zadeh, K., Kis, A., Coleman, J.N. & Strano, M.S. Electronics and optoelectronics of two-dimensional transition metal dichalcogenides. Nature Nanotechnology 7, 699-712 (2012).
  * (2) Splendiani, A., Sun, L., Zhang, Y., Li, T., Kim, J., Chim, C.-Y., Galli, G. & Wang, F. Emerging Photoluminescence in Monolayer MoS2. Nano Lett. 10, 1271-1275 (2010).
  * (3) Li, L., Yu,Y., Ye,G.J., Ge,Q., Ou, X., Wu,H., Feng, D., Chen, X.H. & Zhang, Y. Black phosphorus field-effect transistors. Nature Nanotechnology 9, 372-377 (2014).
  * (4) Hanlon, D. et al. Liquid exfoliation of solvent-stabilised black phosphorus: applications beyond electronics. arXiv:1501.01881 [cond-mat.mes-hall]
  * (5) Das , S., Zhang, W., Demarteau, M., Hoffmann, A., Dubey, M. & Roelofs, A. Tunable Transport Gap in Phosphorene. Nano Lett. 14, 5733-5739 (2014).
  * (6) Xia, F. et al. Rediscovering black phosphorus as an anisotropic layered material for optoelectronics and electronics. Nature Communication 5, 4458 (2014).
  * (7) Khoo, K. H., Mazzoni, M. S. C. & Louie, S. G., Tuning the electronic properties of boron nitride nanotubes with transverse electric fields: A giant dc Stark effect, Phys. Rev. B 69, 201401 (2004).
  * (8) Ramasubramaniam, A., Naveh, D., & Towe, E. Tunable band gaps in bilayer transition-metal dichalcogenides. Phys. Rev. B 84, 205325 (2011).
  * (9) Dolui, K., Pemmaraju, C.D. & Sanvito, S. Electric field effects on armchair MoS2 nanoribbons. Acs Nano 6, 4823-4834 (2012).
  * (10) Park, C.-H., & Louie, S. G.Energy Gaps and Stark Effect in Boron Nitride Nanoribbons. Nano Lett. 8, 2200-2203 (2008).
  * (11) Son, Y.-W., Cohen, M. L. & Louie, S. G. Half-Metallic Graphene Nanoribbons. Nature (London) 444, 347-349 (2006).
  * (12) Gava, P., Lazzeri, M., Saitta, A. M. & Mauri, F. Ab initio study of gap opening and screening effects in gated bilayer graphene. Phys. Rev. B 79, 165431 (2009).
  * (13) Perdew, J. P., Burke, K. & Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 77, 3865-3868 (1996).
  * (14) Giannozzi, P; et al., QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 21, 395502 (2009).
  * (15) Dion, M.; Rydberg, H., Schroder, E., Langreth, D. C. & Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 92, 246401 (2004).
  * (16) Cartz L., Srinivasa, S.R., Riedner, R. J., Jorgensen, J. D. & Worlton, T. G. Effect of pressure on bonding in black phosphorus. J. Chem. Phys. 71, 1718-1721 (1979).
  * (17) Warschauer, D. Electrical and optical properties of crystalline black phosphorus. J. Appl. Phys. 34, 1853-1860 (1963).
  * (18) Qiao,J., Kong,X., Hu, Z.-X., Yang, F.& Ji W. High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus. Nature Communications 5, 4475 (2014).
  * (19) Zheng, F., Liu, Z., Wu, J., Duan, W. & Gu, B.-L. Scaling law of the giant Stark effect in boron nitride nanoribbons and nanotubes. Phys. Rev. B 78, 085423 (2008).
  * (20) Malko, D., Neiss, C., Viñes, F., & Görling, A. Competition for Graphene: Graphynes with Direction-Dependent Dirac Cones. Phys. Rev. Lett. 108, 086804 (2012).
  * (21) Rudenko, A.N. & Katsnelson, M.I. Quasiparticle band structure and tight bonding model for single- and bilayer black phosphorus. Phys. Rev. B 89, 201408 (2014).
  * (22) Hasan, M. Z. & Kane, C. L. Colloquium : Topological insulators. Rev. Mod. Phys. 82, 3045-3067 (2010).
  * (23) Liu, Q., Zhang, X., Abdalla,L. B., Fazzio, A. & Alex Zunger, A. Switching a Normal Insulator into a Topological Insulator via Electric Field with Application to Phosphorene. Nano Lett. 15, 1222-1228 (2015).
  * (24) Liang, X., Chang, A.S.P., Zhang, Y., harteneck, B.D., Choo, H., Olynick, D.L.& Cabrini, S. Electrostatic Force Assisted Exfoliation od Prepatterrned Few-Layer Graphene into Device Sites. Nano. Lett. 9, 467-472 (2009).
  * (25) Hart, C. F. Exact k-q Solution for a Bloch Electron in a Constant Electric Field. Phys. Rev. B 38, 2158-2161 (1988).
  * (26) Bengtsson, L. Dipole correction for surface supercell calculations. Phys. Rev. B 59, 12301 (1999).

Supplementary Information for ”Quantum-confinement and Structural Anisotropy
result in Electrically-Tunable Dirac Cone in Few-layer Black Phosphorous”

Figure 9: Supporting Figure S1: The band-decomposed projected density of
states (PDOS) on different orbitals for 4 layer black phosphorus. The
densities of states are projected on 3$s$, 3$p_{x}$, 3$p_{y}$ and 3$p_{z}$
orbitals. Upper, middle and lower panels show the PDOS at $E_{\rm ext}$ =
0.00, 0.21, and 0.72 V/Å, respectively. Note that the critical field is
$E_{c}$ = 0.48 V/Å. At $E_{\rm ext}$ = 0, the PDOS is summed over all the
atoms in the unit-cell. For finite fields, for visualization purposes, we use
the negative colour scale for PDOS that comes from the two layers kept at
positive potential, and the positive colour scale for PDOS coming from the
other two layers kept at negative potential, as indicated in Figure 1 in the
main text. Figure 10: Supporting Figure S2: (a) Electric field induced charge
density distribution across the layers, $\Delta\rho=\rho(E_{\rm
ext})-\rho(E_{\rm ext}=0)$, for 4L-BP and different values of $E_{\rm ext}$.
$\rho(E_{\rm ext})$ is the charge density for an external field, $E_{\rm
ext}$. $\rho(E_{\rm ext})$ is averaged over the transverse direction
($xy$-plane) and plotted along the out-of-plane direction ($z$-axis). The
dashed vertical lines denote the positions of each atomic plane in 4L-BP. (b)
Effective electrostatic potential across the layers, $\Delta U_{\rm
eff}=U(E_{\rm ext})-U(E_{\rm ext}=0)$, for 4L-BP and different values of
$E_{\rm ext}$.

Figure 11: Supporting Figure S3: (a) The polarization (dipole moment per unit
area) , $P$, is plotted against the applied external field strength for 2L
(black), 3L (red), 4L (green), 5L (blue), 6L (yellow), 7L (brown) and 8L
(grey)-BP. (b) ($E_{g}/E_{c}$ is plotted as a function of the number of
layers. Figure 12: Supporting Figure S4: (a) Spin-polarized bandstructure and
(b) density of states of 4L-BP. The blue dashed line denotes the Fermi level.
Red and black lines denote the spin-up and spin-down component of the DOS,
respectively. Figure 13: Supporting Figure S5: The local density of states
(LDOS) of highest occupied states (HOS) left) and lowest unoccupied states
(LUS) (right) at the different k-points around the Dirac point, $\Lambda$, as
indicated in the bandstucture (top). The figure shows that band inversion
occurs at the $\Lambda$ point. Figure 14: Supporting Figure S6: The
bandstructure of 12L-BP above the critical field, i.e., $E_{\rm ext}$= 0.26
V/Å. Figure 15: Supporting Figure S7: Comparison of electronic structure for
different functionals. (a) The variation of the band-gaps calculated within
PBE(black line) and vdW-corrected(red line) exchange correlation functionals,
plotted as a function of the number of layers. Note that the PBE band gap of
bulk BP is 0.12 eV, while the vdW-corrected band gap is 0.34 eV. The dashed
line denotes the value of the experimental band gap (0.34 eV) of bulk BP. (b)
The variation of the band gap of bilayer BP, calculated within PBE(black line)
and vdW-corrected(red line) functionals, plotted against the applied field
strength. Figure 16: Supporting Figure S8: The bandgap of 4L-BP as a
functional of $E_{\rm ext}$ using the unrelaxed (black) and relaxed (red)
coordinates under an applied field $E_{\rm ext}$.