Interplay between Electronic, Magnetic and Structural Degrees of Freedom in Crystalline Transition Metal Compounds

TESE DE DOUTORAMENTO INTERPLAY BETWEEN ELECTRONIC, MAGNETIC AND STRUCTURAL DEGREES OF FREEDOM IN CRYSTALLINE TRANSITION METAL COMPOUNDS Adolfo Otero Fumega ESCOLA DE DOUTORAMENTO INTERNACIONAL DA UNIVERSIDADE DE SANTIAGO DE COMPOSTELA PROGRAMA DE DOUTORAMENTO EN CIENCIA DE MATERIAIS SANTIAGO DE COMPOSTELA 2021 DECLARACIÓN DO AUTOR/A DA TESE Interplay between Electronic, Magnetic and Structural Degrees of Freedom in Crystalline Transition Metal Compounds D./Dna. Adolfo Otero Fumega Presento a miña tese, seguindo o procedemento axeitado ao Regulamento, e declaro que: 1) A tese abarca os resultados da elaboración do meu traballo. 2) De selo caso, na tese faise referencia ás colaboracións que tivo este traballo. 3) A tese é a versión definitiva presentada para a súa defensa e coincide coa versión enviada en formato electrónico. 4) Confirmo que a tese non incorre en ningún tipo de plaxio doutros autores nin de traballos presentados por min para a obtención doutros títulos. En Santiago de Compostela, 13 de Maio de 2021 Asdo. ........... AUTORIZACIÓN DOS DIRECTORES / TITOR DA TESE Interplay between Electronic, Magnetic and Structural Degrees of Freedom in Crystalline Transition Metal Compounds D./Dna. Víctor Pardo Castro D./Dna. José Francisco Rivadulla Fernández INFORMA/N: Que a presente tese, correspóndese co traballo realizado por D/Dna. Adolfo Otero Fumega, baixo a nosa dirección, e a utorizamos a súa presentación , considerando que reúne os r equisitos esixidos no R egulamento de Estudos de Doutoramento da USC, e que como directores desta non incorre nas causas de abstención establecidas na Lei 40/2015. De acordo co indicado no Regulamento de Estudos de Doutoramento, declaramos tamén que a presente tese de doutoramento é idónea para ser defendida en base á modalidade Monográfica con reproducción de publicaciones, nos que a participación do/a doutorando/a foi decisiva para a súa elaboración e as publicacións se axustan ao Plan de Investigación. En Santiago de Compostela, 13 de Maio de 2021 Abstract Interplay between Electronic, Magnetic and Structural Degrees of Freedom in Crystalline Transition Metal Compounds Adolfo Otero Fumega i Adolfo Otero Fumega ii Crystalline Transition Metal Compounds (TMC) show a diverse amount of physical properties. They can be good metals, semiconductors or insulators, and they can develop intriguing ordered phases such as ferromagnetism or superconductivity. The partially filled d electrons provided by the transition metal atom play a fundamental role to determine these. It is the subtle interplay between their electronic, magnetic and structural degrees of freedom that leads to different ground state properties and phase transitions. This heterogeneous physical behaviour makes TMC highly interesting from a technological point of view, but also from a more fundamental perspective. Therefore, this kind of compounds represent an ideal platform to study emergent phenomena in Condensed Matter Physics. In this thesis we have studied from a theoretical point of view the interplay between the electronic, magnetic and structural degrees of freedom in TMC. The first three chapters are devoted to explain the theoretical aspects and methods required to perform the different studies in the subsequent chapters. In what follows, we will summarize the content and results of the chapters presented in this thesis. In Chapter 1 we provide the theoretical aspects to describe materials at the nanoscale. Nuclei and electrons are the fundamentals particles that determine the properties of materials. An unsolvable many-body interacting Hamiltonian is presented as the corner stone to describe a real material. Therefore, different approximations are required to make the problem tractable. We start with the Born-Oppenheimer and adiabatic approximations that allow to separate the analysis for nuclei and electrons. The multielectronic problems is discussed in the framework of the Density Functional Theory (DFT). This theory is developed upon the Hohenberg-Kohn theorems, where the electronic density is demonstrated to uniquely determine the properties of a material. DFT is made practical through the Kohn-Sham approach. An independent electron system is used to mimic the interacting one. However, in this formulation the electron correlation, introduced via the so-called exchange-correlation term, remains elusive. Hence, approximations to this term are required. We provide definitions of different exchange-correlation functionals. After that, we present the methods used to solve the Kohn-Sham equations, from full-electron localized-orbital approaches to pseudopotential plane-wave methods. Moreover, simplified models, like spin or tight-binding Hamiltonians, are introduced at this step. Regarding the nuclear motion, a classical lattice dynamic approach is introduced. We describe first the harmonic approximation to obtain the phonon spectrum of a crystal. After that, we see how to compute anharmonic effects. This will be essential to correctly describe the lattice heat transport. Different methods to compute the lattice dynamics are exposed. Therefore, in this first chapter we pave the way to study from a theoretical point of view any material. In Chapter 2, the fundamental concept of symmetry is analyzed. We start providing a definition of symmetry group. After that we describe the different symmetry groups that are encountered in a crystal, and in particular when dealing with TMC. A detailed exposition of different point groups and their effect on delectrons is given. Translational symmetry is presented as a fundamental symmetry in crystals. The notion of (electronic and phonon) energy bands is shown to be determined by this symmetry. Here, we provide iii Adolfo Otero Fumega x Os compostos cristalinos con metais de transición, TMC das súas siglas en inglés, amosan propiedades físicas moi diversas. Son dende bos metais ata semicondutores ou aillantes, e moitos de eles desenvolven fases ordenadas como o ferromagnetismo ou a supercondutividade. Os metais de transición achegan electróns dparcialmente ocupados, serán estes os que xoguen un papel fundamental á hora de determinar as propiedades destes materiais. A sutil interacción entre os graos de liberdade electrónicos, estruturais e magnéticos é a responsable de establecer diferentes estados fundamentais e de ocasionar transicións entre fases ordenadas. Este comportamento tan heteroxéneo fai que os TMC sexan altamente interesantes dende un punto de vista tecnolóxico, pero tamén dende unha perspectiva puramente teórica. Estes compostos representan unha excelente plataforma para estudar fenómenos emerxentes en Física da Materia Condensada. Nesta tese estudamos dende un punto de vista teórico a interacción entre os graos de liberdade electrónicos, estruturais e magnéticos nos TMC. Os primeiros tres capítulos están adicados a explicar os fundamentos e métodos teóricos necesarios para levar a cabo os estudos da segunda parte. No que segue resumiremos o contido e os resultados de cada un dos capítulos. No capítulo 1 amosamos como se describen os materiais na nanoescala. Electróns e núcleos son as partículas fundamentais que determinan as propiedades dos materiais a través das súas interaccións. Estás veñen descritas por un Hamiltoniano irresoluble de moitos corpos. Polo tanto, será preciso realizar numerosas aproximación para facer o problema tratable. Comezamos introducindo as aproximacións de Born-Oppenheimer e a adiabática. Estas permiten separar a análise para electróns e núcleos. O problema multielectrónico discutirémolo no marco da teoría do funcional da densidade, DFT do inglés. Está teoría baséase nos teoremas de Hohenberg e Kohn que demostran que a densidade electrónica determina de xeito unívoco as propiedades do material. O xeito de facer práctica esta teoría é por medio da aproximación de Kohn e Sham, que considera que o problema de electróns interaccionantes se pode mapear nun de electróns independentes. Non obstante, nesta formulación a correlación entre os electróns, incluída a través do termo de correlación e intercambio, é descoñecida e polo tanto serán precisas aproximacións a este termo. Amosaremos diferentes aproximacións ao funcional de correlación e intercambio. Despois presentaremos os métodos empregados para resolver as ecuación de Kohn-Sham, dende aproximacións con orbitais localizados e tódolos electróns ata ondas planas e pseudopotenciais. Ademáis, introduciremos neste punto un par de modelos simplificados: Hamiltonianos de espín e modelos tipo tight-binding. Con respecto á dinámica nuclear esta será aproximada dun xeito clásico. Describiremos primeiro a aproximación harmónica que permite obter o espectro de fonóns do cristal. Despois veremos como introducir efectos anharmónicos. Estes serán esenciais para describir axeitadamente o transporte de calor por fonóns. Amosaremos diferentes métodos para calcular a dinámica nuclear. Polo tanto, este primeiro capítulo sentará as bases teóricas no estudo de materiais. No capítulo 2 analizamos o concepto fundamental de simetría. Comezaremos dando unha definición de grupo de simetría para despois describir os diferentes grupos de simetría que se atopan nun cristal e, máis en particular, nos TMC. Proporcionaremos unha exposición detallada do efecto causado nos electróns dpor diferentes grupos puntuais de xi Adolfo Otero Fumega simetría. Presentaremos logo o simetría traslacional como unha das simetrías máis fundamentais dun cristal. Veremos que a noción de bandas de enerxía (electrónicas e fonónicas) emerxe desta simetría. Amosaremos a vista de páxaro diferentes características que podemos atopar nas bandas dun cristal. Seguidamente, definimos a simetría de inversión temporal conectándoa con diferentes ordes magnéticos. Por último, presentaremos unha breve introdución ao renovado tema de invariantes topolóxicas en materia condensada. No capítulo 3 estableceremos a conexión entre a nanoescala e a macro escala dun material por medio da mecánica estatística e as propiedades de transporte. Comezaremos definindo cantidades especias no equilibrio térmico como a densidade de estados ou a densidade de estados proxectada. As distribucións de Bose-Einstein e Fermi-Dirac para fonóns i electróns respectivamente serán definidas. Estas dan o número de estados ocupados por unidade de volume a unha enerxía dada para un sistema en equilibrio. A estabilidade dun sistema será discutida en termos da densidade de estados e as funcións de distribución. De seguido levaremos o sistema fóra do equilibrio para estudar as súas propiedades de transporte. Esto farémolo no marco da ecuación de transporte de Boltzmann, onde campos externos modificarán as funcións de distribución tratando de restrablecer o equilibrio nun tempo de relaxación. Este tempo relacionase cos diferentes procesos de collisións entre as partículas que tratan de restaurar o equilibrio. Para transporte electrónico aproximaremos o tempo de relaxación a unha constante, mentres que no caso de fonóns o tempo de relaxación será calculado por medio de termos anharmónicos. Todo isto permitiranos definir os seguintes coeficientes de transporte: condutividade eléctrica, condutividade térmica electrónica, coeficiente Seebeck e condutividade térmica fonónica. Tamén discutiremos a interacción que existe entre estes coeficientes para determinar como de bo termoeléctrico é un material. Despois destes tres capítulos nos que expoñemos a teoría e os métodos para resolver o problema de moitos corpos xa estamos en disposición de comezar os nosos estudos en TMC. En particular, centrarémonos en óxidos complexos e materias laminares de van der Waals. Estas clases de compostos presenta unha gran variedade de propiedades físicas e poden ser sintetizados en multitude de estruturas de forma sinxela. No caso dos óxidos complexos é salientable a súa gran estabilidade en condicións normais de presión e temperatura. No tocante aos materiais laminares de van der Waals, o feito de estar formados por capas permite exfolialos ata chegar ao limite en dúas dimensións. Polo tanto, estas familias de compostos presentanse como prometedores á hora de estudar novos fenómenos emerxentes a temperatura ambiente. Aparte de estes materiais, outro tipo de compostos serán analizados para mostrar o importante rol dos electróns d. Comezaremos no capítulo 4 estudando as propiedades de transporte no SrTiO3. Este composto representa un prototipo para o estudo de diferentes fenómenos físicos. Dexei- to remarcable, o SrTiO3dopado atópase dentro dos mellores materiais termoeléctricos. Este tipo de materiais son capaces de converter de forma eficiente un gradiente térmico nunha diferenza de potencial e viceversa, polo que son moi interesantes dende o punto de vista da innovación. Neste tipo de materiais é importante que o coeficiente Seebeck sexa maximizado coa finalidade de optimizar a diferenza de potencial xerada por un gradiente xii térmico dado. A maiores, tamén é importante que a condutividade térmica fonónica sexa a mínima posible para que así se manteña o gradiente de temperatura. Como recollemos no capítulo 3, estas propiedades de transporte veñen determinadas pola estrutura electrónica e cristalina. Polo tanto, o obxectivo deste capítulo será estudar dende unha perspectiva ab initio cales son as características fundamentais destes coeficientes de transporte e como poden ser alteradas por medio de distorsións estruturais. A primeira parte deste estudo estará adicada ao coeficiente Seebeck. Analizaremos o efecto da dexeneración orbital no coeficiente Seebeck como función de distorsión (compresión e tensión) da celda unidade e do nivel de dopaxe. Para facer isto construiremos un modelo simplificado para o coeficiente Seebeck fundamentado nas ecuacións de Mott e Heikes. Introduciremos un factor dependente da tensión que modelizará o efecto causado pola dexeneración dos orbitais t2g no coeficiente Seebeck. Poderemos calcular dende DFT este factor de dexeneración, desacoplando deste modo o efecto da tensión e a dopaxe debidas á entropía da dexeneración electrónica ou debidas a modificacións da estrutura de bandas. Atoparemos que as masas electrónicas efectivas e a dexeneración electrónica son os actores principais que determinan o coeficiente Seebeck como función da distorsión. Aplicando o modelo ao caso do SrTiO3 atopamos dous réximes de dopaxe. A baixas dopaxes calquer tipo de distorsión diminúe o valor do coeficiente Seebeck debido a que o termo dominante é o das dexeneracións, polo que calquera distorsión reduce o factor de dexeneración. A altas dopaxes, o réxime interesante en aplicacións, o termo enerxético domina e atopamos que aplicar unha tensión aumenta o coeficiente Seebeck. Isto é debido a que a tensión produce un aumento do volume da celda unidade que modifica de xeito favorable o ancho de banda e polo tanto a masa efectiva. Na segunda parte deste capítulo analizaremos a condutividade termica fonónica do SrTiO3. Faremos un escrutinio das principais causas do transporte térmico nas dúas fases estruturais do SrTiO3e como modificalas aplicado distorsións na estrutura. A altas temperaturas teremos unha fase cúbica e a baixas unha tetragonal con rotacións dos octaedros. Na fase cúbica, determinamos que as bandas acústicas e o modo polar triplemente dexenerado son os principais responsables do transporte de calor. Unha redución do volume da celda unidade aumenta a dispersión dos fonóns acústicos e a enerxía dos modos polares. Deste xeito incrementase a velocidade de grupo e redúcense as colisións fonónicas aumentando a condutividade térmica. O comportamento oposto ocorre se aumentamos o volume da celda unidade. No caso de aplicar unha distosión tetragonal a dexeneración dos modos polares rómpese baixando en enerxía o modo na dirección ce modificando a dispersión das bandas acústicas. Todo isto causa unha diminución do transporte térmico. Na fase a baixa temperatura a aparición de rotacións octaédricas proporciona novos graos de liberdade a considerar cando analizamos a condutividade térmica. Do mesmo xeito que ocorre na fase cúbica, os fonóns acústicos e polares son os responsables do transporte. Analizando de xeito independente o efecto do volume da celda unidade, o ratio c/a e o ángulo das rotacións octaédricas observamos que: i) un incremento do volume causa unha diminución da dispersión das bandas acústicas e baixa a enerxía dos modos polares, diminuindo a condutiviade térmica; ii) incrementar o ratio c/a reduce a dispersión das bandas acústicas e con elo a condutividade térmica; iii) atopamos que o ángulo das xiii Adolfo Otero Fumega rotación octaédricas está fortemente acoplado cos modos polares, un aumento do ángulo aumenta a enerxía dos modos polares diminuíndo as interaccións fonónicas e aumentando a condutividade térmica. Se consideramos conxuntamente as dúas análises para os coeficientes de transporte, a resposta termoeléctrica pode optimizarse incrementando o volume da celda unidade con tensión e reducindo as rotacións octaédricas. Os resultados deste capítulo proporcionan un coñecemento fundamental das propiedades de transporte de óxidos similares e guían a procura de mellores materiais termoeléctricos. No capítulo 5 analizeremos a fase ferromagnética e aillante que emerxe no LaCoO3 cando se medra sobre SrTiO3na dirección (001). Os cálculos ab initio que realizaremos pertimirán estudar a estabilidade termodinámica de estruturas con distintas estequiometrías nas que introduciremos vacantes de oxíxeno de acordo coas observacións experimentais. Seremos capaces de determinar o estado fundamental e a estrutura electrónica e propiedades magnéticas do LaCoO3medrado sobre SrTiO3. Atopamos que a estequiometría máis estable é LaCoO2,83; é dicir, un 6% de vacantes de oxíxeno. Así mesmo, determinamos que as vacantes forman cadeas perpendiculares á dirección (001); e que o parámetro de rede c e as distancias entre as capas de átomos de La están en consonancia coas medidas experimentais. Dende o punto de vista electrónico, veremos que os átomos de Co que se encontran no plano das vacantes son unha mistura de Co2+ con alto espín e Co2+ con baixo espín con acoplamentos ferromagnéticos. Os restantes cobaltos manteñense como Co3+ non magnéticos no entorno octaédrico. A configuración que atopamos da lugar a un momento magnético que concorda co experimental, así como a un orde ferromagnético. Con respecto ao comportamento aillante, observamos que un gap enerxético emerxe como consecuencia do desdobramento do campo cristalino inducido polas vancante de oxíxeno a un valor razoable de U. Polo tanto, neste traballo explicamos a fase ferromagnética e aillante dende unha perspectiva iónica moi sinxela. Estes resultados ofrecen unha maior comprensión da natureza dos aillantes ferromagnéticos e guían a procura de novos óxidos con esta inusual fase nos que as vancantes de oxíxeno puidesen xogar un papel relevante. No capítulo 6 analizaremos o efecto da dopaxe con átomos de Sr na estrutura electrónica e propiedades magnéticas dos niquelatos de capa infinita RNiO2(R= La, Nd). Esta familia de compostos levaba anos proposta de xeito teórico como análogo natural dos cupratos supercondutores a alta temperatura. Non obstante, a confirmación experimental deste feito non ocorreu ata vai un par de anos. Os chamados compostos pais; é dicir, sen dopar, dos cupratos e niquelatos amosan un comportamento electrónico e magnético radicalmente diferente. Mentres que os cupratos son aillantes antiferromagnéticos, os niquelatos son metais e non magnéticos. Sen embargo, ao introducir dopaxe de ocos (dopaxe con átomos de Sr) a supercondutividade emerxe nas dúas familias. Polo que a natureza dos ocos preséntase como un estudo a realizar. Neste traballo realizaremos cálculos DFT en RNiO2como función do catión R (R=La, Nd) e da dopaxe con Sr (11% e 25%). Empregaremos dous esquemas diferentes de contaxe dobre nos cálculos LDA+U. Os resultados amosarán cambios moi substanciais na estrutura electrónica destes niquelatos ao seren dopados con Sr. Por unha parte, o Sr reduce o efecto da chamada autodopaxe desprazando as bandas dda terra rara por riba do nivel de Fermi. Esto da lugar a unha descripción xiv cunha sóa banda onde o orbital dx2−y2do átomo de Ni é dominante. Ademáis, redúcese a enerxía transferencia de carga entre os átomos de Ni e O. Con respecto ás propiedades magnéticas, observamos que o estado de baixo espín para o Ni2+ é favorable. Polo tanto, todos estes resultados demostrans que dopando con átomos de Sr os niquelatos adquiren un comportamento semellante ao dos cupratos, suxerindo que a supercondutividade dos niquelatos presentará unha descripción como a dos cupratos. No capítulo 7 estudaremos a estrutura electrónica do Cs2CuCl4, un composto que se ten discutido no marco dos líquidos de espín xeometricamente frustrados. Empregaremos de xeito combinado técnicas experimentais, como a dispersión inelástica por resonancia de raios X (RIXS do inglés) e a espectroscopía por absorción de raios X (XAS do inglés), e cálculos para determinar o espectro deste comoposto. A análise mostra que a estrutura electrónica está deteminada na súa maioría polo desdobramento dos electróns ddo Cu2+ a causa do campo cristalino tetraédrico fortemente distorsionado xerado polos ións de Cl1−. Os cálculos DFT permiten determinar os diferentes desdobramentos enerxéticos observados por RIXS. A maiores tamén identificamos o gap de transferencia de carga entre o cluster de [CuCl4]2−e os estados sdo Cs1+. Este traballo amosa o potencial de combinar experimento e teoría e os resultados sentan as bases para un escrutinio máis detallado das exóticas propiedades deste composto a baixas temperaturas. Nos seguintes capítulos achegarémonos ao estudo de materiais laminares de van der Waals. Este tipo de compostos representa unha plataforma excelente para estudar propiedades en materials con baixa dimensionalidade. En particular, o magnetismo neste límite de baixa dimensionalidade resulta moi interesante polas prometedoras aplicacións que podería ter. No capítulo 8 presentaremos un estudo comparativo entre a estrutura electrónica e as propiedades magnéticas de dous compostos de van der Waals ferromagnéticos. Estudaremos o efecto da dimensionalidade en CrBr3e Cr2Ge2Te6dende un punto de vista experimental e teórico. Para esto, no canto de ir ao límite puramente bidimensional estudaremos a evolución aplicando presión; isto é, índonos a un límite máis tridimensional. Atoparemos que o composto de CrBr3amosa un comportamento de tipo Mott cun gap d−de se mantén aillante ao aplicar presión. De modo diferente, o Cr2Ge2Te6presenta un gap enerxético tipo p−de polo tanto unha descripción de aillante por transferencia de carga. Os nosos cálculos predín unha transición aillante-metal arredor de 6 GPa, que podería estar relacionada coa transición structural observada nos experimentos de raios X a alta presión. O diferente carácter da estrutura electrónica e conseguinte evolución coa presión destes compostos implica unha evolución diferente das súas propiedades magnéticas coa presión. Os nosos experimentos amosas unha redución da temperatura de Curie nos dous compostos ao aplicar presión, non obstante a súa tendencia é diferente. A nosa análise teórica complúe que a diferente evolución dos acoplamentos magnéticos dentro e fóra do plano para cada composto determina a diferente tendencia das temperaturas de Curie. Polo que este traballo resalta o importante rol da estrutura electrónica para explicar as interaccións entre espíns que causan a orde ferromagnética. En particular, salientamos o efecto da presión, e por tanto da dimensionalidade, na estrutura electrónica e nas propiedades magnéticas en materias de van der Waals baseados en Cr. xv Adolfo Otero Fumega No capítulo 9 seguiremos a estudar materiais de van der Waals que se presentan como prometedores candidatos para atopar orde ferromagnético de largo alcance no límite bidimensional. En concreto centrarémonos na análise de dous dicalcoxenuros d metais de transición (TMDs do inglés), CrTe2e VSe2. Moitos destes compostos desenvolven unha fase de densidade de carga ondulatoria, ou fase charge density wave (CDW) en inglés, a baixas temperaturas que interacciona co procurado orde ferromagnético. Polo tanto o obxectivo principal deste capítulo será estudar a interacción entre os ordes ferromagneticos e CDW nos TMDs. Na primeira parte centrarémonos no CrTe2con infinitas capas, un ferromagnético a temperatura ambiente cos momentos apuntando dentro dos planos de van der Waals. Despois, pasamos ao estudo do límite monocapa, é dicir bidimensional, e atopamos que unha fase CDW emerxe a baixas temperaturas. Proporcionamos unha análise detalla da nova simetría cristalina e como recoñecela experimentalmente. O fonón responsable desta transición é identificado o que abre a posibilidade para saltar dunha fase á outra activando o devandito modo. Na fase CDW ábrense pseudogaps enerxéticos arredor do nivel de Fermi e os orbitais d localízanse máis. Isto fai que os momentos magnéticos apunten fóra do plano, superando as restricións impostas polo teorema de Mermin e Wagner e dando lugar a orde ferromagnético de largo alcance na monocapa. Tamén atopamos que unha distorsión tensionante da monocapa causa un efecto similar, localizando os electróns d e levando os momentos a apuntar fóra do plano. Na segunda parte deste capítulo centrámonos no estudo do VSe2. Experimentalmente demóstrase que a forma multicapa deste composto desenvolve unha fase CDW por debaixo de 110 K e non se reporta ningún tipo de comportamento magnético. O noso estudo ab initio amosará que a fase CDW destrúe o magnetismo no VSe2multicapa. Esta minuciosa análise permite reconciliar o desacordo entre os experimentos e cálculos DFT anteriores nos que o efecto da CDW non se estaba a considerar; e que polo tanto predicían de xeito erróneo un comportamento ferromagnético neste material. A maiores, a estrutura CDW que obtivemos é capaz de explicar o salto atopado experimentalmente no coeficiente Seebeck á temperatura de transición. O estudo da monocapa continúa a ser un campo de investigación moi activo; sen embargo, neste traballo amosamos que a presenza dunha fase CDW afecta de xeito negativo ao ferromagnetismo. Polo tanto, neste capítulo salientamos o papel fundamental que xoga a estrutura CDW para determinar as propiedades magnéticas nos TMDs. Por último, no capítulo 10 presentamos un estudo teórico no BaSn2. Este composto laminar é un semimetal cunha liña nodal topolóxica cando se despreza o acoplamento espín orbital. Os materiais topolóxicos están a atraer moito interese debido a que presentan estados superficiais robustos que poderían dar lugar a fenómenos emerxentes e a prometedoras aplicacións. Neste estudo empregaremos ao BaSn2como prototipo para investigar transicións topolóxicas inducidas pola presión. Para elo, utilizaremos unha combinación de cálculos DFT e un modelo tight-binding. Atopamos que a presión uniaxial produce unha transitión de fase entre topoloxías non-triviais arredor de 4 GPa. Na fase a alta presión o número de liñas nodais topolóxicas aumenta. Nesta transición os electróns d xogan un papel fundamental, xa que é o aumento coa presión das enerxías de hopping xvi entre os orbitais dxz edyz do Ba cos pzdo Sn o que fai emerxer as novas liñas nodais. Ademáis, amosamos a forma de detectar experimentalmente este tipo de transicións mediante medidas de magnetotransporte. Os aspectos e métodos teóricos presentados nesta tese, xunto cos estudos nos distintos TMC, proporcionan unha mostra representativa do coñecemento fundamental que se está a levar a cabo na área de Física da Materia Condensada. xvii Adolfo Otero Fumega xviii Contents Introduction and Objectives xxiii I Theory and Methods 1 1 The Physical System 3 1.1 A Many-Body Quantum System . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 Nuclei are heavier than Electrons . . . . . . . . . . . . . . . . . . . . . . . 6 1.3 InteractingElectrons.............................. 9 1.3.1 Density Functional Theory . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.2 What about Electron Independence? . . . . . . . . . . . . . . . . . 13 1.3.3 Spin and Relativistic Effects . . . . . . . . . . . . . . . . . . . . . . 18 1.3.4 Electronic Structure Methods . . . . . . . . . . . . . . . . . . . . . 19 1.3.5 SimplifiedModels............................ 24 1.4 NucleiandVibrations ............................. 26 1.4.1 HarmonicApproach .......................... 27 1.4.2 Anharmonic Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 1.4.3 Nuclear Dynamics Methods . . . . . . . . . . . . . . . . . . . . . . 30 2 Symmetry 35 2.1 Definition of a Symmetry Group . . . . . . . . . . . . . . . . . . . . . . . . 35 2.2 PointGroupSymmetry............................. 36 2.3 Translational Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4 Time-Reversal Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 2.5 Topological Invariants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3 Statistics and Transport 49 3.1 ThermalEquilibrium.............................. 49 3.2 TransportProperties .............................. 52 xix Adolfo Otero Fumega in different structures. Importantly, this kind of compounds tend to be extremely stable in normal conditions due to the oxidized atmosphere where we live. In fact, many of the most intriguing oxides occur naturally. On the other hand, layered van der Waals materials are formed by layers that are bonded via van der Waals interactions. This weak bonding allows to grow and manipulate this kind of compounds in the few-layer limit. Even more, layers can be exfoliated to obtain a monolayer. This dimensionality control prompts to emergent phenomena. Considering all that, the main body of this thesis will be devoted to the study of this kind of materials. However, other classes of compounds will also be invited to exemplify the important role that d electrons play in nature. Part II will start with studies on oxides. Chapter 4 will be devoted to the analysis of the thermoelectric properties of SrTiO3, in particular we will see how the thermopower and the lattice thermal conductivity can be tuned by structural modifications. In Chapter 5 we will see how LaCoO3becomes an insulating ferromagnet when grown on top of SrTiO3. The effect of doping in infinite-layer nickelates will be analyzed in Chapter 6. Recently, this class of compounds has been found to be high-temperature-superconducting- cuprate analogs. In Chapter 7 we make a break before entering in the study of layered van der Waals materials. We show how the experimental spectrum of Cs2CuCl4, a compound discussed in the framework of spin liquids, can be explained in combination with the methodology presented in this thesis. A comparative study between the electronic and magnetic properties of two similar Cr-based van der Waals ferromagnets, CrBr3and Cr2Ge2Te6, is provided in Chapter 8. In Chapter 9 the interplay between ferromagnetic and charge density wave phases is discussed. Two different transition metal dichalcogenides are studied, CrTe2and VSe2. Finally, in Chapter 10 we analyzed how d electrons can play a role in topological phase transitions by studying BaSn2under pressure. The articles where these works have been published are referred in the List of Publications section, at the end of the thesis. All the theoretical aspects and methods presented in Part I, as well as the works described in Part II, will provide many of the basics to understand the state-of-the-art research in Solid State Physics. This thesis will delve into the fundamentals of Materials Science, and hopefully will serve to achieve some of the global goals that we have discussed in this introduction. xxvi Part I Theory and Methods 1 1 The Physical System “The more we simplify our material needs the more we are free to think of other things." —Eleanor Roosevelt In this chapter we will describe from a fundamental point of view the physical system that describes a material. We will show that different approximations are needed to deal with that system. Moreover, we will present the methods that we will use in the analysis of the second part. 1.1 A Many-Body Quantum System In the Introduction we have been talking about the prominent role materials have for our modern societies. However, we have not provided yet a formal definition of what a material is. In terms of its constituents a material is formed by atoms, which were first postulated by Democritus in Ancient Greece and finally discovered by the fathers of modern Chemistry Antoine Lavoisier and John Dalton. Thus, the building blocks of any material are the atoms of the different elements collected in the Periodic table. The way in which these elements are grouped leads to a plethora of materials with completely different properties that can be used for technological applications. It was at the beginning of the 20th century when it was discovered that atoms are formed by nuclei and electrons. This defines a length scale of nanometers to describe the structure of materials, the nanoscale. Electrons have negative electric charge and surround the nucleus, and this, in turn, is formed by neutral neutrons and positive protons. The number of protons in the nucleus, called the atomic number, determines the particular element of an atom. Many models were proposed to explain the atomic structure, i.e., how electrons and nuclei are arranged together. The Quantum Theory was developed with that purpose since electrons are a kind of object that cannot be interpreted in terms of the Classical theory. It was Schrödinger’s equation that first treated the electrons in a satisfactory way in terms of a wavefunction [7]. Then, Dirac provided a more natural description of spin by writing a relativistic equation [8]. Thus, materials at the nanoscale 3 Adolfo Otero Fumega must be described quantum mechanically. We have established that nuclei and electrons are the fundamental particles that constitute materials.1In a quantum mechanical system like this, physical states for n electrons and Nnuclei are codified by a vector of a complex Hilbert space, usually called wavefunction Ψ(r1,r2, ..., rn,R1,R2, ..., RN). This wavefunction depends on the coordinates riand RI(including position and spin) of each electron and nucleus respectively, and can be expressed as a sum of a complete set of sstates Ψs(r1,r2, ..., rn,R1,R2, ..., RN) weighted by a probability in that Hilbert space. The measurable quantities will be given by Hermitian operators. The results of a measurement and their probability to occur can be computed from the wavefunction and the eigenvectors of such operator. An essential case to consider is the Hamiltonian ˆ H, which corresponds to the total energy operator of a system. The eigenvalues and eigenvectors of the Hamiltonian provide the possible values for the energy Esand the states |Ψsithat the system can take: ˆ H|Ψsi=Es|Ψsi(1.1) The Hamiltonian is given by the sum of the kinetic and potential energies of all the particles. Therefore, in general, we can write the Hamiltonian of any material with n electrons and Nnuclei as ˆ H=− N X I=1 ~2 2MI∇2 I+1 2 N X I=1 N X I6=J ZIZJe2 |RI−RJ| −~2 2me n X i=1 ∇2 i+1 2 n X i=1 n X i6=j e2 |ri−rj|− n X i=1 N X I=1 ZIe2 |ri−RI|, (1.2) where we encounter 5 terms well differentiated. The first is the kinetic energy of the nuclei of mass MI, the second is the Coulomb interaction between the different nuclei of charge ZI, the third one is the kinetic energy of electrons with mass me, the fourth term is the Coulomb interaction between electrons and the final term is the Coulomb interaction between electrons and nuclei. Note that spin terms have been avoided in the many-body Hamiltonian written in eq. (1.2). These terms will be of capital importance for the correct description of many systems such as magnetic and topological ones. We will come back to them in the following sections. The important idea that we want to highlight at this point is that materials are described by many-body Hamiltonians, i.e., Hamiltonians that contain interactions between many different particles (electrons and nuclei). Since the wavefunction Ψsdepends on the set of coordinates {r}and {R}of nelectrons and Nnuclei respectively, we have an enormous amount of degrees of freedom coupled by the different terms in 1Note that the word “fundamental" is usually referred to the elementary particles described by the Standard model in High Energy Physics. However, the energy scale needed to describe materials is much lower, and hence it is fine to consider nuclei and electrons as the fundamental particles of our theory. 4 1 The Physical System Figure 1.1: Scheme with the theoretical approach to the many-body problem. The different parts will be explained in the main text. the many body-Hamiltonian of eq. (1.2). This represents a highly difficult problem to solve as soon as we have a system with more than two particles (like the hydrogen atom). A real material is a macroscopic system formed by a number of atoms of the order of Avogadro’s number (NAv ≃1023). This makes the problem impossible to solve. Hence, 5 Adolfo Otero Fumega we will need to find ways to simplify the many-body Hamiltonian. This will be the task of the following sections. In order to have at a glance all the approximations that we are going to performed, we have summarized them in the sketch of Fig. 1.1. 1.2 Nuclei are heavier than Electrons The first approximation that we will perform on the many-body Hamiltonian of eq. (1.2) is the so-called Born-Oppenheimer approximation [9]. Taking a close look to that Hamiltonian, we can observe that the nuclear kinetic energy terms are very small compared to the others, as long as the masses MIof the nuclei are much heavier than the electron mass me. This allows to treat the nuclear kinetic energy terms as a perturbation upon the Hamiltonian of eq. (1.2). The result leads to the decoupling of the electronic ({r})and nuclear ({R})degrees of freedom in the many-body wavefunction Ψs({r},{R}) = X s0 χss0({R})ψs0({r};{R}),(1.3) where we have expressed the coupled states in terms of a complete set of s0electronic states ψs0({r};{R})for each nuclear configuration in the set {R}. The χss0({R})are functions of the nuclear coordinates and the coefficients of the electronic states ψs0({r};{R}). Therefore, we can reduce the many-body (nuclear and electronic) Hamiltonian of eq. (1.2) to a purely electronic Hamiltonian ˆ He. We have to neglect the kinetic energy of the nuclei and consider the nuclear coordinates as parameters instead of variables. Then, the Hamiltonian reduces to ˆ He=ˆ Te+ˆ Vext +ˆ Vint +ENuc (1.4) where ˆ Teis the kinetic energy operator for electrons ˆ Te=−~2 2me n X i=1 ∇2 i,(1.5) ˆ Vext is the potential that acts on electrons due to the nuclei, in which positions are treated now as parameters ˆ Vext =− n X i=1 N X I=1 ZIe2 |ri−RI|,(1.6) ˆ Vint is the Coulomb interaction between electrons ˆ Vint =1 2 n X i=1 n X i6=j e2 |ri−rj|,(1.7) 6 1 The Physical System and finally ENuc is the Coulomb interaction between nuclei. It contributes to the total energy but is not relevant to describe the electronic states. ENuc =1 2 N X I=1 N X I6=J ZIZJe2 |RI−RJ|(1.8) Diagonalizing the electronic Hamiltonian of eq. (1.4) we can obtain the energy Es0({R})and the electronic states ψs0({r};{R}). The coefficients χss0({R})that couple the nuclear-electron system can be thought as the nuclear wavefunctions. They can be obtained inserting the expression (1.3) into eq. (1.1) and integrating over the electron variables {r}leading to the following equations hˆ TNuc +Es0({R})−Esi|χss0({R})i=−X i0 Csi0|χsi0({R})i(1.9) where ˆ TNuc is the nuclear kinetic energy and the matrix elements Csi0are given by Csi0= N X J=1 me MJhψs({r};{R})|∇J|ψi0({r};{R})i∇J N X J=1 me 2MJhψs({r};{R})|∇2 J|ψi0({r};{R})i (1.10) we can apply now the adiabatic approximation by ignoring the off-diagonal terms of Csi0matrix elements.2,3 Thus, the electrons are considered to remain in their state s0 as the nuclei move. In spite of the change in the energy Es0({R})and the wavefunction ψs0({r};{R}), the s0state does not change and the nuclear motion does not cause any electronic excitation. Therefore, in the adiabatic approximation, the nuclear dynamics are described by a purely electronic potential and the full set of states sis a product of nuclear and electronic states. If all the Csi0elements are ignored,4the equation for the nuclear motion simplifies to hˆ TNuc +Es0({R})i|χn0s0({R})i=En0s0|χn0s0({R})i(1.11) where n0here labels nuclear states. The nuclear potential energy, that we have seen that corresponds to the electronic energy Es0({R}), is usually termed as the Born- Oppenheimer energy surface. A sketch of a potential energy surface is shown in Fig. 1.2. We can see there several minima, they correspond to different structural phases in the 2In most books the Born-Oppenheimer and adiabatic approximations are referred as the same approach even though we have seen they are not strictly the same. 3Off-diagonal terms are responsible of the electron-phonon coupling. 4Diagonal terms are neglected for the calculation of phonon energies in the adiabatic approximation in the frozen-phonon and perturbation methods that we will see below. 7 Adolfo Otero Fumega set of possible nuclear configurations {R}. In principle, the ground state of the system will be the configuration {R0}that corresponds to the absolute minimum of the potential energy surface. Phase transitions between different structural phases can be driven by external fields such as pressure, electromagnetic fields or temperature. Figure 1.2: Sketch of a Born-Oppenheimer energy surface Es0({R}). Different minima have been drawn, they correspond to different structural phases in the set of possible nuclear configurations {R}. Phase transitions between these phases can be driven by external fields such as pressure, electromagnetic fields or temperature. The ground state configuration {R0}, i.e., where the nuclei are located in the ground state, can be computed using the force or Hellmann-Feynman theorem [10]. It allows to compute the nuclear forces as FI=−∂Es0 ∂RI =−hψs0|∂ˆ He ∂RI|ψs0i,(1.12) where RIis the position of nucleus Iand FIis the force on it. In the ground state configuration the forces on all atoms must be zero, since the {R0}configuration is a minimum of the potential energy surface. The force theorem can be generalized to any parameter λof the Hamiltonian as ∂Es0 ∂λ =hψλ|∂ˆ He ∂λ |ψλi,(1.13) and a finite energy difference between two different λstates, considered to be ground states for each λvalue, can be obtained as an integral over a continuous variation of the Hamiltonian from one state to the other. This is called the adiabatic connection [11]. Finally, note that the adiabatic approximation that we have considered in this section is well justified as long as the gap in the electronic excitation spectrum is larger than the energies involved in the nuclear dynamics. In metals or systems where the electronic states are degenerate the approximation might fail. 8 1 The Physical System In this section we have split the many-body problem into two: on one side electrons and on the other nuclei. We will see in the following two sections how to deal with each of them separately. 1.3 Interacting Electrons The Born-Oppenheimer and adiabatic approximations left us with an electronic problem to solve. The interacting electron Hamiltonian that describes this system is given by eq. (1.4). The total energy Eof the electronic system is the expectation value of that Hamiltonian ˆ He E=hΨ|ˆ He|Ψi hΨ|Ψi,(1.14) where Ψis the many-body electronic wavefunction Ψ({r};{R}). At this point, we stop using explicit notation for the dependence on the nuclear coordinates {R}since they are just parameters in the electronic problem. If we consider Ψas a variable trial function, the eigenstates of the Hamiltonian correspond to the stationary points of the energy E. As we have seen in the first section of this chapter, the eigenvalues and eigenstates of the Hamiltonian can be obtained by solving the time-independent Schrödinger’s equation (eq. (1.1)). The ground state wavefunction Ψ0will be the state with the lowest energy. We can obtain it by minimizing the total energy with respect to all variables {r}in Ψ. All the ground state properties will be determined by the ground state wavefunction Ψ0. Therefore, in order to characterize the electronic structure of a material, at least its intrinsic electronic properties, we need to obtain the ground state wavefunction.5 A widely known method to obtain the ground state wavefunction Ψ0is the Hartree- Fock method [12]. This approach assumes that the many-electron wavefunction Ψ({r}) can be expressed as an antisymmetrized determinantal wavefunction for a fixed number n of electrons. The requirement to be antisymmetric comes from the fact that electrons are fermions. In case that no spin-orbit interaction is included in the electronic Hamiltonian ˆ He, the wavefunction takes the form of a Slater determinant [13] Ψ({r}) = 1 (n!)1/2 φ1(r1)φ1(r2). . . φ1(rn) φ2(r1)φ2(r2). . . φ2(rn) . . .. . ..... . . φn(r1)φn(r2). . . φn(rn)  (1.15) 5The other stationary points in the energy with respect to Ψcorrespond to excited states. In this thesis we will focus on the ground state and small perturbations to it, such as the ones produced by small nuclei displacements. 9 Adolfo Otero Fumega ELDA xc [ρ] = Zd3rρ(r)εhom xc (ρ(r)).(1.31) The exchange energy of a homogeneous electron gas as a function of the density can be computed analytically. On the other hand, the correlation term has been computed with great accuracy with quantum Monte Carlo methods. The LDA exchange-correlation functional provides in general a good approximation for solids. It gives bond-lengths with an accuracy around the 1% error. Thus, it is capable of making highly accurate structural predictions. However, the energy gaps tend to be underestimated within this approach. •Semilocal Generalized Gradient Approximation (GGA): this approach is the next natural step to the LDA. The gradient of the density ∇ρ(r)is introduced in the exchange-correlation energy as EGGA xc [ρ] = Zd3rρ(r)εhom xc (ρ(r),|∇ρ(r)|).(1.32) There is some freedom in the way that the gradient is incorporated, so different GGA approaches can be defined. We will use the one given by Perdew, Burke and Enzerhof (PBE) [19]. This approach has been used in a huge number of systems for benchmarking most of the codes that we will use. It improves LDA for many structural predictions. However, as happens with LDA, band gaps are also underestimated, and sometimes a metallic behaviour (no energy gap) is predicted for correlated insulators. This behaviour has been termed the band gap problem of LDA and GGA. •Localized-orbital approximation (LDA+U): localized orbital approaches have been developed to overcome the band gap problem in systems where electrons tend to show a localized behaviour, like transition metal oxides and rare-earth elements [20]. These systems have partially filled dand fstates in which electron-electron interactions become strong, and so does the correlation energy. One of such localizedorbital approximations is the LDA+U, which is based on the LDA and includes the correlation term U given by a Hubbard model.8The method divides the system in two kinds of electrons: –Delocalized electrons that are treated with the usual LDA approach. 8The Hubbard model is a widely used model Hamiltonian in the field of highly-correlated electron systems. It leads to useful results for Mott insulators, in which correlations play a fundamental role to explain the insulating behaviour. It has also been used as the building block to treat more exotic systems like high-temperature superconductors or spin-liquids. 16 1 The Physical System –Localized electrons that include an on-site Coulomb interaction with the form 1 2UPi6=jninj, where niand njrepresent orbital filling. The orbital energies i of these localized electrons can be computed leading to i=LDA i+U1 2−ni(1.33) where LDA iis the LDA contribution to the energy and the last term take into account the Coulomb interactions. Thus, we have that for occupied orbitals ni= 1, and hence the energy of those orbitals is i=LDA i−1 2U and for unoccupied ni= 0 i=LDA i+1 2U. We can observe that the effect of the U is to produce an energy shift between the occupied and the empty states. Thus, leading to an energy gap due to electron-electron interactions in Mott-Hubbard systems. The U term can be calculated by means of many-body methods that are far beyond the scope of the present thesis. The important point for us is that LDA+U improves substantially (quantitative and qualitative) LDA and GGA approaches for systems in which electron localization happens. Since this thesis is devoted to the study of transition metal systems, we will find in the LDA+U approach a good description for many of such systems. It will upgrade band gaps, magnetic moments or exchange coupling constants. •Modified Becke-Johnson potential (MBJ): this method is able to provide Kohn- Sham gaps close to experimental values. The potential is built upon the Becke- Johnson potential (BJ) [21], whose goal is to reproduce the uniform gas limit and describe shell structures which are characteristic of exact exchange potentials. The kinetic energy density is introduced to treat correctly hydrogenic systems. The BJ potential was then modified by Tran and Blaha [22], and is the one that we will use in this thesis. The MBJ potential is a semilocal potential, what makes MBJ computationally cheaper than non-local hybrid functionals. However, it must be noted that it is a potential and not a functional, so it cannot be used to minimize the total energy. It is employed along another functional such as LDA or GGA used to solve the Kohn-Sham problem. The density and eigenfunctions are then used to build the MBJ potential and to compute the desired band gaps with high accuracy. 17 Adolfo Otero Fumega •Van der Waals (vdW): these non-local functionals are designed to treat dispersive forces that are originated from the long-range interactions of quantum fluctuating electric dipoles. The local and semilocal LDA and GGA are not able to describe any interaction of such characteristics. In our studies on transition metal compounds we will deal with layered compounds, whose layers are slightly bonded via van der Waals interactions. Therefore, we will need to use such a kind of non-local functionals to correctly predict their structural properties. An alternative is to include van der Waals interactions in a semiempirical way. This is done by adding a dispersion energy term Edisp to the total energy: Edisp =−1 2 N X I=1 N X I6=J C6IJ |RI−RJ|6fd(|RI−RJ|)(1.34) where the summations range over the number of atoms N,C6IJ is the dispersion coefficient for the atom pair IJ, and fd(|RI−RJ|)is a damping function that avoids divergences at near distances. It varies smoothly from 0 at close distances to 1 at large ones. Finally, some general remarks about exchange and correlation are: i) the exchange energy dominates over correlation, since the main effect of exchange is to remove the spurious self-interaction introduced by the Hartree energy, and ii) the inclusion of exchange and correlation always lowers the ground state energy of the system. They introduce what has been termed as the exchange-correlation hole, which is a decrease in the probability of finding electrons near each other due to the exclusion principle and the repulsive Coulomb interactions. 1.3.3 Spin and Relativistic Effects So far we have not directly talked about spin. As we know, electrons are Dirac fermions of spin 1/2. Spin arises naturally in the covariant Dirac formulation of quantum mechanics. The Dirac equation is linear in time and momenta in order to be Lorentz invariant. Solutions to this equation involve commutation requirements fulfilled by certain matrices9and a four-component wavefunction, called spinor. We see that this new formulation entails more complexity to the many-electron problem that we have been simplifying. As a result, the approach that is taken is to continue working with the Schrödinger’s equation that we have been using thus far, and include the relativistic effects that come from the Dirac equation on it. The electronic wavefunction is then expressed in terms of a two component wave function: Ψ(r) = Ψ↑(r) Ψ↓(r)(1.35) 9The way in which these matrices are chosen leads to describe Dirac or Majorana fermions. 18 1 The Physical System where Ψ↑(r)and Ψ↓(r)represent each spin wavefunction along a chosen axis.10 New terms due to the magnetic nature of spin and relativistic effects are included in Schrödinger’s equation:11 i~∂Ψ(r) ∂t =ˆ He=π2 2me +V+ˆ HZ+ˆ HSOΨ(r)(1.36) where the effects of a magnetic field Bintroduce the vector potential Ain the momentum operator pas π=p−(e/c)A. Also a Zeeman term ˆ HZ=µBσB, where µB=e~/2mecis the Bohr magneton. Finally, the spin-orbit interaction ˆ HSO, which can be viewed as the interaction of the electron with momentum pand spin σmoving in an electric field, or equivalently as a magnetic field in the rest frame of the electron. This term is the only one that mixes spin ↑and ↓components as ˆ HSO =ˆ HSO(↑↑)ˆ HSO(↓↑) ˆ HSO(↑↓)ˆ HSO(↓↓)(1.37) Equation (1.36) provides a way to introduce spin-orbit and magnetic effects in the non-interacting electron Hamiltonian. Therefore, they can be included in the Kohn- Sham DFT leading to a spin-DFT. The inclusion of a magnetic field in DFT is not trivial, it would require going to a density-current functional theory to treat correctly the first term in eq. (1.36). This term will be responsible for the effect of the magnetic field in the orbital motion.12 However, the Zeeman term can be easily included, since it affects only the spin, leading to an energy shift between both spin channels. Thus the effective potential that appears in the Kohn-Sham equation will have from now on a spin σdependency Vσ eff . Therefore, we will have two different spin channels for the electron density ρ(r) = ρ↑(r) + ρ↓(r). In the absence of an external Zeeman field, it may occur that the lowest-energy solution is spin polarized, and hence ρ↑(r)6=ρ↓(r), thus, breaking the spin degeneracy and leading to a magnetic system. In our calculations, the spin-orbit interaction will be introduced in a second-variational manner using the scalar relativistic approach. This approximation is derived when spinorbit is treated as a perturbative effect [23]. 1.3.4 Electronic Structure Methods In this part we will see different implementations to solve the Kohn-Sham equations that we have been describing above. All the methods that we will see agree when used 10In the calculations presented in this thesis we will consider that spin is quantized along a given axis. This simplifies enormously the treatment of spin effects. However, it must be noted that for the study of certain systems, like the ones that present spin textures, a noncollinear treatment of spin is required. 11Spin-independent scalar terms (Darwin and mass-velocity) are not included. 12In general, we will not deal with the inclusion of magnetic fields in this thesis. The only study in which we have included a magnetic field is the one in Chapter 10. The effect on the orbital motion is included there through a Peierls substitution in a tight-binding Hamiltonian. 19 Adolfo Otero Fumega with enough care and convergence [24]. However, each of them has its advantages and intricacies that must be considered when tackling a particular problem. The methods can be divided by the choice of basis set used to expand the Kohn-Sham orbitals and also by the number of electrons and potential considered. A natural basis set to expand the Kohn-Sham orbitals in crystals is the plane-wave basis set. Crystals are the kind of system that we will study in this thesis and, as we will see in the next Chapter, have a spatial periodicity that can be exploited by the use of plane waves. Plane waves are eigenfunctions of the Schrödinger equation with constant potential, so they are an excellent choice for systems in a nearly-free-electron regime. However, near the nuclei electrons tend to be localized, so their wavefunctions vary rapidly. They are usually called core electrons. The use of plane waves to describe these wavefunctions becomes inappropriate, since an enormous amount of plane waves would be required to expand such localized core wavefunctions. Thus, the methods that we present below correspond to the ones implemented in the codes that we will use. They basically show different ways to deal with core electrons. The APW+lo method The augmented plane wave plus local orbitals (APW+lo) is one of the most accurate and efficient methods to solve the Kohn-Sham equations [25, 26]. The method is based on the original augmented plane wave method APW by Slater in 1937 [27]. The idea is to divide the space in two regions as shown in Fig. 1.4. We will have the muffin-tin spheres centered on the atoms13 and an interstitial region between them. Different basis functions are used in each region: ϕ(r) = (Pncneiknrif r∈Interstitial region Plm Almul(r, El)Ylm if r∈Muffin-tin sphere (1.38) So we see that the interstitial region is expanded by plane waves, while within the muffin-tin spheres a combination of radial functions ul(r, El)times spherical harmonics Ylm is used, since they are solutions to the radial Schrödinger equations. The coefficients of the expansion are cnand Alm. These must fulfill certain requirements in terms of crystal periodicity and boundary matching of the muffin-tin and interstitial regions. In the original APW method the radial basis set was energy dependent ul(r, El), which induces undesired non-linearities in the Kohn-Sham problem. Different schemes were proposed to solve this issue like the linearized augmented plane wave method (LAPW) [28, 29]. In the APW+lo, the energy Elof the radial functions ul(r, El)is fixed in order to keep linear the problem. This reduces the variational flexibility of the radial basis set. In order to recover it, local orbitals ϕlo(r)are included within the muffin-tin region: ϕlo(r) = [almul(r, E1l) + blm ˙ul(r, E1l)] Ylm (1.39) 13The name muffin-tin spheres comes from their similar shape to kitchen molds for muffins. 20 1 The Physical System Figure 1.4: Schematic representation of the space division in terms of muffin-tin spheres around the atoms and the interstitial region between them. where alm and blm are coefficients of the expansion and E1lare the linearization energies where the solutions of the radial Schrödinger equation are evaluated. These local orbitals are required to be zero at the sphere boundary as well as normalized. The potential is expanded in the APW+lo method as V(r) = (PkVkeikr if r∈Interstitial region Plm Vlm(r)Ylm if r∈Muffin-tin sphere (1.40) and the charge density in an analogous way. So, no approximations are made to its shape, which is the reason why this method is a full-potential all-electron method. Thus, compared to other methods is by far the most efficient and accurate. The convergence of the APW+lo basis set is controlled by the parameter RmtKmax, where Rmt is the smallest atomic sphere radius and Kmax is the value of the largest kn vector in the expansion given by eq. (1.38). Values of RmtKmax ∼7provide a good convergence for the transition metals that we have studied in this thesis. The APW+lo method is implemented in the wien2k code [30] that we will use in most of the studies presented in the second part. Pseudopotentials and the PAW method As we have mentioned, our aim is to use plane waves as the basis set. However, the basis set becomes too large for localized electrons near the nuclei, i.e., for the core electrons. Thus, making plane waves impractical. The idea of pseudopotentials was first proposed by Hans Hellman in 1934 [31]. This approach is based on the fact that in a material, the electrons that are going to play a role in its properties are the valence electrons, while the core electrons of a specific atom remain practically unaffected by the other atoms. The bonding between atoms will be determined by the valence electrons, and hence they will be the important electrons to describe. 21 Adolfo Otero Fumega The pseudopotential approach proposes to build an effective potential and replace with it the atomic all-electron full-potential eliminating the core electrons. Valence electrons are then described by pseudowavefunctions ϕi,ps that oscillate much less, i.e., they have less nodes, and hence it is possible to describe them with a reasonable plane-wave basis set. A schematic of the pseudopotential idea is shown in Fig. 1.5. Note the smoothness of the pseudowavefunction which allows to describe it with fewer plane waves. Figure 1.5: Schematic of a pseudopotential and pseudowavefunction. A smaller planewave basis set is required to describe the pseudowavefunctions due to their less oscillating behaviour. The pseudowavefunctions and pseudopotential match the real ones outside the cut-off radius. Pseudopotentials are usually constructed from ab intio atomic calculations. An atomic reference state is considered by defining a cut-off radius Rcfor the core. Outside Rc, the pseudowavefunctions and the real wavefunctions for the valence electrons are required to have the same energies and amplitude (density). Soft pseudopotentials are those that have a large cut-off radius. This helps to reduce the size of the plane wave basis set needed to describe the pseudowavefunctions. However, these pseudopotentials are less transferable, i.e., accuracy is lost depending on the chemical environment. Thus, pseudopotentials deal with a compromise between transferability and basis-set size. The so-called small core approximation is taken by assuming that there is no overlap between core and valence states. Problems might arise when dealing with semicore electrons. The most common forms of pseudopotentials are the norm-conserving and ultrasoft pseudopotentials [32, 33]. In the case of norm-conserving pseudopotentials, the norm of the pseudowavefuntions and the real Kohn-Sham wavefuntion are required to be the same inside Rc. Ultrasoft pseudopotentials relax this constraint, which reduces the size of the 22 1 The Physical System required basis set. In return, a generalized eigenvalue problem is introduced. If we define the non-zero norm difference ñij inside Rcas: ñij =hϕi|ϕji−hϕi,ps|ϕj,psi,(1.41) so now the normalized state of the pseudo Hamiltonian ˆ Hps obeys the equation: ˆ Hpsϕi=iˆ Tupsϕi,ps,(1.42) with the operator ˆ Tups being a transformation ˆ Tups =1+X r<Rc,ij |piiñij hpj|,(1.43) where piare projectors that form a basis with the pseudo reference state inside Rc. We have seen that the idea of pseudopotentials entails neglecting the core electrons. However, based on the projectors introduced for ultrasoft pseudopotentials and the notion of augmentation explained for the APW method, core states can be somehow recovered within the Projector Augmented Wave (PAW) method [34]. The idea is to introduce a transformation ˆ TPAW similar to the one in eq. (1.43) for ultrasoft pseudopotentials. However, in this new transfomation, the full-electron Kohn-Sham wavefunctions ϕicalculated in the atomic reference system are still involved: ˆ TPAW =1+X i [|ϕii−|ϕi,psi]hpi|.(1.44) The PAW method is usually combined with the frozen core approximation, where core wavefunctions are considered to be unaffected by the environment. Therefore, they are not modified through the calculation. Pseudopotentials and the PAW method are implemented in the quantum espresso [35] and VASP [36] codes. We will make use of them in the studies shown in the second part. A clear advantage of the pseudopotentials and PAW method is that they are much faster than the APW+lo method and therefore materials with more atoms can be treated more easily. Forces are also implemented in an easier way. This is the reason why we will use many times pseudopotentials for structural relaxations and related properties. In contrast, the APW+lo method is the most accurate method for solving the Kohn-Sham equations. Analysis of the electronic structure and magnetic properties of transition metal compounds requires a good description of the d electrons. These electrons might be considered as semicore states when they appear localized at low enough energies. The pseudopotential method could be scarce in this case, meanwhile the APW+lo method ensures a fair treatment of the d electrons. 23 Adolfo Otero Fumega 1.3.5 Simplified Models In this final subsection we will see that it is possible to work with simple models that describe the essentials of a particular electronic structure. These models can be built from theoretical assumptions or can be distilled from ab initio calculations. We will discuss two simple models: the so-called tight-binding and Heisenberg models. They are probably the easiest models that we can find. However, they allow to gain insight into the physical problem and to compute certain properties that otherwise could not be calculated. Tight-Binding Model The first tight-binding (TB) model was introduced by Bloch in 1928 [37]. He considered a simple s-symmetry function as the basis for it. In 1934 Jones, Mott and Skinner 1934, extended the model to a basis of different atomic orbitals [38]. The idea of the TB model is to represent the electronic wavefunctions as localized orbitals ai(r−RI), each associated with an atom at position RI, where ilabels all the states in the basis. Hence the matrix elements of the Hamiltonian Hij for a pair of states iand jcentered at RIand RJrespectively, are: Hij =Zdrai(r−RI)ˆ Haj(r−RJ),(1.45) and the overlap between the states Sij is: Sij =Zdrai(r−RI)aj(r−RJ).(1.46) Thus, the problem is represented in a matrix form, where the basis set does not need to be explicitly specified. Only the matrix elements of the overlap and the Hamiltonian are required to solve the electronic structure problem. Nothing but a matrix diagonalization is needed to obtain the eigenvalues of the energy and the eigenfunctions as linear combinations of atomic orbitals. Atomic orbitals are usually taken as the localized basis, since the different bonding types between them are well known.14 In Fig. 1.6 σand πbondings between s and p orbitals are depicted. In 1954 Slater and Koster set the fundamentals of the TB method, where they included the celebrated Slater-Koster tables for the interatomic matrix elements [39]. Orbital and translational symmetries are considered to derive the matrix elements that appear on those tables. The real space approach provided by a TB Hamiltonian enables to compute things like surface states of non-trivial topological materials, or to include the effect of a magnetic field in the orbitals. Scalability is another advantage for this kind of models, thus, allowing to treat systems with a huge number of atoms. 14Other basis sets like the ones obtained from a Wannierization are commonly used. 24 1 The Physical System Figure 1.6: Schematic representation of the different σand πbondings between s and p atomic orbitals used to built a tight-binding model. The TB Hamiltonian can be formulated in the frame of second quantization leading to a summation of kinetic and potential energy terms: ˆ H=X i Eia† iai+X i6=j tij ha† iaj+a† jaii,(1.47) where Eiis the so-called onsite energy for a given orbital, i.e., the sum of the kinetic and potential energies of orbital i, and the tij is the so-called hopping term between the orbital iand j. It represents the kinetic energy an electron needs to hop from one orbital to the other. More complex terms could be included in the simple TB Hamiltonian given by eq. (1.47). Those new terms would lead to very different models like the Hubbard model that can be used to explain metal-insulator transitions, or many others to explain more exotic phenomena like superconductivity or topological insulators. We will work with TB Hamiltonians on the second part. The PythTB [40] and Pybinding [41] packages will be used to build TB Hamiltonians and to compute different properties. Heisenberg Model Finally, we present the Heisenberg model. This is the most representative example of a spin Hamiltonian. In systems that show a magnetic insulator behaviour it is possible to get rid of the charge degrees of freedom and build an effective spin Hamiltonian. The insulator behaviour could also be relaxed when the spin Hamiltonian is built for localized states, like the d electrons found in the transition metal compounds that we analyze in the second part. Thus, the system is simplified to a spin model, where spins lay localized at a given site iand interact between them. A schematic of the model is depicted in Fig. 1.7. In the case of the Heisenberg model, continuous rotations of spins are allowed. The Hamiltonian is given by: 25 Adolfo Otero Fumega linear response to a single perturbation is of the same order as the one for the unperturbed density. One particular advantage of DFPT to compute the harmonic phonon spectra compared to other methods is that the response to perturbations of different wavelengths is decoupled. In return it must be noted that it is a perturbative method, which entails that the approximation breaks for high-energy anharmonic corrections to the phonons or when dealing with unstable structures. DFTP is implemented in the quantum espresso and VASP codes. We will make used of it in the studies of the second part to compute harmonic properties. Frozen Phonon method A different approach to obtain the coefficients of the Taylor expansion of the potential energy surface is the so-called frozen-phonon method or finite-displacement approach [47]. The idea of this method is very simple, a set of structures is generated where different atomic displacements ∆uIα are performed. The interatomic force constants DIα,I0α0are then approximated as: DIα,I0α0≈ −FI0α0+∆uIα −FI0α0 ∆uIα (1.69) Thus, for the set of structures with the displacements we will have to compute the atomic forces FI0α0+∆uIα with respect to the undistorted structure. Usually FI0α0= 0 since the undistorted structure is the one where the atomic positions are in equilibrium. The interatomic force constants can be obtained from eq. (1.69); the Parlinski-Li- Kawazoe method [48], which is just a numerical fitting approach to obtain the force constants from a set of forces and displacements. Therefore, only the computation of the forces from a set of structures with different displacements is required. The size of the set can be cleverly limited by introducing the symmetries of the structure. An analogous procedure can be developed to obtain the third-order force constants. One advantage of this method is that it allows to include anharmonic effects that would be difficult to study through DFPT. Moreover, it can be implemented in any standard DFT code able to compute total energies and forces. However, the frozen-phonon approach is much more computationally demanding, in the case of crystals it might require unaffordable supercells to compute the phonon spectra at a desired wavelength vector. The frozen-phonon approach to compute the second and third-order force constants is implemented in the Phonopy [49] and Phono3py [50] packages. We will use them in the studies of Part II. 32 1 The Physical System Ab initio Molecular Dynamics and the Temperature Dependence Effective Potential method Finally, we will briefly present how to include the effect of finite temperature in a phonon calculation via the temperature effective potential (TDEP) method [51]. Until now, we have been working with the potential energy surface or Born-Oppenheimer energy surface. The energy minimum of such surface gives the equilibrium positions of the nuclei. This is appropriate when the effect of temperature can be neglected, i.e., if there are no temperature-dependent dynamical instabilities or strong anharmonic effects that lead to structural transitions. In such cases, where finite temperature or pressure is introduced, the potential energy surface is not suitable, and the Gibbs free energy becomes the proper thermodynamic quantity to minimize. We will be interested in the situation where a structural phase transition occurs at a transition temperature Ts, i.e., equilibrium positions {R0}will be different in each phase. Thus, we will have two different regimes. On one hand, the low temperature one that can be analyzed by the methods that we have been describing for the potential energy surface. On the other hand, the high temperature phase that requires a different treatment. Figure 1.8: Schematic representation of the free energy surface as a function of temperature. Equilibrium positions {R0}({R0 0}) are represented in the low (high) temperature phase T < Ts(T > Ts) as the minimum of the solid (dashed) line. The TDEP is based on the idea of generating a set of thermalized structures at a finite temperature and extracting the force constants for a model Hamiltonian that matches the configurations. The way to explore the phase space at a given temperature is by means of an ab initio molecular dynamics (aiMD) simulation in the canonical ensemble. At a certain number of time steps the system will be equilibrated at the chosen temperature. Then, the interatomic force constants at a given order can be extracted by minimizing 33 Adolfo Otero Fumega the difference in the forces obtained in the aiMD and those from the model Hamiltonian introduced by the TDEP method.18 We end this chapter in which we have described how to tackle the many-body problem of nuclei and electrons. The Born-Oppenheimer and adiabatic approximations allowed to split the analysis in nuclei and electrons. Density functional theory was described for the treatment of electrons and lattice dynamics for nuclei. Their practical implementation was presented along with methods to solve them. All the approximations and the methods that we have been describing so far are summarized in Fig. 1.1. 18In an aiMD atomic positions are moved through a period of time. A small time step is defined, a DFT calculation is performed at each of these steps to determine the forces and produce the atomic displacements that lead from one step to the next one. Thus, producing the dynamics. In the canonical ensemble the set of possible configurations is restricted by the given temperature. The way to introduce it in the calculations is by means of a heat bath that enters the Hamiltonian constraining the phase space. Different thermostats can be found. 34 2 Symmetry “Beauty is rather a light that plays over the symmetry of things than that symmetry itself." — Plotinus In this chapter we will talk about the concept of symmetry and how it is applied to the many-body problem. Symmetry is probably the most fundamental idea in Physics. According to Emmy Noether’s theorem, every conservation law is a consequence of the presence of a particular symmetry in the system [52]. For instance, translational symmetry of space entails linear momentum conservation, or angular momentum is related to rotational invariance. In a quantum system, the eigenvalues’ degeneracy or the eigenfunctions’ properties are determined by symmetry, and their transformation due to phase transitions can be understood as changes in the symmetry. More than that, many of the methods that we have already been describing in Chapter 1, are based on symmetry considerations, such us the basis sets needed to describe core and valence electrons. Therefore, such a fundamental concept deserves a closer look. The aim of this chapter will be to provide the basics to analyze a system of nuclei and electrons through the eyes of symmetry groups. We will see in the second part how symmetry guides all the analysis that we perform. 2.1 Definition of a Symmetry Group The way in which symmetry enters in the mathematical description of a system is through transformations. If a transformation Uis performed to the wavefunction Ψthat codifies the state of a quantum system as UΨ = Ψ,(2.1) and we obtain that the transformation keeps the system in the same state Ψ, we say that the system is invariant under such transformation U. A set of transformations {U}forms a group if all its elements fulfill the following four requirements: 35 Adolfo Otero Fumega •The product of any two elements of the group is itself an element of the group. •The associative law applies. •The unit (or identity) element Eis present in any group and its product with any element of the group leaves the element unaltered. •For every element of the group there is an inverse element such that the product of the element and its inverse leads to the identity element. The group will be a symmetry group if its elements leave the system invariant as given by eq. (2.1). Moreover, the group will be called Abelian if all the elements of the group commute. In general, this case does not happen. Non-Abelian groups are an active area of research due to their promising applications. We have stated in the previous chapter that the Hamiltonian of a system describes all the interactions between nuclei and electrons. The ground state of the system is obtained as the lowest energy eigenstate of the Hamiltonian. For this reason, the symmetry analysis is usually performed on the Hamiltonian. If a system is invariant under a transformation of a given symmetry group, the Hamiltonian must commute with the elements of that group. For instance, for the case of eigenstates with the same energy (degenerate), we can obtain one of the eigenfunctions and then obtain the others by applying to it the symmetry operations that leave the Hamiltonian invariant. Since this thesis is focused on the analysis of crystalline transition metal compounds, in the following sections we will analyze the kind of transformations that are encountered in those systems. Such transformations define the different symmetry groups that characterize the properties of the crystal. 2.2 Point Group Symmetry Rotations, inversions, reflections and their different combinations that leave a system invariant are called point symmetries. The operation of a point transformation Rpon a function f(r)(such as the density of a system) is: Rpf(r) = f(Rpr).(2.2) If a system has inversion symmetry, it means that it is invariant under the change of sign of all the coordinates, for instance f(x, y, z) = f(−x, −y, −z). In the case of mirror symmetry (or reflection), the system is invariant under sign inversion with respect to the mirror plane. For example, if the mirror plane is the xy plane, the system is invariant as f(x, y, z) = f(x, y, −z). Let us focus now on the case of rotations. For a free atom, all rotations about any axis are symmetry operations, i.e., they commute with the Hamiltonian. Thus, the free atom has full rotation symmetry. The basis functions for the full rotation group are the 36 2 Symmetry spherical harmonics Ylm. This is the reason why the atomic orbitals s,p,dand so on, are expressed as a product of a radial function Rnl and a linear combination of spherical harmonics. Thus, for the sorbitals we have s=Rn0Y00,(2.3) for the porbitals pz=Rn1Y10 px=Rn1 √2(Y11 +Y1−1) py=Rn1 i√2(Y11 −Y1−1), (2.4) and for the dorbitals dz2=Rn2Y20 dx2−y2=Rn2 √2(Y22 +Y2−2) dxy =Rn2 i√2(Y22 −Y2−2) dxz =Rn2 √2(Y21 −Y2−1) dxy =Rn2 i√2(Y21 +Y2−1), (2.5) where the subscript nhere denotes the principal quantum number. We see that these orbitals expand the electronic states of a free atom. More than that, core electrons hardly feel the presence of neighbouring atoms in a material, i.e., the interactions of core electrons with the electrons of neighbouring atoms are tiny. For this reason, atomic-like orbitals are a good basis to describe them. This is the idea behind the APW method explained in the previous chapter, where core electrons inside the muffin-tin sphere are expanded as a product of radial functions and spherical harmonics (eq. (1.38)). In the free atom (or nearly free atom, like the core states that we mention), the energy levels of a given nand lquantum numbers are degenerate. For instance, the five 3dorbitals are degenerate in energy.1This situation changes when breaking the full rotational symmetry by placing the free atom within a crystal. The environment created by the neighbouring atoms breaks the full rotational symmetry, thus, driving the infinite symmetry group to a finite one, in which only rotations of a certain angle leave the system invariant. 1At this point we are not considering spin degeneracy explicitly. However, note that for the 3dorbital we have 10 degenerate energy levels if we count spin. 37 Adolfo Otero Fumega Figure 2.1: Crystal-field splittings of the d-orbital energy levels in different environments. The central atom is depicted in blue and ligands in red. 38 2 Symmetry The finite rotation groups that appear by the placement of nuclei in a crystal break the degeneracy according to their symmetry operations. This effect is what is usually known as crystal-field splitting of the orbitals. Since our aim is to study crystalline transition metal compounds, in Fig. 2.1 we summarize the crystal-field splitting on dorbitals for different kinds of environments (ligand atoms surrounding a central atom). Octahedral, square planar, tetrahedral environments and different distortions to them are analyzed there. We show how the energy levels corresponding to the dorbitals of the central atom are split. Of course, in the case of sorbitals there is no crystal-field splitting, since the sorbital is non-degenerate. The porbitals, which are triply degenerate in the free atom, might be split by a crystal field. These splittings can be explained in terms of spatial anisotropies. If there is no anysotropy, the three porbitals remain degenerate. If there is anisotropy along one axis, the orbital that corresponds to that axis gets split from the other two. Finally, if the three axes are inequivalent, all porbitals become non-degenerate. We have been analyzing point group symmetry in terms of the electronic wavefunctions. However, point group symmetry affects or, better to say, constrains the energy levels of the phonon modes too. The polarization vectors of the modes can be analyzed in an analogous way to the electronic states by point group symmetry. Even more, anisotropies in macroscopic properties such as conductivity are also related to the point group symmetries of the nanostructure. 2.3 Translational Symmetry In this thesis we are concerned with the study of crystals, in which translational symmetry plays a fundamental role. Crystals are an ordered state of matter in which the equilibrium positions of the nuclei are repeated periodically in space. Therefore, all the properties of a crystal can be specified regarding a unit cell (the spatial periodicity) and the atoms that it contains. In Fig. 2.2 a crystal lattice formed by a two-atom unit cell is illustrated. A unit cell is characterized by the set of lattice vectors aαand by the position of the atoms inside the unit cell given by their vectors τI, usually known as the basis of the unit cell. Therefore, any translation Tnproduced by a linear combination of integer multiples nαof the lattice vectors Tn=n1a1+n2a2+n3a3(2.6) will leave the system invariant, i.e., any function of the system f(r)(such as the equilibrium positions of the nuclei or the electron density) is invariant under such translations as f(r) = f(r+Tn).(2.7) 39 Adolfo Otero Fumega This periodicity encountered in crystals allows the representation of those periodic functions in terms of Fourier components at wavevectors k, which are defined in the Fourier transform of the real space, know as the reciprocal space. Figure 2.2: Schematic of a crystal lattice with two atoms (grey and white) per unit cell. Six unit cells are depicted, the lattice vectors a1and a2, and the atomic basis vector τ2 for the white atom. If we consider a crystal formed by a large number of unit cells Nc=Nc1×Nc2×Nc3, where Ncα denotes the number of cells in each spatial direction α, we can simplify the problem by restricting to a discrete set of Fourier components periodic in that crystal volume. This implies that each Fourier component must fulfill a set of periodic boundary conditions that, in the limit of large crystal volume, restrict the wavevectors to: kaα= 2π×integer (2.8) where the set of wavevectors that fulfill this condition form the reciprocal lattice. Thus, we can define the set of reciprocal lattice vectors bβas: bβ·aα= 2πδβα.(2.9) The reciprocal lattice vectors bβdetermine a volume Ω1BZ in the reciprocal space known as the first Brillouin zone. The Hamiltonian that describes the many-body system is invariant under lattice translations. Considering the previous chapter, we can consider the Kohn-Sham Hamil- 40 2 Symmetry tonian for electrons and the harmonic Hamiltonian for phonons. They both commute with the translation operator ˆ Tn, thus eigenstates of the Hamiltonian can be chosen to be eigenstates of the translation operators and vice versa. Therefore, the eigenstates of the translation operator can be used to block diagonalize the Hamiltonian. This result is known as the Bloch theorem [37] ˆ Tnφ(r) = φ(r+Tn) = eikTnφ(r).(2.10) where kis a wavevector in reciprocal space: k=X β nβ Nβ bβ,(2.11) whose range can be restricted to the first Brillouin zone by taking nβ< Nβ. Therefore, the eigenstates of the Hamiltonian can be chosen with definite value of kwithin the first Brillouin zone: φk(r) = eikr uk(r).(2.12) These φk(r)are termed Bloch functions and we can see that they are built as the product of a plane wave eikr and a periodic function uk(r). The Bloch theorem is the reason why the methods that we have been describing for solving the Kohn-Sham and the lattice dynamics equations are based on plane waves. Therefore, the Hamiltonians that lead to these equations are block-diagonal in k, and consequently Kohn-Sham equations for independent electrons (eq. (1.30)) and the secular equation for the harmonic phonon spectrum (eq. (1.56)) become independent for each wavevector kin the first Brillouin zone. Then, we can compute the eigenvectors independently for each kand obtain a set of eigenvalues labeled by εi,k. In the macroscopic limit of the crystal the spacing between kpoints, given by eq. (2.12), tends to zero and hence we can take a continuum limit for k. This results in the concept of continuum energy bands εi(k). Figure 2.3 shows a schematic of electronic energy bands for a metal and for a band insulator. In general, we can determine the occupied levels as those that lay below the Fermi energy. The energy bands that lay well below the Fermi energy correspond to core electrons. We can notice that the energy dispersion in kof these bands is practically zero. This means that core levels are well localized that barely suffer the effect of neighbour atoms, being that the reason that justifies the use of pseudopotentials and atomic-like functions to describe them. Near the Fermi level, we found the valence bands below the Fermi level and the conduction bands above. If there is a gap between the highest energy valence band and the lowest energy conduction band the system is an insulator. A band insulator requires the complete filling of the band, with degenerate spin states. A different kind of insulating behaviour, e.g., the Mott insulator, will need a symmetry breaking driven by electronic interactions. A metallic behaviour arises for the gapless case. The energy dispersion of these bands is in general different. However, localized 41 Adolfo Otero Fumega order to grasp the idea of topology. Thus, different signatures to detect the topological character of materials are defined other than the Chern number. Crystal symmetries, like time reversal or inversion play a key role in the Berry Curvature of a system. In this chapter we have analyzed the fundamental concept of symmetry in the frame of condensed matter physics. A definition of symmetry group has been provided. Point groups, translations and time-reversal symmetry groups were described. Finally, a brief introduction to topological invariants was also given. 48 3 Statistics and Transport “Why do they called it rush hour when nothing moves?" — Robin Williams In this small chapter we will present some important quantities that help to reveal the physical behaviour of a system. We will provide a definition for the density of states per unit energy and the distribution functions for electrons and phonons in thermal equilibrium. Amongst other things, the analysis of the density of states allows to determine the stability of a system as well as to compute different equilibrium properties. Moreover, we will define transport coefficients through the Boltzmann transport equation. We will see how to compute them from first principles calculations. Thus, connecting the nanoscale description and the macroscopic properties of a material. 3.1 Thermal Equilibrium One of the most important quantities that we can obtain when solving the many-body problem is the density of states per unit energy DOS(E): DOS(E) = X i 1 Ω1BZ Z1BZ dkδ(εi,k−E)(3.1) where εi,kcan represent the energy of an electron or phonon. For independent particle approximations, like the Kohn-Sham for electrons or the harmonic one for phonons, eq. (3.1) corresponds to the number of states per unit energy. The dimensions of a system, its symmetries and particle interactions play a fundamental role determining the shape of the DOS a given configuration might present. The number of occupied states per unit volume at a given energy for a system in thermal equilibrium at a temperature Tis given by the product of the density of states and a probability distribution function f(E). Since phonons are bosons, a Bose-Einstein distribution function fBE(E)is applied: 49 Adolfo Otero Fumega fBE(E) = 1 eE−µ kBT−1 ,(3.2) where kBis the Boltzmann constant and µis the chemical potential. In a different way, electrons are fermions and hence they will be ruled by a Fermi-Dirac distribution function fFD(E): fFD(E) = 1 eE−µ kBT+ 1 ,(3.3) where the chemical potential µis called the Fermi level EFwhen T= 0. Figure 3.1: Schematic for different electronic DOS cases. The Fermi level is depicted as a dashed line. 50 3 Statistics and Transport Figure 3.1 shows a schematic of different DOS situations. A semiconductor or insulator is shown in the top panel. Note the presence of an energy gap above the Fermi level. In the middle panel a metallic case is exemplified. In the lowest panel an unstable configuration is depicted. When a high density of states occurs at the Fermi level for a given configuration, the system might be unstable. This high DOS at the Fermi level can be seen as a hint that prompts to a phase transition to a lower energy ground state of the system. A small perturbation will lead the system to a different ground state. In the framework of a mean field theory, it can be shown that the system gets stable by breaking a symmetry. An order parameter is associated to that transition. Note that a peak in the Fermi level could occur, but in such cases electron correlations or topological features would be required. As we have stated in Chapter 1, we are dealing with a many-body problem in which different degrees of freedom compete. Therefore, the path the system takes to get stabilized depends on energy considerations. For instance, the system might break timereversal symmetry and become a ferromagnet. In the Kohn-Sham description the spin degeneracy is broken and one of the spin channels gets more occupied than the other.1 This is the case of a Stoner ferromagnet [55]. Alternatively, the system might undergo a structural transition by breaking spatial symmetry. Jahn-Teller distortions [56] and charge density waves [57] are examples of structural transitions. Electronic interactions might also play a fundamental role driving the system to an insulating scenario. This is known as Mott insulator behaviour [58]. In all these phenomena, and in many others, a symmetry is spontaneously broken lowering the DOS at the Fermi level, and thus lowering the energy of the system. The emergence of gaps or pseudogaps2around the Fermi level is associated to that transitions. Apart from that, note that temperature is introduced through the Fermi-Dirac distribution function (eq. (3.3)), and it will play a major role in the stability of a given configuration. The DOS can be used to compute different physical quantities like the internal energy or the specific heat. Finally, let us define the so called projected density of states (PDOSj(E)). It provides the contribution of a particular atom or orbital jto the total DOS. In the case of phonons the projection is performed along the polarization vector of one of the modes and atoms. Meanwhile, in the the case of electrons the projection can be performed on a set of orthonormal states ϕjthat can represent an atomic orbital or a localized state of a given atom. PDOSj(E) = X i |hϕj|φii|2 Ω1BZ Z1BZ dkδ(εi,k−E),(3.4) where φican be seen as the Kohn-Sham states determined in the electronic problem. Thus, the total DOS is given by adding all the contributions of the projected PDOSj(E). Analysis of the DOS and PDOS will be crucial in the studies presented in the second part. 1The DOS of each spin channel is usually represented with different sign. Positive for majority spin and negative for the minority channel. 2Pseudogaps are energy gaps that occur at certain kpoints of the Brillouin zone. 51 Adolfo Otero Fumega 3.2 Transport Properties In the previous section we have introduced the concept of density of states and the distribution functions for phonons and electrons. As we have seen, these functions provide the number of occupied states for a system in thermal equilibrium. In this part we will analyze the situation when an electric field Eor thermal gradient ∇Tis present in the crystalline solid. In such cases, a steady-state flow of charge or heat is established as a consequence of the external fields and the internal scattering processes trying to restore equilibrium. This phenomenon is known as transport, and therefore a non-equilibrium statistical theory is required to study it. We will make use of the Boltzmann transport equation to fulfill this task. In the presence of external fields, the Boltzmann formulation considers that the system undergoes a small deviation from the equilibrium distribution. Hence, a linear expansion on time tcan be performed on the distribution function fthat describes the steady-state that is established. This leads to the following equation: df dt =∂f ∂t force +∂f ∂t diff +∂f ∂t scatt = 0,(3.5) where we can see that a steady flow is achieved by the competition between the deviations from the equilibrium distribution caused by the external fields (first term), the diffusion of the particles (second term) and the scattering processes trying to restore the equilibrium (third term). We will analyze separately electronic and phonon transport. Different approximations, and hence different levels of theory, will be used for each kind of particles. In order to underline that crystals are the kind of systems that we are going to study, from now on we will explicitly include the k-dependence in the development of those transport approximations. The relaxation time approximation is a widely used approach to solve the Boltzmann transport equation for electrons. It is based on the idea that the scattering processes can be related to a relaxation time τ(k): ∂f ∂t scatt =−f(k)−f0(k) τ(k),(3.6) where the numerator describes the deviation of the distribution f(k)from the equilibrium f0(k), and the relaxation time τ(k)provides how the system returns to equilibrium. Introducing this last equation on eq. (3.5), and linearizing it by considering that the external fields are weak and hence the deviation of the distribution from equilibrium is small, we arrive to the following equation for the perturbed distribution function f(k): 52 3 Statistics and Transport f(k) =f0(k)−τ(k)v(k)∂f0(k) ∂T ∇T+∂f0(k) ∂E(k)eE =f0(k)−τ(k)v(k)f2 0(k) kBTe E(k)−µ kBT(E(k)−µ) T+eE, (3.7) where E(k)is the electronic energy and v(k)is the electronic group velocity: v(k) = 2π ~ ∂E(k) ∂k.(3.8) Therefore, eq. (3.7) provides the electron population perturbed by a weak electric field and a small temperature gradient. Note that in an independent electron theory like the Kohn-Sham, the distribution is generalized for each electronic state with energy εi(k). Now that we know the perturbed population of electrons produced by the external fields, we can analyze the currents that emerge leading the system to the steady state. The electron current density Jeand the heat current density JQcan be defined as: Je=2e 8π3Zv(k)f(k)dk JQ=2 8π3Z(E(k)−µ)v(k)f(k)dk, (3.9) and introducing eq. (3.7) we arrive at the following relationships: Je=e2K0E−eK1 T∇T JQ=eK1E−K2 T∇T, (3.10) where the Knterms are integrals defined as: Kn=1 4π3~ZZ τ(k)(E(k)−µ)nv(k)v(k) |v(k)|−∂f0(k) ∂E(k)dSdE, (3.11) where the integral over volume in khas been converted into an integral over surfaces of constant energy. This is quite convenient since only states near the Fermi level (or chemical potential) contribute to the transport. Moreover, the relaxation time is usually approximated as a constant τ. This is known as the constant relaxation time approximation. Apart from that, in ab initio calculations, the so-called rigid-band approximation is commonly introduced. It assumes that the electronic band structure is independent to changes in temperature or chemical potential. Note that eq. (3.10) shows that either an electric field or thermal gradient can generate both an electron or heat current. Such relationship allows us to define the usual transport 53 Adolfo Otero Fumega coefficients (electrical conductivity, electronic thermal conductivity, Seebeck coefficient, etc.) that can be measured in experiments.3Thus, they provide a clear connection between macroscopic properties of materials and their nanoscale description. The isothermal electrical conductivity σprovides the relation Je=σEbetween the applied electric field and the current generated in the material. It is given by: σij =e2τ 4π~ZE=µ vi(k)vj(k) |v(k)|dSE.(3.12) For metals the electronic thermal conductivity κecan be well approximated by the Wiedemann-Franz law through the electrical conductivity as: κe=(πkB)2 3e2Tσ, (3.13) where the ratio κe/σ is known as the Lorenz number. For semiconductors and insulators there is no substantial charge flow, hence heat transport is governed by phonons. We will analyze it below. The relationship between the electric field generated in a material by a thermal gradient E=Ss∇Tis ruled by the Seebeck coefficient Ss. It is given by: Ss=(πkB)2eT 3σ∂K0 ∂E E=µ .(3.14) Note that the Seebeck coefficient, also known as the thermopower, is independent of the scattering rate in the constant relaxation time approximation. Thus, it will be one of the most accurate coefficients that we can compute in that approach. The thermopower will play a fundamental role for thermoelectric applications. In Part II we will make use of the BoltzTrap2 code [59] to solve the Boltzmann transport equations in the constant relaxation time approximation and obtain the electronic transport coefficients. In the study of phonon transport properties only thermal gradients will be considered as the external fields. Therefore, in this case, the linearized Boltzmann transport equations lead to a perturbed distribution function f(k) = f0(k)−τ(k)v(k)∂f0(k) ∂T ∇T, (3.15) where now the equilibrium distribution function f0(k)corresponds to the Bose-Einstein distribution given by eq. (3.2). In the harmonic approximation for phonons the perturbed distribution of eq. (3.15) is applied for each phonon mode. Differently from the electronic case where the relaxation time τ(k)was considered a given constant, for the case of phonons we can compute the relaxation time using the 3Note that these coefficients are second rank tensors that provide the relationship between the external fields and the current that they generate in the material. 54 3 Statistics and Transport anharmonic approximation described in Chapter 1. It can be computed as the inverse of the imaginary part of the self-energy (eq. (1.63)). The lattice thermal conductivity κlprovides the relationship JQ=−κl∇Tbetween a thermal gradient and a heat flow. It can be computed as: κij l=1 kBT2ZE(k)τ(k)f0(k) (f0(k) + 1) vi(k)vj(k)dk(3.16) where usually the phonon energy is substituted by the angular phonon frequency as E(k) = ~ω(k). Moreover, it is common to obtain the total κlas a summation of the contributions of each phonon mode. In the studies of Part II we will compute the lattice thermal conductivity using the ShengBTE code [60]. Finally, now that we have defined the different transport coefficients, we can see how they are coupled in a thermoelectric material. A thermoelectric material is characterized by showing a large thermoelectric effect [61], i.e., a thermoelectric material can generate an electric potential ∆Vfrom a temperature gradient ∇T(and vice versa) as: ∆V=−Ss∇T. (3.17) Therefore, this class of materials present potential applicability in the area of efficient energy management. As shown in eq. (3.17), a high Seebeck coefficient Ssis required in order to have a good thermoelectric. However, this is not the only quantity to pay attention when designing a good thermoelectric. When a temperature gradient is imposed to a thermoelectric material, its charge carriers tend to diffuse from the hot end to the cold end. A net charge is formed at the cold side (negative in the case of electrons and positive for holes). This produces an electrostatic potential. Thus, an equilibrium is reached between the chemical potential for diffusion and the electrostatic repulsion of the total net charge at the cold end. A large thermoelectric effect will be developed if the material is able to maintain the thermal gradient, i.e., if it presents a low thermal conductivity κ, and also if it shows a high carrier mobility, i.e., if it presents a high electrical conductivity σ. The coupling between the different transport coefficients that defines a good thermoelectric is introduced through the so-called thermoelectric figure of merit zT as zT =S2 sσ κT. (3.18) Therefore, the goal when designing a good thermoelectric is to maximize the figure of merit. Moreover, note that a good thermoelectric contains only a single type of carrier, it will be either a n-type or a p-type semiconductor. This is due to the fact that if we had electron and hole conduction at the same time, the net charge at the cold end would be zero and the induced electric field would not arise. 55 Adolfo Otero Fumega In this chapter we have introduced the fundamental concept of density of states and how it is used along distribution functions to compute some of the properties that describe a material in equilibrium. Moreover, we have presented the Boltzmann transport equation that allows to define and compute different transport coefficients. The ingredients to find a good thermoelectric material were also exposed. All this provides the connection between the nano and the macroscale of a material. 56 Part II Studies on Crystalline Transition Metal Compounds 57 Adolfo Otero Fumega S(T→∞)=−kB e ∂log g ∂N (4.10) where kBis the Boltzmann constant, eis the electron charge and Nis the number of electrons. In this paper the authors obtain the Seebeck coefficient in different situations. We will show some of the results in order to clarify the equation that we will use. For a system of Nspinless fermions with NAsites to be occupied, the number of possible configurations gis: g=NA! N!(NA−N)! (4.11) using equation (4.10), the so-called Heikes formula is obtained: S(T→∞)=−kB elog (1 −η) η(4.12) where η=N/NA. We want now to find an equation similar to (4.12), but for the case of fermions with spin. We will be interested in the case in which the repulsion energy between fermions (for our purpose these fermions will be electrons) is larger than the thermal energy, so in each site there can be just one particle. The number of possible configurations is the one seen for spinless fermions (4.11) multiplied by a 2Nfactor that takes into account the spin degeneracy. Hence, we get: g=NA! N!(NA−N)!2N(4.13) so the Seebeck coefficient that we obtain is: S(T→∞)=−kB elog 2(1 −η) η(4.14) Now that we have seen how to introduce the spin degeneracy, we are in the position to expand equation (4.13) to account for the t2gorbital degeneracy that occurs at the bottom of the conduction band in electron-doped STO. We will introduce heuristically a factor ft2gin eq. (4.13), that weighs the t2gdegeneracy in a similar way as the factor 2 does for the spin degeneracy: g=NA! N!(NA−N)!2NfN t2g(4.15) so inserting this in equation (4.10) we obtain: S(T→∞)=−kB elog 2ft2g (1 −η) η(4.16) but we have not talked about the shape of this factor ft2gyet. Our goal is to find an expression that reproduces the degeneracy behavior in different strain situations. We 64 4 Transport Properties of SrTiO3 will model it as a function of the number of electrons in the dxy orbital, Nxy and the total number of electrons in the t2gmanifold, Nt2g. Figure 4.1 shows the three limiting cases we will have to consider in order to model the degeneracy in every case. The unstrained case implies all t2gbands are degenerate, and hence there is a triple orbital degeneracy. If tensile (compressive) strain occurs, the dxy (xz/yz) band/(s) lies (lie) lower in energy and then a single (double) degeneracy occurs. Figure 4.1: We show the three limiting cases for the t2gdegeneracy. In the case without strain (middle panel), all the orbitals are degenerate (ft2g= 3). If strain is applied, degeneracy is broken towards having a lower-lying singlet (tensile strain: ft2g= 1, top panel) or a lower-lying doublet (compressive strain: ft2g= 2, bottom panel). So we encounter that the factor ft2gmust be a function of x=Nxy/Nt2g, its value must be contained in the interval [1,3], and the maximum being equal to 3at x= 1/3. Considering all these facts, we are going to fit the points shown in Fig. 4.1 to the following function (we chose a Gaussian because it is the simplest smooth function that one can think of that has only one adjustable extremum in the [1,3] interval): ft2g=γ+α βpπ 2 e−2((x−1/3) β)2(4.17) where the constants α,βand γare completely determined by the three limiting cases shown in Fig. 4.1. A plot of the factor as a function of xcan be seen in Fig. 4.2. In some previous works the weight of the degeneracy in the thermopower has been studied [101]. Considering eq. 4.16 and the range of values that the degeneracy factor ft2gcan have ([1,3]), the maximum variation in Sthat can be introduced due to the degeneracy ∆Smax (T→∞)is: ∆Smax (T→∞)=−kB e(log [3] −log [1]) ∼ −95 µV/K (4.18) However, in any realistic situation, the modification will not be that large. We have found an analytical function for the high temperature thermopower from a statistical calculation. A degeneracy factor for the t2gmanifold was introduced as a 65 Adolfo Otero Fumega Figure 4.2: ft2gfunction that results from eq. 4.17. This degeneracy factor fulfills the three limiting cases analyzed in Fig. 4.1 for the effect of strain in a t2gelectron system. It accounts for the orbital degeneracy of the system as a function of strain. Nxy and Nt2g need to be determined ab initio. function of Nxy/Nt2g. This procedure will make possible to calculate the entropy-related term S(T→∞)from our DFT calculations just by computing Nxy/Nt2gfor each strain introduced in the STO structure. What we have done so far is being able to decouple the influence in thermopower variations with strain that comes from electronic degeneracies and from band structure modifications (such as effective mass changes in a parabolic-band description), that strain will also introduce. 4.2.3 Results of the Model on SrTiO3 We have run calculations in SrTiO3simulating different strain situations by fixing the lattice parameter ato that of the different standard perovskite substrates mentioned in the computational procedures section. For each a, we have optimized the lattice parameter c. We have not considered oxygen octahedral tilts, since this kind of distortions are hardly dependent on the impurity introduced as dopant. Therefore, we have limited our analysis to a cell volume distortion caused by biaxial strain. We have found that tensile (compressive) strain increases (reduces) the unit cell volume (as shown in Fig. 4.3. This has been experimentally reported previously for other oxides [102] and also for this system [103]. In Fig. 4.3 we plot the evolution of the unit cell volume with strain. We observe that, tensile strain leads to an increase in unit cell volume. The volume reduction that occurs for compressive strain leads to a smaller Ti-Ti distance, which increases the hopping integrals and hence leads to an increased bandwidth. This bandwidth increase should be related to a decrease in effective mass. 66 4 Transport Properties of SrTiO3 Figure 4.3: Bandwidth and unit cell volume for different strain situations. We see that bandwidth decreases with volume and volume is increased as positive (tensile) biaxial strain is applied. With these structures we have performed electronic structure calculations including spin-orbit coupling with the TB-mBJ exchange-correlation potential, that allows to give an accurate band gap without additional computational cost. Figure 4.4 shows the band structures (only the conduction bands close to its bottom) for different strain situations together with the corresponding DOS’s on the same energy scale. GSO corresponds to the tensile strain limit and LAO to the compressive one. We can observe that the former leads, as explained above, to a reduction of the t2gbandwidth, together with some shifts in the DOS peaks. The bottom of the conduction band is at Γ. We can see also that the position of the unoccupied t2gbands is higher in energy for compressive strain and lower for the tensile strain case. Figure 4.5 shows the energy bandgap variation with cell volume and hence with strain. It can be seen that an increase in the lattice parameter of the substrate produces a decrease in the energy bandgap. We have calculated the number of electrons in the dxy orbital, Nxy and the total number of electrons in the t2gmanifold, Nt2g. Both values were obtained by integrating inside the Ti muffin-tin sphere. We will assume that the degree of localization will be identical for the three t2gorbitals and hence integrating inside the muffin-tin spheres is sufficient for our purposes since we are only interested in the ratio. This ratio is needed to obtain the degeneracy factor ft2gand consequently the thermopower at high temperature S(T→∞)as a function of doping. 67 Adolfo Otero Fumega Figure 4.4: Band structure representation and density of states (DOS) calculated for different strain values. It can be seen how the bandwidth increases for compressive strain and decreases for tensile strain. Figure 4.5: Energy bandgap variation with cell volume. We see that tensile strain decreases the energy gap. Figure 4.6 shows the degeneracy factor dependence with doping for the different strain situations analyzed. We can see there, as a function of electron doping (in terms 68 4 Transport Properties of SrTiO3 of electrons doped per unit cell, or Nb atoms substituting Ti, thus adding one extra electron to the conduction band per dopant). We can see that at low doping, the three limiting cases we discussed above are obtained (in the pure single-ion picture), with the unstrained case having threefold degeneracy and the compressive or tensile strain cases tending to a twofold degenerate state or a singlet, respectively. As doping increases, the t2gbands become more heavily populated and the actual band splittings become less and less important. Above 5% doping, the degeneracy factor becomes very close to 3, almost independent of strain. This illustrates that, even though in principle strain can tune degeneracies (by removing them), thermopower at high temperatures might not be affected if doping is substantial. This limit will of course depend on the system, but for STO we see that it is obtained (within a rigid band approximation) at values which are on the order of those sought for in the case of TE applications (∼1020 cm−3). The important point to notice here is that the effect of strain on STO will not damage the TE efficiency at high temperature by reducing orbital degeneracies, when working at those high doping levels we are discussing. Figure 4.6: Degeneracy factor as a function of doping. In the low doping regime the degeneracy factors collapse to the three limiting cases exposed in Fig. 4.1, so they become largely dependent on strain. As we increase the doping level the factor becomes strain independent and tends to 3as the t2gbands become more populated. Our calculations show that this occurs between 2-5 % Nb-doping. In the case of unstrained STO we see that the degeneracy factor is completely independent of the doping level since there is always a triple degeneracy in that case. This figure suggests that degeneracy effects vanish above ∼5% doping. Introducing the results for the degeneracy factor ft2gin eq. (4.16), we can obtain the thermopower at high temperature S(T→∞)as a function of doping for the five strain situ- 69 Adolfo Otero Fumega ations that we are analyzing. Once S(T→∞)is computed, we have to solve the Boltzmann transport equations to obtain the parameter A. After that, we are in the position to get the value of T0using eq. (4.4) and consequently determine the three parameters of our model (eq. (4.2)). We must warn the reader that the Boltzmann transport equations used to compute the parameter Awould be less accurate at lower temperatures (in particular the phonon drag term is not included). Consequently, we will restrict our conclusions to analyze the evolution with strain and not focus too much on the actual values predicted for the thermopower. Also, discrepancies with experimental values could be due to the use of the constant relaxation time approximation which is used routinely but could have its limitations [104]. Figure 4.7 shows six representations of the thermopower obtained using the model that we propose (eq. (4.2)) at different doping levels for the five strain cases analyzed, utilizing the values of the parameters that we have just calculated. We can see that the high temperature thermopower decreases in absolute value as we increase the doping level, as one would expect. The parameters A,S(T→∞),T0and also the obtained effective masses are summarized in Table 4.1. Figure 4.8 shows the dependence that effective masses have with strain for the high doping regime. We can observe there that the theoretical effective masses obtained increase their value as doping increases. However, eq. (4.9) proposed to obtain the effective masses from parameter A is just a mechanism to get the effective mass tendency with strain. As we have said, it will only be valid in the high doping regime so we will focus on analyzing the evolution of the effective mass with strain but not so much on its actual doping dependence, since this will be quite dependent on the type of dopant. This feature is analyzed in Ref. [105] concluding that the effective mass will also depend on the kind of defect introduced to dope STO and the kind of distortions these introduce in the lattice. Coming back to Fig. 4.7, we can analyze now each doping regime, we see that for low doping levels the strain tendency is given by the high temperature Seebeck S(T→∞), i.e. we see in Fig. 4.7, for 0.1% doping (∼1.7×1019 cm−3) that the Seebeck evolution with temperature follows the same strain dependence than S(T→∞). The unstrained case (with aequal to that of STO) is the one that has the largest Seebeck value (in absolute value) at any temperature. Meanwhile, at high doping levels, greater than or equal to 2%, according to the results shown in Fig. 4.6, S(T→∞)stops being strain dependent. Consequently, the only dependence on strain in the evolution of the Seebeck with temperature comes from the Aparameter. We can easily check that the curves in the 2%,4% and 50% doping levels (∼3.3×1020,∼6.6×1020 and ∼8.4×1021 cm−3respectively) follow the strain dependence given by parameter Ain Table 4.1, i.e., the inverse behaviour that effective masses have with strain (see Fig. 4.8). This comes about due to the fact that ft2gbecomes very close to 3 and strain independent at those doping levels. Hence, the only changes to the thermopower occur via variations of the energy-transport term, in which (for high doping regimes of 2% or more) the thermopower is enhanced (in absolute value) by tensile strain. 70 4 Transport Properties of SrTiO3 Figure 4.7: Representation of the Seebeck coefficient as a function of temperature using eq. (4.2) for six different doping regimes. The parameters obtained in each case are compiled in Table 4.1. It can be seen that in the low-doping regime (below 2.0%) the entropy-related term dominates the evolution with strain, leading to a reduction in the Seebeck coefficient (in magnitude) caused by any kind of strain. However, for high doping, it is the low temperature (energy-transport term) thermopower the one that dominates the strain tendency because S(T→∞)becomes strain independent (see Fig. 4.6). In that case, we see that tensile strain increases the thermopower (in magnitude). 71 Adolfo Otero Fumega Table 4.1: Model’s parameters for the six representative doping values shown in Fig. 4.7. Values for the effective masses m∗/m0obtained using eq. (4.9) are also shown. Doping 0.1%A(µV )S(T→∞)(µV/K)T0(K)m∗/m0 GSO 163000 -703 232 0.145 DSO 158000 -730 217 0.150 STO 159000 -737 215 0.151 LSAT 166000 -730 228 0.145 LAO 174000 -720 242 0.140 Doping 0.5%A(µV )S(T→∞)(µV/K)T0(K)m∗/m0 GSO 195000 -592 329 0.353 DSO 199000 -605 323 0.348 STO 190000 -609 312 0.368 LSAT 199000 -607 328 0.354 LAO 210000 -596 352 0.341 Doping 1.0%A(µV )S(T→∞)(µV/K)T0(K)m∗/m0 GSO 215000 -533 403 0.509 DSO 215000 -537 400 0.513 STO 219000 -545 403 0.506 LSAT 230000 -543 423 0.489 LAO 240000 -534 450 0.473 Doping 2.0%A(µV )S(T→∞)(µV/K)T0(K)m∗/m0 GSO 241000 -481 501 0.720 DSO 245000 -483 508 0.714 STO 252000 -486 519 0.700 LSAT 265000 -485 546 0.672 LAO 272000 -481 567 0.661 Doping 4.0%A(µV )S(T→∞)(µV/K)T0(K)m∗/m0 GSO 296000 -425 696 0.932 DSO 303000 -426 711 0.917 STO 311000 -427 729 0.899 LSAT 325000 -427 761 0.870 LAO 330000 -425 775 0.868 Doping 50.0%A(µV )S(T→∞)(µV/K)T0(K)m∗/m0 GSO 1490000 -154 9690 0.995 DSO 1560000 -154 10100 0.957 STO 1590000 -154 10300 0.949 LSAT 1670000 -154 10900 0.910 LAO 1750000 -154 11400 0.878 In order to analyze these trends in more detail, we have obtained the effective masses and its strain and doping dependence (see Table 4.1). Despite the fact that we show the effective masses calculated in six doping cases, it is hard to analyze its strain dependence away from the high doping limit. At 2% doping, we can observe that the effective mass is larger for the unit cells with a larger volume (volume increases for tensile strain). Increasing the volume leads to a reduction in the hopping parameter, smaller band widths and hence a larger effective mass [106]. Above 2% doping, such dependence with strain of the effective mass is retained. For lower doping values, the dependence with strain is more complex due to the partial involvement of different bands. The effective mass increases monotonously with doping and it reaches values close to 1.0 m0at about 50% doping. In 72 4 Transport Properties of SrTiO3 that doping region, a single parabolic band picture makes more sense than at low doping, where the sole concept of a single effective mass where all the bands contribute in a similar fashion is more difficult to identify. Figure 4.8: Effective mass as a function of strain in the high-doping regime. Effective mass increases (decreases) with tensile (compressive) strain and barely changes once all the bands become populated. The limit of 50% doping, even though beyond the validity of the rigid band approximation is shown as the limiting case found in LaAlO3/STO interfaces. In Fig. 4.9 we can see a comparison between the thermopower calculated using Bolz- TraP (energy-transport term) and the prediction using our model, in which the entropy term is included. This plot reinforces the idea that low doping regimes (0.1%) the degeneracy has a great effect on the thermopower evolution with strain (degeneracy is different for each strain value). However, the trend with strain obtained for high doping regimes (4%) retains that of the energy-transport term (BoltzTraP) and is independent of the entropy term (degeneracy is independent of strain). In Ref. [98] Janotti et al. analyze the directional effective masses dependency with strain. In order to compare with their data, we can obtain a mean effective mass m∗ ave performing a harmonic mean of the directional effective masses mΓ−X(parallel) and mΓ−X (perpendicular to the strain direction) given in this article and we can compare these m∗ ave with the effective masses we have obtained for 2% doping. In our calculations LSAT corresponds to biaxial strain −1.0%, STO to 0.0% and DSO to +1.0%. The results are summarized in Table 4.2. Their results, when analyzed this way, agree with ours in the strain dependence of the effective mass. However, the actual values are largely doping dependent, as we have discussed here. 73 Adolfo Otero Fumega Figure 4.13: Weighted phase space for the high temperature phase at 300 K as a function of frequency. An increase in the weighted phase space is related to an increase in the anharmonic scattering rate, thus decreasing the thermal conductivity. The arrows indicate the peaks associated to the polar modes that tune the thermal conductivity. (a) Comparison between cubic (black) and tetragonal (green) structures. (b) Evolution with the unit cell volume. of magnitude of the STO thermal conductivity is given by the dispersion of the acoustic bands, since these modes are the main heat carriers, but also tuned by the polar modes. The mean group velocity associated to the acoustic modes at 300 K that we have obtained is 5.24 km/s, which is in good agreement with previous experiments [II]. In Fig. 4.13 we have plotted the weighted phase space, as defined in [128]; it will help us to identify the main sources of scattering that tune the thermal conductivity of STO. The weighted phase space provides the scattering rate as a function of frequency weighted by the harmonic frequencies. Therefore, it is very sensitive to the phonon spectra and hence can be used to identify the modes responsible for the stronger scattering processes. 80 4 Transport Properties of SrTiO3 Figure 4.14: (a)Phonon lineshape for three different volumes in the Pm3mphase at 300 K. From top panel to bottom we present an increase in the volume of the unit cell, V0 being the experimental unit cell volume. At the Γpoint and between 100 −200 cm−1, we identify the triply degenerate polar mode. Increasing the volume decreases the frequency of the polar mode. The band dispersion of the acoustic modes increases when reducing the volume. (b) Phonon mean group velocity of the acoustic bands as a function of the unit cell volume. 81 Adolfo Otero Fumega It also helps to analyze how the scattering rate is modified when different distortions, that change the phonon spectra, are applied to the structure. In the case of cubic STO at 300 K, we can see in Fig. 4.13a (black dots) that the weighted phase space has a peak around 250 cm−1. At this frequency, this peak is associated to the scattering processes of the polar modes. So, we can conclude that the polar modes shown in Fig. 4.11 are the main source of scattering in STO. This triply degenerate polar mode was previously identified, using other methods, as one of the main sources of phonon-phonon scattering [111, 121]. Apart from that, in a recent work, this mode was also reported to be electron-phonon active [129]. The subsequent discussions will be based on how the acoustic and polar modes evolve under different distortions. In particular shifting the polar modes to lower frequency may be expected to scatter more strongly acoustic phonons lowering the thermal conductivity. Figure 4.14a shows the lineshape for three different unit cell volumes at 300 K for a cubic structure. We can observe that a decrease in the volume of the unit cell increases the dispersion of the acoustic bands. This is translated into an increase of the mean group velocity of the acoustic bands (Fig. 4.14b). A reduction of the unit cell volume also increases the frequency of the triply degenerate polar mode due to the reduction of the Ti-O bond length. This increase in the frequency (depicted in Fig. 4.11) reduces the phonon scattering, as can be seen in Fig. 4.13b, where a clear reduction of the phase space peak occurs by reducing the unit cell volume. Hence, both effects on the acoustic and polar modes produce an increase in the thermal conductivity. The calculated values are: 7.6, 13.4 and 18.0 W/mK for the largest unit cell volume to the smallest one, respectively. This increase of the thermal conductivity when the volume of the unit cell is reduced is in agreement with previous experiments on other family of oxides that report an increase of κwhen pressure is applied [130]. Figure 4.15 shows the lineshape computed at 300 K when a tetragonal distortion is applied. This distortion keeps the volume of the unit cell constant and equal to the experimental one, while the c/a ratio was kept fixed at 1.06. We can see that the triple degeneracy of the polar mode is broken as explained in Fig. 4.11. Since the clattice parameter becomes larger it will increase the bond length of the Ti with the apical O atoms of the octahedra, reducing the frequency of the mode along the z-axis as shown in Fig. 4.15. Moreover, due to the increase of the clattice parameter and the decrease of a, the acoustic modes change their dispersion as compared to the cubic case. The calculated lattice thermal conductivity is κxx =κyy = 10.4,κzz = 8.0W/mK. We see that the tetragonal distortion breaks the triple degeneracy of the κ-tensor and reduces the thermal conductivity at 300 K. The degeneracy break of the polar mode is the main responsible for decreasing the thermal conductivity. In Fig. 4.13a we can see how the peak associated to the polar modes is divided in two peaks (green dots) when a tetragonal distortion is applied. The difference that we observe between cubic and tetragonal thermal conductivity would be diminished taking into account the volume reduction that occurs in the cubic to tetragonal phase transition, which produces an increase on the thermal conductivity. Furthermore, we will see in the next subsection that including the octahedral rotations in 82 4 Transport Properties of SrTiO3 Figure 4.15: Phonon lineshape for a tetragonal distortion of the cubic phase at 300 K (c/a = 1.06). The triple degeneracy of the polar mode is broken due to this distortion as explained in Fig. 4.11. There is a softening of the polar mode along the z-axis. the tetragonal phase will increase the thermal conductivity. Thus, the difference obtained between both phases would be reduced by including these effects. The low temperature I4/mcm phase Below 105 K, STO undergoes a phase transition from the cubic Pm3mspace group to the tetragonal I4/mcm [75]. This tetragonal phase is characterized by the stretching of the unit cell and also by a TiO6octahedra antiphase rotation along the c axis [131]. The structure of this phase is depicted in Fig. 4.11. We have optimized the lattice parameters and the atomic positions. We have then applied different structural distortions and study how they affect the phonon dispersion and also the thermal conductivity. We will see that, as obtained for the high temperature cubic phase, the polar modes play a fundamental role in the heat transport. The optimized structure that we have obtained has a0= 3.868 Å, c0= 3.935 Å, lattice parameters and α0= 6.6◦octahedral rotation angle in agreement with previous calculations [124]. Both c0/a0ratio and α0are overestimated by GGA compared to experimental values [132]. Since the optimized structure is a stable structure at 0 K, we can perform the study using DFT calculations. Small distortions to this structure can be applied without introducing imaginary frequencies in the phonon spectrum that would invalidate our analysis. In this subsection we will study the effect on 83 Adolfo Otero Fumega Figure 4.16: Results for the variation of the unit cell volume in the I4/mcm phase. Octahedral rotation angle and c/a ratio are kept constant. (a) Phonon band structure for three different volumes, in red the largest and in blue the smallest, V0is the optimized unit cell volume. An increase in the volume reduces the frequency of the polar modes and the phonon dispersion of the acoustic modes increases. (b) Thermal conductivity as a function of temperature. An increase in the volume reduces the thermal conductivity. the spectrum and the thermal conductivity of unit cell volume, c/a ratio and octahedral rotation angle. We will keep two of these parameters constant while varying the third one to try to decouple their effects. Note that we will plot the thermal conductivity of this low temperature phase for temperatures above 105 K, the temperature of the structural phase transition. We do this to better visualize the variations in κwith each distortion. Moreover, below 80 K the thermal conductivity is not plotted since it is dominated by the sample dependent boundary scattering. Therefore, the realistic temperature range for the computed κis between 80 and 105 K. 84 4 Transport Properties of SrTiO3 Figure 4.17: Mean group velocity of the acoustic modes for the low temperature phase. Its increase produces an increase of the thermal conductivity. (a) Evolution with the unit cell volume, (b) c/a ratio and (c) octahedral rotation angle. Figure 4.16a shows the phonon dispersion for three different volumes of the unit cell. We can see that a reduction of the volume increases the dispersion of the acoustic bands, thus increasing the group velocity of the compound (Fig. 4.17a). We can also see that the polar modes increase their frequency when the volume is reduced. This causes a reduction of the phonon-phonon interaction. Figure 4.18a shows how the peak of the weighted phase space associated to the polar modes decreases when the volume is reduced. These two consequences in the phonon spectrum act in a positive way increasing the thermal conductivity when reducing the volume (see Fig. 4.16b). This trend is in agreement with the result obtained for the high temperature cubic phase. The thermal conductivity at 100 K for the optimized structure is κzz = 16.1W/mK and κxx =κyy = 14.2W/mK. Considering that V0is overestimated by GGA, these are in pretty good agreement with the experimental value [II] from the List of Publications and also with previous calculations 85 Adolfo Otero Fumega Figure 4.18: Weighted phase space for the low temperature phase at 80 K as a function of frequency. The high peak around 200 cm−1is associated to the polar modes. Its increase produces an increase of the anharmonic scattering rate and hence a reduction of the thermal conductivity (a) evolution with the unit cell volume, (b) c/a ratio and (c) octahedral rotation angle. [111]. We see that a decrease in 1% in the unit cell volume produces an increase of around 30% on the thermal conductivity. Figure 4.19a shows the phonon dispersion for three different c/a ratios. We can see that a reduction of the c/a ratio increases the dispersion of the acoustic bands, thus increasing the group velocity of the compound (Fig. 4.17b). We can also see that the 86 4 Transport Properties of SrTiO3 Figure 4.19: Results for variation of the c/a ratio in the I4/mcm phase. Unit cell volume and octahedral rotation angle are kept constant. (a) Phonon band structure for three different c/a ratios, in red the largest and in blue the smallest, c0/a0is the optimized c/a ratio. The acoustic bands undergo a softening, mostly at the X point, when the c/a ratio is increased. At the Γpoint, we see that an increase of the c/a ratio increases the gap between the xy and z-axis polar modes. (b) Thermal conductivity as a function of temperature. Increasing the c/a ratio lowers the thermal conductivity. z-axis polar mode softens when the c/a ratio is increased. Moreover, the energy gap between this polar mode and the other doubly degenerate polar mode increases with c/a (see Fig. 4.11). Figure 4.19b shows the increase of the thermal conductivity when c/a is reduced. We observe that a decrease of 1% in the c/a ratio produces an increase of around 14% on the thermal conductivity. The effect of the c/a ratio on the thermal conductivity is relatively smaller compared to the effect of the unit cell volume. The reason for this is that one polar mode softens and the doubly degenerate ones raise their frequency in the 87 Adolfo Otero Fumega Figure 4.20: Results for the variation of the octahedral rotation angle in the I4/mcm phase. Unit cell volume and c/a ratio are kept constant. (a) Phonon band structure for three different angles, in red the largest and in blue the smallest, α0is the optimized rotation angle. An increase in the angle increases the frequency of the polar modes. (b) Thermal conductivity as a function of temperature. Increasing the angle raises the thermal conductivity. case of the c/a distortion, while the three polar modes soften when modifying the volume (see Fig. 4.11), increasing the phonon-phonon scattering and having a larger effect on κ. Figure 4.20a shows the phonon dispersion for three different octahedral rotation angles α. We can see that increasing αhas not a big effect on the dispersion of the acoustic bands, compared to the aforementioned volume distortions or varying the c/a ratio. Mean group velocities of the acoustic modes increase with the angle, but barely compared to the other distortions. We can also see that the z-axis polar mode softens when αis reduced (see Fig. 4.11). Figure 4.20b shows the increase of the thermal conductivity when αis 88 4 Transport Properties of SrTiO3 increased. This can be understood in terms of the evolution of the weighted phase space with the octahedral angle. Its increase causes a reduction of the peak. We observe that an increase of 1◦in αproduces an increase of around 14% on κzz while only 3% on κxx =κyy. The effect of varying αon the thermal conductivity is lower compared to the effect of the two other distortions we have considered and analyzed. The main reason for this is that the acoustic bands are barely modified by α, since there is no change in the lattice parameters. The competition between the antiferrodistortive (octahedral rotations) and ferroelectric phases was already studied in literature, by suppressing the rotations one can in principle achieve a ferroelectric phase [133–135]. This could be achieved via strain engineering since different octahedral rotation patterns and angles can be obtained [136]. We have shown here the coupling between the octahedral rotations and the polar mode responsible for ferrolectricity and also the effect that this has on the thermal conductivity of STO. 4.4 Conclusions In this chapter we have described two independent works on transport properties of SrTiO3. They are presented together due to the applicability of STO as a thermoelectric material. In the first work we have analyzed the thermopower of STO and its strain dependence for different n-doping levels, being able to decouple the different origins of these dependencies: those coming from electronic degeneracies and entropy-related terms from those coming from the actual modifications of the band structure. We have proposed a simplified model for the Seebeck coefficient and applied it to analyze the thermopower strain dependency of STO. We have obtained that in the low-doping regime the electron degeneracy dominates the Seebeck so any kind of strain decreases the thermopower, while at high doping the cell volume effect dominates, i.e., increasing volume (tensile strain) leads to an increase in the Seebeck. In the second work we have analyzed the thermal conductivity of STO. We have studied the main sources of thermal transport and how these are modified when a distortion is introduced in the structure. We have analyzed both the high temperature cubic phase and the low temperature tetragonal one. We have identified the acoustic bands and the low-energy triply degenerate polar mode as the main heat carriers. Reduction of the unit cell volume or the c/a ratio (i.e., approaching a cubic structure) are found to increase the lattice thermal conductivity. The octahedral rotation angle was found to be strongly coupled to the phonon-active polar modes. We report that an increase of this angle raises the thermal conductivity. Considering together both works, we can figure out some conclusions about the thermoelectric response of this material. We have seen that in order to optimize this response we need to increase the thermopower and reduce the lattice thermal conductivity as much 89 Adolfo Otero Fumega charge, supposing they become fully ionized, as is found in other oxides [136]. Finally, E0_vac is the energy of the configuration without vacancies 0_vac, this term represents the zero of energies when obtaining Efrom eq. (5.1). We have assumed that all the vacancies form O2and we have also neglected the electrostatic interaction between the homogeneous background charge and the charged vacancies. This assumption was made considering previous works [165–167] and the fact that the most stable compound for the oxygen element is the O2molecule. All terms in eq. (5.1) were calculated using the GGA scheme, but Eewas computed using the LDA+Uscheme for various Uvalues, since GGA leads to metallic solutions and the LDA+Umethod allows to compare the band alignment between two realistic insulating solutions, both for the stoichiometric and non-stoichiometric systems. The results for the energy of each oxygen-vacancy configuration using eq. (5.1) are shown in Fig. 5.3. We observe that only the configuration 2_vac_D is more stable than the one without vacancies. This corresponds to chains of oxygen vacancies perpendicular to the (001) and contained in the plane of vacancies (see Fig. 5.2). Figure 5.3: Energetics for the different oxygen-vacancy configurations computed using eq. (5.1). We can see that the most stable configuration is the 2_vac_D. This configuration corresponds to vacancy chains perpendicular to the (001) direction and contained in the vacancy plane. Only results with U= 7 eV are presented in this figure. Before analyzing in detail configuration 2_vac_D, which is the ground state, we can explore the other configurations. For the case without vacancies, the system is a nonmagnetic insulator. We also see that this configuration is quite stable compared to the others that present oxygen vacancies (2_vac_D is the only configuration with vacancies 96 5 Ferromagnetic and Insulating Behavior of LaCoO3 that is more stable than the stoichiometric solution). If we put an eye now on the two configurations with one vacancy, we will realize that one of them is much more stable than the other. The difference in energy between them is a sizable 0.37 eV/Co. This result suggests us that in the case of adding a second vacancy, it could be energetically favorable to include it in an equivalent position to configuration 1_vac_B, i.e, forming chains perpendicular to the (001) direction. We have demonstrated this statement calculating the energy for the 7 possible 2-vacancy configurations. As we have said 2_vac_D gives a stable solution. All the other n= 2 cases are also FMI solutions, but higher in energy, except 2_vac_C, which is a ferromagnetic metal. We have tried several values of Uin a wide range, but an energy gap was not opened for that particular configuration. Table 5.1: Energies (in eV ) for the different configurations at different values of U. The zero of energies is represented by the configuration without vacancies. It can be observed that configuration 2_vac_D is a robust ground state in a wide range of values of U. It is the only configuration that is stable for U≥5.5eV . U(eV ) 4.1 5.5 6.8 8.2 1_vac_A 5.78 5.40 5.08 5.16 1_vac_B 1.33 0.96 0.66 0.70 2_vac_A 2.29 1.44 0.72 0.87 2_vac_B 2.61 1.60 1.23 1.43 2_vac_C 1.00 -0.16 0.08 0.47 2_vac_D 1.93 -0.24 -0.64 -0.41 2_vac_E 3.13 2.41 2.45 2.90 2_vac_F 1.38 1.10 0.99 1.56 2_vac_G 2.61 1.62 1.33 1.75 Table 5.1 shows a summary of all the energetics of the configurations considered as a function of Ufor four selected values of U. We can see that the only stable solution that includes O-vacancies is the one we have analyzed thoroughly in the paper: configuration 2_vac_D. We see that it is the ground state for a wide range of Uvalues, U≥5.5eV . The typical Uvalues utilized for this kind of cobaltates in similar calculations are in the range of 6−8eV [159, 168, 169]. Thus, our calculations show a robust trend indicating that configuration 2_vac_D is indeed the ground state of the system. Figure 5.4 shows that there is a range of Uvalues where 2_vac_D becomes stable, and such Uis within that range. We have also seen that configuration 2_vac_D is a robust ground state in a wide range of Uvalues above 5.5eV . We will determine the origin of this stabilization in the next subsection as a change in the electronic structure produced by the electronic correlations. If we analyze in more detail the structure of the 2_vac_D configuration that we have obtained from our calculations, we find the following characteristics that can be 97 Adolfo Otero Fumega Figure 5.4: Energy of the 2_vac_D configuration as a function of Ucomputed using eq. (5.1). The zero of energies would be the stoichiometric solution without vacancies. We observe that for values of U > 5.5eV the configuration labeled 2_vac_D becomes stable. compared with the experimental results reported in Ref. [156]. The lattice parameter along the (001) direction undergoes a contraction compared to the bulk value of LCO. We have obtained c= 3.810 Å, which is in reasonable agreement with the experimental value c= 3.767 Å[156]. The distance between LaO planes parallel to the plane of oxygen vacancies depends on whether one considers the distance between two LaO planes that contain the plane of vacancies or two stoichiometric planes. In the former case, with the plane of oxygen vacancies inside, we found that the distance is 4.02 Å. In the other case we found it is 3.85 Å. Experimentally, the same trend is observed but with a larger ratio [160]. Therefore, we have found in this subsection a configuration of oxygen vacancies that agrees with the experimental microscopy data obtained for LCO films grown on STO. 5.3.2 Electronic Structure and Magnetic behavior We analyze now the electronic structure of configuration 2_vac_D. As we have already said, this presents a ferromagnetic-insulator behavior and is found to be a ground state. We will try first to explain the simple ionic model that accounts for the magnetic moments observed and its consistency with the ab initio calculations of the partial density of states (DOS) of the Co atoms shown in Fig. 5.5. In it, we can see the 3 types of Co atoms that appear in the converged solution: a non-magnetic Co atom six-fold coordinated and two magnetic Co atoms with 5-fold oxygen coordination, one with a higher 98 5 Ferromagnetic and Insulating Behavior of LaCoO3 value of the magnetic moment than the other. Figure 5.5: Partial DOS for Co2+ HS, Co2+ LS and Co3+ LS computed in the LDA+U scheme, for U= 7 eV . The Fermi level is set to the zero of energies. Upper panel, Co2+ HS: the t2gand the egmajority-spin states are occupied, and also the dxz, dyz minority orbitals. Middle panel, Co2+ LS: the t2gmajority and minority orbitals are occupied and so is the dz2majority. Lower panel, Co3+ LS: the t2gmajority and minority orbitals are occupied. From an ionic picture, if we consider the unit cell shown in Fig. 5.1, we have La3+, O2−and Co3+ for the stoichiometric compound without vacancies. The t2gorbitals of Co3+, which are in an octahedral environment, are occupied in the low-spin (LS) state, see Fig. 5.6 panel c), which gives rise to a zero magnetic moment and an insulating state caused by crystal-field splitting between the Co t2gand eglevels that leads to a gap opening around the Fermi level in the simple diamagnetic configuration. If we remove now two oxygens from the unit cell (Fig. 5.1), we have the stoichiometry La12Co12O34. Structurally, there are four Co atoms that are five-fold coordinated and the remaining eight Co cations are six-fold octahedrally coordinated and remain nonmagnetic, suggesting they retain the original 3+ valence. A simple electron count would naively imply that the other 4 Co atoms are Co2+ cations. These are in a five-oxygen environment, that could be analyzed in an ionic picture as an octahedral environment distorted along the z-axis. Such distortion breaks the degeneracy of the t2gand the eg levels as sketched in Fig. 5.6. The coordinate system that we use to analyze the electronic structure is local to each Co atom with z lying along the O-Co-vacancy direction. Our 99 Adolfo Otero Fumega Figure 5.6: Ionic model for Co2+ HS, Co2+ LS and Co3+ LS. a) Co2+ HS: the t2gand the egmajority states are occupied, and also the dxz, dyz minority orbitals. b) Co2+ LS: the t2gmajority and minority orbitals are occupied and also the dz2majority. c) Co3+ LS: the t2gmajority and minority orbitals are occupied. calculations show that the magnetic moment values are consistent with an ionic picture where half of these Co2+ cations are in the LS state (in the 2_vac_D panel of Fig. 5.2 the Co atoms on the left side), and the other half in the HS state (the ones on the right side). The actual values for the magnetic moments inside the muffin-tin spheres obtained from our calculations are approximately 0.8 and 2.4 µB, which account for the ionic values plus substantial hybridizations. Taking into account Fig. 5.6: for the case of Co2+ LS (Fig. 5.6 b)), we observe that the t2gbands are fully occupied and there is one extra electron occupying the lower-lying egstate, which is of dz2parentage due to the missing apical oxygen. Thus, the LS Co2+ cations are in an S= 1/2state. In the case of Co2+ HS (Fig. 5.6 a)), the majority channel is fully occupied, i.e., the t2gand the egmajority states are 100 5 Ferromagnetic and Insulating Behavior of LaCoO3 occupied, there are two minority t2gelectrons occupying the lower-lying doublet xz/yz, which is split from the higher-lying xy orbital due to the missing apical oxygen in the fivefold-coordinated environment. This leads to S= 3/2HS state. This simple ionic picture can be observed to be approximately reproduced (with substantial hybridizations, as is common in transition metal oxides) in the partial DOS for each inequivalent Co cation shown in Fig. 5.5. In this figure we show that the DFT calculations can be understood with the simplistic ionic model sketched in Fig. 5.6. The labels in Fig. 5.5 help to see the correspondence with the ionic model. Spectroscopy data confirm that most of the sample is Co3+ non-magnetic and that some magnetic Co2+ exist [156]. The data do not discard the existence of Co3+ HS yet they cannot confirm it. Our calculations show that Co3+ HS does not appear. No analysis of the possible existence of Co2+ LS that our calculations predict is provided in Ref. [156], thus our calculations are in principle consistent with spectroscopic evidences. Figure 5.7: Partial DOS Co2+ LS for two different values of U,4and 7eV respectively. Upper panel: we observe that the dz2orbital (highlighted in yellow) is not fully occupied which results in a solution that is not the ground state. Lower panel: we observe in this case that the dz2orbital is fully occupied, which stabilizes the 2_vac_D configuration to become the ground state. The total spin magnetic moment in this solution is 0.67 µB/Co. If we now include spin-orbit coupling in the calculations, the total ordered moment obtained as the sum of lz+2szis about 0.8µB/Co, since the magnetic Co cations present a non-negligible orbital 101 Adolfo Otero Fumega angular momentum that becomes partially unquenched. This is in reasonable agreement with the saturation magnetization of the LCO thin films obtained experimentally of about 0.85 µB/Co [156]. We also observe in Fig. 5.5 that the energy gap is between the dz2 orbital of the Co2+ LS (top of the valence band) and the egof the Co3+ LS (bottom of the conduction band). Of course the particular value of the gap will depend on the value of Uchosen for the LDA+Ucalculations. In particular, for the DOS presented in Fig. 5.5, the value used was 7eV . In the previous subsection we stated that for values of Uless than 5.5eV the 2_vac_D configuration was not stable. This can be analyzed through the change in the electronic structure that occurs when comparing the low-Uand high-Usolutions. Plots of the partial DOS of Co2+ LS for two different values of U, one for an unstable case and another one for the stable case are shown in Fig. 5.7. We can see that the dz2orbital of the Co2+ LS is not fully occupied for the unstable case, while it becomes totally occupied for the stable case at larger Uvalues. This change in the electronic structure is correlated with the stabilization of the 2_vac_D configuration. Therefore, according to our results, the insulating properties of these films occur naturally due to the additional band splittings introduced by the non-octahedral environments and the correlation energy of the d electrons. Table 5.2: Energy of the whole unit cell (12 Co atoms) for the different magnetic configurations. They are all referred to EF M . Energy (eV) EFM 0.0 EAFM12.1 EAFM20.5 EAFM32.2 EAFM41.9 EAFM52.2 Finally, let us analyze the magnetic order that occurs in LCO. In order to do so, we have performed calculations on different magnetic configurations in the ground state structural configuration 2_vac_D. These configurations are shown in Fig. 5.8. The corresponding magnetic exchange coupling constants we have introduced into a Heisenberg model of the system are depicted as the various J’s. The energies of each magnetic configuration are shown in Table 5.2. They are all referred to the energy of the ferromagnetic (FM) configuration EF M , which is taken as zero. Total energy is given for the whole unit cell (12 Co atoms in total). The arrangement of the Co2+ HS and LS is preserved for all the magnetic configurations and the geometry of the structure is not changed. We see that the FM configuration is the most stable. 102 5 Ferromagnetic and Insulating Behavior of LaCoO3 Figure 5.8: Ferromagnetic and antiferromagnetic configurations of the 2_vac_D in the plane of O-vacancies. The spin direction of each Co is depicted with arrows. The exchange interactions are labeled with J1, J2, J3and J4. The numbers inside the red circles represent the modulus of the magnetization for each Co atom (1 and 3 µB) in the corresponding ionic limit. Moreover, we can obtain the value of the different exchange interactions as a function of the energies of the magnetic configurations: J1=EAFM2+EAFM3−EAFM1−EF M 24 = 300K J2=EAFM1+EAFM3+ 2(EAFM4−EAF M5)−EAFM2−EFM 8= 4600K J3=EAFM1+EAFM3−2(EAFM4−EAF M5)−EAFM2−EFM 72 = 730K J4=EAFM1+EAFM2−EAFM3−EF M 24 = 150K (5.2) We observe from eqs. (8.1) that all the the exchange interactions (J’s) are positive, which entails a ferromagnetic order according to our sign convention. 103 Adolfo Otero Fumega We have characterized in this section the electronic structure of the 2_vac_D configuration. We have related the calculated DOS to a simple ionic model that explains both the ferromagnetic ordering with a consistent value of the total magnetic moment and the insulating behavior. 5.4 Conclusions In this chapter we have analyzed the ferromagnetic-insulating behavior of LCO when it is grown on (001) STO. We have used DFT calculations to analyze the electronic structure properties of LCO under tensile strain and for various stoichiometries including different oxygen vacancy configurations and concentrations. We have found that the ground state of LCO when it is grown on STO is given by an off-stoichiometry of the form LaCoO2.83 produced by about 6% oxygen vacancies. We have shown that the vacancies form chains perpendicular to the (001) direction and are contained in the plane perpendicular to the film/substrate interface, consistent with experimental findings. Other structural features like the distance between La layers or the lattice parameter in the (001) direction agree with the experimental measurements. We found that the Co atoms that lie in the plane of vacancies are a mixture of Co2+ LS and Co2+ HS. The total magnetic moment is in close agreement with experimental measurements. The ferromagnetic insulating behavior of the film is readily obtained as a ground state, the gap opening occurring naturally in that electron count due to the appearing crystal field splittings together with the addition of a reasonable Uvalue. Our ab initio total energy calculations confirm that all the nearest-neighbor magnetic exchange couplings are indeed ferromagnetic. In conclusion, we can say that the strain introduced by STO on LCO favors the presence of O-vacancy chains perpendicular to the (001) direction. Ferromagnetism arises from the inclusion of those vacancy chains and can be explained with an ionic model. The results presented in this work could lead to a better understanding and the eventual design of other ferromagnetic-insulator oxides in which oxygen vacancies could play an important role. 104 6 Sr-doping on Infinite-layer Nickelates Just a year ago before the pandemic, Sr-doped NdNiO2was found to be a high temperature superconductor. This so-pursued result has opened the possibility to study nickelates as cuprates analogs. In this chapter we study the effect of Sr-doping in the electronic and magnetic properties of infinite-layer nickelates as well as the nature of the holes. Our results show that doping induces a cuprate-like character on infinite-layer nickelates. The work presented here can be found in Ref. [V] from the List of Publications and it is part of a collaboration with Antía Botana’s group. 6.1 Introduction One of the major achievements in the field of oxides was the discovery of hightemperature superconductors (HTS) in 1986 [170]. Since then, many different strategies have been taken to try to decipher the origin of HTS [171]. Cuprate nanostructuring in which CuO2layers are arranged with different dopants in multiple ways has been probably the most studied route. However, a more illuminating approach could be to substitute Cu2+ with isoelectronic Ni1+: d9[172]. This nickel oxidation state is reported to occur in infinite layered nickelates of the form RNiO2(where R= La, Nd) [173–177]. It took more than 30 years to experimentally achieve a doped infinite-layer nickelate [178]. In 2019, Srdoped NdNiO2was reported to be a superconductor with Tc∼15 K and a dome-shaped doping dependence as occurs for cuprates [179]. In this chapter we will analyze from first principles how the electronic structure and magnetism of infinite-layer nickelates is affected by Sr doping. The parent phase of 112 nickelates (at d9filling), i.e. the stoichiometric phase, is quite different from that of cuprates. Experimental measurements show that RNiO2family is metallic and no antiferromagnetic order is found [173–176]. Apart from that, and unlike cuprates, previous electronic-structure calculations on RNiO2compounds report the presence of low-lying R-5d states crossing the Fermi level. Specifically, the Ni 3dz2 and R 5dz2hybridization creates a small spherical electron pocket centered at the Γpoint, and the Ni 3dz2and R dxy hybridization creates another electron pocket at the corners of the Brillouin zone. Electrons occupying these pockets originate from the otherwise filled 105 Adolfo Otero Fumega Figure 6.5: Evolution of the orbital resolved DOS for Ni-dz2and Ni-dx2−y2states upon increasing Sr-doping in LaNiO2, obtained using the AMF scheme with U= 5 eV. weight of the Ni dz2band around the Fermi level as Sr-doping is introduced, leaving a dominant dx2−y2contribution that makes Sr-doped 112 nickelates a more cuprate-like, single-band system (see Fig. 6.5). All of the above described trends imply that, as Sr dopants are introduced in RNiO2materials, some of their electronic-structure features become closer to those of the cuprates: low-spin dopant states, reduced charge-transfer energy, and a single Ni dx2−y2band around the Fermi level. Finally, the 4 ×4 supercells with an average d8.75 filling given by a 25% Sr substitution that allow, via clustering of all the Sr dopants, for one Ni in the cell to be nominally d8 (the Ni atom completely surrounded by first-neighbor Sr cations in Fig. 6.2). In this scenario, we find that a LS (S= 0) state is preferred in both RNiO2(R= La, Nd), for the Ni ion surrounded by Sr atoms even within FLL. In the AMF scheme the low-spin solution is not only the lowest in energy, but a high-spin solution does not even exist, as attempts to start the self-consistency procedure with a S=1 state of Ni2+ ion lead to a vanishing magnetic moment. The rest of the Ni atoms preserve the expected magnetic moments. Figure 6.6 contrasts the orbital resolved density of states for the nominally d8 Ni cation (that is surrounded by Sr atoms) in both a LS state and a HS state in this 4×4 supercell within FLL. We choose once again to show calculations for Sr-doped LaNiO2but the situation is identical in the Nd-material. In the low-spin state the t6 2g d2 z2configuration is clear with the two dx2−y2orbitals remaining unoccupied for both spin channels. In the high-spin case, the t2gorbitals are also completely occupied, but now one electron occupies the majority spin dz2and dx2−y2orbitals. The LS state for the Ni2+ is strongly connected to a reduction in the dz2character around the Fermi level. This can be seen in Fig. 6.6 where the Ni dz2band crossing the Fermi level for the HS state becomes fully occupied, well below the Fermi level for the LS state. This effect is concomitant with the reduction of the La-d self-doping effect upon increasing Sr-doping 112 6 Sr-doping on Infinite-layer Nickelates described above. Overall, the stable LS state solution we find gives then rise to an explicit cuprate-like scenario with planes of S=1/2 ions that are lightly doped with mobile lowspin S=0 ions, a configuration that is directly analogous to the low-spin S=0 Cu3+ ion situation, mediated by O-p holes as explained above. Figure 6.6: DOS of the nominally Ni2+ dopant in the 4×4 supercell with 25% Sr-doping in the two different spin states studied. HS (LS) configuration is shown in the left (right) panel. It can be noticed that the lower-energy LS configuration leads to a depletion of Ni dz2states around the Fermi level. We note that in other intensively studied nickelates (such as in La2NiO4) the high-spin (S=1) configuration of Ni2+ is favored [199, 200] and not the low-spin (S=0) configuration as we find here. La2NiO4is structurally different to RNiO2as it preserves apical oxygen atoms and has an octahedral environment for its Ni2+ cations. We show here that if the Ni2+ ions are forced into a square planar local environment as happens in the 112 materials, they prefer a low-spin d8(S=0) cuprate-like state instead, in agreement with recent experimental reports [201]. 6.4 Conclusions In this chapter we have analyzed how Sr-doping affects the electronic structure and magnetic state of infinite-layer nickelates RNiO2(R= La, Nd). We have found that including Sr reduces the self-doping effect by shifting the rareearth d bands up in energy away from the Fermi level. This leads to a more single-band- like picture, with the Ni dx2−y2band dominating. We have also analyzed the magnetic spin configuration as a function of doping (11% or 25%) and rare-earth cation in the framework of LDA+Ucalculations with two different double-counting methods. In any 113 Adolfo Otero Fumega case, a low-spin Ni2+ state is found to be the lower energy solution. Moreover, Sr-doping highly reduces the charge transfer energy between the Ni and O atoms. Our study shows substantial changes in the electronic structure of these nickelates upon Sr-doping. A cuprate-like electronic structure and magnetic state is obtained when doping these compounds. Therefore, this result suggests a cuprate-like description of the superconducting state that arises in doped infinite-layer nickelates. 114 7 Spectroscopy of Cs2CuCl4 In this chapter we study the electronic structure of Cs2CuCl4. This material has been discussed in the framework of frustrated antiferromagnet quantum spin liquids. A combined set of experimental techniques and Density Functional Theory has allowed us to determine the spectrum of this compound. This is mainly determined by the crystal field splitting that Cu2+ ions feel when they are placed in the strongly distorted tetrahedral environment of this compound. This unpublished work is part of a collaboration led by Dr. Santiago Blanco-Canosa. 7.1 Introduction In the previous chapter we have analyzed the electronic structure and magnetic properties of infinite-layer nickelates. We found that when doped with Sr their electronic structure can be understood to resemble that of the high-temperature superconducting cuprates. In this chapter we continue analyzing materials with promising exotic behaviour. However, we abandon at this point the study of oxides just before focusing on van der Waals layered materials in the next chapter. Therefore, this work will serve as an interlude between oxides and van der Waals transition metal compounds. It will also highlight the ability to combine experiment and theory to determine the electronic structure of a material. Quantum spin liquids (QSL) were hypothesized to occur in nature by Philip Anderson in 1973 [202]. This is a state of matter in which the electrons’ spins remain fluid-like, i.e., fluctuating without entering in a long-range order phase, even at 0 K [203]. Interestingly, the excitations of a QSL, known as spinons, are fractional, leading to a characteristic transport behaviour. Geometrically frustrated spin textures, like triangular, Kagomè and honeycomb lattices, are considered to be the most promising systems to find a QSL [204]. The magnetic order and its excitations will be determined by the competition between the different sign of the exchange interactions on the frustrated spin lattice. Among the solids showing frustrated magnetism, significant attention has been drown to Cs2CuCl4due to the experimental observations that report a two-dimensional QSL phase by means of specific heat measurements [205], electron spin resonance (ESR) [206] 115 Adolfo Otero Fumega Figure 7.1: (a) Orthorhombic unit cell of Cs2CuCl4. Green, red and blue balls denote Cs, Cu and Cl atoms, respectively. The tetrahedral coordination of Cu is highlighted. The magnetic exchange interactions commented in the text have been depicted. (b) Crystal field splitting of the Cu2+ ions in the distorted tetrahedral environment, where the eg orbitals lower their energy with respect to the t2g. The Jahn-Teller distortion further splits the energy levels with b1g(dx2−y2), eg(dxz,dyz), a1g(dz2) and b2g(dxy) symmetries. and neutron scattering [207], which uncovered an extensive two-spinon continuum. The ESR and neutron scattering measurements revealed a weaker interchain exchange coupling, J’/kB= 1.4 K, between Cu atoms than the intrachain interaction, J/kB=4.7 K, along baxis [208] (see Fig. 7.1a). Moreover, interlayer coupling in Cs2CuCl4is smaller than Jand J’ by more than one order of magnitude, J”= 0.13 K, in good agreement with the theoretical description of quasi-1D weakly coupled S=1/2 Heisenberg chains [209].1 1Be aware that in the literature antiferromagnetic exchange constants are commonly defined as positive when dealing with QSL models. 116 7 Spectroscopy of Cs2CuCl4 Besides, central to the study of spin exchange interactions in frustrated magnets is the precise knowledge of the crystal-field ground state symmetry and hopping energies of the S= 1/2 Cu2+ ions. In fact, the experimental determination of the crystal field and Jahn- Teller splittings, ∆CF and ∆JT has remained elusive. Moreover, the orbital occupation is commonly inferred from the oxidation state in an ionic picture, but is usually limited in covalent systems, where the strong overlap between dand pbands gives rise to charge transfer effects [210]. Indeed, the most noticeable effect of the different orbital occupation is highlighted in the manganese perovskites, RMnO3(R= rare earth) [211] and cobaltates [146], where the ground state properties are strongly influenced by the site symmetry of the transition metal ion and the electronic rearrangement within the d-shell [212]. Experimentally, the crystal field splitting, especially in 4fsystems, has been measured by inelastic neutron scattering [213, 214], but the information is sometimes hampered by phonons or small amount of crystals. In transition metal compounds, X-ray absorption spectroscopy (XAS) is traditionally used to probe the magnitude of the crystalline electric field splittings in the 3dshell [215] and the degree of the 3dhybridization in the ground state. Nevertheless, the multiplet effects are often hindered by the metal-ligand hybridization, which smears out the splittings of the 3dstates. Therefore, a comprehensive determination of the different energy scales is usually a tedious task and relies on several complementary techniques. On the other hand, band structure calculations have reported the electronic properties of Cs2CuCl4and predicted a strong dependence on the Jahn-Teller distortion, ∆JT , as a function of the correlation exchange functional [216]. However, a direct comparison between theory and experiment is still missing. In this work, we have circumvented these drawbacks by combining resonant inelastic X-ray scattering (RIXS), Density Functional Theory (DFT) and cluster calculations, which allows us to provide a satisfactory and comprehensive description of the Cs2CuCl4 electronic spectrum. 7.2 Experimental and Computational Methods Single crystals of Cs2CuCl4were grown by our experimental colleagues Dr. F. Rodríguez and Dr. S. Blanco-Canosa. The quality of the single crystals was checked by X-ray diffraction, resulting in lattices parameters a= 9.77 Å, b= 7.61 Å, c= 12.41 Å, and inelastic Raman scattering [217]. DFT calculations [17, 18] were performed using the all-electron, full-potential wien2k code [30] based on the augmented plane wave plus local orbital (APW+lo) basis set. The generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof [19] scheme was used for the exchange correlation functional, with a fully converged k-mesh of RmtKmax=7.0 and muffin-tin radii of 2.5, 2.22, 1.91 a.u. for Cs, Cu and Cl, respectively. RIXS experiments were performed by Dr. S Blanco-Canosa at the U41-PEAXIS beamline at BESSY II at 20 K and combined energy resolution, ∆E≈150 meV. 117 Adolfo Otero Fumega 7.3 Results and Discussion The crystal structure of Cs2CuCl4is orthorhombic with space group Pnma. Each Cu2+ atom is surrounded by 4 Cl−ions in a distorted tetrahedral coordination (Fig. 7.1a) which splits the d-levels into low energy 2-fold egand high energy 3-fold degenerate t2gorbitals (Fig. 7.1b). The isolated CuCl2− 4tetrahedra units are further Jahn-Teller distorted towards a lower symmetry D2dpoint group [218]. This environment provides a spatially anisotropic spin-1/2 triangular antiferromagnet of Cu-Cl chains on a geometrical bc plane, bonded by Cl−ions along the adirection. Our magnetic susceptibility measurements show no trace of magnetic order down to 5 K, which is in good agreement with a weakly frustrated magnet [219]. Figure 7.2 summarizes the DFT results of Cs2CuCl4for the simple unit cell depicted in Fig. 7.1a. Only the ferromagnetic state was considered, since there is only one inequivalent Cu atom. Calculations assuming a non-magnetic ground state directly give a non physical metallic solution. We have found that calculations on a 1×2×1supercell provide an antiferromagnetic order as the lower energy solution for the ground state, in agreement with the report of Foyevtsova et al. [216]. Nevertheless, we have opted for the ferromagnetic ground state in the simple cell, since our aim is not the determination of the magnetic order and does not have a substantial effect on the electronic energy levels. Therefore, the ferromagnetic solution allows for an easier visualization of the electronic structure and the direct comparison between DFT and the RIXS spectra. Between -4 and 0 eV (Fig. 7.2a) the DOS shows a clear contribution of the Cu 3dand Cl 3pbands, while the Cs bands do not contribute to the DOS at the Fermi level and, therefore, do not hybridize with Cu. The Cu and Cl bands are separated from the next unoccupied states by a gap of 4.5 eV, having mostly Cs character. We call this charge transfer gap and corresponds to the energy difference between the [CuCl4]2−cluster and the Cs+ions. Note that it must be differentiated from the traditional charge transfer gap defined between the transition metal dand the ligand pstates in oxides [220]. On the other hand, the d-d gap appears below 1 eV and represents the energy gap between the Fermi level and the empty dxy band. Figure 7.2b delves into the band structure of Cs2CuCl4. Bonding and antibonding bands show up between -4 and -3 eV and between -1 and 1 eV, respectively. The untangled electronic band structure gives a clue about the breaking of the orbital degeneracy sketched in Fig. 7.1b. In a tetrahedral coordination, the 3d9electronic configuration splits the crystal field generated by Cl−ions surrounding a Cu2+ ion into the energetically lower Cu eg(dx2−y2and dz2) doublet and higher Cu t2gtriplet (dxy,dxz, and dyz). Owing to the Jahn-Teller-like uniaxial distortion of the tetrahedron, the t2gtriplet is further split into the doubly degenerated dxz,dyz states (egirreducible representation) and the half-filled dxy states (b2g), and the dx2−y2(b1g) and dz2(a1g) become also energetically nonequivalent (Fig. 7.1b). Having carried out a theoretical description of the electronic structure of Cs2CuCl4, we proceed with the experimental analysis of the RIXS spectra. 118 7 Spectroscopy of Cs2CuCl4 Figure 7.2: (a) Projected DFT density of states (DOS) for the ferromagnetic structure of Cs2CuCl4, Cs (Cu, Cl) atom in green (red, blue). Both charge transfer and d-d energy gaps have been highlighted. (b) Energy bands for the majority (minority) spin channel on the left (right) panel. The d-orbital character of the Cu atom is depicted as red circles. 119 Adolfo Otero Fumega The fast improvement of the RIXS instrumentation in energy resolution of soft X-ray RIXS has allowed to study the low energy electronic properties of correlated oxides, giving detailed information of collective magnetic, charge and orbital excitations in oxides [221]. Figure 7.3b illustrates the case of the Cu2+ ion with a 3d9electronic configuration in a tetrahedral crystal field (Td), with a hole in an |xyistate. In the initial step of the RIXS process, a photon resonant at the Cu L3edge (2p→3dtransition) is excited from the ground state, |ii, into the 3dshell (intermediate state, |ni) filling the |xyiorbital and creating an excited core-hole state. In the final step (|fi), the core hole is annihilated via decay towards the ground state (elastic scattering, Eloss=0 eV) or an excited state (magnons, phonons, d-d transitions) [222]. (a) (c) (b) (d) (e) Figure 7.3: (a) Experimental (black) and calculated (red) X-ray absorption (XAS) of Cs2CuCl4. (b) Schematics of the RIXS process showing the initial, intermediate and final steps. (c) RIXS map plotting the energy loss, Eloss, against the incoming energy, Ein, for Cs2CuCl4at 20 K. The maximum intensity of the inelastic features appears at Ein= 931 eV. (d) Close up view of the RIXS scan at Ein=931 eV, highlighting the elastic, d-d excitations and the charge gap in the inset. (e) In-plane momentum dependence of the d-d excitations, showing no orbital dispersion. Figure 7.3c displays the incident energy (Ein) vs energy loss (Eloss) RIXS map. The maximum intensity of the resonant features is observed at Ein= 931 eV (Fig. 7.3d), 0.5 eV below the maximum of the L3absorption edge (Figure 7.3a), which consists on a featureless absorption band. In nice agreement with the DFT calculations, the charge transfer 120 7 Spectroscopy of Cs2CuCl4 excitations resulting from the [CuCl4]2−cluster to the Cs 6sstates are observed as a broad band at Eloss= 4 eV, matching the optical gap observed by absorption spectroscopy [223]. As shown in the inset of Fig. 7.3d, this charge transfer gap displays 2 broad bands at 3.7 and 4.2 eV, corresponding to transitions from the Fermi level and the dxy orbital above the Fermi level (Fig. 7.2b) to the upper Cs 6sbands. The region between Eloss=0.5-1.2 eV (matching the energy difference between the 2 broad charge transfer bands) corresponds to the so called optically forbidden crystal field d-d excitations, as widely reported in superconducting cuprates [221], and identified in the DFT calculations below ≈1 eV. The orbital assignment of these excitations was verified by comparing their energy with the DFT calculations. Considering that the 3dand 3pbands of Cu and Cl atoms are disentangled in energy from the rest of the bands (Fig. 7.2b), an expansion of the Bloch manifold can be performed in real space in terms of localized Wannier functions. The 68 Bloch bands in the energy window between -5 and 1 eV have been used to generate the localized and atom centered Wannier functions spanning such Bloch manifold [224]. A symmetry analysis of those Wannier functions allows to identify how the crystal field splitting affects the Cu 3dorbitals. As shown in Fig. 7.1b, we found that the dx2−y2(b1g) and the dxy (b2g) orbitals correspond to the lowest and highest energy levels, respectively. Apart from that, the strong orthorhombic distortion of the tetragonal environment inverts the energy levels of the doubly degenerate dxz and dyz (eg) and dz2(a1g). Therefore, d-d transitions originate from the decay of an electron from the dx2−y2,dz2,dxz and dyz orbitals due to the broken degeneracy of the 3dstates. To better understand the d-d excitations in the RIXS spectra, we have adopted the hole language, where the ground state represents a hole in the dxy orbital, hence, adyz orbital excitation corresponds to moving a hole from the dxy to the dyz orbital. Within the energy resolution of our experimental setup, we can discriminate the 3 orbital intratomic transitions; 2 sharp excitations at 0.67 and 1.02 eV, respectively and a shoulder at 1.21 eV (see Fig. 7.4b). Having considered only the tetrahedral point group,Td, this would result in a 10Dqvalue of 0.45 eV. Further, we see no orbital dispersion (Fig. 7.3e) indicating highly localized d-d excitations, as expected due to the strongly ionic character of Cs2CuCl4compound. Since these orbital d-d excitations are intra-atomic and well localized, they can be simulated within a full-multiplet calculation considering a single site of a Cu2+ ion of the D4hpoint group, isomorphic with the D2d[225], which is exemplified by a regular tetrahedron elongated along one of its C2axes. The RIXS simulations were carried out by Dr. S. Blanco-Canosa with the Quanty code [226, 227] including the Coulomb interactions and multiplet effects leaving free the radial integrals Dq(crystal field splitting), Dsand Dτ(distortions of the tetrahedra), which are the splitting terms for the Y0 2and Y0 4spherical harmonics, as input parameters in the calculation. Figure 7.4a displays the calculated RIXS map for Cu2+ within the D4hpoint group. The 3 orbital transitions corresponding to the dx2−y2,dz2and doubly degenerate dxz,dyz orbitals are clearly identified. As shown in Fig. 7.4b, the experimental and theoretical RIXS spectra is fairly well reproduced after normalization to the same height of the d-d features, with an additional instrumental and experimental broadening of 0.1 and 0.25 eV, respectively. 121