Eric Pante
email: epante@univ-lr.fr

Vanessa Becquet

Amélia Viricel

Pascale Garcia

Investigation of the molecular signatures of selection on ATP synthase genes in the marine bivalve Limecola balthica

Keywords: Molecular evolution, local adaptation, genetic incompatibilities, selection, ATP synthase, Macoma balthica

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

It is now well-accepted that marine organisms can show local adaptation to environmental conditions despite their high dispersal potential, which was long thought to prevent adaptive differentiation (e.g. review of [START_REF] Sanford | Local Adaptation in Marine Invertebrates[END_REF]. In fact, recent research has shown that adaptations can emerge and be maintained in spite of gene flow [START_REF] Tigano | Genomics of local adaptation with gene flow[END_REF]. For instance, focusing on the high-gene flow sea urchin Strongylocentrotus purpuratus, [START_REF] Pespeni | Signals of selection in outlier loci in a widely dispersing species across an environmental mosaic[END_REF] found evidence for local adaptation to temperature in an environmental mosaic, at multiple protein-coding genes that were previously detected using a genome scan.

Limecola balthica is a highly dispersive marine bivalve, broadly distributed in the northwestern Europe [START_REF] Väinölä | Biosystematics of Macoma balthica in northwestern Europe[END_REF]. Its pelagic larval phase is estimated to last 2 to 5 weeks [START_REF] Caddy | Development of mantle organs, feeding and locomotion in postlarval Macoma balthica (L.) (Lamellibranchiata)[END_REF]). In the Northern Atlantic Ocean, two lineages of L. balthica occur, due to multiple events of trans-arctic dispersal from the northern Pacific Ocean [START_REF] Väinölä | Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the Macoma balthica complex[END_REF][START_REF] Nikula | Diversity and trans-arctic invasion history of mitochondrial lineages in the North Atlantic Macoma balthica complex (Bivalvia: Tellinidae)[END_REF]. A Pacific lineage (L. balthica balthica) occurs in the White and Baltic Seas, and an Atlantic lineage (L. balthica rubra) is found along the Atlantic coasts from Norway to France and around the British Isles [START_REF] Väinölä | Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the Macoma balthica complex[END_REF]Luttikhuisen et al 2003;[START_REF] Nikula | Diversity and trans-arctic invasion history of mitochondrial lineages in the North Atlantic Macoma balthica complex (Bivalvia: Tellinidae)[END_REF][START_REF] Becquet | Glacial refugium versus range limit: conservation genetics of Macoma balthica, a key species in the Bay of Biscay (France)[END_REF]. The rubra lineage extends to the Gironde Estuary, France (with sparse populations down to Arcachon Basin), which corresponds to the present-day range limit for the species [START_REF] Bachelet | Growth and recruitment of the tellinid bivalve Macoma balthica at the southern limit of its geographical distribution, the Gironde Estuary (SW France)[END_REF]. A hybrid zone between rubra and balthica was detected at the entrance of the Baltic Sea (Kattegat Strait; [START_REF] Nikula | A broad transition zone between an inner Baltic hybrid swarm and a pure North Sea subspecies of Macoma balthica (Mollusca, Bivalvia)[END_REF].

Within rubra, populations north and south of the Finistère Peninsula (France) show significant mitochondrial and nuclear (microsatellite loci) genetic differentiation (Figure 1b; [START_REF] Becquet | Glacial refugium versus range limit: conservation genetics of Macoma balthica, a key species in the Bay of Biscay (France)[END_REF].

An FST -based genome scan was performed on transcriptomic data from pooled individuals, across the Kattegat and Finistère transition zones, in a preliminary effort to look for loci associated with local adaptation [START_REF] Pante | SNP detection from De Novo transcriptome sequencing in the bivalve Macoma balthica: marker development for evolutionary studies[END_REF]. Three geographically-disjunct populations were sampled (Figure 1a). The sites of Aytré (Bay of Biscay, France) and Somme Bay (English Channel, France) correspond to two discrete populations of the rubra lineage [START_REF] Becquet | Glacial refugium versus range limit: conservation genetics of Macoma balthica, a key species in the Bay of Biscay (France)[END_REF]. Aytré is located near the southern limit of the species distributional range and is characterized by warmer sea surface temperatures than Somme Bay (Figure 1b). Furthermore, the Bay of Biscay is subjected to warming surface water temperatures [START_REF] Goikoetxea | Trends and anomalies in sea surface temperature, observed over the last 60 years, within the southeastern Bay of Biscay[END_REF]. A third population, corresponding to the balthica lineage, was sampled in Gdańsk Bay (Baltic Sea, Poland). This FST scan revealed multiple genes involved in the oxidative phosphorylation (OXPHO) system (including genes coding for subunits of the F O /F 1 ATP synthase and NADH dehydrogenase complexes, and an ADP/ATP carrier). In particular, the nuclear gene atp5c1 encoding the gamma subunit of the ATP synthase F 1 rotor, may be under strong selection: 23/41 mapped SNPs were detected as FST outliers (maximum FST = 0.838). This gene is paramount to the good functioning of the F O /F 1 complex, as its rotary action promotes ATP synthesis (reviewed in [START_REF] Sielaff | Twisting and subunit rotation in single FOF1-ATP synthase[END_REF].

In this contribution, we examine molecular signatures of selection at atp5c1 by further analysing the transcriptomic data of [START_REF] Pante | SNP detection from De Novo transcriptome sequencing in the bivalve Macoma balthica: marker development for evolutionary studies[END_REF]. Our goals were to determine whether SNPs that are highly differentiated among populations cause modifications of the protein structure, and whether these changes may have functional repercussions. Given that the F O /F 1 complex is generally encoded by 12 nuclear and 2 mitochondrial genes (reviewed in [START_REF] Rand | Cytonuclear coevolution: the genomics of cooperation[END_REF], we asked whether high-FST, non-synonymous mutations fall within predicted sites of inter-protein interactions, suggesting the implication of genetic incompatibilities in maintaining barriers to gene flow between southern and northern populations of L. balthica.

Material and Methods

Specimen sampling, and preparation of genetic data cDNA libraries were prepared from three pools of 10 individuals (one pool per site) and sequenced on a Roche 454 GS-FLX. Details on sequence quality control and assembly, read mapping, polymorphism detection, and gene annotation can be found in [START_REF] Pante | SNP detection from De Novo transcriptome sequencing in the bivalve Macoma balthica: marker development for evolutionary studies[END_REF]. Here we focus on contig G_c113, detected in this latter study, and identified as coding for atp5c1 (characterized in [START_REF] Matsuda | Gene structure of human mitochondrial ATP synthase gamma-subunit. Tissue specificity produced by alternative RNA splicing[END_REF].

Sequence evolution and predicted changes in protein function

Differences in allele frequencies between population pairs were tested with Fisher's Exact tests, using Popoolation2 [START_REF] Kofler | PoPoolation2: Identifying differentiation between populations using sequencing of pooled DNA samples (Pool-Seq)[END_REF]. The molecular signature of selection was investigated by looking at the ratio between non-synonymous (dN) and synonymous (dS) substitutions (e.g. [START_REF] Kimura | Preponderance of synonymous changes as evidence for the neutral theory of molecular evolution[END_REF], calculated using the R package seqinr [START_REF] Charif | SeqinR 1.0-2: a contributed package to the R project for statistical computing devoted to biological sequences retrieval and analysis[END_REF]. A sliding window (150 nt wide, sliding every 3 nt) was used to detect variation in this ratio along the gene. dN/dS ratios were calculated based on population sequence consensuses, as in [START_REF] Barreto | Interpopulation patterns of divergence and selection across the transcriptome of the copepod Tigriopus californicus[END_REF]. The programs SIFT [START_REF] Kumar | Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm[END_REF][START_REF] Ng | Predicting the effects of amino acid substitutions on protein function[END_REF] and SNAP [START_REF] Bromberg | SNAP: predict effect of non-synonymous polymorphisms on function[END_REF] were used to test whether variations in amino-acid (AA) compositions may have an impact on protein function.

SIFT predicts whether AA substitutions can be tolerated by aligning and comparing homologous sequences retrieved by PSI-Blast; protein structure is not taken into account. SNAP uses sequence information, but can incorporate functional and structural annotations if such data are available [START_REF] Bromberg | SNAP: predict effect of non-synonymous polymorphisms on function[END_REF]. In SIFT, median conservation of sequences was set at 3.00 and sequences ≥ 90 % identical to the query were removed. The hypothesis that local environmental conditions, such as sediment temperature and pH, affect the stability of different alleles of atp5c1 was tested using I-Mutant 2.0 [START_REF] Capriotti | I-Mutant2.0: predicting stability changes upon mutation from the protein sequence or structure[END_REF]. Protein stability (∆∆G) was estimated at 5, 15, and 30°C, and at pH 6, 7, and 8 (conditions based on observed data; Ekebom, 1999 and Lavergne and Barnett, personal communication). In order to further evaluate the impact of AA mutations on protein function, the structure of the gamma subunit was predicted using the I-TASSER server [START_REF] Roy | I-TASSER: a unified platform for automated protein structure and function prediction[END_REF][START_REF] Zhang | I-TASSER server for protein 3D structure prediction[END_REF] and visualized with the UGENE toolkit [START_REF] Okonechnikov | Unipro UGENE: a unified bioinformatics toolkit[END_REF].

Finally, we looked for non-synonymous mutations in the other nuclear-encoded FO/F1 ATP synthase subunits that could be detected in the transcriptome.

Results

Replacement mutations on atp5c1

Amino acid sequences of atp5c1 from populations sampled some 2800 km apart are, overall, highly conserved, with a global dN/dS ratio of 0.08. The sliding window analysis (150 nt wide, sliding every codon; Figure 2a) revealed that 65% of the gene is under negative selection (dN/dS ≤ 0.1, min = 0), while 6% of the gene is under positive selection (dN/dS ≥ 1, max = 2.03). Five AA changes occur between residues 126 and 171, meaning that non-synonymous changes are clustered within a 46-residue fragment spanning 15% of the sequence. This region was predicted by I-TASSER (best model C-score = 0.04; Figure 2b) to overlap with two known peptide binding sites (Protein Data Bank IDs 2W6I and 2HLD; C-scores of 0.09 and 0.07, respectively). CDD annotation (Conserved Domains Database, Marchler-Bauer et al., 2013) following a BLASTX [START_REF] Altschul | Basic local alignment search tool[END_REF] search for G_c113 suggest that these mutations occur at the interface with the core domain (interaction with the alpha-beta hexamere) and with the epsilon subunit. These non-synonymous changes were predicted to be evolutionarily tolerated based on SIFT and SNAP (expected SNAP accuracy 85-94%).

The first four AA changes (L126A, N131D, D136N, H156D; nomenclature: Aytré considered as mutant) discriminate Aytré from Gdańsk Bay and Somme Bay, while the fifth one (T171S) distinguishes Gdańsk from the two French populations. N131D and H156D have a stabilizing effect on the protein, while all others have a destabilizing effect (I-Mutant tests between 15 and 30°C at pH7; Figure 3). In all cases, a relative increase in ∆∆G was observed with increasing temperature and pH (i.e., in warmer conditions, a destabilizing mutation was less destabilizing). In one case (N131D), ∆∆G changed sign (from negative to positive) with increased temperature (in all pH conditions), indicating that the expected effect of the mutation shifted from destabilizing to stabilizing. Considering the 4 AA mutations described above, the Aytré population contained a mix of two alleles, the major allele (85.5% of the reads) being unique to Aytré and one corresponding to the type found at Somme and Gdańsk (14.5% of the reads). The median depth of coverage for contig G_c113 was 55, 45 and 18 for Aytré, Gdańsk and Somme, respectively (Table 1).

Potential for incompatibilities with other ATP synthase subunits

Contig G c1077 was identified as coding for the alpha subunit and bear one outlier FST when comparing the samples from Aytré and Somme (FST = 0.672, [START_REF] Pante | SNP detection from De Novo transcriptome sequencing in the bivalve Macoma balthica: marker development for evolutionary studies[END_REF]. This contig could be reliably placed in ORF, which revealed that all nine mutations were synonymous. All mutations segregated Aytré from Somme and Gdańsk. Several contigs were identified as coding for the beta subunit; one (A c2900) was included in our initial FST scan but did not stand out as bearing outlier SNPs (in addition, this contig could not be reliably placed in ORF). Unfortunately, the epsilon subunit could not be detected in our transcriptome dataset.

Discussion

Adaptive divergence at OXPHO loci Mitochondrial function can be significantly influenced by temperature (e.g. [START_REF] Dahlhoff | Effects of temperature on mitochondria from abalone (genus Haliotis): adaptive plasticity and its limits[END_REF]. In Dutch populations of L. balthica, respiration rate was shown to increase with temperature intra-seasonally [START_REF] Hummel | The respiratory performance and survival of the bivalve Macoma balthica (L.) at the southern limit of its distribution area: a translocation experiment[END_REF]. In addition, Dutch specimens transplanted in the Spanish estuary of Bidasoa, 200 km south of the known species range limit (i.e. the Gironde Estuary, some 80 km from our sampling site in Aytré, France) showed respiration rates significantly higher than in their native populations [START_REF] Hummel | The respiratory performance and survival of the bivalve Macoma balthica (L.) at the southern limit of its distribution area: a translocation experiment[END_REF]. Our analyses on the effect of temperature on the stability of subunit gamma suggest a small, positive effect of each replacement mutation in the southern population. While only five replacement mutations were detected, previous studies have shown that even few differences in amino acid composition can have a significant impact on mitochondrial performance. In the marine copepod Tigriopus californicus, three amino acid changes at the nuclear cytochrome c, separating populations characterized by different thermal regimes, correlated with significant differences in cytochrome c oxidase activity [START_REF] Rawson | Functional coadaptation between cytochrome c and cytochrome c oxidase within allopatric populations of a marine copepod[END_REF].

In the seed beetle Callosobruchus maculatus, metabolic rates were detected between cytotypes that differed by one amino acid at COI and one amino acid at Cyt-B [START_REF] Arnqvist | Genetic architecture of metabolic rate: environment specific epistasis between mitochondrial and nuclear genes in an insect[END_REF]

see below).

In killer whales (Orcinus orca), single amino acid changes at Cyt-B are associated with Antarctic ecotypes, and are therefore possibly implicated in improved mitochondrial performance in polar waters [START_REF] Foote | Positive selection on the killer whale mitogenome[END_REF]. Cyt-B was also recently implicated in thermal adaptation in the European anchovy [START_REF] Silva | Thermal adaptation and clinal mitochondrial DNA variation of European anchovy[END_REF]. Finally, the synergistic impact of temperature and few replacement mutations on protein stability and function was recorded at the nuclear locus encoding cytoplasmic malate dehydrogenase (cMDH) in limpets [START_REF] Dong | Temperature adaptation of cytosolic malate dehydrogenases of limpets (genus Lottia): differences in stability and function due to minor changes in sequence correlate with biogeographic and vertical distributions[END_REF]. As the surface water of the Bay of Biscay is known to be warming at rate of 0.26°C/decade (1977-2007 time period, [START_REF] Goikoetxea | Trends and anomalies in sea surface temperature, observed over the last 60 years, within the southeastern Bay of Biscay[END_REF] and the geographical range [START_REF] Jansen | Macoma balthica in Spain, a few decades back in climate history[END_REF] of southern L.

balthica is thought to be receding, it becomes increasingly important from a conservation standpoint to understand how genetic diversity at OXPHO genes relates to adaptive potential. Indeed, as gene flow is limited between southern and northern rubra populations at atp5c1, a narrowing geographical range could result in the loss of adaptive alleles if the southern population does not persist.

Implications for intrinsic genetic incompatibilities among ATP synthase subunits

Our data suggest that five replacement mutations may interact with other subunits, likely epsilon (encoded by atp5e in humans, [START_REF] Tu | Cloning, characterization and mapping of the human ATP5E gene, identification of pseudogene ATP5EP1, and definition of the ATP5E motif[END_REF]. This subunit sits on the F O stator, which is composed of nuclear-and mitochondrial-encoded genes and is embedded in the mitochondrial membrane.

While atp5e was not detected in our dataset, characterizing genetic variation at atp5e and other ATP synthase subunits (including mitochondrial and nuclear ones) across the Finistère transition zone seems central to understanding how intrinsic genetic incompatibilities are involved in maintaining barriers to gene flow in L. balthica. The sequencing of mitochondrial genomes from individuals sampled on either sides of the Baltic and Finistère hybrid zones [START_REF] Saunier | Mitochondrial genomes of the Baltic clam Macoma balthica (Bivalvia: Tellinidae): setting the stage for studying mito-nuclear incompatibilities[END_REF] will help us shed some light on possible incompatibilities between the nuclear and mitochondrial genes coding for ATP synthase subunits. In addition, the recent discovery of sex-linked heteroplasmy in L. balthica suggests that mitochondria in this species are characterised by Doubly Uniparental Inheritance (DUI; [START_REF] Pante | Putative sex-linked heteroplasmy in the tellinid bivalve Limecola balthica (linnaeus, 1758)[END_REF], and unpublished male mitogenome draft). In DUI species, the somatic tissues of males are characterized by a 'female' mitotype (as in all tissues of females), while the male germ line is characterized by a 'male' mitotype that is passed on from fathers to sons (reviewed in Zouros, 2013). Female and male mitotypes sampled from a single individual can be highly divergent, reaching up to 52% in freshwater unionoid mussels (Doucet-Beaupré et al 2010).

In L. balthica, genetic incompatibilities could therefore occur in multiple ways, as OXPHO epistatic interactions among nuclear genes, nuclear and female mitochondrial genes, and among nuclear and male mitochondrial genes in the sperm of interpopulational hybrids.

Interplay between intrinsic genetic incompatibilities and temperature

One fascinating avenue for research is the characterization of interactions between intrinsic genetic incompatibilities and extrinsic environmental forces. [START_REF] Arnqvist | Genetic architecture of metabolic rate: environment specific epistasis between mitochondrial and nuclear genes in an insect[END_REF] were the first to demonstrate environmental effects on mitonuclear epistatic interactions by crossing mitochondrial and nuclear genomes of C. maculatus. In this very elegant experiment, the authors showed that the negative effects of genetic incompatibilities on metabolic rate were only detectable when hybrids were exposed to different temperatures [START_REF] Arnqvist | Genetic architecture of metabolic rate: environment specific epistasis between mitochondrial and nuclear genes in an insect[END_REF]. In L. balthica, the replacement mutations mapped on atp5c1 (i) are located in an inter-peptide interaction zone (most likely in the area were the gamma subunit interacts with epsilon), and (ii) seem to influence the stability of the gamma subunit depending on temperature and pH. This preliminary study is based on very limited data, as we used few, pooled individuals from three populations to investigate the molecular signatures of selection on atp5c1. Future research should therefore focus on determining, with larger population sampling, if subunits of the F O /F 1 ATP synthase complex are indeed involved in mitonuclear genetic incompatibilities enforcing genetic barriers, and whether these putative incompatibilities are affected by the environment.

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. (a) Sliding window analysis of the dN/dS ratio along atp5c1 (window 50 AA wide, sliding every codon). The position of replacement mutations is provided on the right (grey dots), as well as the position of conserved interaction sites detected by CDD (Marchler-Bauer et al., 2013) (black, blue and red: interactions with core domain, delta and epsilon subunits, respectively). (b)Threedimensional model of the ATP synthase gamma subunit of Limecola balthica (I-TASSER prediction based on reference contig from Gdańsk). The yellow segment represents the 46 residue-long fragment bearing replacement mutations (grey dots).

Figure 3 .

 3 Figure 3. Effect of temperature and pH on the stability (represented as ∆∆G) of the ATP synthase gamma subunit. Error bars represent one standard deviation of three measures of ∆∆G per temperature, using different pH values (see methods). The name of replacement mutations is provided on the right.

Table 1 .

 1 Read counts at each position of codons containing non-synonymous mutations, for each population. The significance of differences in allele frequencies are presented as -log10 of p-values from Fisher's Exact tests (significance for at least one pair of populations: *: -log10(p)>1; **: -log10(p)>10).

					Fisher's Exact
		Aytré (A)	Somme (S)	Gdansk (G)	Test
	Locus Pos	Major Minor	Major Minor	Major Minor	A-S A-G S-G
	L126A 376	G (47) T (14)	T (25)	T (50)	11.3 18.4 0.0	**
	377	C (51) T (15)	T (26)	T (51)	12.0 19.2 0.0	**
	378	A (49) G (15)	G (26)	G (52)	11.7 18.8 0.0	**
	N131D 391	G (50) A (13)	A (24)	A (48)	11.7 19.0 0.0	**
	392	A (58)	A (24)	A (47)	0.0	0.0 0.0
	393	C (55)	C (24)	C (44)	0.0	0.0 0.0
	D136N 406	A (45) G (13)	G (27)	G (51)	12.0 18.4 0.0	**
	407	A (58)	A (27)	A (53)	0.0	0.0 0.0
	408	C (58)	C (25)	C (53)	0.0	0.0 0.0
	H156D 466	G (36) C (14)	C (26)	C (48)	9.8	14.9 0.0	**
	467	A (50) G (1)	A (25) G (1)	A (49)	0.0	0.0 0.5
	468	C (43) T (5)	C (23)	C (44)	0.8	1.2 0.0	*
	T171S 511	A (53)	A (21)	A (42)	0.0	0.0 0.0
	512	G (51) C (2)	G (21)	C (35) G (7)	0.0	15.9 10.4 **
	513	C (56)	C (21)	C (40)	0.0	0.0 0.0

Acknowledgements

We thank Patrick Triplet, Antoine Meirland (Association GEMEL, Picardie) and Rafal Lasota for specimen collections, and the Molecular Core Facility at the University of La Rochelle. We also thank two anonymous reviewers for improving the manuscript. This work was funded by the French Agence Nationale de la Recherche (Hi-Flo project ANR-08-BLAN-0334; HySea project ANR-12-BSV7-0011); salary for EP was covered by a "Contrat de Projet Etat-Région." Funding from the Université de La Rochelle (ACI JCJC) allowed the authors to present at PhysioMar17.


Bhuvanesh Srinivasan

Sylvain Le Tonquesse

Alain Gellé

Cédric Bourgès

Leo Monier

Isao Ohkubo

Jean-François Halet

David Berthebaud
email: david.berthebaud@cnrs.fr

Takao Mori
email: mori.takao@nims.go.jp

Screening of transition (Y

Keywords: Thermoelectrics, GeTe, Zr-doping, defects, band engineering

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Energy sustainability is one of the major key technological challenges in this 21 st century. Growing awareness and alarming concern for the environment and renewable energy supplies has revitalized advances in materials engineering and technologies for energy conversion in recent years. Thermoelectric (TE) devices with their potential to reversibly convert waste heat into useful electricity propound the likelihood of an all-solid-state technology for power generation, refrigeration, temperature stability, and control. [1][2][3][4] The TE material's performance is evaluated in terms of a dimensionless figure of merit, zT = S 2 σT/κ where S, σ, T and κ are the Seebeck coefficient, electrical conductivity, absolute temperature and total thermal conductivity (sum of the electronic contribution, κe, and the lattice contribution, κlatt), respectively. Boosting the TE power factor (S 2 σ) and/or suppressing the thermal transport () are the key paradigm to realize high zT. Utilizing phononic scattering mechanism, that is achieved via nanostructuring, which includes the engineering of nanoscale defects (dislocations, stacking faults, point defects, nanopores), and multiscale hierarchal architecting, [5][6][7][8][9] imparting rattling impurities, 10 intrinsic bond / strong lattice anharmonicities, [11][12][13] etc., to suppress thermal transport (κlatt) is a well-known strategy to enhance zT. However, the electrical transport properties, S and σ are interrelated via carrier concentration, but in a contrary way (as per Pisarenko's relation) and confront a bigger challenge in improving the power factor, paramount for better energy conversion efficiency. In the last decade, this issue has been masterly addressed with the revitalization of the 'band structure engineering' approach, which encompasses electronic band convergence, 14,15 fostering resonance states or local distortion of the density of states (DOS) near the Fermi level (by impurity atoms), 16 quantum confinement of charge carriers, 17 nestification, 18 electron filtering effect, 19,20 dimensionality reduction, 21 modulating light bands with low effective mass, 22 deformation potential coefficient, 23 blocking of minority charge carriers, 24,25 which effectively helped in decoupling S and σ to a certain extent. Even the ideas such as the utilization of magnetic effect, [26][27][28][29] partial crystallization of semiconducting chalcogenide glasses (based on phonon-glass electron-crystal approach) [30][31][32][33][34] and processing techniques like Hybrid Flash-SPS techniques [35][36][37] were explored with notable degrees of success. This concept of band engineering has been effectively used to tune the TE properties in various p-and n-type materials such as PbTe, 38,39 SnTe, [40][41][42][43] Mg2Si, 44 half-Heusler, 45 etc. In this category of midtemperature range TE candidates, GeTe has gained prominence in the past few years. It exists in a rhombohedral crystal structure (r-GeTe) and turns into a cubic structure (c-GeTe) at ~700 K. 46 Due to higher valley degeneracy, the TE performance is known to maximize at the cubic high symmetric domain 35,47 and hence, reducing the phase transition temperature can boost/stabilize the TE performance over a broad temperature range. Though this has been the general norm/convention with GeTe, there have also been recent studies where high TE performance were reported even for the low temperature r-GeTe, because of the convergence of split band valleys derived from the rhombohedral distortion. 48- 50 Generally, the high carrier concentration (in the order as high as ~10 21 cm -3 ) due to the low formation energy of Ge vacancies in pristine GeTe affects the TE performance in this p-type highly degenerated semiconductor, 46,51 thus it is essential to optimize the charge carrier level to improve the TE performance in GeTe. Besides the high carrier concentration and phase transition, the presence of multiple valence bands in GeTe, provide these p-type intriguing materials with diverse degrees of freedom to play and tune their transport properties. 52 Some of the earlier reported approaches for GeTebased materials to improve their electrical transport and/or to reduce the thermal transport have been adopted on single dopant containing compositions such as Ge1-xPbxTe, 53 Ge1-xBixTe, 54 Ge1-xInxTe, 55 Ge1-xGaxTe, 35 Ge1-xCrxTe, 56 Ge1-xTixTe, 57 Ge1-xAlxTe, 58 GeTe1-xSe, 59 Ge1-xSbxTe, 60,61 Ge1-xAgxTe, 62 Ge1-xMnxTe, 63,64 LiI-GeTe, 65 etc., where the zTmax value generally ranged between 1.0 -1.5 in most of these singly doped cases. Similarly, many co-doped compositions such as GeTe-AgSbTe2 (TAGS), 66 GeTe-AgInTe2, 67 GeTe-LiSbTe2, 68 (GeTe)nSb2Te3, 69 Ge9Sb2Te12-x, 70 GeTe-AgSbSe2, 71 (Bi2Te3)nGe1-xPbxTe, 72 Ge1-xSbxTe1-ySey, 73 Ge1-x-ySnxPbyTe, 74 GeTe-GeSe-GeS, 75 Ge1-x-yBixInyTe, 76,77 Ge1-x-yBixSbyTe, 47 Ge1-x-yPbxBiyTe, 48 (Ge1-xDyxTe)y(AgSbTe2)1-y, 78 Ge1-x-yCdxBiyTe, 79 Ge1-x-yGaxSbyTe, 35 Ge1-x-yPbxSbyTe, 49 (Ge1-xYbxTe)y(AgSbTe2)1-y, 80 Ge1-x-ySbxZnyTe, 81 Ge1-x-yMnxBiyTe, 63,82 Ge1-x-yMnxSbyTe, 63 Ge1-x-yPxSbyTe, 83 (Ge1-xDyxTe)y(AgSbTe2)1-y, 80 and more recently Ge1-x-yCrxBiyTe 56 and Ge1-x-yTixBiyTe 57 show zTmax values generally ranged between 1.5 -2.0 or even more for some compositions.

The interesting TE results obtained from the aforementioned compositions, especially the couple of recent findings on transition metal (Ti and Cr) doped GeTe, 56,57 trigger us to further explore some of the other rarely studied transition metals like Y, Zr, Hf, V, Nb, Mo, and Ru, and rare-earth metals like La and Pr as potential dopant candidates for GeTe. A thorough experimental screening of these elements is accompanied with theoretical computations to understand how the addition of these transition and rare-earth elements can influence the electronic band structure (and thus the resulting TE properties) of GeTe. In most general cases, pristine GeTe exhibits zT ~0.85 at 700 K. 58,65 In this work, the strategy was that, if any of these transition or rare-earth elements doping to GeTe showed promising TE properties (i.e., zT > 1), then they were optimized for their composition, and the resulting optimized compounds were finally codoped with Sb, which is widely reported to improve the band degeneracy in GeTe by reducing the phase transition temperature to produce more hole pockets. 35,47 In a nutshell, this work investigates the effect of doping a series of transition and rare-earth elements and their best performing optimized codoping (Zr-Sb) on the structural, electronic and thermoelectric properties of GeTe.

Materials & Methods

Reagents

Ge (Ingot, Rare metallic, 99.999 %), Te (Shot, Rare metallic, 99.999 %), Mo (Powder, Rare metallic, 99.9%), Nb (Powder, Wako, 99.9%), V (Turnings, Wako, 99.5%), Y (Ingot, Aldrich, 99.9%), Ru (Powder, Wako, 99.9%), Zr (Sponge, Aldrich, 99%), Hf (Turning, Wako, 99.9%), Pr (Turning, Wako, 99%), La (Turning, Wako, 99%) and Sb (Shots, Aldrich, 99.999%) were used for synthesis without any further purification.

Synthesis

The samples were synthesized using a vacuum-sealed tube melting process. Appropriate stoichiometric amounts of the starting elements were introduced into a quartz tube, and the ampoules that were sealed under a vacuum of ~2 x 10 -3 mbar were then placed in a melting furnace and slowly heated to 1223 K for 8 hours, and then held at that temperature for 10 hours and finally cooled down to room temperature (~6 hours). Obtained ingots were hand crushed into powders, which were then consolidated by Spark Plasma Sintering (SPS-1080 -SPS Syntex Inc.) at 750 K (heating rate ~80 K/min) for 5 mins (holding time) with an axial pressure of ~60 MPa. After SPS processing, highly dense disk-shaped pellets were obtained with theoretical densities of > 98%. The sintered discs were then cut and polished to the required shapes and dimensions for various thermoelectric measurements.

Powder X-ray diffraction

Powder X-ray diffraction (PXRD) patterns were recorded at room temperature in the 2θ range of 10 o -90 o using Rigaku Smart Lab 3 diffractometer (Cu radiation, λKα1 = 1.5418 Å, and λKα2 = 1.5444 Å in a ratio λKα2 / λKα1 = 0.5) with a step size of 0.02 o and a scan speed of 2 o /min. Le Bail refinements were performed using the WinPLOTR software and the FullProf algorithm. [84][85][86] The shape of the diffraction peaks was modeled using a Thompson-Cox-Hastings pseudo-Voigt profile function. 87 Zero-point shift, peak shape, and lattice parameters were systematically refined, and the background contribution was fitted using a 6-order polynomial. More realistic estimated standard deviations were calculated using the Bérar method. 88

Electrical and thermal transport

The Seebeck coefficient and electrical resistivity were measured simultaneously on ~2×2×9 mm 3 bars using a commercial instrument (ZEM-2, ADVANCE RIKO Inc.) under partial pressure of He.

The temperature-dependent thermal diffusivity, D, was measured using the laser flash diffusivity method in LFA-467 Hyperflash (Netzsch) on graphite coated disc-shaped samples of ~10 mm diameter and ~2 mm thickness. The temperature-dependent heat capacity, Cp, was derived using a standard sample (pyroceram-9060) in LFA-467, which is in good agreement with the Dulong-Petit Cp value. The total thermal conductivity was calculated using the formula κ = D×Cp×ρ, where ρ is the density of the sample (measured using Archimedes' principle). In some cases, to better understand the thermal transport properties, the contributions from electronic and lattice parts were calculated. The lattice thermal conductivity (latt) was estimated from  by subtracting the electronic contribution (e) via the Wiedemann-Franz law, as in the following equation,

𝜅 𝑒 = 𝐿𝜎𝑇 (1)
where L is the Lorenz number computed by the condensed version of Single Parabolic Band model with acoustic phonon scattering (SPB-APS) 89 , as in the following equation,

𝐿 = 1.5 + exp [- |𝑆| 116 ] (2)
where the Lorenz number (L) is in 10 -8 WΩK -2 and the Seebeck coefficient (S) is in µVK -1 .

The uncertainty in the results for the values of electrical and thermal transport properties was 5% and 7%, respectively, and for the overall zT was  12%. 90 To increase the readability of the figures, the error bars are not shown in the figures.

Computational procedures

First-principles calculations using Density Functional Theory (DFT) were performed to understand the effect of transition metal and/or rare-earth dopants on the electronic states of GeTe. The projectoraugmented-wave (PAW) approach 91 implemented in the Vienna ab initio simulation package (VASP) 92 was used. Computations were performed using the generalized gradient approximation (GGA) for the exchange-correlation term parameterized in the Perdew-Burke-Ernzerhof functional. 93 Spin-orbit coupling (SOC) was included in the calculations.

Calculations were performed on cubic structural models, as we were interested in the high-temperature behavior of doped GeTe compounds. Dopants were substituted for Ge atom in a 4 × 4 × 4 super-cell (MGe63Te64 cubic models; M = transition / rare-earth element). Similar to the work of Hoang et al. on the impurity clustering in GeTe, 94 we chose a cluster made of one M atom surrounded by 6 Sb atoms (Sb and M atoms were placed far from each other in the supercell) for the M-Sb codoped GeTe composition (Sb6MGe57Te64). In all the calculations, the positions were fully relaxed. For Ru doping, calculations using the GGA + U + (J) method were also carried out with Hubbard value U = 2 eV and Hund value J = 0.5 eV). 95,96 Almost no differences were found compared to the GGA calculations, hence only the results obtained from GGA calculations will be presented in the article.

Results & Discussion

Firstly, the primary motive is to figure out which among the planned transition elements (Nb, Y, Hf, Zr, Ru, Mo, V) and rare-earth elements (Pr, La) as dopants can potentially improve the TE performance in GeTe. Considering the fact that in most general cases the pristine GeTe exhibits zTmax ~0.85 at 700 K, the threshold that we set here to consider that transition / rare-earth dopant to be of any potential benefit for GeTe thermoelectrics, only if the zT exceeds 1 in those doped compositions. In order to shortlist the potential candidate among Nb, Y, Hf, Zr, Ru, Mo, V, Pr, and La, a small proportion of each of the element was doped to GeTe. Depending on their oxidation states, the doping level was fixed at either 2 mol% (for elements with 2 + or 3 + oxidation state like La, Y, V, Ru) or 0.5 mol% (for elements with 4 + oxidation state like Zr, Hf).

Structural analysis

PXRD results for the pristine GeTe, Ge0.08Y0.02Te, Ge0.995Hf0.005Te, Ge0.995Mo0.005Te, Ge0.98V0.02Te, Ge0.995Zr0.005Te and Ge0.98La0.02Te compositions are presented in Figure 1. In all cases, the sharp reflections from PXRD indicated the high crystallinity of the SPS processed samples. The major reflections from the PXRD pattern could be indexed to a rhombohedral GeTe phase (R3m space group, lattice parameters, a ~ 4.165 Å, and c ~ 10.670 Å), which was confirmed further with the presence of double reflections [(024) and ( 220)] in the range of 2θ values between 41 o to 44 o in Ge1-xMxTe (M = transition / rare-earth element; x = 0.02 or 0.005). Besides the low-temperature rhombohedric GeTe phase (main phase), a small amount of Ge impurities (cubic Ge, 𝐹𝑑3 ̅ 𝑚 space group, a ~ 5.656Å) were also observed for the majority of the samples. These minor cubic Ge reflections could be due to the thermodynamic nature of the GeTe, which is known to preferentially exist as a Te-rich composition. 47,50,62,97 Moreover, the formation energy of Ge vacancies is much lower (about just onethird) when compared to the formation energy of other 3 defects in GeTe, i.e., Ge antisite, Te vacancy, and Te antisite, thus making it the most easily and favorably formed intrinsic defect. 98,99 In thermal equilibrium, Ge vacancies are the dominant defect for any Fermi level (EF), even if Ge is the only carrier reservoir, and are envisaged to form impromptu for EF ~0.25 eV above the valance band. 99 Generally, it is reported that this minor cubic Ge impurity phase in a small proportion in GeTe has no apparent detrimental influence on the TE performance. 50,62 Since rhombohedric GeTe is a faintly distorted rocksalt lattice along the [111] direction, 100,101 the impact of these doping elements on the crystal structure can be captured by the variation in lattice parameter and/or the α angle in the rhombohedral setting. The PXRD patterns of pristine and doped GeTe compounds were fitted by the Le Bail method 86 and their respective refined lattice parameters and α angles can be found in the Supporting Info (SI, Table S1).

All these samples (pristine, Y, La, Zr, Hf, and Mo) have a fairly similar unit cell parameters, with a and c being 4.165 Å and 10.670 Å, respectively, and the hexagonal cell volume varying between ~160.2 -160.5 Å 3 . The only exception is the Ge0.98V0.02Ge compound, which possessed a significantly smaller c = 10.573 Å, thus resulting in a slightly smaller unit cell volume of ~159.0 Å 3 . This difference or stand out behavior of Ge0.98V0.02Ge compound, when compared to other GeTe compounds presented here, is not just in terms of its crystallographic structural point of view, but also in terms of its electrical and thermal transport properties, which are discussed in the forthcoming section. PXRD results for the Ge0.98Ru0.02Te, Ge0.98Pr0.02Te and Ge0.98Nb0.02Te compositions are presented in Figure 2. In addition to the low-temperature rhombohedric GeTe main phase, these compounds that were doped with Ru, Pr and Nb also include a non-negligible amount of high-temperature cubic GeTe phase (space group 𝐹𝑚3 ̅ 𝑚). The presence of the secondary cubic GeTe phase in these compounds (Ru, Pr and Nb doped GeTe) can affect their TE transport properties, as observed earlier in the case of Ge0.94Ag0.06Te. 62 From the PXRD patterns, it is important to notice that Bragg peaks corresponding to the crystallographic plans normal or almost normal to the c direction are characterized by peak shape with larger tails, i.e. more Lorentzian characteristic, than other reflections. It affects the sample's diffraction pattern and is more prominent for the doped compositions, for example, for Ge0.98V0.02Te compound than for undoped GeTe (SI, Figure S1). Such an effect, that is embodied by more Lorentzian feature, can be attributed to the presence of defects, which were already elucidated in detail (using HRTEM imaging) in several GeTe-based materials in the past. 51,102 The high values of the agreement factors (Rp, Rwp, χ 2 ) from the fitted patterns are due to the anisotropic peak broadening arising due to the defects in these GeTe compounds. These defects would significantly influence and play a key role in the thermal transport properties of these compounds.

Effects of V, Nb, Mo, Ru, and Pr on TE properties of GeTe

The TE transport properties of some of the transition-metal and rare-earth-doped GeTe compounds are presented in Figure 3. The short inference would be that the dopants V, Nb, Mo, Ru, and Pr were found to be detrimental to the thermoelectric properties of GeTe. Partially substituting these elements (V, Nb, Mo, Ru, and Pr) for Ge decreased their zT when compared to that of pristine GeTe (Figure 3e). Out of those, Ru was found to be the worst, followed by Nb, V, Pr, and Mo. Most of these dopants neither helped in significantly improving the power factor, PF = S 2 σ (Figure 3c) nor in suppressing the thermal transport (Figure 3d), thus resulting in a poor TE performance. The only exemption is the V-doping, which substantially decreased the thermal transport, however, the PF was adversely affected, thus offsetting the effect of suppressed κ.

To have a more cogent understanding of the impact of these transition and rare earth dopants on the electrical transport properties in GeTe, DFT calculations were performed. These DFT calculations were carried out on 4 x 4 x 4 supercells derived from the cubic structural arrangement of GeTe (hightemperature phase). The electronic densities of states (DOS) computed for the cubic models of V, Nb, Mo, and Ru-doped compounds are presented and compared with undoped c-GeTe in Figure 4. Their electronic band structures were also computed and analyzed, aiming to provide some insight into how the valence band structure is modified upon doping. The band structures are plotted in Figure 5 along some high symmetry lines of the cubic Brillouin zone (BZ). The band gap (Eg) at L-point for pristine c-GeTe was calculated to be 226 meV, in agreement with the previous reports. 35,47,58,103 The DOS calculations (Figure 4a-d), which shows the shift of the Fermi level away from the valence band (i.e., upwards / toward the conduction band) in the case of Nb and Mo doping, clearly suggests the suppression of the p-type charge carriers (holes) with the substitution of these dopants. It must be noted that the Eg mentioned here and throughout this paper, corresponds to the energy difference between the Te and the Ge bands at the L-point. When there is no intruder state (i.e., dopant/impurity state) in the gap, it corresponds to the vertical gap at the L-point. Since the variation of this value reflects the difference in the energy gap between Te and Ge bands, it is more suitable to interpret the changes in the electrical properties. Band structure computations show that for the cases of Ru, Mo, and V doping, intruder states are present in the gap (Figures 5d-f). According to the perturbation theory (PT), the intruder states are defined as correction functions having energy (eigenvalue of H0) close to the reference energy 𝐸 𝑚 (0) , thus causing convergence problem in the perturbation summation, thus being branded as an Achilles' heel in multireference PT. 104,105 In the case of V-doped GeTe, the intruding Vbands are flat and localized within the bandgap (Figure 5f). The band gap at the L-point, without considering the intruder states, has increased from 226 meV to 250 meV when compared to the pristine c-GeTe. This can explain the reduction in its electrical conductivity (Figure 3a). V-doping modifies and distorts the DOS with respect to that of GeTe, as it seems to have induced resonance states (marked by the presence of a sharp hump), adjacent to EF (Figure 4d). The presence of such resonance levels near EF can help to increase S (as per Mott's relationship), like in the present case, as V-doped GeTe stands out in terms of thermopower (at least in the lower temperature domain) when compared to any other compositions presented here. This standout behavior is not just in terms of their transport properties, but also from the crystallographic point of view, where they possessed a significantly smaller c = 10.573 Å, thus resulting in a slightly smaller unit cell volume of ~159.0 Å 3 , and a reduced c/a ratio of 2.537. At the same time, it must be noted that too much distortion of DOS near EF can affect the carrier mobility, and this adds up to the reasons for this V-doped GeTe exhibiting poor σ. Had these resonance states been more close, i.e., at the vicinity of EF (instead of being adjacent), then there could have been a much significant beneficial impact on the thermopower, as in the case of previously reported In-doped GeTe. 55 Such impurity states are also induced for Ru-doped GeTe and Mo-doped GeTe, but they are less likely to make any gainful influence on the thermopower, as they are not present close to the proximity of EF.

For both Mo and Ru-doped GeTe, a part of the impurity states are localized near the top of the valence band and another part expands in the whole valence band. Compared to 3d orbitals of V, the 4d orbitals of Mo and Ru are more coupled to the p orbitals of Te. For Mo-doped GeTe and V-doped GeTe, besides the presence of intruder states, the first-principles calculations performed and reported here reveal some magnetic characters, i.e., DOS calculations yielded non-zero total magnetization. Though the introduction of magnetic ions (albeit with very careful tuning of the composition) have recently been shown to boost the power factor by balancing between improved effective mass (m*) and reduced mobility (µ) in certain specific compounds like CuGaTe2, Bi2Te3, SnSe, 8,29,106 it seems not to be the case here in GeTe. Previous spin-polarized band structure computations have proved that the introduction of magnetic ions in GeTe will induce spin scattering, which in turn will negatively affect the carrier mobility in these compounds. Typically, if the carrier density is optimized, then it is the weighted mobility µW = µ(m*/m0) 3/2 , where mo is the free electron mass, that governs the maximum thermoelectric power factor in a given material. Apart from the variety of reasons explained above (such as the presence of intruder states, large Eg, heavy DOS distortion), this magnetic contribution (and thus the low µ) adds up to explain why V-doped GeTe has exhibited poor PF, which is least among the studied materials here (Figure 3c). The high resistivity resulting from V-doping has resulted in reduced κ (Figure 3d), mainly due to the largely suppressed electronic contribution, κe. Given that the Mo-doping (with some predicted magnetic interactions) did not affect PF as badly as V-doping, in fact showing comparable PF as that of pristine GeTe sample (Figure 3c), there is a possibility to improve PF further in Mo-doped GeTe by playing with this interaction between carriers and magnetic moments via composition tuning to achieve a better trade-off between m* and µ, thus throwing an interesting aspect for potential future work. However, the calculations predicted the Ru atom to be non-magnetic with a d 6 electronic configuration. This is coherent with the crystal field theory in octahedral symmetry, where the t2g orbitals are fully occupied. GeTe band structure manifests a second maximum of the valence band in the  → K () direction. The energy separation between these light and heavy hole valence bands (ELΣ) for undoped c-GeTe was found to be 64 meV, in agreement with the literature report. 103 This energy difference (ΔELΣ) for high-temperature c-GeTe is much lower than that of the lowtemperature r-GeTe (ELΣ = 150 meV). 103 This is beneficial and is consistent with the fact that the temperature increases the band convergence in GeTe. In the case of Nb and Ru doping, the negative value of the energy separation between the two valence bands, ELΣ (Table S2, SI), indicates that the heavy hole band is larger in energy than the light hole band (it can be seen from Figures 5c,d). Nb states are located just at the bottom of the conduction band. The higher κ in Ru, Nb and Pr-doped GeTe samples can be attributed to the presence of a significant proportion of secondary cubic GeTe phase, as mentioned in the previous section and shown from the PXRD patterns given in Figure 2. Such high symmetric cubic GeTe secondary phases are generally known to exhibit higher lattice thermal conductivity than the low symmetric rhombohedral GeTe main phase. 48 It is not a surprise that Rudoped GeTe, with a culmination of aggregate effects (mostly negative on TE transport properties) such as the presence of intruder states, valence band divergence (hence lower S), secondary cubic GeTe phase (higher κlatt), languishes at the bottom as the least thermometrically performing material (least zT) among the transition and rare-earth element doped GeTe compounds that are reported in this work (Figure 5e).

Effects of Y and Hf on TE properties of GeTe

Substitution of Y and Hf for Ge was found to be relatively promising, as they marginally improved zT in GeTe from ~0.85 at 723 K (for pristine case) to ~0.95 (Y-substituted) and ~0.98 (Hf-substituted) at 723 K (Figure 5e). In the case of Y-doping, the marginal improvement in zT is due to the suppressed thermal transport, while in the case of Hf-doping, it is due to the marginally improved power factor. Despite their slightly improved TE performance, Hf and Y-doped GeTe compounds exhibited zT < 1, hence they were not shortlisted for further experimental studies in this work. Of course, there is still some scope for further improvements with these Y and Hf dopants, if the compositions are properly optimized, presenting scope for future studies. DFT results show that Y-doping does not change much the energy separation between the two valence band maxima (Figure 6c) when compared to that of pristine c-GeTe. The DOS results (Figure 6a) also shows that the projections arising from the Y-doping do not impart any change near EF or near the top of the valance band or near the bottom of the conduction band, and instead, the dopant states arising from Y-doping seem to be buried deep within the conduction band. This explains why Y-doping did not influence the power factor values at high temperatures. Hfdoping to GeTe slightly tends to open up the band gap and tends to slightly decrease the energy separation between the two valence band maxima (Figure 6d), promoting band convergence and must explain their improved power factor (and hence the overall marginal improvement in zT). The values of Eg and ΔELΣ for Hf-doped GeTe (Table S2, SI) are computed for the HfGe63Te64 model, which has a much higher Hf-content than the experimental composition (Ge0.995Hf0.005Te), so only the tendency has to be taken into account.

Effects of La and Zr on TE properties of GeTe

In the present study on transition and rare-earth elements as potential dopants for GeTe, La and Zr were found to be effective (Figure 7). Substituting La and Zr for Ge considerably enhanced zTmax of GeTe from ~0.85 at 723 K to ~1.1 at 673 K and 723 K, respectively, for La-doping and Zr-doping, as shown in Figure 7e. In both cases, i.e., La-doping and Zr-doping, the thermal conductivity did not change much when compared to pristine GeTe (Figure 7d), however, the power factor was considerably improved (Figure 7c). Indeed, pristine GeTe exhibited a peak PF ~3.3 x 10 -3 W/m.K 2 at 673 K, while it substantially increased for La-doped GeTe (PFmax ~4 x 10 -3 W/m.K 2 at 673 K) and Zr-doped GeTe (PFmax ~4 x 10 -3 W/m.K 2 at 723 K). From the band structure computations, the band gap at L-point increased from 226 meV for the pristine c-GeTe to 249 meV for the La-doped GeTe (Figure 8c). In the case of La-doping, the energy difference between the two valance band maxima remained unchanged when compared to that of pristine c-GeTe (69 meV). The impurity states arising from La are located far from EF and seem to be buried deep within the conduction band of GeTe. Similar to that of La, doping Zr slightly opens up the band gap from 226 meV to 255 meV (Figure 8d). The energy difference between the two valence band maxima is also reduced to 53 meV (for 0.5 at.% Zr-doping) from 69 meV (undoped GeTe). Refer SI (Table S3) for the data extrapolation for 0.5 at.% Zr and for the band structure computed with a larger super-cell (Figure S3). This marginal convergence of the valence bands explains the improved PF in Zr-doped GeTe. As both La-doped GeTe and Zr-doped GeTe compounds exhibited zT > 1, they were shortlisted for further experimental studies. The upcoming sections in this paper will exclusively focus only on Zr-doped GeTe. The optimizations and codoping strategies with La-doped GeTe will be reported in future communication.

Ge-deficient vs stoichiometric Zr-doped GeTe

The next stage was aimed at optimizing the composition of Zr-doped compounds. 0.5 at% Zr-doping to GeTe (Ge0.995Zr0.005Te) boosted zT up to ~1.1 at 723 K. To find the optimum content of Zr to achieve the best zT in the GeTe-solid solution, experiments were carried out with varying Zr-content in the nominal composition Ge1-xZrxTe (0.005 ≤ x ≤ 0.015). Indeed, the x = 0.005 composition that exhibited the best zT in that series, as the addition of more Zr content did not have any significant improvement in the TE properties of GeTe. From the band structure computations, it is found that Zr-doping progressively increased the optical band gap and, at the same time, considerably reduced ELΣ (Table S3, SI). However, the effect of a large increase in Eg (thus lower σ) with higher Zr-content has negated the benefits of the converging trend of the valence bands. As altering the Zr-content did not yield any TE improvement, the next strategy was formulated with an aim at tuning the Ge-content to optimize the composition. By controlling the Ge-vacancy levels (by nearly making it vacancy-free), recent reports have shown some significant improvement in the TE performance in GeTe, thanks to the suppression of hole concentration and recovery of high carrier mobility. 50,51,107 In particular, Dong et al.,50 has effectively demonstrated this strategy in Bi-doped Ge-excess Ge1+xTe, where these excess Ge acted as a source to subdue the formation of Ge vacancies. In this work, a reverse strategy of Ge-deficient Ge1-x+yZryTe compositions, where the aliovalent Zr can balance out or modulate the effects caused by Gedeficiency was explored. For this purpose, the self-compensated samples of Ge1-xZr0.005Te (x = 0.01 -0.03) were prepared. The self-compensated Ge0.98Zr0.005Te (x = 0.02 in Ge1-xZr0.005Te) was found to be the optimized composition. The TE transport properties of this Ge-deficient Zr-doped GeTe composition (Ge0.98Zr0.005Te) is compared with that of its stoichiometric counterpart (Ge0.995Zr0.005Te), as shown in Figure 9. Creating some small Ge-deficiency in the Zr-doped GeTe compound did not alter much the electrical transport properties, σ and S (Figure 9a-c) when compared to stoichiometric Ge0.995Zr0.005Te compound. To understand this feature, the electronic structure and defect calculations (based on DFT) on GeTe, both on rhombohedral and cubic models that were meticulously computed and reported by Edwards et al.,99 must be carefully analyzed. The simple model presented there for a Ge vacancy predicted one lesser s band and four fewer electrons in the valence bands. 99 The removal of each Ge atom will result in one empty valence band. The reported projected DOS near the Ge vacancy (for both rhombohedral and cubic arrangements) showed no new features in the gap. The new features that are rather present were buried beneath the valence band edge in the Te s and p bands. Since these defect states are found to be hyper-deep, the Ge-deficiency induced vacancies in GeTe are electrically and paramagnetically dormant. This can explain the fairly similar electrical transport properties between the Ge0.98Zr0.005Te and Ge0.995Zr0.005Te compounds. Moreover, the empty states resulting due to these Ge vacancies are not localized and they are virtually unperturbed states from the top of the valence band. The inconsequential defect relaxation energy was also reported to agree with the trivial charge localization around the vacancy, 99 and those empty states at the top of the valence band arising due to the inherently present Ge-vacancies, i.e., intrinsic defects in GeTe (due to their thermodynamic nature), justifies the p-type metallic conduction in GeTe. This is further corroborated by the DFT calculations reported by Bayikandi et al.,51 where the Fermi level moved down progressively into the valence band with increasing Ge vacancy levels; and in those Ge-vacancy models, the DOS contributions to the valence band were dominated by Te-5p 4 and Ge-4s 2 orbitals, while the majority contribution to the conduction band minima came from the Ge-4p 2 orbital. Imparting some Ge-deficiency in the Zr-doped GeTe composition, though did not affect much the overall electrical transport properties, i.e., PF, it has influenced their thermal transport by notably decreasing the thermal conductivity (Figure 9d), especially the lattice contribution (Figure 9e). Such similar reduction in thermal transport with the introduction of cation deficiency has also been observed in the past for n-type Pb-deficient Pb0.98-xSbxTe compounds. 108,109 To explain this phenomenon in reduction in κ with the Ge-deficient Zr-doped Ge0.98Zr0.005Te when compared to the stoichiometric Ge0.995Zr0.005Te, their structural results must be carefully introspected and compared. Le Bail fitted PXRD patterns of Ge-deficient Ge0.98Zr0.005Te and stoichiometric Ge0.995Zr0.005Te compounds are shown in Figure 10 and Figure S3 (SI), respectively, and their refined lattice parameters are shown in Table S4 (SI). The a and c lattice parameters of the two compositions are very similar. The major difference between the two diffraction patterns come from the peak shape of the Bragg reflections normal or almost normal to the [001] direction. In order to investigate this effect, several well-defined and non-overlapping Bragg peaks were fitted individually with a pseudo-Voigt function, Vg(x). This function consists of a linear combination, Vg(x) = η L(x) + (1-η) G(x), where L(x) and G(x) are the Lorentzian and Gaussian functions, while η is the refined mixing parameter and it shifts the profile more towards pure Gaussian or pure Lorentzian when approaching 0 or 1, respectively. Figure . 11 shows η and the full width at half maximum (FWHM) for several Bragg peaks of Ge0.98Zr0.005Te and Ge0.995Zr0.005Te patterns. In the case of stoichiometric Ge0.995Zr0.005Te, all the Bragg peaks were well fitted with η close to 1, which corresponds to pure Lorentzian peak shape.

The result is similar in the case of Ge-deficient Ge0.98Zr0.005Te, except for the reflections ( 003), ( 104) and ( 105) which are better fitted with η abnormally larger than 1 corresponding to a 'super-Lorentzian' profile. Such index-dependent super-Lorentzian peak shapes are a clear indication of coherent planar defects in crystalline solids 110 . In the current scenario, other than the deformation of the peak shape, no notable dependence of the FWHM with the Miller indices is observed. In any case, the 'super-Lorentzian' profile for Ge-deficient Ge0.98Zr0.005Te can be attributed to the presence of the planar-defects in GeTe, as already mentioned above. Lee et al. 102 showed using extensive TEM analyses that planar defects in GeTe are caused by the phase transition from the low-temperature rhombohedric to the hightemperature cubic phase. The phase transition can be described as an elongation of cubic unit cell, which can occur along with one of the four possible and equiprobable [111] directions, which will result in a complex 3D herringbone microstructure. 102 The diffraction pattern of such a complex crystal with a variety of defects is difficult to be simulated, as the current advanced programs can only handle 1D planar displacement [FAULTS program]. 110,111 The deformation of Bragg peaks on the Ge-deficient Ge0.98Zr0.005Te diffraction pattern seems to indicate a high density of defects for this composition. As the synthesis conditions were the same between these two samples (Ge-deficient Ge0.98Zr0.005Te and stoichiometric Ge0.995Zr0.005Te), it seems likely that the higher density of defects in Ge0.98Zr0.005Te is certainly related to the formation of Ge vacancies. Figure 11 clearly evidences the presence of the 'super-Lorentzian' profile only for the Ge-deficient Ge0.98Zr0.005Te composition, pointing that Gedeficiency has readily favored the formation of the larger density of planer defects in GeTe, than it would normally occur in the stoichiometric cases. This explains the lower lattice thermal conductivity in the Ge-deficient Zr0.005Ge0.98Te when compared to that of its stoichiometric counterpart over the whole measured temperature range. The Ge-deficiency induced vacancy domains can create a barrier or hinder the flow of heat-carrying phonons of different mean free paths (MFP). The contrary approach adopted by Dong et al.,50 where they used excess Ge (Ge1+xTe) to suppress these vacancies resulted in an increasing tendency of κlatt with excess Ge-content, thus the phonon scattering was reportedly weakened due to the lack of Ge vacancies and the resulting additional Ge precipitates also affected κ (as κ of pure Ge is much higher). The defect/vacancy induced a reduction in thermal conductivity has helped in conspicuously increasing zT in the Ge-deficient composition (Figure 9f). The Ge-deficient Ge0.98Zr0.005Te compound exhibited a peak zT ~1.3 at 673 K, which is significantly higher when compared to that of the stoichiometric Ge0.995Zr0.005Te composition (zTmax ~1.1 at 723 K) and to that of pristine GeTe (zTmax ~0.85 at 723 K). This strategy of intentionally creating more Ge vacancies may contradict the conventional approach (which is rather aimed at suppressing those intrinsic defects), it seems to work well for such transition metal-doped GeTe. These results further reinforce the very recent findings by Shuai et al.,56 where such an opposite directional approach by decreasing the Ge vacancy formation energy (i.e., the formation of more Ge vacancies and Ge precipitates) via another transitional metal (Cr) substitution at Ge site has also resulted in an improvement in zT.

TE properties of Zr and Sb codoped GeTe

To the optimized Ge-deficient Ge0.98Zr0.005Te composition, 10 at.% Sb was codoped to form the Gedeficient Zr-Sb codoped compound (i.e., Ge0.88Zr0.005Sb0.1Te). An improved band degeneracy has been widely reported in 8-10 at.% Sb doping in GeTe. 35,47,81,103 Hence Sb was codoped to the optimized Zrdoped GeTe intending to realize coadjutant synergistic band effects that would play a key role to ultimately boost the TE performance. The PXRD pattern of the Zr-Sb codoped Ge0.88 Zr0.005Sb0.1Te compound is presented in Figure 12. Le Bail lattice refined parameters (Table S5, SI) reveal that, with the insertion of Sb, the lattice parameter a increases from 4.166 Å for Ge0.98Zr0.005Te to 4.199 Å for Ge0.88Zr0.005Sb0.1Te. The other lattice parameter, c, however decreases from 10.661 Å for Ge0.98Zr0.005Te to 10.530 Å for Ge0.88Zr0.005Sb0.1Te. These changes in the lattice parameters have resulted in an expansion of the hexagonal unit-cell volume from 160.3 Å 3 for Ge0.98Zr0.005Te to 160.8 Å 3 for Ge0.88Zr0.005Sb0.1Te, which can be explained by Ge 2+ being substituted by larger Sb 3+ ions. It must be noted that the c/a ratio for pristine GeTe and Ge0.98Zr0.005Te are ~2.56 and ~2.55, respectively, while it significantly reduces for Ge0.88Zr0.005Sb0.1Te (c/a ratio ~2.507). Such variations in the crystal field effect can notably influence the TE properties (by tuning the electronic band structures), as recently reported by some of us in the case of Ti-Bi codoped GeTe. 57 The characteristic rhombohedral double reflections [( 024) and ( 220)] in the range of 2θ values between 41 o to 44 o , got closer and almost merged when Sb was partially substituted for Ge in Ge0.98Zr0.005Te, as can be seen from the PXRD pattern in Figure 12. The evolution of the unit-cell parameters can be visually confirmed in the inset of Figure 12, which shows the displacement of the ( 105) and ( 006) reflections toward higher angles and ( 113) and ( 202) toward lower angles for Ge0.88Zr0.005Sb0.1Te when compared to the Ge0.98Zr0.005Te. Besides, an important broadening of the ( 006) and ( 105) peak shapes can also be observed (inset of Figure 12). This would indicate the formation of a much larger density of planar defects with the insertion of Sb in Ge0.98Zr0.005Te, which could play a critical role in influencing their thermal transport properties. However, no further in-depth peak profile analysis could be performed because of the overlapping of the important reflections.

The TE transport properties of the Zr-Sb codoped Ge0.88 Zr0.005Sb0.1Te are presented and compared with Ge-deficient Ge0.98Zr0.005Te and pristine GeTe in Figure 13. The addition of Sb to Zr-GeTe extensively decreased the electrical conductivity, which is due to the suppression of the charge carriers arising due to the aliovalent donor dopant effect of Sb in GeTe. However, this reduction in σ is recompensated by the drastic rise in the S values (Figure 13b). The Zr-Sb codoped Ge0.88Zr0.005Sb0.1Te sample exhibited a high thermopower of 125 µV/K at room temperature and a maximum of 240 µV/K at higher temperatures, a notable improvement by 75% when compared to that of pristine GeTe and/or Zr-doped GeTe. First-principles calculations (DFT) can corroborate the reasons for this impressive rise in thermopower with Zr-Sb codoping. Computations with possible arrangements for the distribution of Sb atoms were carried out. However, the differences in the band structures between those structural arrangements were found to be small, and they showed consistent results (Figure 14 and Figure S5, SI). Quite significantly, Zr-Sb codoping in GeTe has resulted in the convergence of electronic band valleys by markedly shrinking the energy separation between the light hole and heavy hole valence bands, ELΣ, to 21 meV, a sharp reduction by 70% when compared to the ELΣ of pristine c-GeTe, as shown in the band structure sketched in Figure 14. Codoping of Zr and Sb modifies the electronic band structure of GeTe by activating hole pockets at the bottom of the conduction band (refer to the DOS in Figure 14a). Zr-Sb codoping has resulted in adding donor states just below the conduction band. The exact position of EF in these calculations should be slightly considered with caution, as these computations were based on stoichiometric GeTe (i.e., without any Ge-vacancy). In any case, the band convergence would favor the raise in effective mass, m* owing to the contributions from both the light and heavy hole bands (L and Σ bands) from the heavier effective mass near the EF. 51,52 This phenomenon in a reduction in energy separation ELΣ by tuning the crystal field effect with the decrease of c/a ratio for Zr-Sb codoped is something that has also been observed for Ti-Bi codoped GeTe. 57 There is a slight change of trend or tendency with the electrical and thermal transport properties in Zr-Sb codoped Ge0.88 Zr0.005Sb0.1Te at ~573 K, which is an indication of the structural transition 𝑅3𝑚 → 𝐹𝑚3 ̅ 𝑚. It must be noted that this structural transition for pristine GeTe occurs at around 673 K. Hence, codoping of Zr and Sb has stupendously reduced the phase transition temperature and readily increased the symmetric nature, which in turn will facilitate higher density of states effective mass (m*DOS) and weighted mobility. The Zr-Sb codoped Ge0.88Zr0.005Sb0.10Te sample exhibited  1.5 W/m.K at room temperature, which accounted for a massive reduction by 360% when compared to that of pristine GeTe (  7 W/m.K at room temperature) and reduction by 300 % when compared to Ge-deficient Ge0.98Zr0.005Te composition ( 6 W/m.K at room temperature). More importantly, the thermal conductivity remained low (between 1.5 -1.2 W/m.K) throughout most of the measured temperature range. This ultra-low  can be attributed to the synergistic suppression of both electrical (e) and lattice (latt) contribution to the thermal transport. The increased electrical resistance with the substitution of aliovalent Sb 3+ for Ge 2+ (and hence decreased carrier density), resulted in suppressed electronic contribution. The large density of planar defects, predicted by the PXRD results (Figure 12 and its inset), must have suppressed the transport of thermal phonons. This is consistent with a recent report on Sb-doping in GeTe, 107 where the presented HRTEM results showed a plethora of such planar defects, ranging from twinned and inter-twined planar domains, thickened herringbone domain boundaries, presence of stacking faults, interstitials all of which contribute to the scattering of phonons of a wider spectrum of mean-free paths. Hong et al. have established that Sb alloying in GeTe will alter the phonon dispersions, including the closing of the acoustic-optical phonon band gap, reducing the phonon group velocities and increasing the phononphonon scattering rates. 81 Importantly, Lee et al. observed 'colonies' of stacked herringbone domain structures in Sb-doped GeTe, 102 which are coherent with the large density of planar defects predicted from our XRD results (based on the anisotropic peak broadening and peak shifts of some particular reflections) on Zr-Sb codoped GeTe. The commingled effects of band structure modulations to markedly improve the thermopower and the suppression of both electronic and lattice contribution to stupendously reduce the thermal conductivity in Zr-Sb codoped Ge0.88Zr0.005Sb0.10Te sample have enabled it to manifest a high zTmax ~1.8 at 723 K. This zT value is on par with or better than some of best reported GeTe materials that were reportedly codoped with Bi-Sb, 47 Mn-Bi, 82 Mn-Sb, 63 Ga-Sb, 35 Ti-Bi, 57 etc. However, it must be mentioned that the recent improvisations (especially the symmetrybreaking approach) have reportedly yielded zT ≥ 2 in some of the codoped GeTe compounds. 48,77,81,103,107 According to the DFT calculations, doping of Sb reportedly decreases the energy of mixing in GeTe (i.e., the formation energy of the alloy), 35,47,81,103 and this provides scope for increased solubility of Zr in such Sb-codoped GeTe compositions. This renders further avenues for future work to explore more and optimize the Zr-content in these Zr-Sb codoped GeTe compositions to reach zT > 2.

Conclusion

To summarize, combined experimental and theoretical methods allowed us to extensively screen a family of transition and rare-earth elements as potential dopants for GeTe compounds for thermoelectric application. Each of their effects is summed up as follows,

(1) The presence of the cubic GeTe secondary phase affected the thermal conductivity in Nb-doped GeTe.

(2) Ru-doping produced several intruder states in the band gap of GeTe and profoundly affected the electrical transport properties. The presence of the secondary GeTe cubic phase also severely affected the thermal conductivity. The combined detrimental effect to both electrical and thermal transport properties with Ru-doping resulted in the least TE performance (least zT) among the transition and rareearth materials that are reported in this work.

(3) Mo-doping produced intruder states in the band gap of GeTe and also exhibited some magnetic character. Overall, PF was not affected much when compared to pristine GeTe, however, κ was increased.

(4) Pr-doping, though, did not influence much the electrical transport properties, adversely affected κ due to the presence of a large proportion of secondary GeTe cubic phase.

(5) V-doping produced flat and localized intruder states within the gap, and also distorted the density of states (induced resonance states adjacent to the Fermi level). These band modifications coupled with some magnetic characterall of which are known to suppress the charge carrier mobilityresulted in high resistivity and hence low PF. Unit cell volume also shrank with the substitution of V for Ge.

(6) Y-doping did not change much the electronic band structure of GeTe, and hence comparable PF are measured. But the κ was reduced, thus resulting in a marginal improvement in zT.

(7) Hf-doping tends to marginally decrease the energy separation between the two valence bands, thus slightly boosting PF and overall zT.

(8) La-doping neither changed κ nor the energy separation between the two valence bands, but the band gap was slightly opened up, and PF (hence overall zT) was slightly improved.

(9) Zr-doping did not change κ, but opened up the band gap and decreased the energy separation between the two valence bands, resulting in improvement in PF and overall zT.

(10) Creating Ge-deficiencies (vacancies) and modulating it with Zr-doping resulted in a large density of planar defects (deduced by the super-Lorentzian profile from XRD), leading to a significant reduction in κ and a marked improvement in zT from 0.85 (for pristine GeTe) to 1.3 at 673 K for Ge-deficient Zrdoped GeTe. Moving forward, this novel unconventional method of producing more electrically dormant Ge-vacancies and balancing it with an aliovalent metallic dopant as a strategic route for producing highly efficient GeTe-based TE materials via suppression of thermal transport properties without noticeably affecting the electrical transport properties.

(11) Zr-Sb codoping in GeTe resulted in an improved band degeneracy and facilitated the convergence of electronic band valleys (thus improving the Seebeck) with the tuning of the crystal field effect (reduction in c/a ratio) and suppressed the thermal transport, thus resulting in an exceptionally high zT ~ 1.8 at 723 K, comparable with some of the best reported high-performance codoped GeTe materials. and5, apply also here.

Table of Contents / Graphical Abstract

Reports an improvement in thermoelectric performance in GeTe via an intentional creation of more electrically dormant Ge-vacancies (in contrast to the classical approaches) and modulating/balancing it with optimized Zr-doping.

Figure 1 .

 1 Figure 1. Room temperature PXRD patterns for pristine GeTe, Ge0.08Y0.02Te, Ge0.995Hf0.005Te, Ge0.995Mo0.005Te, Ge0.98V0.02Te, Ge0.995Zr0.005Te, and Ge0.98La0.02Te. Simulated patterns for rhombohedral GeTe (main phase) and cubic Ge (impurity) are shown at the bottom of the plots.

Figure 2 .

 2 Figure 2. Room temperature PXRD patterns for Ge0.98Ru0.02Te, Ge0.98Pr0.02Te and Ge0.98Nb0.02Te. Simulated patterns for GeTe (cubic and rhombohedral) and Ge are shown at the bottom of the plots. The additional 2θ peaks around 66 o for Nb-GeTe came from the sample holder (insufficient powder during XRD sample preparation).

Figure 3 .

 3 Figure 3. Temperature-dependent (a) electrical conductivity (σ), (b) Seebeck coefficient (S), (c) power factor (PF = S 2 σ), (d) total thermal conductivity (), and (e) figure of merit (zT) for Ge1-xMxTe (x = 0.00, 0.02, 0.005; M = Nb, Y, Hf, Ru, Mo, Pr, and V) compounds.

Figure 4 .

 4 Figure 4. Calculated DOS for (a) NbGe63Te64, (b) RuGe63Te64, (c) MoGe63Te64, and (d) VGe63Te64 cubic models. The DOS of each model is compared with that of the pristine c-GeTe (c-Ge64Te64). The Fermi level (EF) of pristine GeTe is set arbitrarily at 0 eV. The dashed line represents the shifted EF for the doped compositions. Additional Gaussian smearing of 25 meV was applied and the projected DOS (for dopants) were magnified for better readability of the curves.

Figure 5 .

 5 Figure 5. (a) Brillouin zone of c-GeTe. Band structures of cubic models using a 4 × 4 × 4 supercell showing band folding in the → K () direction for (b) pristine c-Ge64Te64, (c) NbGe63Te64, (d) RuGe63Te64, (e) MoGe63Te64, and (f) VGe63Te64. The dopant (M = Nb, Ru, Mo, and V) projections are highlighted in green. The line thickness is proportional to the projection of the wave function on the orbitals (in green) of the respective dopants (M = Nb, Ru, Mo, and V).

Figure 6 .

 6 Figure 6. First-principles computed DOS for (a) YGe63Te64 and (b) HfGe63Te64, and computed band structures for (c) YGe63Te64 and (d) HfGe63Te64 cubic models. Additional figure captions (regarding EF, additional Gaussian smearing, and projection of the wave function) given in Figures4 and 5, apply also here.

Figure 7 .

 7 Figure 7. Temperature-dependent (a) electrical conductivity (σ), (b) Seebeck coefficient (S), (c) power factor (PF = S 2 σ), (d) total thermal conductivity (), (e) figure of merit (zT) for Ge0.995Zr0.005Te and Ge0.98La0.02Te compounds.

Figure 8 .

 8 Figure 8. First-principles computed DOS for (a) LaGe63Te64 and (b) ZrGe63Te64, and computed band structures for (c) LaGe63Te64 and (d) ZrGe63Te64 cubic models. Additional figure captions (regarding EF, additional Gaussian smearing, and projection of the wave function) given in Figures4 and 5, apply also here.

Figure 9 .

 9 Figure 9. Temperature-dependent (a) electrical conductivity (σ), (b) Seebeck coefficient (S), (c) power factor (PF = S 2 σ), (d) total thermal conductivity (), (e) lattice thermal conductivity (latt), and (f) figure of merit (zT) for Ge-deficient Ge0.98Zr0.005Te in comparison with stochiometric Ge0.995Zr0.005Te and pristine GeTe. The Ge deficient Ge0.98Zr0.005Te compound is shown exhibiting superior a zTmax ~1.3 at 673 K.

Figure 10 .

 10 Figure 10. Le Bail refined room temperature XRD pattern of Ge0.98Zr0.005Te. The experimental pattern is plotted in red (symbol), the calculated pattern in black (line), and the difference in blue (line). The vertical ticks indicate the Bragg positions of rhombohedral GeTe phase (black) and secondary cubic Ge (red) phase.

Figure 11 .

 11 Figure 11. Refined η and FWHM obtained by individually fitting several Bragg peaks of stoichiometric (a) Ge0.995Zr0.005Te and (b) Ge-deficient Ge0.9Zr0.005Te diffraction patterns with a pseudo-Voigt function. The 'super-Lorentzian' profile with η significantly larger than 1 signifies the presence of a large density of planar defects in the Ge-deficient Ge0.98Zr0.005Te compound.

Figure 12 :

 12 Figure 12: PXRD pattern of Zr-Sb codoped Ge0.88Zr0.005Sb0.1Te compound. The inset is the magnified view of the Ge0.88Zr0.005Sb0.1Te (red) and Ge0.98Zr0.005Te (black) patterns showing the peak shifts and broadening with the insertion of Sb in Zr0.005Ge0.995Te.

Figure 13 .

 13 Figure 13. Temperature-dependent (a) electrical conductivity (σ), (b) Seebeck coefficient (S), (c) power factor (PF = S 2 σ), (d) total thermal conductivity (), and (e) figure of merit (zT) for Zr-Sb codoped Ge0.98Zr0.005Sb0.1Te in comparison with Ge-deficient Ge0.98Zr0.005Te and pristine GeTe.

Figure 14 .

 14 Figure 14. First-principles computed DOS (a) and band structures highlighting (b) Zr projections (in green) and (c) Sb projections (in red) for the Sb6ZrGe67Te64 cubic model. Additional figure captions (regarding EF, additional Gaussian smearing, and projection of the wave function) given in Figures 4 and 5, apply also here.

Acknowledgment

B.S and S.L.T acknowledge the JSPS postdoctoral fellowships P19720 and PE19749, respectively. Support from JSPS KAKENHI 19F19720, JP16H06441 and JST-Mirai JPMJMI19A1 project grants are also gratefully acknowledged.