# Magnetoexcitons in phosphorene monolayers, bilayers, and van der Waals
heterostructures

Roman Ya. Kezerashvili1,2 and Anastasia Spiridonova1,2 E-mail contact: 1New
York City College of Technology, The City University of New York, Brooklyn, NY
11201, USA
2The Graduate School and University Center, The City University of New York,
New York, NY 10016, USA

###### Abstract

We study direct and indirect excitons in Rydberg states in phosphorene
monolayers, bilayer and van der Waals (vdW) heterostructure in an external
magnetic field, applied perpendicular to the monolayer or heterostructure
within the framework of the effective mass approximation. Binding energies of
magnetoexcitons are calculated by a numerical integration of the Schrödinger
equation using the Rytova-Keldysh potential for direct magnetoexcitons and
both the Rytova-Keldysh and Coulomb potentials for indirect one. The latter
allows to understand the role of screening in phosphorene. We report the
magnetic field energy contribution to the binding energies and diamagnetic
coefficients (DMCs) for magnetoexcitons that strongly depend on the effective
mass of electron and hole and their anisotropy and can be tuned by the
external magnetic field. We demonstrate theoretically that the vdW phosphorene
heterostructure is a novel category of 2D semiconductor offering a tunability
of the binding energies of magnetoexcitons by mean of external magnetic field
and control the binding energies and DMCs by the number of hBN layers
separated two phosphorene sheets. Such tunability is potentially useful for
the devices design.

##  I Introduction

The most recent addition to the growing family of 2D material is phosphorene
that is the monolayer of black phosphorus (BP)_._ The BP is composed of
individual phosphorene layers. The blossoming interest in phosphorene derives
partially from its direct gap that is retained in both monolayer and bulk
structure, high mobility, and high on-off current ratio in field-effective
transistors [1, 2, 3, 4]. In contrast, these characteristics simultaneously
are not present in graphine, transitional metal dichalcogenides (TMDCs): WSe2,
WS2, MoSe2, MoS2, or Xenes: silicene, germanene, stanene, making phosphorene a
promising material for electronic and optical applications. Black phosphorus
is the most stable of phosphorus allotropes. It was first synthesized in 1914
by Bridgman [5]. While bulk black phosphorus has interesting properties and
was extensively studied in twentieth century [6, 7, 8, 9, 10, 11], it was not
used in the design of electronic devices.

However, everything has changed after the discovery of graphene in 2004 [12]
which paved the way for two-dimensional materials. Initially, the research was
focused on graphene [13, 14] and TMDCs [15, 16, 17, 18, 19, 20, 21, 22, 23]
then silicene and other group VI elements have followed [24, 25, 26, 27].
There are some drawbacks in the above materials. Graphene has high mobility
and on-off ratio but a lack of a gap impedes its use in the design of
electronic devises [4, 2]. TMDCs have a gap [28, 29, 30], but its carrier
mobility is orders of magnitudes lower than in graphene. The phosphorene that
has been synthesized in 2014 has advantages over graphene, TMDCs, and Xenes.
Its most remarkable properties include thickness-dependent band gap, strong
in-plane anisotropy, and high carrier mobility.

BP has a structure where monolayer appears to be composed of two distinct
planes made of puckered honeycomb structure [31, 32, 33]. This results in an
anisotropic electronic structure. Due to phosphorene unique topological
structure and differences between the armchair (AC) and zigzag (ZZ)
directions, it displays strong in-plane anisotropy. Many properties of
phosphorene in these two principal directions are drastically different. The
anisotropic structure is strongly reflected in effective masses of charge
carriers and leads to anisotropic effective masses of electrons and holes
along AC and ZZ ($x$ and $y$) directions. Along the AC direction the effective
electron and hole masses are smaller than along the ZZ direction [9, 31, 32,
33]. The distinct features of BP is that the direct gap located at the
$\Gamma$ point is preserved in the monolayer and bulk structure. This material
is a direct band gap semiconductor with a strongly anisotropic dispersion in
the vicinity of the gap. It is worth noting that the unusual structure of
phosphorene sets it aside from graphene and other widely studied 2D
semiconductors.

Since the synthesis of phosphorene, it has been extensively studied. For
example, optical and thermal properties have been considered in Refs.[34, 35,
33, 36, 37], in plane electric field has been studied in Refs. [38, 39, 40,
41], the effects of the strain on different properties of phosphorene have
been addressed in Refs. [42, 43, 44], the Landau levels have been reported in
Refs. [45, 46, 47, 48, 49, 50].

Phosphorene hosts tightly bound excitons [51, 52, 31, 53, 54]. Similarly to
other monolayer semiconductors, reduced dimensionality and reduced screening
of the Coulomb attraction lead to a high exciton binding energy in phosphorene
_._ Moreover,__ the exciton binding energy in phosphorene is larger than the
one in other 2D materials. Excitons in semiconductors in the presence of the
external magnetic field have been studied for the past sixty years. Elliot and
Loudon [55] and Hasegawa and Howard [56] developed the theory of the Mott
exciton in the strong magnetic field. Authors in Refs. [57, 58, 59] addressed
Mott excitons. Excitons in TMDCs monolayers, bilayer, and the double-layer
heterostructure in the presence of the external magnetic field have been
extensively studied. The diamagnetic shifts have been reported in Refs. [60,
61, 62, 63, 64, 65, 66, 67, 23], and Zeeman shifts have been reported in Refs
[68, 69, 70, 71, 72, 73]. We cite these works, but the recent literature on
the subject is not limited by them. However, there is a lack of similar
research on excitons in phosphorene. This motivates us to study the effect of
the external magnetic field on the binding energies of Rydberg states of
magnetoexcitons in monolayer, bilayer, and the double-layer heterostructure
composed of phosphorene and to calculate diamagnetic coefficients (DMCs).

In this paper, we study the dependence of the magnetoexciton binding energy of
Rydberg states: 1$s$, 2$s$, 3$s$, and 4$s$, on the external magnetic field
perpendicular to a monolayer or heterostructure that is varied between 0 and
60 T and report the diamagnetic coefficients (i) for the direct
magnetoexcitons in a freestanding (FS) and encapsulated phosphorene
monolayers, (ii) for the indirect magnetoexciton in bilayer composed of two
phosphorene monolayers, (iii) for the indirect magnetoexcitons in
heterostructure formed by two phosphorene monolayers and separated by $N$ hBN
monolayers. The latter van der Waals heterostructure is denoted as vdW. In our
approach, we numerically solve the Schrödinger equation for the magnetoexciton
and obtain eigenfunctions and eigenvalues. Then we obtain the energy
contribution from the magnetic field to the binding energies and use it to
calculate diamagnetic coefficients. For the direct exciton, we solve the
Schrödinger equation with the Rytova-Keldysh potential [74, 75], and for the
indirect exciton, we solve the Schrödinger equation with the Rytova-Keldysh
and the Coulomb potentials to understand the role of the screening in
phosphorene.

The remainder of this paper is organized in the following way. In section II
the theoretical model for the description of an electron-hole system in the
external magnetic field with the charge carriers effective mass anisotropy is
presented. Results of calculations of binding energies of Rydberg states and
DMCs for direct magnetoexcitons in a freestanding and encapsulated phosphorene
and indirect magnetoexcitons in FS bilayer and vdW phosphorene
heterostructures are presented in Sec. III. Finally, conclusions follow in
Sec. IV.

##  II Theoretical Model

It is known that electrostatically-bound electrons and holes in the external
magnetic field form magnetoexcitons. In this section, following Refs. [58, 67,
27, 23], we introduce briefly the theoretical model for the description of the
Mott-Wannier magnetoexciton in phosphorene. We consider the energy
contribution from the external magnetic field to the Rydberg states binding
energies of magnetoexcitons and diamagnetic coefficients (DMCs). In the system
under consideration, excitons are confined in a 2D freestanding and
encapsulated phosphorene monolayer, FS bilayer phosphorene and van der Waals
heterostructure, where $N$ layers of hBN monolayers separate two phosphorene
monolayers. In the latter two cases equal number of electrons and holes are
located in parallel phosphorene monolayers at a distance $D$ away. The
corresponding schematic illustrations of these systems are shown in Fig. 1.

It is worth mentioning that we are considering the monolayers of hBN as an
insulator. Phosphorene encapsulated between hBN layers is robust to oxidation
and exhibits high mobilities. Moreover, hBN also has a high dielectric
constant, resulting in a strong damping of the electrostatic repulsion by
charged impurities which are responsible for the decrease of carriers
mobility.

Figure 1: (Color online) Schematic illustration of magnetoexcitons in
phosphorene monolayers and heterostructures. $(a)$ A direct magnetoexciton in
a freestanding phosphorene monolayer. $(b)$ A direct magnetoexciton in an
encapsulated phosphorene monolayer. $(c)$ An indirect magnetoexciton in a
freestanding bilayer phosphorene heterostructure. $(d)$ An indirect
magnetoexciton in phosphorene van der Waals heterostructure.

Let us introduce the coordinate vectors of the electron and hole for the Mott-
Wannier exciton in the phosphorene layer. The following in-plane coordinates
$\mathbf{r}_{1}=(x_{1},y_{1})$ and $\mathbf{r}_{2}=(x_{2},y_{2})$ for an
electron and hole, respectively, will be used in our description. We assume
that at low momentum $\mathbf{p}=(p_{x},p_{y})$, i.e., near the $\Gamma$
point, the single electron and hole energy spectrum
$\varepsilon_{l}^{(0)}(\mathbf{p})$ is given by

$\displaystyle\varepsilon_{l}^{(0)}(\mathbf{p})=\frac{p_{x}^{2}}{2m_{x}^{l}}+\frac{p_{y}^{2}}{2m_{y}^{l}},\
\ l=e,\ h,$ (1)

where $m_{x}^{l}$ and $m_{y}^{l}$ are the electron and hole effective masses
along the $x$ and $y$ directions, respectively. We assume that $OX$ and $OY$
axes correspond to the armchair and zigzag directions in a phosphorene
monolayer, respectively.

The anisotropic nature of the 2D phosphorene atomic semiconductor, in contrary
to other 2D isotropic materials such as graphene and TMDC semiconductors,
breaks the central symmetry and requires for the description of excitons the
use of the Cartesian coordinates. The asymmetry of the electron and hole
dispersion in phosphorene is reflected in the Hamiltonian for the Mott-Wannier
magnetoexciton, and within the framework of the effective mass approximation
the Hamiltonian for an interacting electron-hole pair in phosphorene in the
external magnetic field reads ($\hbar=c=1$)

$\displaystyle\hat{H}=$
$\displaystyle\frac{1}{2m_{e}^{x}}\Big{(}i\nabla_{x_{e}}-eA_{x}(r_{e})\Big{)}^{2}+\frac{1}{2m_{e}^{y}}\Big{(}i\nabla_{y_{e}}-eA_{y}(r_{e})\Big{)}^{2}+\frac{1}{2m_{h}^{x}}\Big{(}i\nabla_{x_{h}}+eA_{x}(r_{h})\Big{)}^{2}+\frac{1}{2m_{h}^{y}}\Big{(}i\nabla_{y_{h}}+eA_{y}(r_{h})\Big{)}^{2}$
(2) $\displaystyle+V\Big{(}|\mathbf{r}_{e}-\mathbf{r}_{h}|\Big{)},$

where the $m_{i}^{j},~{}j=x,y,~{}i=e,h$ correspond to the effective mass of
the electron or hole in the $x$ or $y$ direction, respectively, and
$V\left(|\mathbf{r}_{e}-\mathbf{r_{h}}|\right)$ describes the electrostatic
interaction between the electron and hole. When the electron and hole are
located in 2D plane, we use the Rytova-Keldysh (RK) potential [74, 75] that is
widely used for the description of charge carriers interaction in 2D
materials. The RK potential is a central potential, and the interaction
between the electron and hole for direct excitons in a monolayer has the form
[74, 75]:

$V_{RK}(r)=-\frac{\pi
ke^{2}}{2\kappa\rho_{0}}\left[H_{0}\left(\frac{r}{\rho_{0}}\right)-Y_{0}\left(\frac{r}{\rho_{0}}\right)\right],$
(3)

where $r=r_{e}-r_{h}$ is the relative coordinate between the electron and
hole. In Eq. (3), $e$ is the charge of the electron,
$\kappa=(\epsilon_{1}+\epsilon_{2})/2$ describes the surrounding dielectric
environment, $\epsilon_{1}$ and $\epsilon_{2}$ are the dielectric constants
below and above the monolayer, $H_{0}$ and $Y_{0}$ are the Struve and Bessel
functions of the second kind, respectively, and $\rho_{0}$ is the screening
length. The screening length $\rho_{0}$ can be written as
$\rho_{0}=2\pi\chi_{2D}/\kappa$ [76], where $\chi_{2D}$ is the 2D
polarizability, which in turn is given in [75]. At the long-range distances
$V_{RK}$ retains $1/r$ behavior as the Coulomb potential, for smaller
distances the potential well is logarithmic.

Following the standard procedure [77] for the separation of the relative
motion of the electron-hole pair from their center-of-mass motion one
introduces variables for the center-of-mass of an electron-hole pair
$\mathbf{R}=(X,Y)$ and the relative motion of an electron and a hole
$\mathbf{r}=(x,y)$, as
$X=(m_{x}^{e}x_{1}+m_{x}^{h}x_{2})/(m_{x}^{e}+m_{x}^{h})$,
$Y=(m_{x}^{e}y_{1}+m_{x}^{h}y_{2})/(m_{x}^{e}+m_{x}^{h})$, $x=x_{1}-x_{2}$,
$y=y_{1}-y_{2}$, $r^{2}=x^{2}+y^{2}$. The Schrödinger equation with
Hamiltonian (2) has the form:
$\hat{H}\Psi(\mathbf{r}_{1},\mathbf{r}_{2})=\mathcal{E}\Psi(\mathbf{r}_{1},\mathbf{r}_{2})$,
where $\Psi(\mathbf{r}_{1},\mathbf{r}_{2})$ and $\mathcal{E}$ are the
eigenfunction and eigenenergy, respectively. One can write
$\Psi(\mathbf{r}_{1},\mathbf{r}_{2})$ in the form
$\Psi(\mathbf{r}_{1},\mathbf{r}_{2})=\Psi(\mathbf{R},\mathbf{r})=e^{i\mathbf{P}\cdot\mathbf{R}/\hbar}\Phi(x,y)$,
where $\mathbf{P}=(P_{x},P_{y})$ is the momentum of the center-of-mass of the
electron-hole pair and $\Phi(x,y)$ is the wave function of the relative motion
of the electron-hole pair. Following Refs. [78, 67, 27, 23], after lengthy
calculations one obtains the equation that describes the Mott–Wannier exciton
in Rydberg optical states in the external magnetic field perpendicular to the
phosphorene monolayer. Finally, the equation for the relative motion of the
electron and hole in the phosphorene monolayer with zero center-of-mass
momentum reads

$\left[-\frac{1}{2\mu_{x}}\frac{\partial^{2}}{\partial
x^{2}}-\frac{1}{2\mu_{y}}\frac{\partial^{2}}{\partial
y^{2}}+\frac{e^{2}}{8\mu_{x}}B^{2}x^{2}+\frac{e^{2}}{8\mu_{y}}B^{2}y^{2}+V(x,y)\right]\Phi(x,y)=E\Phi(x,y).$
(4)

where $\mu_{x}=\frac{m_{x}^{e}m_{x}^{h}}{m_{x}^{e}+m_{x}^{h}}$ and
$\mu_{y}=\frac{m_{y}^{e}m_{y}^{h}}{m_{y}^{e}+m_{y}^{h}}$ are the reduced
masses, related to the relative motion of an electron-hole pair in the $x$ and
$y$ directions, respectively. In Eq. (4) the anisotropy is present in the
kinetic and magnetic terms, while the potential term has isotropic form.

Following Refs. [79, 80] one can obtain from (4) the expectation value for the
ground state energy. For the case of $V_{RK}(x,y)$ potential we have
$E=\left\langle-\frac{1}{2\mu_{x}}\frac{\partial^{2}}{\partial
x^{2}}\right\rangle+\left\langle-\frac{1}{2\mu_{y}}\frac{\partial^{2}}{\partial
y^{2}}\right\rangle+\left\langle\frac{e^{2}}{8\mu_{x}}B^{2}x^{2}\right\rangle+\left\langle\frac{e^{2}}{8\mu_{y}}B^{2}y^{2}\right\rangle+\left\langle
V_{RK}(x,y)\right\rangle$. The later expression could be viewed as the sum of
the average values of the operators of kinetic and potential energies in 2D
space obtained by using the corresponding eigenfunction $\Phi(x,y)$ of the
exciton. Phosphorene exhibits charge carriers effective mass anisotropy with
lighter effective masses along the armchair direction and heavier effective
masses along the zigzag direction that leads to the reduced mass anisotropy:__
$\mu_{y}>\mu_{x}$. For the set of masses presented in Table 1 the ratio
$\mu_{y}/\mu_{x}$ varies from $\sim$ 7 to 16. Therefore, from the latter
expression one can conclude that contributions from the kinetic and magnetic
energies to the total ground state energy can be by the order of magnitude
different for AC and ZZ directions. In particular, the kinetic energy along
the AC direction is much larger than that along the ZZ direction. Thus, the
anisotropic nature of the 2D phosphorene atomic semiconductors, in contrast to
other 2D materials such as graphene, Xenes, and TMDC semiconductors, allows
excitons to be confined in a quasi-one-dimensional space predicted in theory
[81, 1], leading to remarkable phenomena arising from the reduced
dimensionality and screening. In Fig. 2 the total potential
$W(x,y)=V_{RK}(x,y)+\frac{e^{2}}{8\mu_{x}}B^{2}x^{2}+\frac{e^{2}}{8\mu_{y}}B^{2}y^{2}$
that acts on the electron-hole system in phosphorene is presented. One can
observe the anisotropic structure of this potential and its asymmetry with
respect to the AC and ZZ directions.

Let us now to consider the indirect magnetoexcitons formed by electrons and
holes located in two different phosphorene monolayers in the bilayer or vdW
heterostructure. In the latter case the phosphorene layers are separated by
$N$ layers of hBN monolayers. Such magnetoexcitons have a longer lifetime than
the direct excitons due to longer recombination time. Equation (4) still
describes the indirect exciton. However, for indirect excitons, the
expressions for the interaction between the electron and hole in Eq. (4) can
be written as:

$V_{RK}(\sqrt{\rho^{2}+D^{2}})=-\frac{\pi
ke^{2}}{2\kappa\rho_{0}}\left[H_{0}\left(\frac{\sqrt{\rho^{2}+D^{2}}}{\rho_{0}}\right)-Y_{0}\left(\frac{\sqrt{\rho^{2}+D^{2}}}{\rho_{0}}\right)\right],$
(5)

for the RK potential, and

$V_{C}\left(\sqrt{\rho^{2}+D^{2}}\right)=-\frac{ke^{2}}{\kappa\left(\sqrt{\rho^{2}+D^{2}}\right)}$
(6)

for the Coulomb potential and where $\rho^{2}=x^{2}+y^{2}$. Equations (5) and
(6) describe the interaction between the electron and hole that are located in
different parallel phosphorene monolayers separated by a distance
$D=h+Nl_{\text{hBN}}$, where $l_{\text{hBN}}=0.333$ nm is the thickness of the
hBN layer and $h$ is the phosphorene thickness given in Table 1. Therefore,
one can obtain the eigenfunctions and eigenenergies of magnetoexcitons by
solving Eq. (4) using the potential (3) for direct magnetoexcitons, or using
either potential (5) or (6) for indirect magnetoexcitons.

Figure 2: Total potential
$W(x,y)=V_{RK}(x,y)+\frac{e^{2}}{8\mu_{x}}B^{2}x^{2}+\frac{e^{2}}{8\mu_{y}}B^{2}y^{2}$
is plotted as a function of $x$ and $y$. $W(x,y)$ and $x$ and $y$ are given in
eV and nm, respectively. The potential is calculated at $B=30$ T with the
parameters of the reduced masses of the set 1. Table 1: Parameters for
phosphorene. Four sets of reduced masses for phosphorene that are used in
calculations as input parameters. The anisotropic reduced masses $\mu_{x}$ and
$\mu_{y}$ are obtained based on effective masses of the electron and hole
given in Ref. [82] for the set 1, in Ref. [83] for the set 2, in Ref. [84] for
the set 3, and in Ref.[1] for the set 4. $\mu_{x}$ and $\mu_{y}$ are in units
of the electron mass, $m_{0}$. $\chi_{2D}$ and $h$ are the polarizability and
thickness of phosphorene, respectively, in nm.

[b]  $\mu_{x}$ $\mu_{y}$ $\mu_{y}$/$\mu_{x}$ $\chi_{2D}$ (nm) $h$ (nm) set 1
0.06296 0.96774 15.37 0.41a 0.541b set 2 0.06667 0.88780 13.32 set 3 0.09122
0.6599 7.23 set 4 0.07969 0.9498 11.92

  * a

Reference [85].

  * b

Reference [86].

##  III Results of calculations and discussion

We report the dependence of the energy contribution from the external magnetic
field to the binding energies of magnetoexcitons in Rydberg states 1$s$, 2$s$,
3$s$, and 4$s$ in the FS (Fig. 1($a$)) and encapsulated by hBN (Fig. 1($b$))
phosphorene, the FS phosphorene bilayer (Fig. 1($c$)), and the vdW
heterostructure (Fig. 1($d$)) on the magnetic field. The diamagnetic
coefficients for the monolayers, bilayer, and vdW phosphorene heterostructure
are reported for the first time.

In calculations, we use effective masses for electron and hole found in
literature Refs. [82, 83, 84, 1], which are obtained by using the first
principle calculations. The lattice constants in Refs. [82, 83, 84, 1] do not
coincide with each other and different functionals for the correlation energy
and setting parameters for the hopping lead to some difference in results for
anisotropic masses. It is obvious that these can cause the difference in the
band curvatures and, therefore, effective masses along the armchair and the
zigzag directions. The latter motivate us to use in calculations of binding
energies of magnetoexcitons and DMCs the different sets of masses from Refs.
[82, 83, 84, 1] that allows to understand the dependence of the binding energy
of magnetoexcitons and DMCs on effective masses of electrons and holes and the
role of anisotropic masses. The corresponding reduced masses $\mu_{x}$ and
$\mu_{y}$ along the armchair and zigzag directions, respectively, are given in
Table 1. We use these masses as input parameters and refer to them as the set
1, the set 2, the set 3, and the set 4. Below, we report results for the sets
of masses 1 and 3 since they give an upper and lower bounds on the binding
energies and the energy contribution from the magnetic field to the binding
energies of the Rydberg states. Results from the sets of masses 2 and 4 fall
within results of the sets 1 and 3.

A numerical solution of the Schrödinger equation (4), using the aforementioned
interaction potentials and input parameters from Table 1, is performed using
the finite element method, which yields pairs of eigenenergies and
eigenfunctions which are solutions to the Schrödinger equation, corresponding
to the most-strongly-bound states. The method is based on using the finite
element method implemented in Wolfram Mathematica in the NDEigensystem
function. We modify the code successfully implemented in Ref. [36] in a way
that it explicitly contains $\frac{e^{2}}{8\mu_{x}}B^{2}x^{2}$ and
$\frac{e^{2}}{8\mu_{y}}B^{2}y^{2}$ terms. Our results for the binding energies
of direct excitons in FS and encapsulated hBN phosphorene monolayers and
indirect excitons in FS phosphrene bilayer for the masses sets 1 and 3 are
given in Table 2. One can see that the binding energy of excitons in FS and
encapsulated monolayers significantly depend on the anisotropic reduced
masses. Moreover, a much lower screening of excitons in the freestanding
monolayer makes their binding energies, for example, in the ground state over
three fold bigger than when the monolayer is encapsulated by hBN.

To check the validity of the code, we use the input parameters from respective
papers listed below and calculate the binding energies of direct and indirect
excitons. The code reproduces the theoretical binding energies of direct
excitons in phosporene and of indirect excitons in bilayer reported in Ref.
[3], where binding energies are calculated using Wannier effective-mass theory
in 2D space, within 8%. The binding energies for 1$s$ state obtained using a
semi-analytical perturbation theory approach [54] are reproduced within 4%.
The code reproduces binding energies for $n=0,...,4$ states for FS
phosphorene, phosphorene on SiO2 substrate, and encapsulated by hBN
phosphorene when energies obtained using numerical procedure where the
Schrödinger equation is solved via a finite element representation of the
exciton wave equation [87] within 1% and using $k\cdot p$ theory [88] within
4%.

Table 2: Binding energies of direct magnetoexcitons in FS and encapsulated hBN phosphorene monolayers and indirect magnetoexcitons in FS phosphorene bilayer for the sets of masses 1 and 3. For indirect magnetoexcitons in the bilayer phosphorene binding energies are calculated using $V_{RK}$ and $V_{C}$ potentials . Energies are given in meV. State | Monolayer | Bilayer
---|---|---
| Set 1 | Set 3 | Set 1 | Set 3
|  hBN |  FS |  hBN |  FS |  $V_{RK}$ |  $V_{C}$ |  $V_{RK}$ |  $V_{C}$
1$s$ |  187.94 |  718.71 |  199.66 |  746.06 |  388.88 |  827.17 |  395.54 |  854.03
2$s$ |  78.62 |  488.84 |  74.37 |  478.75 |  298.56 |  653.05 |  294.50 |  644.97
3$s$ |  53.83 |  394.04 |  50.35 |  377.72 |  248.52 |  525.71 |  239.53 |  501.07
4$s$ |  35.67 |  319.55 |  32.03 |  311.22 |  211.61 |  430.68 |  209.07 |  429.89

###  III.1 Contribution from the external magnetic field to binding energies
of magnetoexcitons in a monolayer

Figure 3: Dependencies of the energy contribution from the magnetic field to
the binding energies of magnetoexcitons in states 1$s$, 2$s$, 3$s$, and 4$s$
for FS ($a$) and encapsulated by hBN ($b$) phosphorene monolayers on the
squared magnetic field. The boundary dashed and solid curves correspond to the
sets of masses 1 and set 3, respectively. The contributions for the sets of
masses 2 and 4 fall within the shaded region. The inset in Fig. 3($b$) shows
the dependence of the energy contribution from the magnetic field on the ratio
$\frac{\mu_{y}}{\mu_{x}}$.

In Fig. 3 ($a$)-($b$) we report the energy contribution from the external
magnetic field to the binding energies of Rydberg states for direct
magnetoexcitons in FS and encapsulated by hBN phosphorene, respectively, as a
function of $B^{2}$. We report results for the sets of masses 1 and 3. The
dashed and solid boundary curves in Fig. 3 correspond to the results obtained
with the sets of masses 1 and 3, respectively. The results from the sets of
masses 2 and 4 fall in the shaded region.

Interestingly enough, in phosphorene, we observe the same tendencies for
magnetoexcitons as in TMDCs [67, 23] and Xenes [27] monolayers. The binding
energies and the energy contribution from the magnetic field to the binding
energies of Rydberg states of magnetoexcitons strongly depend on the exciton
reduced mass, even though phosphorene has anisotropic mass along $x$ and $y$
directions. The reduced mass of the exciton along the AC direction, depending
on the set, is between seven and fifteen times smaller than the reduced mass
along the ZZ direction. Since according to Eq. (4), the Schrödinger equation
contains terms that are proportional to $\frac{1}{\mu_{x}}$ and
$\frac{1}{\mu_{y}}$, the terms with significantly smaller $\mu_{x}$ dominate
in Eq. (4). This is consistent with the anisotropic behavior of the potential
$W(x,y)$ shown in Fig. 2 that exhibits the anisotropy with respect to the
armchair and zigzag directions. Since the set 1 has the lowest $\mu_{x}$ and
the highest $\mu_{y}$, the set 3 is vice verse, we report results for the sets
1 and 3 as representative cases.

Here should be noted that in contrast to binding energies of direct excitons
in TMDCs and Xenes, the binding energies of direct excitons in phosphorene
have more complicated dependence on the anisotropic reduced masses. In
phosphorene the central symmetry is strongly broken due to the large
anisotropy of the electron and hole effective masses. This can observed, for
example, in the encapsulated by hBN phosphorene. The insert in Fig. 3$(b)$
shows the dependence of the magnetic field contribution to the binding energy
of the direct magnetoexciton on the ratio $\frac{\mu_{y}}{\mu_{x}}$. Here we
plot $\Delta E$ using the values $\frac{\mu_{y}}{\mu_{x}}$ from Table 1.
$\Delta E$ for $\frac{\mu_{y}}{\mu_{x}}=8,9,10,11,12.5$ and 14 are calculated
by extrapolation of $\mu_{x}$ from Refs. [82, 83, 84, 1] and obtaining
corresponding $\mu_{y}$ using the value of $\frac{\mu_{y}}{\mu_{x}}$. The
bigger $\frac{\mu_{y}}{\mu_{x}}$ ratio corresponds to the more anisotropic
system that becomes quasi-one-dimensional [82], and the increase of
$\frac{\mu_{y}}{\mu_{x}}$ leads to the increase of the contribution due to the
magnetic field to the binding energy of the magnetoexciton. The character of
these magnetoexcitons is that their spatial distribution of wave functions is
strongly anisotropic. These magnetoexcitons form striped-like patterns,
similar to those in nanotubes [89] or nanowires [90] that are 1D systems.

Interestingly, the set 1 with the lowest $\mu_{x}$ and the highest $\mu_{y}$
has the lowest 1$s$ state binding energy, but for states 2$s$, 3$s$, and 4$s$
the set 3 with the highest $\mu_{x}$ and the lowest $\mu_{y}$ gives the lowest
binding energies. The binding energy percent difference between the sets 1 and
3 in states 2$s$, 3$s$, and 4$s$ is 6%, 7%, and 11%, respectively. Also, based
on Fig. 3($a$), for states 1$s$-3$s$ the set 1 gives the highest energy
contribution to the binding energy and the set 3 gives the lowest $\Delta E$.
However, for the state 4$s$ the trend is reversed, and the set 3 gives the
highest $\Delta E$, and the set 1 gives the lowest $\Delta E$. In addition, in
contrast to magnetoexcitons that dissociate in the TMDCs and Xenes monolayers
encapsulated by hBN [27], in states 3$s$ and 4$s$ the magnetoexcitons in
phosphorene monolayers stay bound when the magnetic field is varied in the
range from 0 to 60 T. This feature is also the result of quasi-1D character of
magnetoexcitons in phosphorene.

###  III.2 Contribution from the external magnetic field to the binding
energies of magnetoexcitons in bilayer and vdW heterostructure

Figure 4: The energy contribution from the magnetic field to the binding
energy of Rydberg states 1$s$-4$s$ for the FS bilayer obtained using $V_{RK}$
and $V_{C}$ potentials. The boundary dashed and solid curves correspond to the
sets of masses 1 and 3, respectively. The yellow and green shaded areas
correspond to the calculations with $V_{RK}$ and $V_{C}$ potentials,
respectively.

For indirect magnetoexcitons in the FS bilayer phosphorene and the vdW
heterostructure, we examine how interaction potential affects Rydberg states
binding energies and the energy contribution from the external magnetic field
to the binding energies. In addition, for the vdW heterostructure, we consider
an additional degree of freedom to tune binding energies by number of hBN
layers that separate two phosphorene monolayer.

In Fig. 4, we report results for FS bilayer for the sets 1 and 3. The set 1
gives an upper bound on the energy contribution and is denoted by the dashed
line. The set 3 gives a lower bound on $\Delta E$ and is denoted by the solid
line. As was observed in the monolayer, in the state 4$s$ in bilayer there is
a flip of boundaries: the set 1 gives a lower boundary and the set 3 gives an
upper boundary. Similar to results reported for magnetoexcitons in Xenes and
TMDCs bilayer systems in Ref. [27, 23], in phosphorene bilayer, the following
relation always holds $\Delta E_{RK}>\Delta E_{C}$. In addition, it is worth
mentioning that energy contribution from the magnetic field to the binding
energy is higher for indirect magnetoexcitons in the FS bilayer system than
for direct magnetoexcitons in FS phosphorene. Based on the results in Table 2
for the bilayer phosphorene, Rydberg states binding energies calculated with
the Coulomb potential are twice as big as the binding energies calculated with
the Rytoval-Keldysh potential. Similar to the monolayer system,
magnetoexcitons in the bilayer in states 2$s$, 3$s$, and 4$s$ are more bound
than magnetoexcitons in TMDCs bilayer system.

(a) 1$s$ |  (b) 2$s$
---|---
|
(c) 3$s$ |  (d) 4$s$
|
Figure 5: The energy contribution for indirect magnetoexcitons in the vdW
heterostructures as a function of the magnetic field and the number of hBN
layers. The data are plotted for the set of masses 3.

For the vdW heterostructure, in Fig. 5 we compare $\Delta E_{RK}$ and $\Delta
E_{C}$. Results are given for the set 3 that is taken to be the representative
case. We plot $\Delta E$ as a function of the external magnetic field and the
number of hBN layers. As can be seen from Fig. 5 as the number of hBN layers
increases $\Delta E_{RK}$ and $\Delta E_{C}$ converge. Once again, in contrast
to indirect magnetoexcitons in TMDCs and Xenes double-layer heterostructure,
magnetoexcitons in the vdW phosphorene heterostructure stay bound in states
3$s$ and 4$s$ while the magnetic field is varied between 0 and 60 T. This can
be explained by the fact that the binding energies of states 3$s$ and 4$s$ are
higher for indirect magnetoexcitons in phosphorene due to their effectively
quasi-one-dimensional nature. The quasi-one-dimensionality of the
magnetoexcitons in phosphorene is also demonstrated by the comparison of the
contribution of the external magnetic field to binding energies of the vdW
phosphorene heterostructure shown in Fig. 5 and TMDC heterostructures: $\Delta
E$ due to the magnetic field in the vdW phosphorene heterostructure is always
higher than the one in TMDC heterostructures. We report in Table 3 the binding
energies of Rydberg states 1$s$, 2$s$, 3$s$, and 4$s$ for $N$=1,2,3,4,5,6 when
Eq. 4 is solved with $V_{RK}$ and $V_{C}$. The binding energies calculated
using the Coulomb potential are always larger than the one obtained with the
RK potential. The increase of the phosphorene interlayer distance in the vdW
heterostructure leads to the deduction of electron-hole interaction and
results in the decrease of the binding energy of indirect magnetoexcitons.

Table 3: Binding energies for indirect excitons in the van der Waals heterostructures calculated using $V_{RK}$ and $V_{C}$ potentials for the set of masses 3. Energies are given in meV. $N$ | $V_{RK}$ | $V_{C}$
---|---|---
|  1$s$ |  2$s$ |  3$s$ |  4$s$ |  1$s$ |  2$s$ |  3$s$ |  4$s$
1 |  98.50 |  55.13 |  37.45 |  28.11 |  128.38 |  66.00 |  43.02 |  30.91
2 |  86.39 |  51.29 |  35.27 |  27.15 |  104.90 |  58.98 |  39.30 |  29.39
3 |  76.91 |  47.76 |  33.31 |  26.15 |  89.28 |  53.38 |  36.34 |  27.96
4 |  69.32 |  44.60 |  31.55 |  25.16 |  78.03 |  48.84 |  33.90 |  26.64
5 |  60.10 |  41.80 |  29.96 |  24.30 |  69.49 |  45.06 |  31.82 |  25.43
6 |  57.93 |  39.31 |  28.53 |  23.29 |  62.75 |  41.88 |  30.03 |  24.32

##  IV Diamagnetic coefficients

In this paper, we calculate the diamagnetic coefficients in the same fashion
as was done in Refs. [67, 27, 23]. For phosphorene the diamagnetic
coefficients for magnetoexcitons Rydberg states are reported for the first
time. The magnetic field range between 0 and 30 T has been used to calculated
the DMCs. In this range of the magnetic field one can observe the linear
dependence of energy on $B^{2}$. Here, we follow the notation in literature
and denote the diamagnetic coefficient as $\sigma$. The diamagnetic
coefficients for magnetoexcitons in monolayer and bilayer systems are reported
in Table 4. $\sigma$ for the indirect magnetoexcitons in the vdW
heterostructure are reported in Table 5.

Table 4: The diamagnetic coefficients, $\sigma$, of direct magnetoexcitons in FS and encapsulated by hBN monolayers for the sets of masses 1 and 3. $\sigma$ is given in $\mu$eV/B2. DMCs are obtained for the range of the magnetic field between 0 and 30 T and correspond to $R^{2}$ = 0.9998 for the linear regression model.  State | Monolayer | Bilayer
---|---|---
| Set 1 | Set 3 | Set 1 | Set 3
|  hBN |  FS |  hBN |  FS |  $V_{RK}$ |  $V_{C}$ |  $V_{RK}$ |  $V_{C}$
1$s$ |  0.94 |  0.41 |  0.46 |  0.21 |  0.97 |  0.44 |  0.52 |  0.23
2$s$ |  2.31 |  0.68 |  1.41 |  0.40 |  1.34 |  0.55 |  0.80 |  0.32
3$s$ |  4.13 |  0.99 |  2.96 |  0.67 |  1.81 |  0.71 |  1.22 |  0.47
4$s$ |  |  1.34 |  |  1.88 |  2.32 |  0.93 |  3.30 |  1.30
Table 5: The diamagnetic coefficients, $\sigma$, of indirect magnetoexcitons in the vdW heterostructure for the sets of masses 1 and 3. $\sigma$ is given in $\mu$eV/B2. DMCs are obtained for the range of the magnetic field between 0 and 30 T and correspond to $R^{2}$ = 0.9998 for the linear regression model. State |  $N$ | Set 1 | Set 3
---|---|---|---
|  | $V_{RK}$ | $V_{C}$ | $V_{RK}$ | $V_{C}$
1$s$ |  1 |  2.36 |  1.72 |  1.24 |  0.88
|  2 |  2.80 |  2.22 |  1.49 |  1.16
|  3 |  3.23 |  2.72 |  1.74 |  1.45
|  4 |  3.67 |  3.20 |  2.00 |  1.73
|  5 |  4.09 |  3.68 |  2.27 |  2.02
|  6 |  4.51 |  4.14 |  2.53 |  2.30
2$s$ |  1 |  3.64 |  2.85 |  2.24 |  1.73
|  2 |  4.04 |  3.37 |  2.50 |  2.06
|  3 |  4.46 |  3.88 |  2.77 |  2.39
|  4 |  4.88 |  4.38 |  3.05 |  2.71
|  5 |  |  |  3.34 |  3.03
|  6 |  |  |  3.62 |  3.35
3$s$ |  1 |  |  |  |  3.39

We report the diamagnetic coefficients in Table 4 for two sets of masses.
$\sigma$ for the FS and encapsulated by hBN monolayers phosphorene is obtained
when the magnetoexcitons interact through the $V_{RK}$ interaction, while for
the bilayer $\sigma$ is given when Eq. (4) is solved using both $V_{RK}$ and
$V_{C}$ potential. In the monolayer $\sigma_{1}>\sigma_{3}$, where
$\sigma_{1}$ and $\sigma_{3}$ are DMCs obtained with sets of masses 1 and 3,
respectively. For phosphorene encapsulated by hBN, $\sigma$ can be extracted
in states 1$s$, 2$s$, and 3$s$. In the case of FS phosphorene $\sigma$ can
also be extracted in the state 4$s$. In the case of bilayer phosphorene, the
diamagnetic coefficients can be extracted for 1$s$, 2$s$, 3$s$, and 4$s$
states. This is different to the diamagnetic coefficients of bilayer composed
of Xenes [27] and TMDCs [23], where $\sigma$ cannot be extracted for all
examined states.

The DMCs for the vdW heterostructure are reported in Table 5 for two sets of
masses when Eq. (4) is solved using $V_{RK}$ and $V_{C}$ potentials. $\sigma$
can be extracted for the state 1$s$ for all values of examined $N$, in the
state 2$s$ for the set 3 $\sigma$ can be extracted for $N=1,...,6$. But for
the set 1 in the state 2$s$, $\sigma$ can be extracted only for $N=1,...,4$.
From data listed in Tables 4 and 5, it follows that $\sigma_{RK}>\sigma_{C}$,
where $\sigma_{RK}$ and $\sigma_{C}$ are the DMCs obtained using the RK and
Coulomb potentials, respectively. Finally, our calculations show that the
magnetoexcitons binding energies and DMCs in phosphorene vdW heterostructure
can all be broadly tuned by changing the number of stacked hBN layers. This
serves as a convenient and efficient method for engineering the materials
properties.

##  V Conclusions

In the framework of the excitonic Mott-Wannier model, we study direct and
indirect excitons in Rydberg states in phosphorene monolayers, bilayer, and
vdW heterostructure in the external magnetic field, applied perpendicular to
the monolayer or heterostructure. We calculated the binding energies and DMCs
for the Rydberg states 1$s$, 2$s$, 3$s$, and 4$s$ for magnetoexcitons formed
in these systems.

The magnetic field contribution to the binding energy of magnetoexcitons
strongly depends on the effective masses of electron and hole and their
anisotropy. Interestingly, the binding energies of magnetoexcitons are
strongly correlated with the reduced mass anisotropy that makes
magnetoexcitons in phosphorene effectively quasi-1D quasiparticles. An
overall, the unique anisotropic character of these magnetoexcitons in
phosphorene leads to the larger contribution to the binding energy than in
TMDC materials.

The DMCs demonstrate the strong dependence on the effective electron-hole
masses: DMCs in phosphorene monolayers, bilayer, and vdW heterostructure for
the set 1 are about two fold bigger than the one for the set 2. The other
distinct feature for the diamagnetic coefficients in the bilayer and vdW
heterostructures is that $\sigma$ calculated using the Coulomb potential is
always smaller than DMCs obtained with the Rytova-Keldysh potential. In other
word, the reduced dimensionality and screening increase the diamagnetic
coefficients.

We show that the binding energy of direct and indirect magnetoexcitons can be
tuned by the external magnetic field. Also, the binding energy of indirect
excitons and DMC could be efficiently tuned by the stacking hBN layers. Such
tunability of binding energies and DMCs is potentially useful for the devices
design.

Finally, in summary, we have shown theoretically that the vdW phosphorene
heterostructure is a novel category of 2D semiconductor offering a tunability
of the binding energies of magnetoexcitons by mean of the external magnetic
field and control by the number of hBN layers separating two phosphorene
sheets. Also, the DMCs in the vdW heterostructure are tunable by controlling
the number of hBN layers that separate two phosphorene layers. Thus,
phosphorene provides a unique platform for novel optoelectronic applications
and the exploration of the role of the symmetry breaking in anisotropic
exciton physics.

Acknowledgments. This work is supported by the U.S. Department of Defense
under Grant No. W911NF1810433 and PSC-CUNY Award No. 62261-00 50.

## References

  * Qiao  _et al._ [2014] J. Qiao, X. Kong, Z.-X. Hu, F. Yang, and W. Ji, High-mobility transport anisotropy and linear dichroism in few-layer black phosphorus, Nat. Commun. 5, 4475 (2014).
  * Liu  _et al._ [2014] H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tománek, and P. D. Ye, Phosphorene: An unexplored 2D semiconductor with a high hole mobility, ACS Nano 8, 4033 (2014).
  * Castellanos-Gomez  _et al._ [2014] A. Castellanos-Gomez, L. Vicarelli, E. Prada, J. O. Island, K. L. Narasimha-Acharya, S. I. Blanter, D. J. Groenendijk, M. Buscema, G. A. Steele, J. V. Alvarez, _et al._ , Isolation and characterization of few-layer black phosphorus, 2D Mater. 1, 025001 (2014).
  * Li  _et al._ [2014a] L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Black phosphorus field-effect transistors, Nat. Nanotechnol. 9, 372 (2014a).
  * Bridgman [1914] P. W. Bridgman, Two new modifications of phosphorus, J. Am. Chem. Soc. 36, 1344 (1914).
  * Keyes [1953] R. W. Keyes, The electrical properties of black phosphorus, Phys. Rev. 92, 580 (1953).
  * Warschauer [1963] D. Warschauer, Electrical and optical properties of crystalline black phosphorus, J. Appl. Phys. 34, 1853 (1963).
  * Cartz  _et al._ [1979] L. Cartz, S. R. Srinivasa, R. J. Riedner, J. D. Jorgensen, and T. G. Worlton, Effect of pressure on bonding in black phosphorus, J. Chem. Phys. 71, 1718 (1979).
  * Takao and Morita [1981] Y. Takao and A. Morita, Electronic structure of black phosphorus: Tight binding approach, Physica B+C 105, 93 (1981).
  * Morita [1986] A. Morita, Semiconducting black phosphorus, Appl. Phys. A 39, 227 (1986).
  * Noläng  _et al._ [1990] B. Noläng, O. Eriksson, and B. Johansson, The cohesive energy and band structure of black phosphorus, J. Phys. Chem. Solids 51, 1025 (1990).
  * Novoselov  _et al._ [2004] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, Science 306, 666 (2004).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Avouris  _et al._ [2017] P. Avouris, T. F. Heinz, and T. Low,  _2D materials: Properties and devices_ (Cambridge University Press, 2017).
  * Cheng  _et al._ [2014] R. Cheng, D. Li, H. Zhou, C. Wang, A. Yin, S. Jiang, Y. Liu, Y. Chen, Y. Huang, and X. Duan, Electroluminescence and photocurrent generation from atomically sharp WSe2/MoS2 heterojunction p–n diodes, Nano Lett. 14, 5590 (2014).
  * Lee  _et al._ [2014] C.-H. Lee, G.-H. Lee, A. Zande, W. Chen, Y. Li, M. Han, X. Cui, G. Arefe, C. Nuckolls, T. Heinz, _et al._ , Atomically thin p–n junctions with van der Waals heterointerfaces, Nat. Nanotechnol. 9, 676–681 (2014).
  * Cotleţ  _et al._ [2016] O. Cotleţ, S. Zeytinoǧlu, M. Sigrist, E. Demler, and A. Imamoǧlu, Superconductivity and other collective phenomena in a hybrid Bose-Fermi mixture formed by a polariton condensate and an electron system in two dimensions, Phys. Rev. B 93, 054510 (2016).
  * Manzeli  _et al._ [2017] S. Manzeli, D. Ovchinnikov, D. Pasquier, O. V. Yazyev, and A. Kis, 2D transition metal dichalcogenides, Nat. Rev. Mater. 2, 17033 (2017).
  * Berman and Ya. Kezerashvili [2016] O. L. Berman and R. Ya. Kezerashvili, High-temperature superfluidity of the two-component Bose gas in a transition metal dichalcogenide bilayer, Phys. Rev. B 93, 245410 (2016).
  * Berman and Ya. Kezerashvili [2017] O. L. Berman and R. Ya. Kezerashvili, Superfluidity of dipolar excitons in a transition metal dichalcogenide double layer, Phys. Rev. B 96, 094502 (2017).
  * Wang  _et al._ [2018] G. Wang, A. Chernikov, M. M. Glazov, T. F. Heinz, X. Marie, T. Amand, and B. Urbaszek, Colloquium: Excitons in atomically thin transition metal dichalcogenides, Rev. Mod. Phys. 90, 021001 (2018).
  * Brunetti  _et al._ [2018] M. N. Brunetti, O. L. Berman, and R. Ya. Kezerashvili, Optical absorption by indirect excitons in a transition metal dichalcogenide/hexagonal boron nitride heterostructure, J. Phys.: Condens. Matter 30, 225001 (2018).
  * Kezerashvili and Spiridonova [2021a] R. Y. Kezerashvili and A. Spiridonova, Magnetoexcitons in transition metal dichalcogenides monolayers, bilayers, and van der Waals heterostructures, Phys. Rev. Res. 3, 033078 (2021a).
  * Matthes  _et al._ [2013] L. Matthes, O. Pulci, and F. Bechstedt, Massive dirac quasiparticles in the optical absorbance of graphene, silicene, germanene, and tinene, J. Phys.: Condens. Matter 25, 395305 (2013).
  * Molle  _et al._ [2017] A. Molle, J. Goldberger, M. Houssa, Y. Xu, S.-C. Zhang, and D. Akinwande, Buckled two-dimensional Xene sheets, Nat. Mater. 16, 163 (2017).
  * Zheng  _et al._ [2020] J. Zheng, Y. Xiang, C. Li, R. Yuan, F. Chi, and Y. Guo, All-optically controlled topological transistor based on Xenes, Phys. Rev. Applied 14, 034027 (2020).
  * Kezerashvili and Spiridonova [2021b] R. Y. Kezerashvili and A. Spiridonova, Effects of parallel electric and magnetic fields on Rydberg excitons in buckled two-dimensional materials, Phys. Rev. B 103, 165410 (2021b).
  * Mak  _et al._ [2010] K. F. Mak, C. Lee, J. Hone, J. Shan, and T. F. Heinz, Atomically thin MoS2: A new direct-gap semiconductor, Phys. Rev. Lett. 105, 136805 (2010).
  * Splendiani  _et al._ [2010] A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y. Chim, G. Galli, and F. Wang, Emerging photoluminescence in monolayer MoS2, Nano Lett. 10, 1271 (2010).
  * Tongay  _et al._ [2012] S. Tongay, J. Zhou, C. Ataca, K. Lo, T. S. Matthews, J. Li, J. C. Grossman, and J. Wu, Thermally driven crossover from indirect toward direct bandgap in 2D semiconductors: MoSe2 versus MoS2, Nano Lett. 12, 5576 (2012).
  * Carvalho  _et al._ [2016] A. Carvalho, M. Wang, X. Zhu, A. S. Rodin, H. Su, and A. H. Castro Neto, Phosphorene: from theory to applications, Nat. Rev. Mater. 1, 16061 (2016).
  * Batmunkh  _et al._ [2016] M. Batmunkh, M. Bat-Erdene, and J. G. Shapter, Phosphorene and phosphorene-based materials – prospects for future applications, Adv. Mater. 28, 8586 (2016).
  * Qiu  _et al._ [2017] D. Y. Qiu, F. H. da Jornada, and S. G. Louie, Environmental screening effects in 2D materials: Renormalization of the bandgap, electronic structure, and optical spectra of few-layer black phosphorus, Nano Lett. 17, 4706 (2017).
  * Low  _et al._ [2014] T. Low, A. S. Rodin, A. Carvalho, Y. Jiang, H. Wang, F. Xia, and A. H. Castro Neto, Tunable optical properties of multilayer black phosphorus thin films, Phys. Rev. B 90, 075434 (2014).
  * Ezawa [2014] M. Ezawa, Topological origin of quasi-flat edge band in phosphorene, New J. Phys. 16, 115004 (2014).
  * Brunetti  _et al._ [2019] M. N. Brunetti, O. L. Berman, and R. Y. Kezerashvili, Optical properties of anisotropic excitons in phosphorene, Phys. Rev. B 100, 155433 (2019).
  * Yoon  _et al._ [2021] S. Yoon, T. Kim, S.-Y. Seo, S.-H. Shin, S.-B. Song, B. J. Kim, K. Watanabe, T. Taniguchi, G.-H. Lee, M.-H. Jo, _et al._ , Electrical control of anisotropic and tightly bound excitons in bilayer phosphorene, Phys. Rev. B 103 (2021).
  * Koenig  _et al._ [2014] S. P. Koenig, R. A. Doganov, H. Schmidt, A. H. Castro Neto, and B. Özyilmaz, Electric field effect in ultrathin black phosphorus, Appl. Phys. Lett. 104, 103106 (2014).
  * Chaves  _et al._ [2015] A. Chaves, T. Low, P. Avouris, D. Çak ır, and F. M. Peeters, Anisotropic exciton stark shift in black phosphorus, Phys. Rev. B 91, 155311 (2015).
  * Le and Yarmohammadi [2019] P. Le and M. Yarmohammadi, Tuning thermoelectric transport in phosphorene through a perpendicular magnetic field, Chem. Phys. 519, 1 (2019).
  * Kamban and Pedersen [2020] H. Kamban and T. Pedersen, Interlayer excitons in van der Waals heterostructures: Binding energy, Stark shift, and field-induced dissociation, Sci. Rep. 10 (2020).
  * Rodin  _et al._ [2014a] A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Strain-induced gap modification in black phosphorus, Phys. Rev. Lett. 112, 176801 (2014a).
  * Li  _et al._ [2014b] Y. Li, S. Yang, and J. Li, Modulation of the electronic properties of ultrathin black phosphorus by strain and electrical field, J. Phys. Chem. C 118, 23970 (2014b).
  * Lv  _et al._ [2014] H. Y. Lv, W. J. Lu, D. F. Shao, and Y. P. Sun, Enhanced thermoelectric performance of phosphorene by strain-induced band convergence, Phys. Rev. B 90, 085433 (2014).
  * Zhou  _et al._ [2015a] X. Zhou, W.-K. Lou, F. Zhai, and K. Chang, Anomalous magneto-optical response of black phosphorus thin films, Phys. Rev. B 92, 165405 (2015a).
  * Zhou  _et al._ [2015b] X. Y. Zhou, R. Zhang, J. P. Sun, Y. L. Zou, D. Zhang, W. K. Lou, F. Cheng, G. H. Zhou, F. Zhai, and K. Chang, Landau levels and magneto-transport property of monolayer phosphorene, Sci. Rep. 5, 12295 (2015b).
  * Tahir  _et al._ [2015] M. Tahir, P. Vasilopoulos, and F. M. Peeters, Magneto-optical transport properties of monolayer phosphorene, Phys. Rev. B 92, 045420 (2015).
  * Faria Junior  _et al._ [2019a] P. E. Faria Junior, M. Kurpas, M. Gmitra, and J. Fabian, $k\cdot{}p$ theory for phosphorene: Effective $g$-factors, Landau levels, and excitons, Phys. Rev. B 100, 115203 (2019a).
  * Zhou  _et al._ [2020] J. Zhou, T.-Y. Cai, and S. Ju, Anisotropic exciton excitations and optical properties of hittorf’s phosphorene, Phys. Rev. Res. 2, 033288 (2020).
  * Wu  _et al._ [2020] J.-Y. Wu, W.-P. Su, and G. Gumbs, Anomalous magneto-transport properties of bilayer phosphorene, Sci. Rep. 10, 7674 (2020).
  * Tran  _et al._ [2015] V. Tran, R. Fei, and L. Yang, Quasiparticle energies, excitons, and optical spectra of few-layer black phosphorus, 2D Mater. 2, 044014 (2015).
  * Wang  _et al._ [2015] X. Wang, A. M. Jones, K. L. Seyler, V. Tran, Y. Jia, H. Zhao, H. Wang, L. Yang, X. Xu, and F. Xia, Highly anisotropic and robust excitons in monolayer black phosphorus, Nat. Nanotechnol. 10, 517 (2015).
  * Tian  _et al._ [2020] R. Tian, R. Fei, S. Hu, T. Li, B. Zheng, Y. Shi, J. Zhao, L. Zhang, X. Gan, and X. Wang, Observation of excitonic series in monolayer and few-layer black phosphorus, Phys. Rev. B 101, 235407 (2020).
  * Henriques and Peres [2020] J. C. G. Henriques and N. M. R. Peres, Excitons in phosphorene: A semi-analytical perturbative approach, Phys. Rev. B 101, 035406 (2020).
  * Elliott and Loudon [1960] R. Elliott and R. Loudon, Theory of the absorption edge in semiconductors in a high magnetic field, J. Phys. Chem. Solids 15, 196 (1960).
  * Hasegawa and Howard [1961] H. Hasegawa and R. Howard, Optical absorption spectrum of hydrogenic atoms in a strong magnetic field, J. Phys. Chem. Solids 21, 179 (1961).
  * Shinada and Sugano [1965] M. Shinada and S. Sugano, Optical absorption edge in layer-type semiconductors, J. Phys. Soc. Jpn. 20, 1274 (1965).
  * GorKov and Dzualoshinskii [1967] L. P. GorKov and I. E. Dzualoshinskii, Contribution to the theory of the Mott exciton in a strong magnetic field, Zh. Eksp. Teor. Fiz. 53, 717 (1967).
  * Akimoto and Hasegawa [1967] O. Akimoto and H. Hasegawa, Interband optical transitions in extremely anisotropic semiconductors. II. Coexistence of exciton and the Landau levels, J. Phys. Soc. Jpn. 22, 181 (1967).
  * Walck and Reinecke [1998] S. N. Walck and T. L. Reinecke, Exciton diamagnetic shift in semiconductor nanostructures, Phys. Rev. B 57, 9088 (1998).
  * Luckert  _et al._ [2010] F. Luckert, M. Yakushev, C. Faugeras, A. Karotki, A. Mudryi, and R. Martin, Diamagnetic shift of the A free exciton in CuGaSe2 single crystals, Appl. Phys. Lett. 97, 162101 (2010).
  * Choi  _et al._ [2015] B. Choi, Y. Kim, and J. Song, Diamagnetic shift of a InGaP-AlInGaP semiconductor single quantum well under pulsed-magnetic fields, Appl. Sci. Converg. Technol 24, 156 (2015).
  * Van der Donck  _et al._ [2018a] M. Van der Donck, M. Zarenia, and F. M. Peeters, Excitons, trions, and biexcitons in transition-metal dichalcogenides: Magnetic-field dependence, Phys. Rev. B 97, 195408 (2018a).
  * Van der Donck  _et al._ [2018b] M. Van der Donck, M. Zarenia, and F. M. Peeters, Strong valley Zeeman effect of dark excitons in monolayer transition metal dichalcogenides in a tilted magnetic field, Phys. Rev. B 97, 081109 (2018b).
  * Stier  _et al._ [2018] A. V. Stier, N. P. Wilson, K. A. Velizhanin, J. Kono, X. Xu, and S. A. Crooker, Magnetooptics of exciton Rydberg states in a monolayer semiconductor, Phys. Rev. Lett. 120, 057405 (2018).
  * Han  _et al._ [2018] B. Han, C. Robert, E. Courtade, M. Manca, S. Shree, T. Amand, P. Renucci, T. Taniguchi, K. Watanabe, X. Marie, _et al._ , Exciton states in monolayer MoSe2 and MoTe2 probed by upconversion spectroscopy, Phys. Rev. X 8, 031073 (2018).
  * Spiridonova [2020] A. Spiridonova, Magnetoexcitons in monolayer transition-metal dichalcogenides, Phys. Lett. A 384, 126850 (2020).
  * Li  _et al._ [2014c] Y. Li, J. Ludwig, T. Low, A. Chernikov, X. Cui, G. Arefe, Y. D. Kim, A. M. van der Zande, A. Rigosi, H. M. Hill, _et al._ , Valley splitting and polarization by the Zeeman effect in monolayer MoSe2, Phys. Rev. Lett. 113, 266804 (2014c).
  * Srivastava  _et al._ [2014] A. Srivastava, M. Sidler, A. Allain, D. Lembke, A. Kis, and A. Imamoglu, Valley Zeeman effect in elementary optical excitations of a monolayer WSe2, Nat. Phys. 11 (2014).
  * Aivazian  _et al._ [2015] G. Aivazian, Z. Gong, A. Jones, R.-L. Chu, J.-Q. Yan, D. Mandrus, C. Zhang, D. Cobden, W. Yao, and X. Xu, Magnetic control of valley pseudospin in monolayer WSe2, Nat. Phys. 11, 148 (2015).
  * MacNeill  _et al._ [2015] D. MacNeill, C. Heikes, K. F. Mak, Z. Anderson, A. Kormányos, V. Zólyomi, J. Park, and D. C. Ralph, Breaking of valley degeneracy by magnetic field in monolayer MoSe2, Phys. Rev. Lett. 114, 037401 (2015).
  * Plechinger  _et al._ [2016] G. Plechinger, P. Nagler, A. Arora, A. Aguila, M. Ballottin, T. Frank, P. Steinleitner, M. Gmitra, J. Fabian, P. Christianen, _et al._ , Excitonic valley effects in monolayer WS2 under high magnetic fields, Nano Lett. 16 (2016).
  * Stier  _et al._ [2016] A. Stier, K. Mccreary, B. Jonker, J. Kono, and S. Crooker, Exciton diamagnetic shifts and valley Zeeman effects in monolayer WS2 and MoS2 to 65 Tesla, Nat. Commun. 7, 10643 (2016).
  * Rytova [1967] N. Rytova, Screened potential of a point charge in a thin film, Proc. Moscow State University, Phys. Astron. 3, 30 (1967).
  * Keldysh [1979] L. V. Keldysh, Coulomb interaction in thin semiconductor and semimetal films, JETP Lett 29, 658 (1979).
  * Berkelbach  _et al._ [2013] T. C. Berkelbach, M. S. Hybertsen, and D. R. Reichman, Theory of neutral and charged excitons in monolayer transition metal dichalcogenides, Phys. Rev. B 88, 045318 (2013).
  * Landau and Lifshitz [1977] L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory_ , third edition, revised and enlarged. ed. (Pergamon Press, Oxford, 1977).
  * MacDonald and Ritchie [1986] A. H. MacDonald and D. S. Ritchie, Hydrogenic energy levels in two dimensions at arbitrary magnetic fields, Phys. Rev. B 33, 8336 (1986).
  * Ya. Kezerashvili [2019] R. Ya. Kezerashvili, Few-body systems in condensed matter physics, Few-Body Syst. 60, 52 (2019).
  * Kezerashvili [2020] R. Ya. Kezerashvili,  _Recent Progress in Few-Body Physics_ , edited by N. A. Orr, M. Ploszajczak, F. M. Marqués, and J. Carbonell (Springer International Publishing, Cham, 2020) pp. 827–831.
  * Tran  _et al._ [2014] V. Tran, R. Soklaski, Y. Liang, and L. Yang, Layer-controlled band gap and anisotropic excitons in few-layer black phosphorus, Phys. Rev. B 89, 235319 (2014).
  * Peng  _et al._ [2014] X. Peng, Q. Wei, and A. Copple, Strain-engineered direct-indirect band gap transition and its mechanism in two-dimensional phosphorene, Phys. Rev. B 90, 085402 (2014).
  * Tran and Yang [2014] V. Tran and L. Yang, Scaling laws for the band gap and optical response of phosphorene nanoribbons, Phys. Rev. B 89, 245407 (2014).
  * Páez  _et al._ [2016] C. J. Páez, K. DeLello, D. Le, A. L. C. Pereira, and E. R. Mucciolo, Disorder effect on the anisotropic resistivity of phosphorene determined by a tight-binding model, Phys. Rev. B 94, 165419 (2016).
  * Rodin  _et al._ [2014b] A. S. Rodin, A. Carvalho, and A. H. Castro Neto, Excitons in anisotropic two-dimensional semiconducting crystals, Phys. Rev. B 90, 075429 (2014b).
  * Kumar  _et al._ [2016] P. Kumar, B. S. Bhadoria, S. Kumar, S. Bhowmick, Y. S. Chauhan, and A. Agarwal, Thickness and electric-field-dependent polarizability and dielectric constant in phosphorene, Phys. Rev. B 93, 195428 (2016).
  * Kamban  _et al._ [2020] H. C. Kamban, T. G. Pedersen, and N. M. R. Peres, Anisotropic Stark shift, field-induced dissociation, and electroabsorption of excitons in phosphorene, Phys. Rev. B 102, 115305 (2020).
  * Faria Junior  _et al._ [2019b] P. E. Faria Junior, M. Kurpas, M. Gmitra, and J. Fabian, $k\cdot{}p$ theory for phosphorene: Effective $g$-factors, Landau levels, and excitons, Phys. Rev. B 100, 115203 (2019b).
  * Bondarev [2011] I. V. Bondarev, Asymptotic exchange coupling of quasi-one-dimensional excitons in carbon nanotubes, Phys. Rev. B 83, 153409 (2011).
  * Ya. Kezerashvili  _et al._ [2019] R. Ya. Kezerashvili, Z. Machavariani, B. Beradze, and T. Tchelidze, Trions and biexcitons in Zno/ZnMgO, CdSe/ZnS and CdSe/CdS core/shell nanowires, Phys. E 109, 228 (2019).