# Enhanced thermoelectric efficiency of zigzag bilayer phosphorene nanoribbon;
edge states engineering

Shima Sodagar  Hossein Karbaschi<EMAIL_ADDRESS>Morteza Soltani
<EMAIL_ADDRESS>Mohsen Amini  Department of Physics, University of
Isfahan, Isfahan 81746-73441, Iran

###### Abstract

We theoretically investigate the thermoelectric properties of zigzag bilayer
phosphorene nanoribbons (ZBPNR). We first, draw an analogy between the
extended Su-Schrieffer-Heeger (SSH) ladder and ZBPNR edge states and obtain
their corresponding band structure and wave functions analytically. Then, by
applying the energy filtering method, we show that the electric power and
thermoelectric efficiency of the ZBPNRs can be improved remarkably in the
presence of mid-gap edge states. We also argue how to engineer the edge modes
to further optimize thermoelectric power and efficiency of the system by
applying periodic point potentials at the boundaries.

##  I Introduction

With the increase in world population and the advancement of technology, the
need for energy resources is increasing day by day, but fossil fuel resources
continue to decline. More than $60\%$ of the input energy of industrial
processes is lost as waste heat. The thermal energy is also wasted from
vehicles, electronic devices like computers and mobile phones and even the
human body. If this huge amount of wasted energy is converted back into usable
energy, the problems of energy crisis will be solved to some extent.
Thermoelectric power generation is considered as one of the best solutions for
remarkably reducing the use of fossil fuels and preventing an energy crisis.
Thermoelectric power generators (TEGs) capture thermal energy from the
environment and convert it directly into electricity Rowe (2006); Goldsmid
(2009). TEGs are reliable and maintenance-free solid-state devices that cause
no noise and environmental pollution. TEGs have been studied for many years,
but their use was practically very limited until the last decade. Despite the
many advantages of TEG over conventional power generators, TEGs suffer from
inadequate output power and low efficiency. Many studies have shown that using
low-dimensional materials for thermoelectric purposes is an efficient solution
to overcome the problem of low efficiency in TEGs Hicks and Dresselhaus
(1993a, b); Balandin and Lazarenkova (2003); Murphy et al. (2008); Esposito et
al. (2009); Leijnse et al. (2010); Karbaschi et al. (2016); Zhang et al.
(2016); Hung et al. (2016); Karbaschi et al. (2020); Rezaei et al. (2021). The
Improved efficiency in low-dimensional thermoelectric modules has two main
reasons. First is the existence of novel quantum effects related to the
confinement of carriers. The next reason is the increase of boundaries in the
system, which increases the scattering of phonons and reduces the thermal
conductivity.

Two-dimensional (2D) materials refer to classes of low-dimensional crystalline
solids involving a single or few layers with weak inter-layer bondings and
strong intra-layer interactions in a layer. 2D materials have been widely
investigated because of their special structures and unique properties which
are very different from their corresponding 3D materials. This sort of
materials have extensive applications due to their novel transport and optical
properties. Recently, the few-layer form of black phosphorus, namely
phosphorene, has been exfoliated successfully Li et al. (2014); Lu et al.
(2014), and attracted huge interest due to its outstanding optical,
mechanical, electronic, and thermoelectric properties. Unlike graphene, which
has a zero-energy bandgap, phosphorene is a semiconductor with a large energy
bandgap that can be adjusted by defects, electric field, and strain. The
energy bandgap can also be modulated by the number of layers in a stack that
makes bilayer phosphorene more interesting than its monolayer Takao and Morita
(1981); Dai and Zeng (2014); Carvalho et al. (2016). In addition to the high
carrier mobility and large band-gap, bilayer phosphorene is easy to fabricate,
which makes it more attractive to study.

The presence of edges in the ribbons considerably affects the electronic
structure of phosphorene. It has been pointed out a prominent feature of
phosphorene nanoribbons with zigzag edges is the presence of quasi-flat edge
modes in the middle of the bandgap. Mid-gap states have distinctive properties
from the bulk band and play an important role in conduction. Each phosphorene
layer grants a band of two-fold degenerate edge modes. So, two and four
branches are appearing in the bandgap of energy for monolayer and bilayer
zigzag phosphorene nanoribbons correspondingly.

One way to improve thermoelectric properties is to filter electrons at
specific energies by means of a boxcar transmission function. Previous studies
have shown that such a transmission function can block electrons that are
responsible for the reduction of electrical output power and then increase
thermoelectric efficiency Whitney (2014, 2015); Karbaschi et al. (2016).

The mid-gap edge states transmission function of some nanoribbons probably is
the ideal boxcar function. There are usually two issues with using this
transmission function. The first drawback is the small energy gap and the
overlap of the edge states with the valence or conduction bands. The second
issue for most cases is that the width of this transmission function is much
larger than that which is suitable for use for thermoelectric
applicationsKarbaschi et al. (2020); Rezaei et al. (2021). The appropriate
width for the transmission function is determined in proportion to the
operating temperatures of the leads. In this study, we have proposed that the
application of periodic on-site potentials at the edges of zigzag bilayer
phosphorene nanoribbons can tune the width of edge states’ transmission
function and introduce this structure as a worthy candidate for thermoelectric
purposes.

This article is organized as follows. In section II, we generally introduce
the bilayer phosphorene nanoribbons crystal structure with inter-layer and
intra-layer hopping parameters and introduce the TB Hamiltonian in order to
study band structure. The rest of section II is devoted to discussion on
zigzag bilayer phosphorene nanoribbons wave functions using the analogy
between two-dimensional S-S-H ladder and bilayer phosphorene nanoribbons. In
the following, we deal with the effect of all hopping parameters including
intra-layer and inter-layer on wave function and dispersion of edge modes. In
section III, methods for calculating thermoelectric properties are introduced
and in section IV, we presente the results of numerical study of the
thermoelectric properties of zigzag bilayer phosphorene in the presence of
periodic on-site potentials at both edges using the
Landauer-B$\mathrm{\ddot{u}}$ttiker formalism. Lastly, we finish the article
with the conclusion in section V.

Figure 1:  (a) five inter-layer and ten intra-layer hopping parameters in top
view of bilayer phosphorene crystal structure. The blue rectangle indicates
the unit cell. (b) Bilayer phosphorene crystal structure from th side view.
(c) Schematic of considered system. (d) Band structure of zigzag bilayer
phosphorene.

##  II Edge states in bilayer Phosphorene

###  II.1 Bilayer phosphorene structure

Bilayer Phosphorene can be generally described as a two-dimensional material
in which two puckered structure monolayer of Phosphorene are stacked in
parallel. The layers of covalently bonded phosphorus atoms are stuck weakly to
each other by Van der Waals interactions. Bilayer phosphorene can be
considered with different stacking orders. In the A-B stacking which is
considered in this study, the top layer is moved by half of a unit cell along
x or y with respect to the bottom layer, figures 1(a) and (b).

In this study, we focus on nanoribbons with zigzag edges. A monolayer of
phosphorene includes two puckered-structure sublayers. The unit cell of the
phosphorene monolayer contains four atoms and for an extra layer, four
additional atoms are added to the unit cell. The tight-binding Hamiltonian
model is a proper technique to investigate the qualitative electronic
properties of phosphorus based systems Li et al. (2017); Rudenko et al.
(2015):

$\displaystyle H={}$ $\displaystyle\sum_{\langle i,j\rangle\langle
i^{\prime},j^{\prime}\rangle}t^{\parallel}_{\mathbf{i}}(\hat{c}^{{\dagger}A_{1(2)}(B_{1(2)})}_{i,j}\hat{c}^{B_{1(2)}(A_{1(2)})}_{i^{\prime},j^{\prime}}$
(1) $\displaystyle+\sum_{\langle i,j\rangle\langle
i^{\prime},j^{\prime}\rangle}t^{\perp}_{\mathbf{i}}(\hat{c}^{{\dagger}A_{1}(B_{2})}_{i,j}\hat{c}^{B_{2}(A_{1})}_{i^{\prime},j^{\prime}}$
$\displaystyle+\sum_{\langle
i,j\rangle}\epsilon_{i,j}(\hat{c}^{{\dagger}A_{1(2)}(B_{1(2)})}_{i,j}\hat{c}^{A_{1(2)}(B_{1(2)})}_{i,j}),$

where the first and second terms describe the intra-layer and inter-layer
Hamiltonian respectively. In order to avoid confusion let us label each site
by two indexes $(i,j)$ where $i$ denotes the row number along the
$x$-direction and $j$ denotes the position of the site along the
$y$-direction. In Hamiltonian of Eq. (1), $t_{\mathbf{i}}$ describe hopping
parameters between the $\langle i,j\rangle$ and $\langle
i^{\prime},j^{\prime}\rangle$ sites which are shown in figure 1(a).
$\epsilon_{i,j}$ is the on-site energy and there are set to zero for all
lattice sites when no defect and external field exist in the lattice. The band
structure of a ZBPNR can be well described by considering 15 adequate tight-
binding hopping integrals (Table1), 5 inter-layer and 10 intra-layer
Soleimanikahnoj and Knezevic (2017).

In what follows we consider a ZBPNR with periodic-boundary condition in the
x-direction and open boundary condition in the y-direction. Therefore, we can
use the wave number in the $x$-direction, $k_{x}$, as a good quantum number to
represent the band structure of the system. Figure 1 (d) illustrates the
energy band structure of the ZBPNR with $N_{w}=50$ zigzag chains along the
width of ribbon using the above-mentioned hopping parameters.

Table 1: 10 intra-layer ($t^{\parallel}$) and 5 inter-layer ($t^{\perp}$) hopping parameters of bilayer phosphorene, shown in figure 1(a). $\mathbf{i}$ | $t^{\parallel}_{\mathbf{i}}$ | $d_{\mathbf{i}}(\mathrm{\AA})$ | $\mathbf{i}$ | $t^{\parallel}_{\mathbf{i}}$ | $d_{\mathbf{i}}(\mathrm{\AA})$ | $\mathbf{i}$ | $t^{\perp}_{\mathbf{i}}$ | $d_{\mathbf{i}}(\mathrm{\AA})$
---|---|---|---|---|---|---|---|---
1 | -1.486 | 2.22 | 6 | 0.186 | 4.23 | 11 | 0.524 | 3.60
2 | 3.729 | 2.24 | 7 | -0.063 | 4.37 | 12 | 0.180 | 3.81
3 | -0.252 | 3.31 | 8 | 0.101 | 5.18 | 13 | -0.123 | 5.05
4 | -0.071 | 3.34 | 9 | -0.042 | 5.37 | 14 | -0.168 | 5.08
5 | 0.019 | 3.47 | 10 | 0.073 | 5.49 | 15 | 0.005 | 5.44

###  II.2 Analogy between effective SSH ladder model and bilayer phosphorene
edge states

In this section, we introduce a correspondence between an effective SSH ladder
model and ZBPNR system to understand the behavior of edge states analytically.
We discuss how the analogy between such systems may be used in order to find
the dispersive properties of the edge bands in ZBPNRs which will be used in
the following sections. Due to translational invariant along the
$x$-direction, we can make the following Fourier transform:

$\hat{c}^{\dagger}_{(i,j)}=\sum_{k_{x}}e^{ik_{x}i}\hat{c}^{\dagger}_{k_{x},j}$
(2) $\hat{c}_{(i,j)}=\sum_{k_{x}}e^{-ik_{x}i}\hat{c}_{k_{x},j}$ (3)

Since our main goal is to obtain the the edge-states behavior, we restrict our
focus on the hopping parameters of the problems which do not break the
particle-hole symmetry of the system. In this regard, the most important
hopping parameters are $t^{\parallel}_{1}$ and $t^{\parallel}_{2}$ and
$t^{\perp}_{11}$ and hence, we can write our effective Hamiltonian of the
ZBPNR system as :

$H=H_{1}+H_{2}+H_{11},$ (4)

in which

$H_{1}=\sum_{<i,j>}t^{\parallel}_{1}\hat{c}^{{\dagger}A}_{(i,j+1)}\hat{c}^{B}_{(i,j)}+t^{\parallel}_{1}\hat{c}^{{\dagger}A}_{(i+1,j+1)}\hat{c}^{B}_{(i,j)},$
(5)
$H_{2}=\sum_{<i,j>}t^{\parallel}_{2}\hat{c}^{{\dagger}A}_{(i,j)}\hat{c}^{B}_{(i,j)},$
(6)

and

$H_{11}=\sum_{<i,j>}t^{\perp}_{11}\hat{c}^{{\dagger}A_{1}}_{(i,j)}\hat{c}^{B_{2}}_{(i-1,j)}+t^{\perp}_{11}\hat{c}^{{\dagger}A_{2}}_{(i,j)}\hat{c}^{B_{1}}_{(i,j)}.$
(7)

After performing the Fourier transformation the resulting Hamiltonian
$H(k_{x})=H_{1}(k_{x})+H_{2}(k_{x})+H_{11}(k_{x})$ will be obtained as:

$H_{1}(k_{x})=\sum_{k_{x}}\sum_{j}\tau^{\parallel}_{1}(k_{x})\hat{c}^{{\dagger}A}_{(k_{x},j+1)}\hat{c}^{B}_{(k_{x},j)}+\tau^{\parallel}_{1}\hat{c}^{{\dagger}A}_{(k_{x},j+1)}\hat{c}^{B}_{(k_{x},j)}e^{ik_{x}}),$
(8)
$H_{2}(k_{x})=\sum_{k_{x}}\sum_{j}\tau^{\parallel}_{2}(k_{x})\hat{c}^{{\dagger}A}_{(k_{x},j)}\hat{c}^{B}_{(k_{x},j)},$
(9)

and

$H_{3}(k_{x})=\sum_{k_{x}}\sum_{j}(\tau^{\perp}_{11}(k_{x})\hat{c}^{{\dagger}A_{1}}_{(k_{x},j)}\hat{c}^{B_{2}}_{(k_{x},j)}e^{-ik_{x}}+\tau^{\perp}_{11}\hat{c}^{{\dagger}A_{2}}_{(k_{x},j)}\hat{c}^{B_{1}}_{(k_{x},j)}),$
(10)

where

$\tau^{\parallel}_{1}(k_{x})=2t^{\parallel}_{1}\cos(\frac{k_{x}}{2}),$ (11)
$\tau^{\parallel}_{2}(k_{x})=t^{\parallel}_{2},$ (12)

and

$\tau^{\perp}_{11}(k_{x})=t^{\perp}_{11}e^{(\frac{-ik_{x}}{2})}.$ (13)

Here we used gauge
transformation$(\hat{c}^{B}_{k_{x},j}\rightarrow\hat{c}^{B}_{k_{x},j}e^{(\frac{-ik_{x}}{2})})$
and our results indicate that new hopping parameters
$\tau^{\parallel}_{1}(k_{x})$ and $\tau^{\perp}_{11}(k_{x})$ are now momentum
dependent.

The most interesting property of the resulting Hamiltonian $H(k_{x})$ is its
equivalence to an effective SSH ladder with the effective hopping parameters
shown in Fig 2 (a). In order to illustrate this analogy, we have plotted the
resulting band structures of both ZBPNR system and its corresponding SSH
ladder in Fig. 2 (b). As it is shown, the resulting band structure derived
from the effective SSH ladder coincides with one obtained for the ZBPNR system
with width $N_{w}=50$. This means that our SSH ladder can be used as an
analogy to represent the edge states of ZBPNR at the middle of its energy
spectrum. Here it is appropriate to mention that we found a similar analysis
to introduce an effective SSH ladder for the case of graphene nanoribbons in
Ref. Tan et al. (2021) when we were writing the paper.

Figure 2:  (a) The schematic representation of the SHH ladder model in which
the momentum-dependent hopping parameters $\tau^{\parallel}_{1}(k_{x})$,
$\tau^{\parallel}_{2}(k_{x})$ and $\tau^{\perp}_{11}(k_{x})$ are shown
explicitly. (b) Comparison between band structure of the ZBPNR and its
corresponding ladder analogy with 2 intra-layer and 1 inter-layer hopping
parameters.

###  II.3 Wave functions of ZBPNR edge states

We can now proceed to obtain an analytical expression for the wave functions
associated with the zero-energy state in ZBPNR. For simplicity, we generalize
the notation of Refs. Amini and Soltani (2019); Amini et al. (2019) for the
ZBPNR to show the components of the the wave functions in real space. As we
know, each unit cell consists of two different atoms, denoted as $A$ and $B$
which can be used to label the wave function of the atom on the outermost
sites by $\psi_{n,A}$ and $\psi_{n,A}$ where $n$ shows the cell index on the
edge of ribbon. We start our analysis by writing the Schrödinger equation as:

$\hat{H}\ket{\psi}=E\ket{\psi}$ (14)

There are $N_{w}$ atoms across the width of ribbon in each layer and total
wave function is labeled as
$\psi=\\{\psi_{1,A_{1}},\psi_{1,A_{2}},\psi_{1,B_{1}},\psi_{1,B_{2}},\ldots\\}$.
We are interested in the edge states at $E=0$ which means that we need to
solve the equation $\hat{H}\ket{\psi}=0$ with the following Hamiltonian:

$\hat{H}(k)=\begin{pmatrix}0&0&\tau^{\parallel}_{1}(k_{x})&0&0&0&\ldots\\\
0&0&0&\tau^{\parallel}_{1}(k_{x})&0&0&\ldots\\\
(\tau^{\parallel}_{1}(k_{x}))^{\ast}&0&0&0&(\tau^{\parallel}_{2}(k_{x}))^{\ast}&(\tau^{\perp}_{11}(k_{x}))^{\ast}&\ldots\\\
0&(\tau^{\parallel}_{1}(k_{x}))^{\ast}&0&0&(\tau^{\perp}_{11}(k_{x}))^{\ast}&(\tau^{\parallel}_{2}(k_{x}))^{\ast}&\ldots\\\
0&0&\tau^{\parallel}_{2}(k_{x})&\tau^{\perp}_{11}(k_{x})&0&0&\ldots\\\
0&0&\tau^{\perp}_{11}(k_{x})&\tau^{\parallel}_{2}(k_{x})&0&0&\ldots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\end{pmatrix}\quad$ (15)

We propose the following ansatz for the edge states wave function:

$\begin{split}\ket{\psi(k_{x})}_{edge}&=\sum_{n}(\psi_{(n,A_{1})(k_{x})}\ket{n,A_{1}}+\psi_{(n,B_{1})(k_{x})}\ket{n,B_{1}}\\\
&+\psi_{(n,A_{2})(k_{x})}\ket{n,A_{2}}+\psi_{(n,B_{2})(k_{x})}\ket{n,B_{2}}),\end{split}$
(16)

and try to solve the eigenvalue problem using the recursive procedure
mentioned in Ref. Amini et al. (2019). It is clear that the amplitudes of the
wave functions which are localized at the ribbon’s top edge are zero on even
sites and similar conditions hold for the bottom edge wave functions and odd
sites yielding two sets of answers. In order to proceed further we need to
consider the following two different conditions. If the amplitudes of wave
function are equal on both layers,
$\psi_{n,A_{1},k_{x}}=\psi_{n,A_{2},k_{x}}$, the solution is symmetric and to
that end:

$\psi^{+}_{(n,A_{1},k_{x})}=(\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k_{x})+(\tau^{\perp}_{11}(k_{x}))^{\ast}})^{n}\psi^{+}_{(0,A_{1},k_{x})}$
(17)

By the same token, the solution is asymmetric if the amplitudes of the wave
function on different layers is opposite to each other,
$\psi_{n,A_{1},k_{x}}=-\psi_{n,A_{2},k_{x}}$, and therefore:

$\psi^{-}_{(n,A_{1},k_{x})}=(\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k_{x})-(\tau^{\perp}_{11}(k_{x}))^{\ast}})^{n}\psi^{-}_{(0,A_{1},k_{x})}.$
(18)

Now, we expect the $\ket{\psi^{+(-)}{(k_{x})}}_{(edge)}$ to be normalized.
This implies that$\langle{\psi^{+(-)}}\ket{\psi^{+(-)}}_{edge}=1$.

###  II.4 The effect of other hopping parameters on wave function and
dispersion of edge-bands

The resulting edge bands obtained in the previous sub-section is not
dispersive. However, in this section we consider the effect of other hopping
parameters which we ignored before and discuss their effects on the dispersion
of our resulting edge bands. Therefore this section starts out with a quick
look at Table1, which shows that except for the main hopping parameters
(represented in figure 2 (a)), there are some considerable hopping parameters
that may have an influence on the wave function and dispersion of the edge
band. Such important hopping parameters are shown in figure 3 (a) by dash/dot
lines. We first, turn our attention to the hopping integral
$t^{\parallel}_{6}$ which has a relatively larget value among the intra-layer
hopping parameters. This hopping integral connects inequivalent atoms at the
upper and lower layers and can not change the dispersion of the edge bands.
This is because as we have already discussed, the amplitudes of wave function
is zero on even sites. However, it is interesting to calculate its effect on
the analytic form of the wave function.

Figure 3:  (a) The same as figure 2 (a), but for the relevant new hopping
parameters $\tau^{\parallel}_{3}(k_{x})$, $\tau^{\parallel}_{6}(k_{x})$ ,
$\tau^{\parallel}_{8}(k_{x})$ and $\tau^{\perp}_{14}(k_{x})$. (b) The same as
figure 2 (b), but in presence of hopping integral $t_{3}$ which is responsible
for the dispersion of the edge bands.

In order to calculate an analytic form of the wave funcion in presence of
$t^{\parallel}_{6}$ , one needs to obtain the effective momentum dependent
parameter $(\tau^{\parallel}_{6}(k_{x})=t^{\parallel}_{6}exp^{-ik_{x}/2})$ and
then by solving the same Schrodinger equation for the edge states can obtain
the following symmetric and asymmetric solutions:

$\psi^{+}_{(1,A_{1},k_{x})}=\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k_{x})+(\tau^{\perp}_{11}(k_{x}))^{\ast}}\psi^{+}_{(0,A_{1},k_{x})}$
(19)
$\begin{split}\psi^{+}_{(n,A_{1},k_{x})}&=\frac{(-\tau^{\parallel}_{6}(k_{x}))^{\ast}}{\tau^{\parallel}_{2}(k_{x})+(\tau^{\perp}_{11}(k_{x}))^{\ast}}\psi^{+}_{(n-2,A_{1},k_{x})}\\\
&+\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k_{x})+(\tau^{\perp}_{11}(k_{x}))^{\ast}}\psi^{+}_{(n-1,A_{1},k_{x})}\end{split}$
(20)
$\psi^{-}_{(1,A_{1},k_{x})}=\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k)-(\tau^{\perp}_{11}(k_{x}))^{\ast}}\psi^{-}_{(0,A_{1},k_{x})}$
(21)
$\begin{split}\psi^{-}_{(n,A_{1},k_{x})}&=\frac{(-\tau^{\parallel}_{6}(k_{x}))^{\ast}}{\tau^{\parallel}_{2}(k)-(\tau^{\perp}_{11}(k_{x}))^{\ast}}\psi^{-}_{(n-2,A_{1},k_{x})}\\\
&+\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k_{x})-(\tau^{\perp}_{11}(k_{x}))^{\ast}}\psi^{-}_{(n-1,A_{1},k_{x})}.\end{split}$
(22)

Another hopping parameters which seems to be significant is $t^{\perp}_{14}$.
It does not eliminate the degeneracy of edge bands and the effect that the
mentioned parameter has on the deformation of the wave function is exactly the
same as the previous recurrence equation considering that
$\tau^{\perp}_{14}(k_{x})=t^{\perp}_{14}e^{(-ik_{x}/2)}$. $t^{\parallel}_{8}$
is another relevant hopping parameter that we would discuss. This is a hopping
parameter that connects inequivalent sub-lattices from the upper layer of one
chain to an atom in lower layer at a neighboring chain. The presence of this
parameter, similar to $t^{\parallel}_{6}$ and $t^{\perp}_{14}$, does not
remove the degeneracy of edge bands and consequently only affects the
$t^{\parallel}_{1}$ parameter by a small phase change.

The last hopping parameter which we consider here is $t^{\parallel}_{3}$. It
is an important parameter since it is responsible for the dispersion of edge
bands. This parameter connects the nearest neighbor sites of a pair of zigzag
chains in the lower or upper layers. To understand its effect on the edge band
dispersion, we need to perform the same analysis as before which results in
the hopping form $\tau^{\parallel}_{3}(k_{x})=2t^{\parallel}_{3}\cos(k_{x})$.
On the other hand, the presence of $t^{\parallel}_{3}$ shifts the edge states
spectrum with $2t^{\parallel}_{3}\cos(k_{x})$ but the electron-hole symmetry
around fermi level$(E=0)$ will n’t be broken as well as the degeneracy of edge
bands. The results which is plotted in figure 3(b) indicates that the flatness
of the edge bands will be completely removed by the presence of $t_{3}$.

###  II.5 Degeneracy removal using the perturbative analysis

In the remainder of this section we will present the derivation of removing
the degeneracy of the edge bands analytically via perturbation theory. Our
investigation indicates that using the above-mentioned wave functions in the
standard first-order perturbation theory and taking the expectation value of
the most relevant parameters like $t^{\perp}_{12}$ and $t^{\perp}_{13}$ , we
can estimate the dispersion energy of the edge bands. In order to understand
its ladder analogy we have presented the parameter $t^{\perp}_{12}$ in figure
4.

Figure 4: Schematic representation of hopping parameter $t^{\perp}_{12}$ in
the ladder analogy of the ZBPNRs.

Each layer of zigzag phosphorene nanoribbon contributes a band of two-fold
degenerate mid-gap modes, making four mid-gap modes for ZBPNRs. In what
follows we employ the degenerate perturbation theory to obtain the dispersion
of edge bands for ZBPNRs. We consider the Hamiltonian $H$ in such a way that
it can be split into two parts, namely $\hat{H}=\hat{H}_{0}+\hat{H^{\prime}}$,
where $\hat{H^{\prime}}$ is known as the perturbation. Thus, we can write

$\displaystyle\hat{H^{\prime}}=\begin{pmatrix}E_{11}&E_{12}\\\
E_{21}&E_{22}\end{pmatrix}$ (23)

and the perturbative Hamiltonian matrix elements are determined by:

$\displaystyle\begin{cases}E_{11}=\bra{\psi^{+}}\hat{H^{\prime}}\ket{\psi^{+}}\\\
E_{12}=\bra{\psi^{+}}\hat{H^{\prime}}\ket{\psi^{-}}\\\
E_{21}=\bra{\psi^{-}}\hat{H^{\prime}}\ket{\psi^{+}}\\\
E_{22}=\bra{\psi^{-}}\hat{H^{\prime}}\ket{\psi^{-}}.\end{cases}$ (24)

Then, the eigenvalues in the presence of the hopping $t^{\perp}_{12}$ and
$t^{\perp}_{13}$ will be given by the following equation:

$\displaystyle\begin{vmatrix}E-\bra{\psi^{+}}\hat{H^{\prime}}\ket{\psi^{+}}&\bra{\psi^{+}}\hat{H^{\prime}}\ket{\psi^{-}}\\\
\bra{\psi^{-}}\hat{H^{\prime}}\ket{\psi^{+}}&E-\bra{\psi^{-}}\hat{H^{\prime}}\ket{\psi^{-}}\end{vmatrix}=0\quad$
(25)

Applying the same procedure we used in the previous section and considering
the gauge transformation, we try to write the perturbative Hamiltonian as:

$\displaystyle\hat{H^{\prime}}_{12}(k_{x})={}$
$\displaystyle\sum_{k_{x}}\sum_{j}(t^{\perp}_{12})\\{\hat{c}^{{\dagger}A_{1}}_{(k_{x},j+1)}\hat{c}^{A_{2}}_{(k_{x},j+2)}$
(26)
$\displaystyle+\hat{c}^{{\dagger}A_{2}}_{(k_{x},j)}\hat{c}^{A_{1}}_{(k_{x},j+1)}+\hat{c}^{{\dagger}B_{1}}_{(k_{x},j)}\hat{c}^{B_{2}}_{(k_{x},j+1)}e^{(i\frac{k_{x}}{2})}$
$\displaystyle+\hat{c}^{{\dagger}B_{2}}_{(k_{x},j-1)}\hat{c}^{B_{1}}_{(k_{x},j)}e^{(-i\frac{k_{x}}{2})}\\}$
$\displaystyle\hat{H^{\prime}}_{13}(k_{x})={}$
$\displaystyle\sum_{k_{x}}\sum_{j}(t^{\perp}_{13})2\cos(k)\\{\hat{c}^{{\dagger}A_{2}}_{(k_{x},j)}\hat{c}^{A_{1}}_{(k_{x},j+1)}$
(27)
$\displaystyle+\hat{c}^{{\dagger}A_{1}}_{(k_{x},j+1)}\hat{c}^{A_{2}}_{(k_{x},j)}+\hat{c}^{{\dagger}B_{2}}_{(k_{x},j-1)}\hat{c}^{B_{1}}_{(k_{x},j-2)}$
$\displaystyle+\hat{c}^{{\dagger}B_{1}}_{(k_{x},j-2)}\hat{c}^{B_{2}}_{(k_{x},j-1)}\\}$

This equation can be solved and the resulting eigen energies are the
following:

$E^{(1)}_{12}=t^{\perp}_{12}(\sum_{j}\psi^{\ast}_{(j+1,A_{1})}\psi_{(j+2,A_{2})})+t^{\perp}_{12}(\sum_{j}\psi^{\ast}_{(j,A_{2})}\psi_{(j+1,A_{1})})$
(28)

and so considering the
$\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k_{x})+(\tau^{\perp}_{11}(k_{x}))^{\ast}}=x$
and
$\frac{-\tau^{\parallel}_{1}(k_{x})}{\tau^{\parallel}_{2}(k_{x})-(\tau^{\perp}_{11}(k_{x}))^{\ast}}=y$,
matrix elements of perturbation hamiltonian consist of:

$E_{11}=t^{\perp}_{12}(x)(\frac{1+x^{2}}{1-x^{2}})$ (29)
$E_{22}=t^{\perp}_{12}(y)(\frac{1+y^{2}}{1-y^{2}})$ (30)
$E_{12}=t^{\perp}_{12}(\frac{y}{(1-x^{\ast})(1-y)})(x^{\ast}y+1)$ (31)

and about another hopping parameter:

$\begin{split}E^{(1)}_{13}&=t^{\perp}_{13}(2\cos(k_{x}))(\sum_{j}\psi^{\ast}_{(j,A_{2})}\psi_{(j+1,A_{1})})\\\
&+t^{\perp}_{13}(2\cos(k_{x}))(\sum_{j}\psi^{\ast}_{(j+1,A_{1})}\psi_{(j,A_{2})})\end{split}$
(32) $E_{11}=2t^{\perp}_{13}\cos(k_{x})(\frac{x+(x^{\ast})}{1-x^{2}})$ (33)
$E_{22}=2t^{\perp}_{13}\cos(k_{x})(\frac{y+(y^{\ast})}{1-y^{2}})$ (34)
$E_{12}=2t^{\perp}_{13}\cos(k_{x})(\frac{1}{(1-x^{\ast})(1-y)})(x^{\ast}+y)$
(35)

This result is consistent with our numerical calculations and our conclusion
is that hopping parameters $t^{\perp}_{12}$ and $t^{\perp}_{13}$ are the
responsible parameters for the removal of the degeneracy between energy bands
associated with each layer.

##  III Thermoelectric properties

In the thermoelectric power generation process, electrical power (output) is
generated in cost of losing heat from hot lead (input). Therefore, it is
desirable to produce the maximum amount of electrical power with the least
heat flow. The efficiency of TEGs is then defined as:

$\eta^{TEG}(\%)=\frac{P^{el}}{Q_{hot}}=\frac{VI}{Q_{hot}}\times 100\%,$ (36)

where $I$ and $Q_{hot}$ are electric and heat currents and can be obtained
using the Landauer-B$\mathrm{\ddot{u}}$ttiker formalism and $V$ denotes the
bias voltage. The electric current in the ballistic transport regime can be
obtained from the following equation:

$I=\frac{2q_{e}}{h}\int dE\;T_{HC}(E)\left(f_{H}-f_{C}\right),$ (37)

where $q_{e}$, $h$ and, $T_{HC}(E)$ denote the elementary charge, the Planck
constant and, total transmission function, respectively. Also, $f_{H}$ and
$f_{C}$ are the Fermi-Dirac distribution functions at hot and cold leads.

The heat current has two electron and lattice parts. In order to calculate the
electron part of the heat current, $e$ should be replaced with the energy
which each electron carries.

$Q_{hot}=\frac{2}{h}\int dE\;T_{LR}(E)(E-\mu_{h})\left(f_{L}-f_{R}\right).$
(38)

Many studies have shown that the lattice part of thermal conductivity of low-
diemensional materials is much less than their corresponding bulk structures
which is mainly due to increased scattering of phonons at the boundaries.
Since the reduction of thermal conductivity leads to an increase in
thermoelectric efficiency, this is the basis for the idea of using low-
dimensional materials for thermoelectric applications.

To further increase the efficiency, several methods have been proposed that
increase the phonon scattering while having no significant destructive effects
on the electronic properties. The use of superlattices Karbaschi et al. (2016)
and heterojunctions Wu et al. (2020); Chen et al. (2020) and introducing
strain Chen and Chen (2020) and defects Rezaei et al. (2021) are the most
relevant approaches to thermal conductivity reduction.

Optimizing the electronic properties of the system is another approach to
enhance electric power and thermoelectric efficiency. As mentioned earlier,
the energy filtration method by means of a boxcar transmission function leads
to optimal output power and efficiency Whitney (2014, 2015); Karbaschi et al.
(2016). A boxcar transmission function can eliminate the destructive current
of electrons that flow from cold to hot lead and cause the enhancement of the
output power.

We select zigzag bilayer phosphorene nanoribbons for this study for two
reasons. The first is that among two-dimensional structures it has a small
thermal conductivity, and the second is the presence of mid-gap edge modes ,
which provide a desired boxcar transmission function. The width of edge states
transmission function in pristine bilayer zigzag phosphorene nanoribbons is
very large and not suitable for room temperature thermoelectric applications.
To regulate the width of this function, on-site potentials have been applied
periodically at the two edges of the nanoribbon. The existence of on-site
potentials also increases the phonon scattering and significantly reduces the
lattice thermal conductivity of nanoribbon. In this study and for the cases of
dense on-site potentials, we neglected the lattice part of thermal
conductivity.

The transmission function can be calculated using the
Landauer-B$\mathrm{\ddot{u}}$ttiker formalism:

$T_{LR}(\epsilon)=Tr(\Gamma_{L}G^{ret}\Gamma_{R}{G^{ret}}^{\dagger}).$ (39)

Here,
$\Gamma_{L(R)}(\epsilon)=i[\Sigma_{L(R)}(\epsilon)-\Sigma_{L(R)}^{\dagger}(\epsilon)]$
is the coupling between device and left(right) lead and $G^{ret}$ is the
retarded Green’s functions:

$G^{ret}(\epsilon)=\frac{1}{\epsilon-
H_{center}-\Sigma_{L}{(\epsilon)}-\Sigma_{R}{(\epsilon)}}$ (40)

where, $H_{center}$ is the device Hamiltonian and $\Sigma_{R}$ ad $\Sigma_{L}$
are the coupling to the leads retarded self energies

$\Sigma_{L}(\epsilon)=\tau_{LD}^{\dagger}~{}G^{ret}_{L}~{}\tau_{LD},\\\
\Sigma_{R}(\epsilon)=\tau_{DR}~{}G^{ret}_{R}~{}\tau_{DR}^{\dagger},$ (41)

in which $\tau_{DR}$ and $\tau_{LD}$ are the coupling Matrices between device
and right and left leads. Also, $G^{ret}_{R(L)}=[\epsilon-H_{L(R)}]^{-1}$ is
the right(left) lead retarded Green’s function and, $H_{L}$ and $H_{R}$ are
the left and right leads Hamiltonian.

The second law of thermodynamics sets a maximum to the heat engine’s
efficiency operating at two temperatures. This ideal limit which is called
Carnot efficiency, depends on the temperatures of the cold and hot leads and
is given by:

$\eta_{Carnot}=1-\left(\frac{T_{c}}{T_{h}}\right)$ (42)

Figure 5:  (a), (d), and (g) band structure plotted over a wide range of
energies for the cases of m=2, 3, and 4, respectively. (b), (e), and (h) show
the band structure and (c) , (f) and (i) demonstrate the transmission function
close to the edge states corresponding to (a), (d), and (g), respectively. The
blue circles indicate the position of on-site potentials for each case.

This ideal efficiency is obtained when electrons of only one certain energy
are transported, which means $T=\delta(E-E_{0})$. Although very narrow
transmission functions lead to high efficiencies, they provide negligible
electrical power output. The use of the boxcar transmission function is the
solution to optimize both electric power output and efficiency.

Such a boxcar transmission function blocks all electrons except low-energy
electrons close to Fermi energy, and since electron-phonon scattering is
significant for high-energy electrons, we omit this electron-phonon
interaction in this study. Also, due to the presence of many impurities in the
system, the phonon thermal conductivity will be significantly reduced. In this
study, we have omitted the heat flow due to phonon thermal conductivity, but
it can later be added to the calculations as a disturbance

Figure 6:  (a), (c), and (e) electrical power output and (b), (d), and (f) the
thermoelectric efficiency as a function of $V$ and $\mu$ for the cases of m=2,
3, and 4, respectively.

##  IV Results

In the nonlinear calculations, the hot and cold lead temperatures are assumed
to be $T_{h}=400~{}\mathrm{K}$ and $T_{c}=300~{}\mathrm{K}$, so the ideal
Carnot efficiency at these temperatures is $25\%$.

As shown in figure 5, periodic on-site potentials with three periodicities are
applied to both edges of the zigzag bilayer phosphorene nanoribbons. In the
first case labeled with m = 2, the on-site potentials are applied to two
adjacent atoms, one atom from the top layer and another atom from the bottom
layer of both edges of ZBPNRs while the next two atoms have no on-site
potentials. In the second and third cases with m=3 and m=4, after two atoms
with on-site potentials, the next 4 and 6 atoms remain intact, respectively.

Figures 5(a), 5(d), and 5(g) show the energy bandstructure of ZBPNRs for the
cases of m=2, 3 and 4, respectively. It can be seen that there are a number of
energy bands corresponding to the edge states in the energy bandgap.

Figures 5(b), 5(e), and 5(h) illustrate a zoom of energy bands and figures
5(c), 5(f), and 5(i) show their corresponding transmission function close to
the edge states.

We see that changing the spacing between the on-site potentials results in a
shift in the position of the transmission peaks. At larger distances between
the on-site potentials, the transmission peaks come closer to each other. This
is similar to the situation that occurs in continuous systems such as
rectangular potential barriers and makes sense because changing the distance
between the on-site potentials leads to a change in resonance between their
wave function. So, more transmission peaks can be observed in the bandgap for
cases with larger spacing between on-site potentials. Also, increasing the
spacing between the on-site potentials reduces the energy bands dispersion and
ultimately leads to narrower boxcar transmission functions. Narrowing of the
transmission windows by increasing the spacing between on-site potentials is
caused by decreasing the overlap of their localized states.

Figures 6(a), 6(c), and 6(e) show the electric power and figures 6(b), 6(d),
and 6(f) demonstrate the thermoelectric efficiency, using the transmission
function of figures 5(c), 5(f), and 5(j), plotted on color scale as functions
of bias voltage $V$ and the average chemical potential $\mu$. We can attain
optimal electric power and thermoelectric efficiency by properly regulating
bias voltage and average chemical potential. As indicated in figure 6, the
maximum thermoelectric efficiency is larger for cases with larger spacing
between on-site potentials which have a narrower transmission functions. These
results show that the output power at the maximum efficiency point is much
smaller than the maximum power, so the maximum power point is probably a
better choice than the maximum efficiency point. The highest maximum electric
power is obtained for the case of $m=2$, which has the widest transmission
function, $P^{m=2}_{max}=4550~{}\mathrm{pW}$ and for the other cases
$P^{m=3}_{max}=4450~{}\mathrm{pW}$ and $P^{m=4}_{max}=2930~{}\mathrm{pW}$. In
contrast to maximum electric power, the highest efficiency at maximum electric
power is is achieved for the case of m=4 which has the narrowest transmission
function, $\eta^{m=4}_{maxP}=12.7\%$ and for the other cases
$\eta^{m=3}_{maxP}=11.0\%$ and $\eta^{m=2}_{maxP}=8.6\%$. Considering the
values obtained for maximum electric power and its corresponding efficiency,
it seems that the optimal electric power and thermoelectric efficiency for the
considered temperatures are achieved for the case of m=3, whose maximum
electric power is very close to the case m=2, but with higher efficiency.

##  V Conclusion

In conclusion, we have investigated the thermoelectric power generation
properties of zigzag bilayer phosphorene nanoribbons in the presence of
periodic on-site potentials at the edges and have calculated the generated
electric power and thermoelectric efficiency. Since the width of the edge
state transmission function in pristine ZBPNRs is very large, we don’t expect
good thermoelectric efficiency at room temperature for this structure. By
periodically applying on-site potentials, we have obtained narrow transmission
functions that are suitable for room temperature thermoelectric applications.
Our results show that excellent electric power can be achieved at high
efficiencies by edge states engineering.

## References

  * Rowe (2006) D. M. Rowe, ed., _Thermoelectrics Handbook: Macro to Nano_ (CRC press, 2006).
  * Goldsmid (2009) H. J. Goldsmid, _Introduction to Thermoelectricity_ (Springer, 2009), Springer Series in Materials Science.
  * Hicks and Dresselhaus (1993a) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 12727 (1993a).
  * Hicks and Dresselhaus (1993b) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 16631 (1993b).
  * Balandin and Lazarenkova (2003) A. A. Balandin and O. L. Lazarenkova, Applied Physics Letters 82, 415 (2003).
  * Murphy et al. (2008) P. Murphy, S. Mukerjee, and J. Moore, Phys. Rev. B 78, 161406(R) (2008).
  * Esposito et al. (2009) M. Esposito, K. Lindenberg, and C. van den Broeck, Eur. Phys. Lett. 85, 60010 (2009).
  * Leijnse et al. (2010) M. Leijnse, M. R. Wegewijs, and K. Flensberg, Phys. Rev. B 82, 045412 (2010).
  * Karbaschi et al. (2016) H. Karbaschi, J. Lovén, K. Courteaut, A. Wacker, and M. Leijnse, Phys. Rev. B 94, 115414 (2016).
  * Zhang et al. (2016) Y. Zhang, X. Wang, M. Cleary, L. Schoensee, N. Kempf, and J. Richardson, Applied Thermal Engineering 96, 83 (2016).
  * Hung et al. (2016) N. T. Hung, E. H. Hasdeo, A. R. T. Nugraha, M. S. Dresselhaus, and R. Saito, Phys. Rev. Lett. 117, 036602 (2016).
  * Karbaschi et al. (2020) H. Karbaschi, N. Nouri, M. Rezaei, and G. Rashedi, Nanotechnology 31, 375403 (2020).
  * Rezaei et al. (2021) M. Rezaei, H. Karbaschi, M. Amini, M. Soltani, and G. Rashedi, Nanotechnology 32, 375704 (2021).
  * Li et al. (2014) L. Li, Y. Yu, G. J. Ye, Q. Ge, X. Ou, H. Wu, D. Feng, X. H. Chen, and Y. Zhang, Nature Nanotechnology 9, 372 (2014).
  * Lu et al. (2014) W. Lu, H. Nan, J. Hong, Y. Chen, C. Zhu, Z. Liang, X. Ma, Z. Ni, C. Jin, and Z. Zhang, Nano Research 7, 853 (2014).
  * Takao and Morita (1981) Y. Takao and A. Morita, Physica B+C 105, 93 (1981).
  * Dai and Zeng (2014) J. Dai and X. C. Zeng, The Journal of Physical Chemistry Letters 5, 1289 (2014).
  * Carvalho et al. (2016) A. Carvalho, M. Wang, X. Zhu, A. S. Rodin, H. Su, and A. H. Castro Neto, Nature Reviews Materials 1, 16061 (2016).
  * Whitney (2014) R. S. Whitney, Phys. Rev. Lett. 112, 130601 (2014).
  * Whitney (2015) R. S. Whitney, Phys. Rev. B 91, 115425 (2015).
  * Li et al. (2017) L. L. Li, D. Moldovan, W. Xu, and F. M. Peeters, Phys. Rev. B 96, 155425 (2017).
  * Rudenko et al. (2015) A. N. Rudenko, S. Yuan, and M. I. Katsnelson, Phys. Rev. B 92, 085419 (2015).
  * Soleimanikahnoj and Knezevic (2017) S. Soleimanikahnoj and I. Knezevic, Journal of Computational Electronics 16, 568 (2017).
  * Tan et al. (2021) T. Tan, C. Li, and W. Yao, Phys. Rev. B 104, 245419 (2021).
  * Amini and Soltani (2019) M. Amini and M. Soltani, Journal of Physics: Condensed Matter 31, 215301 (2019).
  * Amini et al. (2019) M. Amini, M. Soltani, and M. Sharbafiun, Phys. Rev. B 99, 085403 (2019).
  * Wu et al. (2020) D. Wu, X.-H. Cao, P.-Z. Jia, Y.-J. Zeng, Y.-X. Feng, L.-M. Tang, W.-X. Zhou, and K.-Q. Chen, Science China Physics, Mechanics & Astronomy 63, 276811 (2020).
  * Chen et al. (2020) X. Chen, M. Pang, T. Chen, D. Du, and K. Chen, ACS applied materials & interfaces (2020).
  * Chen and Chen (2020) X.-K. Chen and K.-Q. Chen, Journal of Physics: Condensed Matter 32, 153002 (2020).