I. INTRODUCTION
High-performance thermoelectric devices are highly desirable in clean energy applications, such as thermoelectric generators and refrigerators.
Reference Goldsmid1–Reference DiSalvo3
However, the bottleneck of such applications is the energy conversion efficiency, which is closely related to the thermoelectric properties of the material used. In general, the thermoelectric efficiency of a material is characterized by a dimensionless positive quantity called figure of merit, ZT, which is defined as
$ZT = {{{G_{\rm{e}}}{S^2}T} \mathord{\left/ {\vphantom {{{G_{\rm{e}}}{S^2}T} {\left( {{\sigma _{\rm{e}}} + {\sigma _{{\rm{ph}}}}} \right)}}} \right. \kern-\nulldelimiterspace} {\left( {{\sigma _{\rm{e}}} + {\sigma _{{\rm{ph}}}}} \right)}}$
, where G
e is electronic conductance, S is the Seebeck coefficient, T is temperature, σe is the electronic thermal conductance, and σph is the phonon thermal conductance. A larger ZT corresponds to a higher efficiency of energy conversion. For traditional metal materials, ZT is constrained by the Wiedemann–Franz (WF) law,
Reference Ashcroft and Mermin4
which states that the ratio of electronic thermal conductance and electronic conductance is fixed upon to temperature, σe/G
e = LT where L = 2.44 × 10−8 WΩK−2 is the Lorenz number. Unfortunately, the Lorenz number is too large, which caps the upper limit of ZT, and thus limits the practical applications of metals as thermoelectric materials. Interestingly, semiconductors, in particular, semiconducting nanostructures allow the breakdown of the WF law.
Reference Dubi and Di Ventra5
In recent years, two-dimensional (2D) materials were found to have many fascinating physical and chemical properties, and thus they have received a great deal of attention in material science. Studies have shown that 2D materials present many novel properties that are not observed in their bulk counterparts. Reference Bhimanapati, Lin, Meunier, Jung, Cha, Das, Xiao, Son, Strano, Cooper, Liang, Louie, Ringe, Zhou, Kim, Naik, Sumpter, Terrones, Xia, Wang, Zhu, Akinwande, Alem, Schuller, Schaak, Terrones and Robinson6,Reference Butler, Hollen, Cao, Cui, Gupta, Gutiérrez, Heinz, Hong, Huang, Ismach, Johnston-Halperin, Kuno, Plashnitsa, Robinson, Ruoff, Salahuddin, Shan, Shi, Spencer, Terrones, Windl and Goldberger7 However, highly conducting 2D materials such as graphene are of an ultra-high thermal conduction, thus limiting their application in thermoelectric devices. Since some of 2D materials are semiconducting, it is expected that these semiconducting 2D materials are promising for use in high-performance thermoelectric devices. Reference Zhang and Zhang8–Reference Zhang, Liu, Cheng, Wei, Liang, Fan, Shi, Tang and Zhang15
Recently, phosphorene, single-layer black phosphorus, was synthesized experimentally. Reference Li, Yu, Ye, Ge, Ou, Wu, Feng, Chen and Zhang16 It was revealed that phosphorene possesses a direct band gap of 1∼2 eV, a high electron mobility and strain-sensitive properties. Reference Cai, Zhang and Zhang17 Unlike graphene, phosphorene has strong anisotropy such that the transport properties of both electrons and phonons depend strongly on the lattice direction. The electric conductivity in the armchair direction was reported to be substantially larger than that in the zigzag direction. Reference Fei, Faghaninia, Soklaski, Yan, Lo and Yang10,Reference Cai, Zhang and Zhang17 However, the anisotropy of thermal conduction is opposite, and thus the phonon thermal conductance in the zigzag direction is larger than that in the armchair direction. Reference Hong, Zhang, Huang and Zeng18–Reference Zhu, Zhang and Li24
It is noted that previous theoretical studies on the thermoelectric properties of 2D phosphorene sheet and nanoribbons mainly focused on the diffusive transport regime based on the Boltzmann transport theory, Reference Fei, Faghaninia, Soklaski, Yan, Lo and Yang10,Reference Lv, Lu, Shao and Sun12,Reference Qin, Yan, Qin, Yue, Cui, Zheng and Su13,Reference Zhang, Liu, Cheng, Wei, Liang, Fan, Shi, Tang and Zhang15 and various values of maximum ZT were reported, ranging from 0.5 to 7. Reference Fei, Faghaninia, Soklaski, Yan, Lo and Yang10,Reference Lv, Lu, Shao and Sun12,Reference Qin, Yan, Qin, Yue, Cui, Zheng and Su13,Reference Zhang, Liu, Cheng, Wei, Liang, Fan, Shi, Tang and Zhang15 The calculation results; however, are critically dependent on the approximations made in estimating the relaxation times of both electrons and phonons. For example, the electronic thermal conductivity in some of the previous studies was estimated based on the WF law. However, the validity of WF law remains unchecked. Reference Ashcroft and Mermin4 Currently, a systematic investigation of the thermoelectric properties of phosphorene in the same theoretical framework is still lacking. An equal-footing treatment of all the transport quantities, including electric conduction, electronic thermal conduction, Seebeck coefficient, and phonon thermal conduction, is crucial for the accurate evaluation of ZT. Furthermore, when the size of a nanostructure is comparable or shorter than the mean free path of transport carrier (electron or phonon), the transport of electron/phonon becomes ballistic rather than diffusive. Meanwhile, quantum effects also play an important role in such coherent transport regime. Hence, it is both important and necessary to consider the size and quantum effects on the thermoelectric properties of the material. However, systematic examination of these effects is still lacking.
In this work, we employ the nonequilibrium Green's function (NEGF) technique to study the thermoelectric properties of phosphorene in a comprehensive manner. All transport quantities such as electronic conductance, Seebeck coefficient, electronic thermal conductance, and phonon thermal conductance, are evaluated within the same framework to obtain a reliable thermoelectric figure of merit. We exactly solve the quantum transport properties in the ballistic transport regime when nonelastic scattering is not important. Our study shows that the maximum value of figure of merit ZT of monolayer phosphorene is 0.2 at room temperature, and reaches 0.5 at the temperature of 600 K. With the detailed analysis of the causes and effects, possible ways of boosting ZT are suggested.
II. METHODS
A. Theory of electron transport
The atomic structure of phosphorene is shown in Fig. 1(a). The lattice constants in the armchair and zigzag direction are 3.26 Å and 4.34 Å, respectively. The thickness of the phosphorene sheet is 5.5 Å. To evaluate the electronic transport properties, we employ the four-band tight-binding model. This model has been widely used to study the electronic properties of phosphorene and its predictions match well with first-principle calculations near the Γ point, where the band gap is formed. Reference Jiang and Park25–Reference Zhou, Zhang, Sun, Zou, Zhang, Lou, Cheng, Zhou, Zhai and Chang27 In this four-band tight-binding model, the Hamiltonian is modeled by the hopping of electrons between the atomic sites. It is written as:
(a) The crystal structure of phosphorene from both top view and side view. We use the transport in the armchair direction as example, where phosphorene is partitioned into the left, center, and right regime along the armchair direction. A chemical potential or temperature bias is applied between the left and right regime to drive the current. The parameters t 1, t 2, t 3, t 4, and t 5 are hopping energies in tight-binding model. The labeling in the reciprocal lattice is shown in the right-hand side. (b) Electron band structure evaluated from tight-binding model. (c) Phonon dispersion relation evaluated from first-principles calculations.
The hopping constants are determined by the relative position of the atoms as illustrated in Fig. 1(a), they are t 1 = −1.22 eV, t 2 = 3.665 eV, t 3 = −0.205 eV, t 4 = −0.105 eV, and t 5 = −0.055 eV. Reference Rudenko and Katsnelson28 The other long-distance hopping possibilities are small and thus ignored. We set the reference energy E i = −0.415 eV so as to shift the middle of band gap to 0 eV for convenience. By using this tight-binding model, the electronic band structure is shown in Fig. 1(b). The bands have degeneracy at the path Y–S–X. A direct band gap is formed at Γ point with a band gap energy of 1.6 eV.
To evaluate the electronic transport properties, we partition the phosphorene into three regimes along the transport direction. We call them left, center, and right regimes, indicated as the red, pink, and blue atoms in Fig. 1(a). The left and right regimes are regarded as electron source and drain, respectively, so it is assumed to have infinite electrons. In the transverse direction of transport, we use Fourier transformation to fold the infinite-large periodic lattice into Brillouin zone as
where now i, j become the atomic indices inside a unit cell and
$\vec {\bi R}$
is labeling the position of the unit cell. We only apply this transformation to the transverse direction so both
$\vec {\bi R}$
and
$\vec {\bi k}$
are actually one-dimensional vector pointing perpendicular to the transport direction. However, in the transport direction, we still work in the real space, by partitioning
${t_{ij}}(\vec {\bi k})$
into submatrices
$t_{ij}^{00}({\vec {\bi k}})$
,
$t_{ij}^{01}({\vec {\bi k}})$
and
$t_{ij}^{10}({\vec {\bi k}})$
. Here
$t_{ij}^{00}({\vec {\bi k}})$
denotes the hopping constant from site i to site j of the same unit cell,
$t_{ij}^{01}({\vec {\bi k}})$
is the hopping constant from site i to site j of the unit cell in the right-hand side, while
$t_{ij}^{10}({\vec {\bi k}})$
is the hopping constant from site i to site j of unit cell in the left-hand side.
In the NEGF formalism,
Reference Wang, Agarwalla, Li and Thingna29–Reference Wang, Wang and Lü31
the transmission coefficients of electrons with energy E and transverse momentum
$\vec {\bi k}$
are given by:
where G r, (G a) are the retarded (advanced) Green's functions, and ΓL (ΓR) are the spectral density of left (right) regime. The retarded Green's function can be evaluated via:
where S′ is a matrix representing the overlap of electron wavefunctions at different atomic sites. Here we assume it to be identity, whereas the electrons are not overlapping. ηi is small imaginary part with small positive η and ΣL, ΣR are the self-energy of the left, right regime. The advanced Green's function can be obtained from the retarded Green's function as
${G^{\rm{a}}}({E,\vec {\bi k}}) = {[ {{G^{\rm{r}}}({E,\vec {\bi k}})} ]^\dag }$
. The spectral densities can be obtained from self-energies as:
The self-energies is obtained from the surface Green's function of the electron source or drain, via
$\Sigma _\alpha ^{\rm{r}}({\vec {\bi k}}) = {t^{01}}({\vec {\bi k}}){g^{\rm{r}}}({E,\vec {\bi k}}){t^{10}}({\vec {\bi k}})$
. The surface Green's function can be evaluated through iterative methods by using the formula:
However, in this work, we used a generalized eigenvalue method, which is illustrated in Ref. Reference Wang, Wang and Lü31. It is computationally more efficient than simple iteration.
Using the above formulas, we can obtain the mode-based transmission function, and hence we can obtain the energy-based transmission function per cross-sectional area A by integrating all modes:
where d is the thickness of the material. Using the transmission function, we can obtain the transport coefficients via:
where β = 1/k B T is the inverse temperature and µ is the chemical potential. Here we assume that the change of chemical potential will not alter the band structure of phosphorene, but only the distribution of electrons. With the transport coefficients, all the transport quantities can be evaluated. The electronic conductance is
where e is the elementary charge. The Seebeck coefficient is
The electronic thermal conductance is:
For the calculation of ZT, we also need to calculate the phonon thermal conductance.
B. Theory of phonon transport
Same as electron transport, we calculate the phonon thermal conductance using nonequilibrium Green's function formalism. In analogue with the hopping matrix used in electronic NEGF framework, we need the interatomic force constants in the phonon NEGF calculation. Here we obtain the interatomic force constant matrix via density functional perturbation theory by using QUANTUM ESSPRESSO package. Reference Paolo, Stefano, Nicola, Matteo, Roberto, Carlo, Davide, Guido, Matteo, Ismaila, Andrea Dal, Stefano de, Stefano, Guido, Ralph, Uwe, Christos, Anton, Michele, Layla, Nicola, Francesco, Riccardo, Stefano, Alfredo, Lorenzo, Carlo, Sandro, Gabriele, Ari, Alexander, Paolo and Renata32 The local density functional was used as the exchange-correlation functional with a cutoff energy of 50 Ry. The Trouiller-Martins type norm-conserving pseudopotential was adopted for treating the electron-ion interaction. We used 14 × 10 × 1 and 8 × 6 × 1 grids for sampling of Brillouin zone of k and q points, respectively. The structure of monolayer phosphorene was fully relaxed until the forces become smaller than 0.01 eV/Å.
Based on the obtained force constant matrix, we can find the phonon dispersion relation, as shown in Fig. 1(c). A frequency gap is formed nearby 300 cm−1. The cut-off frequency is about 450 cm−1, which is much smaller than that of graphene (∼1500 cm−1).
The phonon transmission function can be evaluated using the same NEGF framework described above. The only differences are to replace the hopping matrix t by force constant matrix K, replace the energy parameter E + ηi with frequency parameter (ω + ηi)2, and replace the electron overlapping matrix to the atomic mass matrix, so that we can get the phonon transmission function T ph[ω]. The phonon thermal conductance can be evaluated according to the Landauer formula Reference Wang, Wang and Lü31
The NEGF formalism described here is a powerful tool to handle ballistic transport, for both electrons and phonons. It is a complete quantum mechanical treatment with quantized thermal conductance predicted. However, in this formalism, nonelastic scattering is not taken into consideration, such that all the electrons and phonons can transport without loss of energy. This is an appropriate assumption when the system is at low temperature, such that the population of electrons and phonons are small. In such a scenario, the possibility of collisions including electron–electron, phonon–phonon, and electron-phonon is low. It is also valid when the size of materials is shorter than the mean free path, therefore, the carriers (electrons or phonons) are able to reach the end before any inelastic scattering.
III. RESULTS AND ANALYSIS
A. The transmission functions
The electron and phonon transmission functions are shown in Fig. 2. We find that the electron transmission function reduces to zero at the regime −0.8 eV to 0.8 eV. This is consistent with the electron band structure since the density of states there vanishes. It also exhibits anisotropic behavior with the transmission in the armchair direction being higher than that in the zigzag direction. The transmission function of phonons also shows a gap, consistent with the phonon dispersion relation. In contrast to the electron transmission function, the phonon transmission function in the zigzag direction is higher than that in the armchair direction.
(a) Plot of electron transmission function per area against electron energy at both armchair direction (AM) and zigzag direction (ZZ). (b) Plot of phonon transmission function per area against phonon frequency at both armchair direction (AM) and zigzag direction (ZZ).

B. The transport properties
Figure 3 demonstrates the results of various transport properties of phosphorene against the applied chemical potential. The chemical potential can be tuned by either adding a gate voltage or by changing the carrier density via doping. Figure 3(a) shows the electronic conductance in both armchair (AM) and zigzag (ZZ) directions. The electronic conductance vanishes when the chemical potential is near the middle of the band gap. This is due to the lack of excited carriers. When the chemical potential approaches both valence band maximum (VBM) or conduction band minimum (CBM), the electronic conductance starts to increase, due to the contribution from the thermally excited carriers. By comparing the curves of 300 K and 600 K, the temperature effect is obvious at the VBM or CBM edge but not obvious when chemical potential increases further. The electronic conductance shows a strong anisotropy: The electronic conductance in the armchair direction is much larger than that in the zigzag direction at all ranges of chemical potentials. This observation is consistent with previous experiments. Reference Li, Yu, Ye, Ge, Ou, Wu, Feng, Chen and Zhang16
(a) Electronic conductance per cross-sectional area against chemical potential at different directions and temperatures. (b) The Seebeck coefficient at different directions and temperatures. (c) Electronic thermal conductance σe and phonon thermal conductance σph against chemical potential. (d) Ratio of electronic thermal conductance and electronic conductance with respect to chemical potential. The Lorenz constant is plotted for comparison.
Figure 3(b) shows the Seebeck coefficient of phosphorene. Interestingly, the curve of Seebeck coefficient in the armchair direction and zigzag direction coincide, indicating that the Seebeck coefficient does not exhibit any anisotropic behavior at both temperatures of 300 K and 600 K. This finding is in strong contrast with the electronic conductance. This means that the electronic current driven by either voltage bias or temperature bias shares the same anisotropic behavior, such that their ratio turns out to be isotropic. However, temperature strongly affects the value of maximum Seebeck coefficient: A higher temperature leads to a reduction of Seebeck coefficient.
The thermal transport properties of phosphorene are shown in Fig. 3(c). The electron contribution to the thermal conduction vanishes when the chemical potential is inside the band gap, since no electron is transferred. However, the electronic thermal conductance increases rapidly when the chemical potential approaches the CBM and VBM, and it eventually surpasses the phonon contribution when the chemical potential approaches ±1 eV. This suggests that phonons are dominant in the thermal conduction of phosphorene in most cases, unless it is highly doped. Importantly, for the electrons, the thermal conductance in the armchair direction is much larger than that in the zigzag direction, but it is opposite for the phonon thermal conductance. Hence, overall, the anisotropy of thermal conductance depends on the dominant type of thermal carriers, and it can be tuned via chemical potential.
With the electronic thermal conductance and electronic conductance, we are now able to check the WF law for phosphorene. Hence we plot the conductance ratio σe/(G e T) against chemical potential in Fig. 3(d). In metals, the ratio is equal to Lorenz constant, L = 2.44 × 10−8 WΩK−2, which is shown as a reference value. We find that when the chemical potential is inside the bands, the conductance ratio approaches the Lorenz constant, indicating the validity of the WF law. However, when the chemical potential is around the CBM (or VBM) edge, the conductance ratio becomes smaller and it approaches another constant of 1.11 × 10−8 WΩK−2. We wipe out the values near µ = 0 since the numbers there are not reliable as both the electronic conductance and electronic thermal conductance approach zero. The change of conductance ratio is due to the absence of density of states at those energy levels such that the electronic/thermal conduction is purely caused by thermal excitation. Reference Zhou, Cai, Zhang and Zhang33 Hence the transport mechanism is fundamentally different from metals, thus violating the WF law. Nevertheless, a reduction of conductance ratio is preferred for good thermoelectric performance.
C. The thermoelectric performance
The thermoelectric performance is characterized by figure of merit ZT. The values of ZT of phosphorene are shown in Fig. 4. With varying the chemical potential, ZT exhibits two peaks near the band edges at around ±0.75 eV. These peaks are the compromise of Seebeck coefficient, electronic conductance and thermal conductance. When
$\left| \mu \right| \lt 0.75\;{\rm{eV}}$
, the electronic conductance is too small to harvest electricity. However, when
$\left| \mu \right| > 0.75\;{\rm{eV}}$
, the Seebeck effect is too weak and the thermal conduction due to electrons becomes significant. As a result, the value of ZT reaches maximum near the band edges. The value of ZT in the armchair direction is significantly higher (4-fold) than that in the zigzag direction, due to higher electronic conductance and lower phonon thermal conductance. The predicted maximum ZT in the armchair direction is 0.22 at room temperature (300 K) and reaches 0.50 at 600 K. Recently, using a first-principles calculation on the electron–phonon interaction based on density functional perturbation theory and Wannier interpolation, Liao et al. studied the thermoelectric properties of phosphorene in the frame of diffusive transport.
Reference Liao, Zhou, Qiu, Dresselhaus and Chen34
The theoretically predicted value of ZT is about 0.06 at room temperature, close to our results although ours are obtained at the ballistic limit.
(a) Figure of merit ZT against chemical potential at different directions and temperatures. (b) Figure of merit ZT as a function of temperature when the chemical potential is optimized (µ = −0.75 eV).

To consider the effect of temperature, we plot the ZT against temperature when the chemical potential is set at −0.75 eV. We find that ZT almost increases linearly with temperature. The temperature will suppress the Seebeck coefficient while enhance the conductance ratio, resulting in an overall effect of constant
$Z = {{{G_{\rm{e}}}{S^2}} \mathord{\left/ {\vphantom {{{G_{\rm{e}}}{S^2}} {({{\sigma _{\rm{e}}} + {\sigma _{{\rm{ph}}}}})}}} \right. \kern-\nulldelimiterspace} {\left( {{\sigma _{\rm{e}}} + {\sigma _{{\rm{ph}}}}} \right)}}$
. At a higher temperature, the difference in the thermoelectric performance between the armchair direction and zigzag direction becomes larger.
IV. DISCUSSION
From the aforementioned analysis, we can conclude that the figure of merit ZT is overall effects of many factors and each of them can be important under certain conditions.
When the chemical potentials is close (but not equal) to the middle of band gap (in-weakly doped condition), the Seebeck coefficient is high. Thus the conductance ratio is the major factor that adversely affects its thermoelectric performance, due to the low electron conduction and high thermal conduction. In this regime, the phonon thermal conductance plays a major role in comparison with electronic thermal conductance. Any suppression of phonon transport will significantly enhance the ZT. Therefore, in this regime, we may introduce various techniques to reduce the phonon thermal transport, such as the creation of superlattice to localize the long-wavelength phonons, introduce edge roughness to increase the boundary scattering, or add mass-disorders to suppress the phonon transmissions. Reference Zhang and Zhang35
When the chemical potential is far away from the middle of band gap (in highly-doped condition), the reduction of Seebeck coefficients and the rapid increase of electronic thermal conduction are the major factors that affect ZT. In this regime, the ratio of electronic thermal conductance and electronic conductance obeys the WF law. Hence the ideas of boosting ZT should focus on the electronic structures. Possible methods include defects/impurities engineering, electromagnetic field modulation, or strain manipulation. Reference Zhang and Zhang8
V. SUMMARY AND OUTLOOK
We presented a comprehensive analysis on the quantum thermoelectric transport properties of phosphorene using nonequilibrium Green's function technique. We investigated the details of electronic conductance, Seebeck coefficient, electronic thermal conductance, and phonon thermal conductance. For the electronic conductance, we found a strong anisotropy, with the electronic conductance in the armchair direction being significantly higher than that in the zigzag direction. The electronic conductance depends sensitively on the chemical potential and it increases rapidly when the chemical potential reaches the edges of band gap. However, in a strong contrast to the electronic thermal conductance, the Seebeck coefficients are precisely isotropic. Furthermore, the Seebeck coefficients are larger when the chemical potential is inside the band gap and it vanishes when the chemical potential goes into either conduction band or valence band. The electronic thermal conductance shows a similar behavior as the electronic conductance as the energies are carried by electrons. It can surpass the phonon thermal conductance when the chemical potential is large. However, when the chemical potential is inside the band gap, the phonon thermal conductance dominates the thermal conductance.
We further checked the validity of the WF law in phosphorene. It was found that the WF law is valid when the chemical potential is inside either conduction band or valence band. However, the conductance ratio σe/(G e T) becomes chemical-potential-dependent when it reaches VBM or CBM. The conductance ratio turns to be around half smaller when the chemical potential is inside the band gap. This violation of the WF law may be used to boost thermoelectric efficiency.
The thermoelectric performance was studied by examining the figure of merit ZT. The value of ZT shows two peaks located near the VBM and CBM. These two peaks are overall effects of the reduction of Seebeck coefficients and increase of electronic (thermal) conductance when the chemical potential is increased. The maximum ZT is 0.2 at room temperature and reaches 0.5 at 600 K. The ZT in the armchair direction is significantly higher than that in the zigzag direction.
Our work also suggests possible ways to enhance the energy conversion efficiency of phosphorene-based thermoelectric devices. For example, for lightly doped phosphorene (small chemical potential), reducing the phonon thermal conduction can effectively boost the thermoelectric performance. For heavily doped phosphorene (high chemical potential), however, the electronic structure engineering may be effective in boosting the Seebeck effect.
ACKNOWLEDGMENTS
This work was supported in part by a grant from the Science and Engineering Research Council (152-70-00017). The authors gratefully acknowledge the financial support from the Agency for Science, Technology and Research (A*STAR), Singapore and the use of computing resources at the A*STAR Computational Resource Center, Singapore.

