Introduction

Thermoelectricity has become fascinating prospect for use in generating green energy especially from the processes of waste heat [1,2,3]. The narrow-band gap (0.3 eV) semiconductor such as Bi2Se3 has been recognized as an outstanding thermoelectric material [4,5,6]. Thermoelectric material’s unique properties is that it can convert waste heat to electrical energy directly [7, 8]. It has contributed to the development of many applications such as waste heat recovery, humidity sensors, solid-state generators and radiation sensors [9, 10]. The most attractive application is a sensor that use for internet of things (IoT). There is considerable demand to power these sensors without replacing the batteries. Thermoelectric device is favorable candidate for wireless sensor networks (WSNs) in IoT due to its capability such as no moving part, no maintenance requirement, high reliability, small size and flexible [11]. Thermoelectric material has advantages such as reliable conversion, scalable and compact compared with other technologies of energy conversion and can be the solution to the environmental crisis because no energy waste produced during operation [12,13,14]. In the several crucial markets, thermoelectric devices are in mass production for heating, temperature-control and cooling applications. For future inquiry, if ZT performance of thermoelectric devices has reached 2, commercial solid-state heating, air-cooling systems and ventilating would become available for public to use at home and office [15].

Thermoelectric performance of thermoelectric material is evaluated by the equation \(ZT={S}^{2}\sigma T/\kappa\), where S is Seebeck Coefficient or thermopower, \(\sigma\) is electrical conductivity, T is temperature, \(\kappa\) is thermal conductivity and ZT is dimensionless figure of merit. From the equation, good thermoelectric material is the one that has higher thermopower and electrical conductivity, and a lower thermal conductivity.

Experimental report by Patil et al. shows one of the highest thermopower that Bi2Se3 has ever attained which is 556 μV/K. Due to the lower electrical conductivity (1.877 S/cm), the performance of Bi2Se3 is low as 0.033 [16]. The experimental report by Kulsi et al. shows the same trend which Bi2Se3 has low performance (0.02) due to the lower electrical conductivity (5.200 S/cm) [17]. Many modifications have been done to increase the performance such as doping for instance. Doping by copper, antimony and nickel have increased the ZT value to 0.54, 0.036 and 0.57, respectively [16,17,18]. There are also modification by nano-structuring such as nanoplates and nanosheets of Bi2Se3 that have increased the performance to 0.14 and 0.17, respectively [19, 20]. The advantage of nano-structuring is that the performance can be improved by reducing thermal conductivity without decreasing power factor greatly [13]. Even though there are many high-performance thermoelectric materials that almost reach unity 1.0 figure of merit, these materials could only operate at higher temperature [21,22,23,24]. Bi2Se3, on the other hand, specialized in operating thermoelectricity at room temperature [25, 26].

Study by Peter et al. has made enhancement through optimization of p-doped Bi2Se3 that has significant advantage over Bi2Te3 on thermoelectric properties at the range temperature from 400 to 600 K. The optimization has produced better ZT performance [27]. Study by Leo et al. made further enhancement from optimized doped Bi2Se3 for both n-type and p-type through a tensile in-plain strain and compressive uniaxial strain. The strain has increased Seebeck coefficient that lead to increase of power factor [28]. There have been many fruitful optimized doping predictions for thermoelectric materials such as done on CuFeS2, Cr-doped Bi2Te3 and ZnSb which shows that analyzing in doped configuration could predict good prediction of good thermoelectric performance [29,30,31].

In this study, we attempt to investigate thermoelectric properties of Bi2Se3 as the function of Fermi level to find its optimal doping level and performance by using first-principles study of density functional theory (DFT). According to Xu et al., Fermi level can be controlled and adjusted in topological insulator [32]. Bi2Se3 also exhibit as topological insulator besides thermoelectric material [4, 32,33,34]. Thus, thermoelectric properties are expected can be significantly tuned and optimized. The study involved the investigation on structural, electronic and thermoelectric properties of Bi2Se3. Structural properties are calculated using Cambridge Serial Total Energy Package (CASTEP) as an introduction of the unique structure of Bi2Se3 especially its Van der Waals interaction that gives influence on thermal conductivity of thermoelectric material [35]. Electronic properties that consist of band structure and density of states (DOS) were calculated using full-potential linearized augmented plane wave (FP-LAPW) that implemented in WIEN2K. The contribution of electronic properties on thermoelectric performance was analyzed. Finally, thermoelectric properties were calculated by employing Boltzmann transport equation that implemented in Boltzmann transport properties (BoltzTraP) computer code.

The investigation of enhancement mechanism of thermoelectric performance of Bi2Se3 is done through the contribution of Seebeck coefficient, electrical conductivity, thermal conductivity, DOS, carrier concentration and effective mass to the maximum value of power factor and figure of merit (ZT). We found that, the power factor and ZT have different kind of contribution to attain their maximum point. However, both contributed by n-type doping at their corresponding optimum level to gain maximum value at 300 K. This study could give insight and guideline to narrow down the necessity parameter and region of doping for experiments to improve the thermoelectric performance.

Computational details

Rhombohedral structure of Bi2Se3 with the space group \({D}_{3d}^{5}(R\stackrel{-}{3}m)\) is used [6]. The initial structural parameter was taken from the experimental study by Nakajima [36]. Structural optimization was calculated using DFT that implemented in CASTEP within Materials Studio. For optimization, on-the-fly generation (OTFG) ultrasoft pseudopotential was used to reduce the inaccuracy with respect to completely converged all-electron DFT calculations [37]. Van der Waals (vdW) was considered within Ortmann–Bechstedt–Schmidt (OBS) schemes since Bi2Se3 has vdW interaction in its structure [38, 39]. Generalize gradient approximation for solid (GGA-PBESol) is chosen as functional exchange–correlation [40]. After performed convergence test, cut-off energy 350 eV and k-points 4 × 4 × 1 are chosen. All forces on atoms were converged to less than 0.03 eV/Å, the maximum ionic displacement was set to 0.001 Å and the total stress tensor was decreased to the 0.05 GPa.

After optimization, electronic properties that consist of band structure and density of states (DOS) were calculated using full-potential linearized augmented plane wave (FP-LAPW) that implemented in WIEN2K. The atomic space in FP-LAPW consists of muffin-tin (MT) region and interstitial region. The former is the space throughout atoms occupied by spheres of radius, RMT that is known as MT spheres. The latter is the remaining space region [41]. GGA-PBE has been employed as the exchange–correlation functional [42]. Spin–orbit coupling (SOC) is not considered in the calculation. We use 500 k-points in the first Brillouin zone for the self-consistent calculation. The calculation was converged with respect to the RMTKmax = 7.0 and Gmax = 12.0.

The calculation of thermoelectric properties was executed by using semi-classical Boltzmann transport equation that implemented in BoltzTraP incorporated into WIEN2K computer code within FP-LAPW [43, 44]. The properties can be computed based on calculated band structure on the rigid band approach [44]. The thermoelectric properties consist of thermopower or Seebeck coefficient S, electrical conductivity \(\sigma\), thermal conductivity \(\kappa\) and figure of merit ZT.

Results and discussion

Thermoelectric material operates under the system that based on the equation \(ZT={S}^{2}\sigma T/\kappa\). From the equation, excellent thermoelectric material receives an advantage from higher Seebeck coefficient and electrical conductivity [45]. In addition, the material must have lower thermal conductivity. Higher Seebeck coefficient can be achieved through larger effective mass and smaller carrier concentration [46, 47]. On the contrary, to attain higher electrical conductivity required high carrier concentration and light effective mass [48, 49]. Besides that, thermal conductivity contributed by electronic thermal conductivity and lattice thermal conductivity [50]. Electronic thermal conductivity has a strong relationship with electrical conductivity based on Wiedemann–Franz law which by decreasing electronic thermal conductivity would also decrease electrical conductivity [51]. These contradictory conditions hindered the development towards a higher performance (ZT) for many years.

This work will set forth on structural properties which is the lattice parameter. Then, the result of electronic properties which are band structure and density of states was analyzed. Then, we will explore the transport properties which consist of Seebeck coefficient S, electrical conductivity \(\sigma\), thermal conductivity κ, power factor, and dimensionless figure of merit ZT. Geometrical optimization is done using GGA-PBESol functional together with the corresponding available experimental and theoretical values are listed in Table 1. The structure consists of five-atom layers known as quintuple layers along z-direction as shown in Fig. 1. The layer of the structure arranged with the sequence system of Se1-Bi-Se2-Bi-Se1 [5]. The structure has a strong bond within the quintuple layer. However, the bond is weak between quintuple layers. The weak bond is known as Van der Waals bond and this the reason we included the OBS Van der Waals scheme in the calculation [52, 53]. As can be seen, the optimized lattice parameters concur very well with the experimental works with the small difference around 0.07–0.17% for lattice a = b and 0.38–0.67% for lattice c. The calculated lattice parameter also agrees very well with the other theoretical works. Unlike local density approximation (LDA) due to the treatment limitation, our GGA-PBESol approach is closer to the available experimental values.

Table 1 Calculated lattice parameters of Bi2Se3 using GGA-PBESol approach along with available experimental and other theoretical works
Fig. 1
figure 1

Rhombohedral structure of Bi2Se3

After structural optimization was performed, electronic properties which are band structure and density of states were calculated. The calculated band structure was performed using GGA-PBE functional. The calculated band structure is shown in Fig. 2 which has a direct band gap of 0.33 eV at P position. This band gap is overestimated by 9.09% when compared with its typical band gap of 0.30 eV as shown in Table 2 [33, 61, 62]. We find our band gap is in agreement with the other experimental work by LaForge et al. [63] and overestimated by 3.08% when compared with the work of Valla et al. [64].

Fig. 2
figure 2

Band structure of Bi2Se3 along selected symmetry direction of X to Z

Table 2 The value of the band gap of Bi2Se3 using GGA-PBE approach along with available experimental and other theoretical works

To understand the origin and orbital contribution to the band structure, the density of states is plotted against the energy as shown in Fig. 3. It shows that total Se orbital is the most contributed at the valence band especially orbital Se-4p when compared with total Bi orbital. The maximum valence band wave vector P is contributed by orbital Se-4p. Orbital Se-4p and Bi-6s are partially hybridized. At the conduction band, both Se and Bi orbital contributed evenly specifically orbital Se-4p and Bi-6p. At the range of energy from 3 to 5 eV, orbital Bi-6p and Se-4p show antibonding states. Minimum conduction band at wave vector Γ and P solely contributed by orbital Bi-6p. Se-4p is fully hybridized with orbital Bi-6p.

Fig. 3
figure 3

Density of states (DOS) of Bi2Se3 against energy from – 10 to 10 eV

Next, thermoelectric properties were calculated using semi-classical Boltzmann transport equation that implemented in BoltzTraP. Under rigid band approximation, the Fermi level can go up or down based on the type of doping (p-type or n-type). However, the fundamental band structure is not allowed to change. Therefore, only one band structure is needed to compute thermoelectric properties at different doping level. This approach has been widely used and helpful for the theoretical investigation of thermoelectric properties based on band structure [44, 65, 66]. The effectiveness of the rigid band model for higher holes or electrons concentration depends the dopant substitution. In our case, Se-4p is the highest contribution in the band structure, specifically at the valence band. Therefore, substituting some of Se atoms at higher concentration would change the band structure, otherwise, the doping does not change the band structure considerably. Thus, for some cases, the rigid band model is effective large concentration [67]. The result of each thermoelectric properties was plotted against the various Fermi level from − 1.0 to 1.0 eV at different temperatures as shown in Fig. 4 to analyze of what gives enhancement to good thermoelectric properties along the Fermi level.

Fig. 4
figure 4

a Density of states, Seebeck coefficient, electrical conductivity and power factor against Fermi level. The dotted vertical line is drawn along the peak value of thermoelectric power factor of each temperature 300 K, 400 K and 500 K. b Seebeck coefficient, electrical conductivity, thermal conductivity and figure of merit (ZT), respectively, against Fermi level. The dotted vertical line is drawn along the peak value of figure of merit (ZT) of each temperature 300 K, 400 K and 500 K

Basically, as the Fermi level reaches band edges, the DOS increases as shown in Fig. 3. From the energy range 10–5 eV, the DOS is very low. The DOS started to increase when the Fermi level reaches 5 eV and band edge at minimum conduction band. Figure 4a shows that the Seebeck coefficient exhibits a high value at the band gap region. It quickly decreases as the Fermi level depart from the band gap and get into band edges. Even though Seebeck coefficient attains high value at the middle gap, for desirable thermoelectric performance, the value of Seebeck coefficient at the band edge is crucial as shown in Fig. 4a. Higher power factor is attained when the Fermi level at the band edge of valence band where Seebeck coefficient is 27.77 μV/K at 300 K, very low when compared with the other reports which has typical values 50–95 μV/K [19, 58, 59, 68,69,70]. However, the power factor attained maximum value mostly through the contribution of electrical conductivity at the Fermi level − 0.42 eV. In addition, large power factor occurred at p-type of doping (corresponding to negative value of Fermi level). The electrical conductivity contributed mainly by the orbital Bi-6s and Se-4p. Large DOS contributed to the large electrical conductivity [71, 72]. The ideal location of Fermi level to attained higher electrical conductivity would be at region − 1.2 to − 0.6 eV. Yet, the Seebeck coefficient is very weak at these regions and that is undesirable.

Besides calculating favorable power factor to attain, we also calculated figure of merit (ZT) to investigate favorable ZT as the function of Fermi level. The results of the calculation are shown in Fig. 4b. Contrary with the highest power factor that occurred at p-type of doping, the highest figure of merit occurred at n-type of doping (corresponding to positive value of Fermi level). In addition, it also occurred exactly at the band edge which is minimum conduction band. Higher figure of merit mostly contributed by the higher Seebeck coefficient with the value − 727 μV/K at 300 K and lower thermal conductivity. However, at 400 K and 500 K, the highest ZT attained occurred at middle gap and mostly contributed by higher Seebeck coefficient and thermal conductivity.

Another particular part of theoretical investigation for thermoelectric material is to discover level of optimal doping at which the power factor and figure of merit accomplish the maximum value. This investigation would deliver a direction in narrowing down doping area in experimental work. To carry out the investigation of discovering level of optimal doping, electrical conductivity with respect to relaxation time and Seebeck coefficient were calculated at wide range of doping concentrations ranging from 1020 to 1021 cm−3 for power factor and included thermal conductivity in the calculation within concentration range from 1016 to 1018 cm−3 for ZT. The result of calculated power factor with respect to the relaxation time and ZT were plotted against carrier concentrations that consist of electron concentrations and hole concentrations as shown in Fig. 5. For each case, the optimal doping was referring to the peak value and all the maximum values are tabulated in Table 3.

Fig. 5
figure 5

a Power factor against hole and electron concentration respectively for each temperature 300 K, 400 K and 500 K. b Figure of merit against hole and electron concentration, respectively, for each temperature 300 K, 400 K and 500 K

Table 3 Level of optimal doping (carrier/u.c) for both p-type and n-type along with the corresponding value of thermoelectric power factor and figure of merit

The power factor attained a maximum value of 39.487 μW/cm K2 s at 0.654 p-type doping level. The value is increased as the temperature increased. It is clear that Bi2Se3 has higher thermoelectric power factor for p-type doping than n-type doping. The optimum doping level of p-type is higher than the optimum n-type doping level which shows holes contributed more to power factor. However, n-type doping has a higher power factor at 300 K when compared with the p-type doping power factor. Most of the high power factors for p-type doping correspond to carrier concentration ranging from 8 × 1020 to 2 × 1021 cm−3. Power factor contributed by the higher charge carrier instead of thermopower. One of the main reasons power factor has lower value for n-type doping is due to the lower DOS at conduction band region which is lower charge concentration as shown in Fig. 3.

The figure of merit (ZT) attained maximum value at 300 K with n-type of doping 0.000027, very small when compared with the optimum level of doping for power factor. The ZT is decreased as the temperature and doping level increased. Besides 300 K, the ZT has a higher value from p-type doping at 400 K and 500 K when compared with n-type doping. It is clear that electrons contributed more to ZT at 300 K. Lower carrier concentration on contribution to the figure of merit indicates electrical conductivity does not contribute much to the accomplishment of higher ZT. Instead, it is lowering electronic thermal conductivity that plays the main role to achieve high ZT.

Finally, we investigate the reason for the higher power factor and figure of merit based on the band structure. Based on Fig. 2, there are two types of bands that occurred on the structure which are flat and curve bands. Flatter band indicates heavy effective mass, while the curve band indicates light effective. Effective mass was calculated as shown in Table 4.

Table 4 The effective mass that calculated based on the curvature of band structure of Bi2Se3 along the specified wave vector

We calculated effective mass at the wave vector P, L and Γ. It can be seen that, conduction band has higher effective mass when compared with the effective mass at the valence band except at the wave vector Γ which the conduction band lower effective mass by 8.78%. Heavy effective mass contributed to the higher Seebeck coefficient as shown its maximum value occurred at the conduction band. The negative value of Seebeck coefficient indicates electron carrier. Light effective mass, on the other hand, contribute to the charge mobility to improve electrical conductivity as shown higher DOS in valence band.

Conclusion

Our density functional theory calculations show that there are two ways to improve the thermoelectric performance of Bi2Se3 which are through enhancing thermoelectric power factor or figure of merit (ZT). The power factor can be enhanced by increasing electrical conductivity with the price of reducing thermopower. The contribution of increased electrical conductivity can be found through large DOS and light effective mass at the valence band of Bi2Se3. The power factor also can be improved through the optimal n-type doping at the temperature 300 K. However, for the higher temperature (400 K and 500 K), p-type is the optimal type of doping. ZT is improved through the higher thermopower and lower thermal conductivity with the price of reducing electrical conductivity. We found that high Seebeck coefficient contributed by heavy effective mass and low thermal conductivity contributed by low DOS. Furthermore, ZT also can be enhanced through the n-type of doping at 300 K and p-type doping for higher temperatures the same case as the power factor. We hope that our theoretical investigation could give help and provide a guideline to improve thermoelectric performance of Bi2Se3.