Giant thermoelectric power factor in charged ferroelectric domain walls of GeTe with Van Hove singularities

Abstract

Increasing the Seebeck coefficient S in thermoelectric materials usually drastically decreases the electrical conductivity σ, making significant enhancement of the thermoelectric power factor σS² extremelly challenging. Here we predict, using first-principles calculations, that the extraordinary properties of charged ferroelectric domain walls (DWs) in GeTe enable a five-fold increase of σS² in the DW plane compared to bulk. The key reasons for this enhancement are the confinement of free charge carriers at the DWs and Van Hove singularities in the DW electronic band structure near the Fermi level. These effects lead to an increased energy dependence of the DW electronic transport properties, resulting in more than a two-fold increase of S with respect to bulk, without considerably degrading the in-plane σ. We propose a design of a nano-thermoelectric device that utilizes the exceptional thermoelectric properties of charged ferroelectric DWs. Our findings should inspire further investigation of ferroelectric DWs as efficient thermoelectric materials.

Introduction

Thermoelectric materials can convert waste heat into electrical power or, in reverse, cool devices using electrical current. Achieving large values of the thermoelectric power factor σS² for efficient energy conversion is the critical barrier to further progress in the field of thermoelectricity^1,2,3. It is hard to realize a simultaneous increase of the Seebeck coefficient S and the electrical conductivity σ due to their interdependence ($S \sim {\left.\frac{1}{\sigma }\frac{{\text{d}}\sigma }{\text{d}}E\right|}_{E\,=\,{{E}}_{\text{F}}}$, where E is the electronic energy and E_F is the Fermi level)¹. So far, thermoelectric power factors have been enhanced using e.g., the concepts of resonant impurities⁴, energy filtering⁵, band convergence⁶, and low dimensionality^7,8. Nevertheless, the reported increases of σS² have rarely been larger than two-fold^{1,2,3,4,6,9,10,11}. Here, we will show that a large enhancement of power factor could be realized in a novel class of thermoelectric materials: ferroelectric domain walls (DWs).

Ferroelectric DWs are the interfaces between regions with differently oriented polarizations (domains) in a ferroelectric material. Polarization discontinuity at DWs can induce electric fields along domains, which in turn confine free charge carriers on these so-called charged DWs (see Fig. 1). The resulting two-dimensional (2D) electron gas can have dramatically different properties than the bulk material^12,13. Several normally insulating perovskite materials show conductive behaviour at DWs^{14,15,16,17,18,19,20}, with the electrical conductivity of DWs up to 13 orders of magnitude larger than that in the domains²¹. These discoveries have lead to a new paradigm of ferroic devices where DWs, rather than domains, are active elements in nanoelectronic circuits¹³, such as diodes²², rectifiers²³, resistive memories^24,25, and memristors²⁶. In contrast, there have been no reported studies of the Seebeck coefficient and the power factor of ferroelectric DWs. Despite over a decade of research of the DW conductive properties, the idea that their thermoelectric properties may be superior to those of bulk has not been considered.

Fig. 1: Ferroelectric domain walls.

Illustration of a ferroelectric material containing charged domain walls (DWs). Coloured arrows denote the polarization orientations. DWs are the interfaces between domains with different polarization orientations. There are head-to-head (H-H) and tail-to-tail (T-T) charged DWs depending of the polarization orientation with respect to the DW plane. The polarization changes at these DWs form bound charge (labelled by plus and minus). The direction of the electric field induced by the bound charge is given with the black arrow. The electric field localizes free charge carriers on the DWs (green—holes, yellow—electrons). In the 39° (11$\bar{1}$) DWs considered here, the angle between the polarizations in neighbouring domains is 39°, and the DWs are in the (11$\bar{1}$) plane.

Here, we investigate the thermoelectric transport properties of charged ferroelectric domain walls in germanium telluride (GeTe). GeTe is an exceptionally good thermoelectric material that also exhibits ferroelectric nature below 700 K^{27,28,29,30,31,32}. Ferroelectric DWs have been observed in GeTe to form a heringbone structure^{33,34,35,36,37,38}. The existence of DW structures in GeTe has been associated with an increase in the thermoelectric figure of merit, but only through a decrease of the lattice thermal conductivity^39,40 and carrier concentration control^41,42. On the other hand, ferroelectric DWs have not been discussed in the context of increasing the power factor in thermoelectric materials.

In this work, we report a prediction that the thermoelectric power factor in the plane of ferroelectric DWs in GeTe can be up to five times larger than that of bulk GeTe. This increase is enabled by more than a two-fold increase of the Seebeck coefficient in the DW plane, which does not significantly deteriorate the electrical conductivity. These enhancements are a direct result of an increased energy dependence of the local density of states (LDOS), group velocities, and relaxation times of the electronic states near the Fermi level that are localized on these DWs and exhibit Van Hove singularities in their band structure. These results demonstrate that domain wall engineering in GeTe, and potentially ferroelectric oxides, could create new opportunities for enhancing thermoelectric power factor at the nanoscale.

Results and discussion

GeTe domain walls

Recently we have identified different types of domain walls that can occur in GeTe and analysed their structural and thermal transport properties³⁹. We note that there is a high energetic cost of creating charged DWs in ferroelectric oxides and most observed DWs are neutral. However, we have predicted that the differences between the DW energies of charged and neutral DWs in GeTe can be smaller than those in perovskites³⁹ due to the smaller band gap of GeTe. This indicates that the formation of charged DWs might be more common in GeTe than in perovskites.

Here we will focus on the electronic properties of charged 39° (11$\bar{1}$) DWs. The angle between the polarization vectors in neighbouring domains in these DWs is close to 39°, and the DW plane lies in the (11$\bar{1}$) crystallographic plane (see Fig. 1). Depending of the polarization orientation with respect to the DW boundary, we have two types of 39° (11$\bar{1}$) DWs: tail-to-tail (T-T) and head-to-head (H-H). We construct GeTe supercells containing both 39° T-T and H-H DWs as explained in the ref. ³⁹. We carry out electronic structure calculations on these supercells using density functional theory, as detailed in the “Methods” section.

Charge confinement

Using Gauss’s law, it can be shown that discontinuity of polarization at a domain wall will lead to the appearance of bound charge on this DW⁴³:

$$\rho =({{\bf{P}}}_{1}-{{\bf{P}}}_{2})\cdot {\bf{n}},$$

(1)

where ρ is the bound charge interface density, P_1,2 are the polarizations in neighbouring domains, and n is the vector normal to the DW plane. Bound charge will create an electric field along domains, as shown in Figs. 1 and 2a, which will be screened by free charge carriers to reduce the electrostatic energy. This will lead to the accumulation of free charge carriers at DWs, see Fig. 1. The charge accumulation on the DWs is evident from Fig. 2b, which shows that the LDOS of both 39° (11$\bar{1}$) T-T and H-H DWs is finite at the Fermi level.

Fig. 2: Confinement of free charge carriers at domain walls.

a Averaged total potential of GeTe containing 39^∘ (11$\bar{1}$) tail-to-tail (T-T) and head-to-head (H-H) domain walls (DWs), represented by the green and yellow shaded regions respectively, as a function of the coordinate normal to the DWs. The potential drop along domains arises from polarization discontinuity at the DWs. b Electronic local density of states in the middle of domains and at the domain walls. The local density of states (LDOS) of domain walls is obtained by summing the LDOS of the atoms in six primitive unit cells near the domain wall plane. E_F is the Fermi energy of the system at 0 K. For the H-H and T-T DW, the Fermi level lies above and below the forbidden gap, respectively.

The accumulation of free charge carriers at DWs can be understood as a self-doping phenomenon. Indeed, Fig. 2b shows that the T-T DW is p-doped and the H-H DW is n-doped. This is in accordance with Fig. 2a, which depicts a drop in the total potential at the H-H DW that attracts electrons, leaving holes on the T-T DW. GeTe DWs are thus intrinsically conductive, similarly to DWs in perovskites^{14,15,16,17,18,19,44}. Self-doping of charged DWs has already been observed in first principles calculations of BaTiO₃ DWs^45,46. It is important to notice that the band gap does not close at these DWs, and they effectively behave as doped semiconductors. However, the value of the band gap does change and becomes somewhat smaller at the H-H DW compared to bulk.

To understand the nature of the electronic states forming on 39° DWs, we calculated the electronic band structure of the GeTe supercell containing 39° H-H and T-T DWs and projected it on the DW plane, see Fig. 3a. Blue lines correspond to the eigenvalues of the states with the reciprocal lattice vector normal to the DW plane (k₃) set to zero. Dispersion of these states in the direction normal to the DW plane is represented by broadening given in red, which illustrates how these eigenvalues vary with k₃. The bands close to the Fermi level E_F do not have any dispersion along the direction normal to the DW plane (i.e., they have no broadening in that direction). Consequently, these states are localized and behave like a 2D electron gas, as we would expect from the potential profile in Fig. 2a. Contrary to that, the bands further from E_F have large broadening and are not localized.

Fig. 3: Localization of electronic states near the Fermi level on domain walls.

a Electronic band structure of GeTe containing 39° (11$\bar{1}$) tail-to-tail (T-T) and head-to-head (H-H) domain walls (DWs) projected on the DW plane. The eigenvalues of the states with the reciprocal lattice vector normal to the DW plane (k₃) set to zero are shown in blue, while the dispersion of these states as k₃ changes is shown in red. The grey area in the graph corresponds to the bulk GeTe band structure projected to the domain wall plane. The horizontal line represents the Fermi level at 0 K. The absence of dispersion for the states near the Fermi level confirms their localization due to electric fields along domains. The side plot shows the local density of states of the T-T DW. b Charge density associated with the states labelled as 1. and 2. in a, as a function of the coordinate normal to the DWs, which are indicated by black vertical lines. These states are clearly localized on the DWs.

To confirm that the electronic states near the Fermi level are localized on the DWs, we computed the average charge densities associated with two representative states close to E_F on planes parallel to the DWs as a function of the coordinate perpendicular to the DWs, see Fig. 3b. These states are indeed confined to the DWs and do not extend into the domains. The valence band states are localized on the T-T DW, while the conduction band states are localized on the H-H DW. Localization of these states on the DWs is induced by the electric field generated by the bound charge on the DWs. In contrast, the absence of electric fields along domains for neutral DWs (e.g., 141° (11$\bar{1}$) head-to-tail DW) produces delocalized states near the Fermi level (see Supplementary Note 1).

2D electronic band structure and Van Hove singularities

The electronic states localized at 39° (11$\bar{1}$) T-T DW have a sharp peak in the LDOS near the top of the valence band, as opposed to the step characteristic for 2D systems (see Fig. 2b). Such distortion of the LDOS resembles those stemming from the resonant impurity levels in PbTe-doped with Tl⁴ and GeTe-doped with In⁴⁷, which lead to an experimentally observed increase of the Seebeck coefficient. We will show that the LDOS peak of the T-T DWs also enhances the Seebeck coefficient, even though its physical origin is very different from that of resonant impurities. The large increase of the LDOS near the Fermi level also occurs in 180° (11$\bar{1}$) T-T and 180° (111) DWs and may be a general feature for these systems (see Supplementary Note 2).

To uncover the origin of the LDOS peak for 39° (11$\bar{1}$) T-T DW, we plot the LDOS of this DW together with the interface projected electronic band structure in Fig. 3a. We can see that the LDOS peak occurs at the similar energy as the valence band maximum on the X − Γ line just below the Fermi level. The 2D dispersion of the electronic states of this particular band is given in Fig. 4. This figure shows that the states on the X − Γ line shown by blue points are actually saddle points in the electronic band structure: the energies of the electronic states decrease along the X − Γ direction away from the saddle points, and increase in the perpendicular direction. These are well-known Van Hove singularities that give a logarithmic divergence for the density of states in 2D systems⁴⁸, which explain the LDOS peak in our calculations.

Fig. 4: 2D electronic band structure and Van Hove singularities of the tail-to-tail domain wall.

Energy of the top valence band for 39° (11$\bar{1}$) tail-to-tail domain wall (in eV) in the interface Brillouin zone. Black arrows represent the directions of the in-plane reciprocal lattice vectors. The blue arrow corresponds to the projection of the polarization change onto the domain wall plane. Van Hove singularities are indicated by blue points.

Similar Van Hove singularities (i.e., saddle points) are also present in the electronic band structure of bulk GeTe (see Supplementary Note 3). However, they are further from E_F than those of the T-T DW and do not lead to the logarithmic divergence of the DOS in three dimensions (3D)⁴⁸. In the T-T DW, the electric field pushes the Van Hove singularities much closer to the Fermi level, thus making them contribute to electronic transport. Furthermore, the 2D nature of DWs makes the effect of the Van Hove singularities on the LDOS much larger than in bulk.

We note that the T-T DW has a very anisotropic valence band surface (see Fig. 4), which is beneficial for thermoelectric transport properties. The direction parallel to the projection of the polarization change along the DW given by the blue arrow (X − Γ) has much lower group velocities compared to the perpendicular direction in the DW plane (Γ − K). This is a consequence of the layered structure of bulk GeTe, because the direction with low group velocities is perpendicular to van der Waals gaps in GeTe, making the states along that direction more confined (see Supplementary Note 4). The 1D cigar like states in the X − Γ direction (red color in Fig. 4) give an additional contribution to the peak of the density of states, which should be beneficial to the Seebeck coefficient as we show later. The Γ − K direction, on the other hand, has higher group velocities, which should result in higher electrical conductivity values for this direction.

Thermoelectric transport properties

We have shown that the electronic states localize at either the T-T or the H-H domain wall, see Fig. 3b. Localization prevents interaction between states localized at different domain walls and allows us to consider T-T and H-H domain walls as separate systems. The Fermi level in the stochiometric GeTe lies within these localized states and far from the domain, bulk-like, states (see Fig. 3a). This means that charge transport in this structure will occur through domain walls only, with DWs acting as independent channels. This can be accomplished by an appropriate placement of electrodes as we show later (see Fig. 7). To calculate the transport properties of the T-T DW (the p-type channel in this system), we consider only the states localized at the T-T DW in the Boltzmann transport equation (BTE) (see Eq. (6) in “Methods” section).

We compute the Seebeck coefficient of 39° T-T DW in the X − Γ and Γ − K directions in the DW plane and in the direction perpendicular to the DW plane. S is plotted in Fig. 5 as a function of the extrinsic free charge density on the DWs. Our calculations show at least a two-fold enhancement of the Seebeck coefficient of 39° T-T DW with respect to bulk GeTe in all directions. To compare the free charge densities of the DW (2D) and bulk (3D), we obtain the volume charge concentration of the T-T DWs by assuming that the separation between the T-T and H-H DWs corresponds to the one in our DFT calculation (5.6 nm) (see “Methods” section). The intrinsic DW charge density due to charge transfer between the H-H and T-T DWs and the related volume charge concentration are indicated by the black vertical line in Fig. 5.

Fig. 5: Seebeck coefficient of the tail-to-tail domain wall.

Calculated Seebeck coefficient for p-type bulk GeTe and 39° (11$\bar{1}$) tail-to-tail (T-T) domain walls (DWs) for various directions at 300 K. Experimental data for bulk GeTe from refs. ^62,63 are shown by red squares and stars, respectively. The top x-axis shows the DW interface free charge density, while the bottom x-axis corresponds to the volume charge concentration for the DW separation of 5.6 nm. Black vertical line shows the intrinsic doping concentration of the T-T DW for this DW separation due to charge transfer between the head-to-head and T-T DWs.

An increase of the domain size (the separation between the domain walls) will lead to an increase in the DW free charge carrier density. We use the domain size for which the electronic states at the head-to-head and tail-to-tail domain walls are decoupled and their band structures are converged. The increased free charge carrier density will more completely screen the bound charge caused by the polarization change at the DWs, leading to a decrease in the electric field along the domains. However, the increase of the domain size does not have a significant effect on the electronic band structure of domain walls (see Supplementary Note 7). Consequently, we expect that the major effect of the increased domain size would be a shift of the black vertical line in Fig. 5 to higher free charge carrier densities.

Our BTE results for 39° T-T DW show that the Seebeck coefficient is the largest in the direction perpendicular to the DW plane (Fig. 5). The enhancement of S is smaller in the in-plane directions, where the electronic states are more dispersive and more conductive. We explain these trends for the Seebeck coefficient using the Mott expression⁴⁹:

$$S=\frac{{\pi }^{2}{k}_{{\rm{B}}}^{2}T}{3e}{\left.\left(\frac{N^{\prime} (E)}{N(E)}+\frac{\langle{v}^{2}(E)\rangle^{\prime} }{\langle{v}^{2}(E)\rangle}+\frac{\tau ^{\prime} (E)}{\tau (E)}\right)\right|}_{E\,=\,{{E}}_{\text{F}}},$$

(2)

where k_B is the Boltzmann constant, T is the temperature, e is the electron charge, N(E) is the density of states, $\langle{v}^{2}(E)\rangle$ is the squared group velocities average, and τ(E) is the relaxation time. Primed quantities represent the first derivatives with respect to energy.

To disentangle the individual contributions of the three terms in Eq. (2), we calculate the Seebeck coefficient using three levels of approximation (see Supplementary Note 5). Firstly, we set group velocities and relaxation times to constant values, and account only for the energy dependence of DOS in the S calculation. This approximation gives an isotropic Seebeck coefficient for the DW that is enhanced compared to bulk (Supplementary Note 5). This increase comes from the LDOS peak near the Fermi level at the DW, which arises due to Van Hove singularities. Secondly, including the energy dependence of group velocities in addition to that of DOS in the S calculation gives rise to an anisotropic Seebeck coefficient for the DW. S increase is larger in the directions of lower group velocities and stronger localization (perpendicular to the DW plane and X − Γ) compared to the direction with higher group velocities (Γ − K) (Supplementary Note 5). Finally, including the energy and momentum dependent relaxation times in the S calculation further increases the Seebeck coefficient of the DW in all three directions with respect to bulk, see Fig. 5. Therefore, the increased energy dependence of the DOS, group velocities and relaxation times near the Fermi level induced by the presence of 2D Van Hove singularities is the key to the S enhancement in the T-T DW.

Next we calculate the electrical conductivity and the power factor of 39° T-T DWs (see Fig. 6). Even though there is a local increase of the free charge concentration at the DWs compared to bulk, the electrical conductivity in the in-plane directions (Γ − K and X − Γ) is somewhat reduced at the DW for some extrinsic doping concentrations. The reason for this decrease are the increased electron scattering rates due to Van Hove singularities. Nevertheless, there is more than a five-fold increase of the in-plane power factor across the whole doping range due to the large S enhancement. Most importantly, the maximum of the power factor for the T-T DW in the Γ − K and X − Γ directions is five times larger than that for bulk. Consequently, the combination of the increased energy dependence of the transport quantities that enhance S and the increased local charge concentration that prevents σ reduction lead to a large increase of the in-plane power factor in the T-T DWs.

Fig. 6: Thermoelectric transport properties of the tail-to-tail domain wall.

Calculated a electrical conductivity and b power factor of p-type bulk GeTe and 39° (11$\bar{1}$) tail-to-tail (T-T) domain walls (DWs) for various directions at 300 K. The top x-axis shows the DW interface free charge density, while the bottom x-axis corresponds to the volume charge concentration for the DW separation of 5.6 nm. Black vertical line shows the intrinsic doping concentration of the T-T DW for this DW separation due to charge transfer between the head-to-head and T-T DWs.

In the direction perpendicular to the T-T DW, the electrical conductivity and the power factor of the DW are both significantly decreased with respect to bulk. This occurs because the electronic states near the Fermi level are localized on the DW, and have low group velocities in the direction perpendicular to the DW. As a result, the thermoelectric transport properties of DWs in the out-of-plane direction are inferior to those in the DW plane.

We predict that the Seebeck coefficient and the power factor of 39° (11$\bar{1}$) H-H DW in the Γ − K direction are also enhanced with respect to bulk (see Supplementary Note 6). The thermoelectric transport properties of the H-H DW (n-type) are not as large as those of the T-T DW (p-type). This is due to the fact that Van Hove singularities in the electronic band structure of the H-H DW do not lead to as prominent peaks in the DOS and are not as near E_F as in the T-T DW (see Fig. 2b). Other types of charged DWs in GeTe have more pronounced DOS peaks for both p-type and n-type DWs e.g., 180° (111) DWs (Supplementary Note 2), and could exhibit comparable thermoelectric properties to those of 39° (11$\bar{1}$) T-T DW.

Nano-thermoelectric device

In GeTe samples with charged DWs, both n-type and p-type DWs will contribute to the Seebeck coefficient, thereby suppressing its overall enhancement. To harness the predicted outstanding thermoelectric properties of ferroelectric DWs, we propose a nano-thermoelectric device consisting of H-H and T-T DWs acting as n-type and p-type legs of the thermoelectric couple⁵⁰, separated by electrically insulating domains (see Fig. 7). The temperature gradient would be applied along DWs, rather than perpendicular to them. The proposed device could be used for nanoscale energy harvesting or cooling. Since DWs already exist in GeTe forming an ordered herringbone structure^{33,34,35,36,37,38}, it should be possible to realize such device with an appropriate placing of contact electrodes. Another important advantage of using DWs as nano-thermoelectric couples is that we can write or erase them by applying an electric field (or light)³⁴. For example, we could grow a pristine GeTe film without domains and then pole it to form desired types of DWs that will act as active device elements³⁵. Manipulating the domain length would also allow us to tune the doping concentrations at the domain walls⁵¹, in addition to vacancies and dopants.

Fig. 7: Proposed design of the nano-thermoelectric device based on domain walls.

White regions represent domains with different polarization orientations, acting as electrical insulators. Head-to-head and tail-to-tail domain walls, shown by red and blue regions, act as an n-type and p-type leg of the thermoelectric device, respectively. Electrodes and insulating parts are shown in black and green, respectively. Heating up the top electrodes causes electron and hole diffusion in n-type and p-type domain walls, respectively, which produces electrical current.

Discussion

In summary, we have shown that the thermoelectric power factor of charged ferroelectric domain walls in GeTe can be enhanced up to a factor of five compared to the bulk values. The bound charge arising from polarization discontinuity on these interfaces causes localization of the states close to the Fermi level on the domain walls. The local density of these two-dimensional states diverges logarithmically near the Fermi level due to the presence of saddle points in their electronic band structure. This significantly increases the Seebeck coefficient in the domain wall plane without considerably reducing the electrical conductivity. Our results clearly show the potential of ferroelectric domain walls for nanoscale thermoelectric applications. Our predictions of exceptional thermoelectric properties of GeTe domain walls could be tested experimentally using scanning thermoelectric microscopy^52,53. To further investigate the thermoelectric properties of GeTe domain walls, it will be important to understand their interaction with intrinsic vacancies⁴⁶ and extrinsic dopants.

Methods

Structure of the domain walls

The detailed description of the construction of supercells containing both 39° ($11\bar{1}$) T-T and H-H DWs is in ref. ³⁹. We construct the supercell containing both head-to-head and tail-to-tail domain walls from the rhombohedral primitive cell of bulk GeTe, which is defined by the lattice vectors:

$$\begin{array}{l}{{\bf{r}}}_{1}=a(b,0,c),\\ {{\bf{r}}}_{2}=a\left(-\frac{b}{2},\frac{b\sqrt{3}}{2},c\right),\\ {{\bf{r}}}_{3}^{\prime}=a\left(-\frac{b}{2},-\frac{b\sqrt{3}}{2},c\right).\end{array}$$

(3)

Here a is the lattice constant, $b=\sqrt{2(1-\cos \theta )/3}$, $c=\sqrt{(1+2\cos \theta )/3}$, and θ is the angle between the primitive lattice vectors. The atomic positions are given as: Ge (0.0,0.0,0.0) and Te (0.5 + τ, 0.5 + τ, 0.5 + τ) in reduced coordinates. The construction of 39° domain walls involves a mirror reflection of the domain across the domain wall plane, here taken to be ($11\bar{1}$) in the conventional pseudocubic cell. The in-plane directions are the vectors that lie in the ($11\bar{1}$) plane, while the out-of-plane direction is the vector perpendicular to this plane.

The relaxation of the domain structure is done in two steps. First we relax Te atoms only, and after reaching a minimum of energy we relax the total structure and all Ge and Te atoms. With this method we reduce the forces on atoms below 10⁻⁶ eV/Å even at the domain walls. Two of the lattice vectors that lie in the domain wall plane are roughly the same as in the unit cell (they change a little during the structural relaxation), and the third lattice vector is perpendicular to them and the domain wall plane. This means that one of the reciprocal lattice vectors (k₃) is colinear with the third lattice vector. The orientation of the other two reciprocal lattice vectors (which are in the domain wall plane) is given in Fig. 4.

Ab initio calculations

We perform electronic structure calculations on these supercells using density functional theory (DFT) and the ABINIT code⁵⁴, employing the Perdew–Burke–Ernzerhof parametrization⁵⁵ for the exchange-correlation potential, and Hartwigsen–Goedecker–Hutter pseudopotentials⁵⁶. We use the energy cutoff of 16 Ha for plane waves and a 12 × 12 × 1k-point grid for the Brillouin zone sampling of the electronic states.

Transport properties calculations

We compute the thermoelectric transport properties of 39° (11$\bar{1}$) DWs and bulk GeTe using the Boltzmann transport equation (BTE). BTE combined with DFT has proven to be an exceptional method for reproducing and predicting transport properties of materials completely from first principles^57,58. In this formalism, the electrical conductivity and the Seebeck coefficient are given as:

$${\sigma }^{{\rm{ij}}}={L}_{0}^{{\rm{ij}}},$$

(4)

$${S}^{{\rm{ij}}}=-{L}_{1}^{{\rm{ij}}}/(eT{L}_{0}^{{\rm{ij}}}),$$

(5)

where σ^ij and S^ij are the ij elements of the electrical conductivity and Seebeck tensors, respectively, e is the electron charge, T is the temperature and:

$${L}_{{\rm{n}}}^{{\rm{ij}}}=\frac{2{e}^{2}}{N{V}_{{\rm{DW}}}}\sum _{{\bf{k}}}\left(-\frac{\partial f}{\partial E}\right){v}_{{\bf{k}}}^{{\rm{i}}}{v}_{{\bf{k}}}^{{\rm{j}}}{\tau }_{{\bf{k}}}{({E}_{{\bf{k}}}-{{E}}_{\text{F}})}^{n}.$$

(6)

Here N is the number of k-points in the reciprocal space, V_DW is the volume of the domain wall, f is the equilibrium Fermi-Dirac distribution function, τ_k is the transport relaxation time of the state k, and ${v}_{{\bf{k}}}^{{\mathrm{i}}}$ is the group velocity of the state k with energy E_k along the i-th direction. We calculate the group velocities using the finite difference method.

An accurate calculation of the scattering rates of the electronic states localized on DWs is very challenging because of the complicated electronic band structure, and the enormous computational cost involved in extracting electron–phonon matrix elements for large supercells containing DWs. Therefore, we have developed a model for the energy and momentum dependent scattering rates of the electronic states localized on DWs, whose all parameters can be obtained from DFT calculations. In our model, we account for scattering due to electron–phonon coupling. Scattering due to ionized impurities is neglected due to the large static dielectric constant of GeTe⁵⁷.

The transport relaxation time of the electronic state with the band number n and the wave vector k at the domain wall due to electron–phonon interaction can be obtained from first order perturbation theory as⁵⁹:

$$\begin{array}{l}\frac{1}{{\tau }_{{\rm{n}},{\bf{k}}}}=\sum\limits _{{\rm{m}},{{\bf{k}}}^{\prime}}\sum\limits _{\lambda ,{\bf{q}}}\frac{2\pi }{\hslash }{\left|\left\langle {\psi }_{{\rm{m}},{{\bf{k}}}^{\prime}}\left|{H}_{\lambda ,{\bf{q}}}\right|{\psi }_{{\rm{n}},{\bf{k}}}\right\rangle \right|}^{2}\left({n}_{\lambda ,{\bf{q}}}+\frac{1}{2}\mp \frac{1}{2}\right)\\ \times \delta ({E}_{{\rm{m}},{{\bf{k}}}^{\prime}}-{E}_{{\rm{n}},{\bf{k}}}\pm \hslash {\omega }_{\lambda ,{\bf{q}}})(1-\cos \theta ),\end{array}$$

(7)

where $\hbar$ is the reduced Planck constant, ψ_n,k is the wave function of the electronic state $\left|n{\bf{k}}\right.\rangle$ with energy E_n,k, n_λ,q, and ω_λ,q are the equilibrium Bose–Einstein distribution and the frequency of phonons for the branch λ and the crystal momentum q, and H_λ,q is the electron–phonon perturbation potential. The term $1-\cos \theta =1-\frac{{{\bf{v}}}_{{{n}},{\bf{k}}}\cdot {{\bf{v}}}_{{{m}},{{\bf{k}}}^{\prime}}}{| {{\bf{v}}}_{{{n}},{\bf{k}}}| | {{\bf{v}}}_{{{m}},{{\bf{k}}}^{\prime}}| },$ represents the scattering angle which accounts for the momentum relaxation rate. We neglect phonon frequencies in the energy conservation condition, but they are partially accounted via Gaussian broadening (standard deviation of 5 meV). We use ${n}_{\lambda ,{\bf{q}}}+\frac{1}{2}\mp \frac{1}{2}\approx \frac{{k}_{{\rm{B}}}T}{\hslash {\omega }_{\lambda ,{\bf{q}}}}$, where k_B is the Boltzmann constant. This is a good approximation at 300 K since the Debye temperature of GeTe is ≈180 K.

We assume that the dominant electron–phonon interactions are non-polar, since we consider systems with high doping concentrations where polar interactions are largely screened by free carriers, as shown in first principles calculations for bulk GeTe⁶⁰. We thus consider the electron–phonon perturbation of the form ${H}_{{\rm{ac}},{\bf{q}}}=\sqrt{\frac{\hslash }{2M^{\prime} {\omega }_{{\bf{q}}}}}{U}_{1}q{{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {\bf{r}}},$ for long-wavelength acoustic phonons with linear dispersions, and ${H}_{{\rm{opt}},{\bf{q}}}=\sqrt{\frac{\hslash }{2\bar{M}\omega }}{U}_{2}{{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {\bf{r}}}$ for long-wavelength optical dispersionless phonons^57,59 ($M^{\prime}$ is the total mass and $\bar{M}$ is the reduced mass of the unit cell). Considering the dispersion relations for acoustic and optical modes, we can further write:

$${H}_{{\rm{ac}},{\bf{q}}}\sqrt{\frac{2\pi {k}_{{\rm{B}}}T}{{\hslash }^{2}{\omega }_{\lambda ,{\bf{q}}}}}=\left(\sqrt{\frac{\pi {k}_{{\rm{B}}}T}{\hslash M^{\prime} }}\frac{1}{v}{U}_{1}\right){{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {\bf{r}}}={\bar{U}}_{1}{{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {\bf{r}}},$$

(8)

$${H}_{{\rm{opt}},{\bf{q}}}\sqrt{\frac{2\pi {k}_{{\rm{B}}}T}{{\hslash }^{2}{\omega }_{\lambda ,{\bf{q}}}}}=\left(\sqrt{\frac{\pi {k}_{{\rm{B}}}T}{\hslash \bar{M}}}\frac{1}{\omega }{U}_{2}\right){{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {\bf{r}}}={\bar{U}}_{2}{{\rm{e}}}^{-{\rm{i}}{\bf{q}}\cdot {\bf{r}}}.$$

(9)

Here v is the speed of sound and ω is the characteristic frequency of optical modes. The final expression for the total relaxation time including both scattering mechanisms can be written as:

$$\frac{1}{{\tau }_{{\rm{n}},{\bf{k}}}}={U}^{2}\sum _{{\rm{m}},{{\bf{k}}}^{\prime}}{I}_{{\rm{n,m}},{\bf{k}},{{\bf{k}}}^{\prime}}^{2}\ \delta ({E}_{{\rm{m}},{{\bf{k}}}^{\prime}}-{E}_{{\rm{n}},{\bf{k}}})(1-\cos \theta ),$$

(10)

where the coefficient ${U}^{2}={\bar{U}}_{1}^{2}+{\bar{U}}_{2}^{2}$ contains all the physical constants and material parameters described above and temperature (kept constant at 300 K in all our calculations). k corresponds to the wave vector component in the DW plane, and the momentum is conserved only in the DW plane⁶¹. The quantity ${I}_{{\rm{n,m}},{\bf{k}},{{\bf{k}}}^{\prime}}$ is defined as⁶¹:

$${I}_{{\rm{n,m}},{\bf{k}},{{\bf{k}}}^{\prime}}=\frac{1}{{N}_{\perp }V}\sum _{{{\bf{q}}}_{\perp }}\int {u}_{{\rm{m}},{{\bf{k}}}^{\prime}}^{* }\ {{\rm{e}}}^{-{\rm{i}}{{\bf{q}}}_{\perp }\cdot {\bf{r}}}\ {u}_{{\rm{n}},{\bf{k}}}\ \,\text{d}{\bf{r}},$$

(11)

where the integration over dr is within the unit cell, q_⊥ is the wave vector component perpendicular to the DW and N_⊥ is the number of q_⊥ vectors. u_n,k is the periodic part of the Bloch wave function of the state $\left|n{\bf{k}}\right.\rangle$, and V is the unit cell volume. Here, we assume that the phonon band structure and the strength of electron–phonon interaction do not change at the DW with respect to bulk.

We can derive the expressions for the scattering rates in bulk GeTe in a similar manner. In this case, the total momentum is conserved in Eq. (10), and ${I}_{{\rm{n,m}},{\bf{k}},{{\bf{k}}}^{\prime}}$ is the overlap of the wave functions of the states $\left|n{\bf{k}}\right.\rangle$ and $\left|m{{\bf{k}}}^{\prime}\right.\rangle$. To compare our conductivity and power factor values for bulk GeTe to experiments, we fit the value of U to the first principles calculations of the bulk electronic conductivity for p-type GeTe at 300 K⁶⁰.

In the case of bulk GeTe, the thermoelectric transport properties are obtained along the trigonal [111] axis and two directions perpendicular to it, and using their average. In the DW case, we calculate the transport quantities along the Γ − X and Γ − K directions in the DW plane and in the direction perpendicular to the DW plane. We have performed convergence studies of the transport properties with respect to the number of k-points. For the bulk calculation, we use a 70 × 70 × 70 grid, while for the DW calculation we use a uniform 2D grid with 1806 k-points in the first Brillouin zone. To make the calculation of the wave function overlaps more computationally tractable, we reduce the energy cutoff for plane waves from 16 to 8 Ha when calculating wave functions. For small k-point grids this method gives a very small error compared to the calculation with the converged energy cutoff of 16 Ha.

We compute the free charge concentration at the T-T domain wall as follows. We assume that the separation between the T-T and H-H DWs is equal to the domain size in our DFT calculation (5.6 nm). We calculate the local density of states (LDOS) of the T-T DW in the supercell containing both H-H and T-T DWs by summing the LDOS contributions from the atoms belonging to the half of the supercell that contains only the T-T DW. We normalize the LDOS of these atoms so that the LDOS integral up to the top of the valence band for each Ge atom is 4 and for each Te atom is 6. We then obtain the surface and volume charge densities of the T-T DW as:

$${p}_{{\rm{S,V}}}=\frac{1}{S,V}{\int\nolimits_{\!-\infty }^{{E}_{\text{V}}}}\left(1-f(E,{E}_{\text{F}},T)\right)N(E){\rm{d}}E,$$

(12)

where S is the DW plane surface, V is the volume of the half of the supercell, E_V is the top of the valence band, and N(E) is the LDOS of the T-T DW structure. We note that the DW volume carrier concentration (p_V) depends on the size of the domain, unlike the surface charge density (p_S).