Iveta Zolotová
email: iveta.zolotova@tuke.sk

Marek Bundzel
email: marek.bundzel@tuke.sk

Tomáš Lojka
email: tomas.lojka@tuke.sk

Industry IoT Gateway for cloud connectivity

Keywords: Cyber-Physical System, Gateway, Internet of Things. 1

New approaches and technologies like Internet of Things (IoT), cloud computing and Big Data are giving rise to another industrial revolution. We propose here an implementation of an industrial gateway architecture adopting the idea of IoT, intelligent methods, Machine-to-Machine and Cyber-Physical Systems. The proposed gateway creates a virtual representation of the physical world scanning the technological layer's devices in real time. It creates a uniform communication interface for the heterogeneous technological layer, enables self-management of devices, diagnostics and self-reconfiguration to improve Quality of Service aided with cloud SCADA and MES services. We have tested the proposed gateway in an experimental setup with a programmable logic controller.

Introduction

Industry is a heterogeneous system consisting of various heterogeneous subsystems [START_REF] Atzori | The Internet of Things: A Survey[END_REF]. The subsystems must be interconnected to reach the complexity needed for an optimal performance of a plant [START_REF] Li | An Agile Governance Method for Multi-tier Industrial[END_REF][START_REF] Breiner | Run-time adaptation of a universal user interface for ambient intelligent production environments[END_REF]. This need for interconnection leads to the next industrial evolution step defined as the Industry 4.0. The term "Industrie 4.0" (German for "industry") refers to the fourth industrial revolution. It originates from a project in the high-tech strategy of the German government, which promotes the computerization of manufacturing. [START_REF] Zuehlke | SmartFactory-Towards a factory-of-things[END_REF][START_REF] Lydon | The 4th Industrial Revolution, Industry 4.0[END_REF]. The first industrial revolution increased production and decreased difficulty of manual work by using steam power. The second revolution did the same with the help of electric power and the third revolution brought electronics and IT to further automate production. The idea of the latest industrial revolution was presented at the Hanover Fair in 2011 [START_REF] Li | An Agile Governance Method for Multi-tier Industrial[END_REF].

Prof. Dr. Ing. Detlef Zühlke (Scientific Director of Innovative Factory Systems at the German Research Center for Artificial Intelligence) is sometimes called the "Father of Industry 4.0". The goal is to create more successful companies quickly capable to endure in the global competition [START_REF] Lydon | The 4th Industrial Revolution, Industry 4.0[END_REF].

Improvement of the interconnection of the industrial subsystems is related to the concept of the Internet of Things (IoT) [START_REF] Lydon | The 4th Industrial Revolution, Industry 4.0[END_REF], especially because of IoT aims to improve the interconnection of the digital and the physical world. Industry 4.0 in the contrary to the typical industry architecture significantly improves management of processes. IoT improves the processes of collecting, analyzing and processing the valuable and now easily accessible information originated in every part of the plant. The industrial IoT changes the machine-machine and machine-human interactions, comprehensive interoperability and intelligence aids to improve analytical description of the environment and better decisions can be taken. [START_REF] Zuehlke | SmartFactory-Towards a factory-of-things[END_REF], [START_REF] Brizzi | Bringing the Internet of Things along the manufacturing line: A case study in controlling industrial robot and monitoring energy consumption remotely[END_REF].

Fig. 1. Chronology of the industrial revolutions

Plants are comprised of different subsystems, modules, devices, machines operating with various communication protocols and interfaces. The problem is that the subsystems do not provide a unified connection to the technological layer and to the higher layer of the industrial hierarchical model. We propose implementation of the IoT in the Supervisory Control and Data Acquisition layer (SCADA). Typically the sensors of the technological layer produce large amounts of data; the SCADA layer collects it and provides the data to the Manufacturing Execution System (MES). This decreases interoperability between the technological layer and MES. The plants typically use strictly defined fixed interfaces for the communication within each layer between the layers of the hierarchy. This constraint does not allow ad-hoc behavior of the plant. Our goal is to create ad-hoc environment connecting technological layer with higher layers. The connections will be not constrained by fixed interfaces between the layers but automatically ad-hoc included and ad-hoc accessible for every layer.

2

Industry 4.0 and IoT.

Industry 4.0 brings greater flexibility, adaptability, self-organization, selfoptimization, self-diagnostics and leverages interaction between business and customer relationship, increases safety, reliability, autonomy and efficiency. The technological part is represented by a cyber-physical system [START_REF] Brizzi | Bringing the Internet of Things along the manufacturing line: A case study in controlling industrial robot and monitoring energy consumption remotely[END_REF]. Every product is individually described by the connection to business and customer layers. Integration of intelligent autonomy and monitoring with increasing automation into industry was important. It is represented by Cyber-Physical Systems (CPS) and IoT. [START_REF] Breiner | Run-time adaptation of a universal user interface for ambient intelligent production environments[END_REF]. Modern factories exceed the boundaries of the traditional ones. The factories are able to react to unpredicted situations autonomously and to select the optimal responses with regard to the plant or the market. The interconnection of the processes and plants is improved and the network security is better. The boundaries of the new industry are exceeded between the regions and between the different plants.

The Concept of Industrial Connection with IoT.

Industry 4.0 connects machines, workers, and factory systems into a network. New technological approaches of IoT, cloud computing, M2M (machine-to-machine), CPS and Big data are implemented into the network. Each of the technologies increases productivity. [START_REF] Spiess | SOA-based integration of the internet of things in enterprise services[END_REF][START_REF] Broy | Cyber-physical systems: Innovation durch softwareintensive eingebettete Systeme[END_REF]. We focus on IoT facilitating the interconnection of heterogeneous factory environments. The communication within the factory and the digital model of the real factory states are improved. The proposed interconnection enables better analysis and decision making. The position of IoT in the hierarchical model is depicted in Fig. 2.

Fig. 2. Industrial connection with IoT, CPS, cloud, and Big data

IoT creates a virtual (digital) representation of the physical world. The software applications access the virtual representations of the devices only. This solution is more suited for the world representation management when using various physical devices [START_REF] Schuh | Cyber-Physical Production Management[END_REF][START_REF] Jammes | Service-oriented paradigms in industrial automation[END_REF]. This layer must implement ad-hoc behavior and it must be manageable by central or distributed nodes [START_REF] Xiong | Multi-sensor management for information fusion: issues and approaches[END_REF][START_REF] Copie | From Cloud Governance to IoT Governance[END_REF].

The IoT layer runs a software surveyor for detection and identification of new devices. It enables adaptation of any physical layer and its integration into the factory architecture. Data/information is accessible to every application inside the plant or inside/outside the plant respectively. The devices detected for the first time are looked up in a devices database and included into the device ontology. [START_REF] Zuehlke | SmartFactory-Towards a factory-of-things[END_REF]. Measurement frequency, measurement range, device position etc. describe the device.

The connection to the technological layer is crucial and it influences the plant's Quality of Service (QoS) significantly. The technological layer of a heterogeneous system will not feature composite communication interfaces in Industry 4.0. The technological layer consists of CPS and smart products. These enable easier but not uniform connection to the technological layer. We propose a solution for the IoT connection to the technological layer here. We have designed a concept of IoT gateway that is described below. The proposed gateway supports communication with a cloud. SCADA and MES services may be implemented in the cloud as we have described in [START_REF] Lojka | Improvement of Human-Plant Interactivity via Industrial Cloud-Based Supervisory Control and Data Acquisition System[END_REF].

Concept Definition of the IoT Gateway

The important part of SCADA is the connection to the technological layer. This connection is crucial for data forwarding to the higher layers of the factory architecture. Therefore, we propose an IoT gateway connecting the physical devices and the higher layers of the Information and Control System (ICS). Regular and stochastic changes occur in a real word environment and the IoT gateway must adapt to the changes. The gateway's main role is the development of a virtual representation of the physical world enabling data collection and forwarding. The gateway's functionality must be robust to fulfill the requirements of all layers in the factory architecture. The IoT gateway distinguishes between the physical device and its virtual representation. That enables implementation of device management in the cases of failure or for saving energy when redundant measurements are taken. Additionally the virtual representation of the devices simplifies the connection between the cloud application and the physical devices (sensors/actuators). The gateway offers a unified connection to the virtual representations of the devices so that any application may use it.

The virtual representation of a device does not simply mirror the state of the device. Functionality for enhanced description of the device, its properties and measurements is added. [START_REF] Schuh | Cyber-Physical Production Management[END_REF][START_REF] Jammes | Service-oriented paradigms in industrial automation[END_REF]. The IoT gateway provides a complex representation of the instances.

We have defined the main tasks the IoT gateway will perform to reach the best QoS:

 Data forwardingthe gateway is responsible for data forwarding from the technological layer to the higher layers of the factory architecture. [START_REF] Xiong | Multi-sensor management for information fusion: issues and approaches[END_REF].  Gateway managementthe gateway is a key integration element in the entire factory architecture. A fixed and no-manageable gateway will cause fewer problems with adaptivity and implementation. [START_REF] Xiong | Multi-sensor management for information fusion: issues and approaches[END_REF].

 Device managementenables runtime configuration of devices, setting the statuses, functionality modes, errors acknowledgements and other. [START_REF] Xiong | Multi-sensor management for information fusion: issues and approaches[END_REF].  Data analysisthe IoT gateway monitors and analyses data in real time using statistics and artificial intelligence. This module cooperates with the device management to maintain QoS, reliability and fault tolerant operation. The faults are detected and the faulty devices are substituted by the available redundant devices. Data analysis with the device management enable self-configuration, adaptive and robust behavior with regard to the technological layer connected to the IoT gateway. The goal is to achieve autonomous management of devices without the need of human intervention aided by M2M communication.  Diagnosticsthe IoT gateway detects errors and faults in the entire technological layer and in the IoT gateway itself.

The operation of the proposed IoT gateway is divided into the management and runtime parts.

Management Part of the IoT Gateway

Management part of the IoT gateway detects, identifies and networks new devices in real-time. The IoT gateway creates a virtual instance of a device after having detected it. A definition of the device is needed to create the virtual instance. The definition may contain measurement ranges, device's location and power management information, lists of fault states, error states and error messages. If the definition of the virtual instance is inferred from the physical device it is uploaded to the device database [START_REF] Jammes | Service-oriented paradigms in industrial automation[END_REF]. The device manager will then create the virtual instance and forward it to the runtime. The virtual instance communicates with the device in the runtime. The data exchanges are based on events to reduce the communication load. However, the devices in Industry 4.0 should be manageable; at least the device's status must be settable.

This industrial solution enables to contact new devices. The contacted devices must provide metadata describing themselves to be correctly recognized. The IoT gateway creates the virtual instances based on the metadata. The virtual instances are frequently refreshed based on the corresponding events. The inactive physical devices are set to idle state by the device manager. Additional functionality is the management of the running virtual devices. This implies additional ability to update the physical and virtual instances in real time. The management part provides additional information to the data access point of the IoT gateway and the access point provides it to the higher layers (to MES, for example). The access point also collects information on configuration of virtual and physical instances from SCADA layer.

The diagnostics module enables management of the physical devices, the virtual instances and of the IoT gateway to ensure reliable control of the devices and data transfers from and to the technological layer. The diagnostics module also enables problem identification and invokes events describing the critical or the alarm states. In the case of a fault occurring on the physical device the virtual instance may link itself to another redundant physical device or self-reconfigure. The functionality of the faulty device is substituted by functionality of another device. The runtime part runs an event manager that keeps the physical devices and the virtual instances synchronized.

5

The Implemented Solution

We have not implemented the full functionality of the above described concept so far. We are try to identify and to use progressive approaches to improve the interconnection of the heterogeneous technological level and the cloud SCADA and MES services. We have designed a clustering method for grouping redundant devices. Substitutes for the faulty devices are easy to find and robustness is increased.

We have developed the management and the runtime parts of the IoT gateway. We have implemented the software agent scanning the defined subnets and identifies all the connected devices. The implemented management part selects the devices with the known communication interfaces and stores them in the device database. After that virtual instances for the devices are created and a synchronization interval is set up.

We have used PLC CompactLogix L23E from Rockwell Automation in our experiment. The PLC has been detected and identified by IoT gateway. An industrial IoT gateway must communicate via various communication protocols. The proposed IoT gateway uses a wide scale of industrial communication protocols; CIP (Common Industrial Protocol) was used for communication with the PLC in our experiment.

We have developed a GUI showing the IoT gateway processes using WPF and MS Visual Studio 2013. A list of the available devices found by the IoT gateway is shown in Fig. 5. The IP address of the PLC is shown in the Device register tab. Its virtual instance properties are displayed in the Instance properties tab. The virtual instance was automatically created upon reception of the metadata downloaded from the PLC. New industrial revolution based on CPS, IoT, cloud computing, and Big Data is on the way. [START_REF] Zuehlke | SmartFactory-Towards a factory-of-things[END_REF][START_REF] Copie | From Cloud Governance to IoT Governance[END_REF]. We have focused CPS and IoT for increased industrial autonomy, flexibility, efficiency, self-configuration, self-adaptiveness and robustness. Important part is communication in industry. Heterogeneity of the technological layer makes communication problematic. We have designed and implemented an IoT gateway creating a unified connection to the technological layer aided by IoT. This proposed gateway allows management of devices, identification of new devices, reading context data of the devices, self-reconfiguration, diagnostics and fault tolerance by creating virtual representations of the physical devices.

The future work includes implementation of artificial intelligence methods for classification of devices and their communication interfaces, for analysis of their behaviors and implementation of M2M to achieve self-managing behavior of the IoT gateways and their cooperation.

Fig. 3 .

 3 Fig. 3. Concept definition of IoT in the industry.

Fig. 4 .

 4 Fig. 4. Management and runtime parts of the IoT gateway

Fig. 5 .

 5 Fig. 5. Graphic representation IoT gateway processes

Acknowledgements. This publication is the result of the Project implementation: University Science Park TECHNICOM for Innovation Applications Supported by Knowledge Technology, ITMS: 26220220182, supported by the Research & Devel-opment Operational Programme funded by the ERDF (50%) and by grant KEGA -001TUKE-4/2015 (50%).


S Poncé

G Antonius

Y Gillet

P Boulanger

J Laflamme Janssen

A Marini

M Côté

X Gonze

The temperature dependence of electronic eigenenergies in the adiabatic harmonic approximation

Keywords: numbers: 63.20.kd, 65.40.-b, 71.15.Mb, 71.38.-k

The renormalization of electronic eigenenergies due to electron-phonon interactions (temperature dependence and zero-point motion effect) is important in many materials. We address it in the adiabatic harmonic approximation, based on first principles (e.g. Density-Functional Theory), from different points of view: directly from atomic position fluctuations or, alternatively, from Janak's theorem generalized to the case where the Helmholtz free energy, including the vibrational entropy, is used. We prove their equivalence, based on the usual form of Janak's theorem and on the dynamical equation. We then also place the Allen-Heine-Cardona (AHC) theory of the renormalization in a first-principle context. The AHC theory relies on the rigid-ion approximation, and naturally leads to a self-energy (Fan) contribution and a Debye-Waller contribution. Such a splitting can also be done for the complete harmonic adiabatic expression, in which the rigid-ion approximation is not required. A numerical study within the Density-Functional Perturbation theory framework allows us to compare the AHC theory with frozen-phonon calculations, with or without the rigidion terms. For the two different numerical approaches without rigid-ion terms, the agreement is better than 7 µeV in the case of diamond, which represent an agreement to 5 significant digits. The magnitude of the non rigid-ion terms in this case is also presented, distinguishing specific phonon modes contributions to different electronic eigenenergies.

I. INTRODUCTION

Electronic properties of solids and nanostructures can be computed from first principles with varying accuracies. In particular, the widespread GW approximation 1 within many-body perturbation theory describes electronic bandgaps with errors in the range of 0.1-0.3 eV with respect to experiment 2 . Excitonic effects can be added based on the Bethe-Salpeter equation (BSE) 3 . However, a crucial ingredient is often disregarded: the electron-phonon interaction. Incidentally, the most advanced self-consistent GW calculations usually lead to overestimated bandgaps 2 , and most of the remaining discrepancy might be due to the electron-phonon interaction. Indeed, the influence of the lattice vibration at 0 Kelvin, known as the zero-point motion renormalization (ZPR), can be as large as 0.37 eV for the indirect bandgap of diamond 4,5 . This correction often leads to a reduction of the bandgap and hence might be crucial to correct the overestimation of self-consistent GW approximation.

The study of the direct effects of the electron-phonon interaction on the electronic structure has a long history. From the fifties to the late eighties, they were investigated and computed in a semi-empirical context. It was first recognized that the temperature dependence of the electronic energies has two different origins: the volume expansion (studied by Shockley and Bardeen 6 ) and the electron-phonon interactions at constant volume. In fact, the effect of the electron-phonon interaction at constant volume is usually the major contribution, and proves to be the most difficult to compute from first principles. This contribution is the focus of this paper [START_REF]Although first-principle studies of thermal expansion are not very frequent[END_REF] . In a pioneering work, Fröhlich introduced a model Hamiltonian that includes these interactions [7][START_REF] Fröhlich | Proc. R. Soc. Lond. A[END_REF][START_REF] Fröhlich | [END_REF] . However, his approach leads to overscreening of phonon frequencies and involves empirical parameters 10 .

In parallel, the self-energy contribution to electronic eigenenergy renormalization due to thermal vibrations was introduced by Fan 11,12 . His theory has no adjustable parameters and is based on the first-order perturbed Hamiltonian. Later, Cohen 13 used the Fan self-energy to compute the temperature dependence of the germanium bandgap. Also, ideas from electron diffraction theories based on thermally averaged nuclear potentials lead Antončík 14 and others [15][16][17] to develop empirical Debye-Waller (DW) corrections to the nuclear potential. As these two lines of thought (Fan and DW) were developed independently of each other, only one of them was usually included in calculations.

At about the same time as Antončík, Brooks [START_REF] Brooks | Theory of the Electrical Properties of Germanium and Silicon[END_REF] emphasized that the electronic gap is actually a free energy difference with respect to varying occupation numbers of electronic and phononic levels. Furthermore, Allen and Hui 19 highlighted the equivalence between the action of the phonon population (atomic position fluctuations) on the electronic eigenenergies and the action of the electronic occupations on the phonon eigenfrequencies. This equivalence, that was later called Brook's theorem, was used in the eighties by several authors to discuss the temperature dependence of eigenenergies 19,20 .

In 1974, Baumann 21 first suggested that both the Fan self-energy and DW terms were needed to describe the influence of lattice vibrations on the electronic eigenenergies. Two years later, Allen and Heine 22 rigorously unified the theory and made the DW term translationally invariant. Their approach, combined with the rigidion approximation (RIA), which is valid for their semiempirical model, allows for a re-writing of the problem in terms of first-order derivatives of the effective potential only. Calculations of electron-phonon renormalization were then led by Cardona and coworkers [23][24][25][26] , including Allen, based on physically motivated models 27,28 or on rigid-ion pseudopotentials approximations and empirical phonon models. The resulting approach is now called the Allen-Heine-Cardona (AHC) theory.

Until then, none of the calculation were based on first principles. In 1989, King-Smith et al 20 computed the temperature-dependent bandgap of Silicon using Density Functional Theory (DFT) [START_REF] Martin | Electronic Structure. Basic Theory and Practical Methods[END_REF] , by evaluating the change of phonon frequencies due to electronic occupations and invoking Brook's theorem.

It took more than one decade before other firstprinciples calculations were performed. Such studies relied on widely varying formalisms and methods, that can be broadly classified in three types, each with their distinct advantages and drawbacks. First, the temperaturedependent eigenenergies can be computed as a time average of the band gap obtained using first-principles molecular dynamics (MD) simulations. Using this method, Franceschetti [START_REF] Franceschetti | [END_REF] studied Si nanocrystals, Kamisaka et al 31 studied CdSe and PbSe quantum dots and Ibrahim 32 the temperature dependence of the optical response of GaAs. Ramírez et al 33,34 simulated the temperature dependence of diamond and 3C-SiC bandgap based on pathintegral molecular dynamics (PIMD). This approach has the interesting characteristics that it includes effects beyond the harmonic approximation. Moreover, while normal molecular dynamics, which involves the harmonic approximation, wrongly delivers a classical Boltzmann statistics for phonons, the more computationally demanding PIMD properly includes nuclear quantum effects and delivers Bose-Einstein statistics for phonons. Thus, with PIMD, zero-point motion effects are observed. However, this MD (or PIMD) method is most suited for finite systems. Indeed, for solids, a supercell has to be used to sample the phonon wavevector space, so that the eigenenergies are not well separated from each others and only the band edges can be clearly identified.

A second method uses frozen phonons (FP): the computation of the change of eigenenergies due to atomic displacements along the normal modes is followed by a Bose-Einstein weighted sum of the contribution of each mode. In 2005, Capaz et al 35 studied the temperature dependence of carbon nanotube bandgap within the framework of a tight-binding method. Patrick et al 36 examined diamondoids, and Han and Bester 37 studied various semiconductor nanoclusters, still using the FP method but this time with DFT simulation. Anharmonic electron-phonon contribution to the temperature dependence of the indirect band gap of diamond were also studied by Monserrat et al 38 with the same methodology. Recently, Antonius et al 39 , still relying on the FP method, computed the renormalization of the diamond bandgap within the GW approximation and observed a large increase of the renormalization with respect to DFT, in better agreement with experimental values. This result is in line with earlier estimations of many-body effects on electron-phonon coupling [40][41][42][43][44][45] .

As a third approach, the diagrammatic method of many-body perturbation theory, from which the AHC approach originates, allowed Giustino et al 46 to compute the ZPR and the temperature dependence of the diamond bandgap with Wannier functions in the Density-Functional perturbation theory (DFPT) 47,48 framework. Marini et al [49][50][51][52] focused on the dynamical effects, beyond the adiabatic approximation, which are absent from the two previous approaches (MD and FP).

There has been some confusion about the (non)equivalence of these three approaches, in the firstprinciples context. Although the first (MD) and second (FP) one are equivalent when considered within the adiabatic harmonic approximation, the third one (AHC) is equivalent to the MD and FP only when the rigidion approximation is valid, which is not the case in the first-principle context. Indeed, Gonze et al 53 pointed out that the FP and AHC formalisms differ by nondiagonal Debye-Waller terms, and computed these for simple diatomic molecules. In some cases, the nondiagonal Debye-Waller terms were as large as the direct terms. They also reformulated the AHC theory using Sternheimer equations instead of summations over empty states, which led to a significant speed-up of their calculations. The difference between AHC and FP was also examined in the above-mentioned study by Antonius et al 39 , where the global effect of the RIA on the ZPR was found to be rather small for the bandgap of diamond. Recently, a thorough validation study by Poncé et al 54 , comparing different first-principles codes, allowed to resolve a persisting disagreement on the value of the ZPR for the direct bandgap of diamond and established a value of -0.41 eV from AHC formalism on top of DFT 54 . The latter study also revealed the very slow convergence with respect to the number of q-points (phonon wavevectors). Such a slow convergence is not restricted only to the AHC method: frozen phonons of many wavevectors should be taken into account for obtaining properly converged FP calculations as well, while huge supercells should be used in the case of MD-based approaches.

To complete this literature overview, the phononinduced lifetime broadening of the electronic states de-rived from the imaginary part of the Fan self-energy was investigated by Lautenschlager 55,56 in a semiempirical context, and, more recently, by Giustino 57 and Restrepo 58,59 with a first-principle implementation.

In this paper, we will clarify or establish links between the different approaches at the first-principle level, for semiconductors and insulators. Strictly speaking, because of the adiabatic approximation, our theory does not apply to metals, as the phonon frequencies cannot be neglected with respect to electronic excitations. The same limitation is also encountered for semiconductors and insulators with a temperature sufficiently high to create a non-negligible population of holes and conduction electrons. So, provided the adiabatic approximation is valid, we establish Brook's theorem, and provide a detailed analysis of the difference between the AHC approach and the FP approach, elaborating on the brief results presented in Ref. 53. Detailed DFT numerical results for diamond will also be provided, going, for DFT, further than Ref. 39.

The structure of the paper is as follows. Section II examines Brooks theorem which relies, at the first-principle level, on Janak's theorem generalized to the case where the free energy including the vibrational entropy is used, instead of a purely electronic expression. Section III links the equations of section II with the AHC theory, in the periodic case (with notations suitable for later practical implementation) and also explores rigorously the rigidion approximation. In particular, one distinguishes the Fan contribution, the diagonal rigid-ion Debye-Waller contribution, the diagonal non-rigid-ion Debye-Waller contribution, and the non-diagonal Debye-Waller contribution. The sum over a large number of bands, present in the AHC theory, can be reduced drastically, by complementing it with an expression based on the projection over high energy bands of the first-order wavefunctions, in the spirit of the Sternheimer equation. Section IV establishes the connection between those equations and an equivalent finite-difference approach. Then, this connection is used in section V to validate the theory of section III as well as its implementation by comparison with finite-difference calculations for the case of diamond. This section also allows one to assess the importance of the rigid-ion approximation in the case of periodic solids.

All equations derived in this work are expressed in Hartree atomic unit where m e = = e = 1.

II. TEMPERATURE-DEPENDENT ELECTRONIC EIGENENERGIES IN THE ADIABATIC HARMONIC APPROXIMATION

Brook's "theorem" states that the electron-phonon interaction can be regarded either as the shift in an electronic eigenenergy (labeled n) when a phonon is added in a phononic mode (labeled m) or, equivalently, as the shift in the m phonon mode eigenfrequency when an electron is placed in the n level. Namely, [START_REF] Brooks | Theory of the Electrical Properties of Germanium and Silicon[END_REF][19][20]

∂ε n ∂n m = ∂ω m ∂f n , (1)
ε n and f n being respectively the eigenenergy and the thermal average of the occupation of the electronic state n, while ω m and n m are the eigenfrequency and the thermal average of the phononic occupation of the phonon mode m, respectively. In this section, we detail this equality in a firstprinciples context: the phonons are obtained in the harmonic adiabatic approximation, from the interatomic force constants, themselves second-order derivatives of the Born-Oppenheimer total energy in which the electronic occupations effects can be traced explicitly.

We decided, for pedagogical reasons, to work in this section with isolated systems but the extension to periodic systems is straightforward using the convention of appendices A 2 and A 3.

A. Fluctuations of atomic positions

The eigenfrequencies ω m and unitless mass-scaled eigendisplacements ξ m,κα for the phonon mode m can be obtained from the dynamical equation 60 κα

D κα κ ′ γ ξ m,κα = ω 2 m ξ m,κ ′ γ , (2)
where κ and κ ′ label atoms in the unit cell, α and γ label the cartesian spatial dimensions, and D κα κ ′ γ are the mass-scaled interatomic force constants. The mass-scaled quantities are expressed as follows

ξ m,κα = M κ U m,κα , D κα κ ′ γ = 1 √ M κ C κα κ ′ γ 1 √ M κ ′ , (3)
where M κ is the mass of atom κ, U m,κα is the phonon eigendisplacements, and

C κα κ ′ γ = ∂ 2 E BO ∂R κα ∂R κ ′ γ , (4)
are the inter-atomic force constants (IFC), where E is the Born-Oppenheimer energy of the system (which excludes the kinetic energy of the nuclei), computed from first principles. Using the fact that the mass-scaled eigendisplacements ξ m,κα are orthonormal and complete, the phonon frequencies can be obtained from

ω 2 m = κα κ ′ γ U * m,κ ′ γ C κα κ ′ γ U m,κα (5)
and the normalisation of phonon eigendisplacements is

κα M κ U * m,κα U m ′ ,κα = δ mm ′ . (6)
Within the harmonic approximation, where normal modes are decoupled from each others, we express the static temperature dependence of the eigenenergies ε n (T ) as a thermal average of the value of the temperatureindependent, position-dependent eigenenergies εn [zU], where zU denotes generically a displacement from the equilibrium atomic positions

ε n (T ) εn [zU] (T ) (7) = 3N m 1 Z m sm e -(sm+1/2)ωm k B T χ * m,sm (z)ε n [zU m,κ ]χ m,sm (z)dz, ( 8
)
with N the number of atom, T the temperature, s m the integer occupation of phonon mode m,

Z m = sm e -(sm+1/2)ωm k B T
the mode-partition function, k B the Boltzmann's constant, χ m,sm (z) the phonon eigenfunctions, z the spatial coordinate associated with a phonon mode and where bold symbols like U m,κ denote cartesian vectors.

We can expand the eigenenergies εn [zU m,κ ] of Eq. ( 7) in a Taylor series up to second order in z, since we are working within the harmonic approximation

εn [zU m,κ ] = εn [0] + z d dz εn [zU m,κ ] z=0 + 1 2 z 2 d 2 dz 2 εn [zU m,κ ] z=0 , (9)
and insert it in Eq. (7). Using the properties of phononic wavefunctions in an harmonic potential (see appendix A 1 for more details), we obtain

ε n (T ) = εn [0] + 1 2 3N m 1 Z m ω m d 2 dz 2 εn [zU m,κ ] z=0 sm e -(sm+1/2)αm 2s m + 1 2 (10)
where we have defined α m ωm kB T . For sake of brevity, we introduce the notation ∆ε n (T )

ε n (T ) -εn [0]
. Using the properties of geometrical series to evaluate the sum over s m , we get

∆ε n (T ) = 1 2 3N m 1 e -αm /2 1-e -αm ω m d 2 dz 2 εn [zU m,κ ] z=0 e -αm/2 e -αm (1 -e -αm ) 2 + 1 2(1 -e -αm ) , (11)
which simplifies to

∆ε n (T ) = 3N m 1 2ω m d 2 dz 2 εn [zU m,κ ] z=0 n m (T ) + 1 2 , (12)
where we have introduced the Bose-Einstein distribution

n(T ) = 1 e ωm k B T -1 . ( 13
)
We can consider the phonon occupation numbers n m (T ) as independent variables in this expression, in which case

∂ε n ∂n m ∂ε n (T ) ∂n m (T ) = 1 2ω m d 2 dz 2 εn [zU m,κ ] z=0 , (14)
where we have defined the short-hand notation on the left hand side of the equation because each term in the central part of the equation depends on the temperature but their ratio does not. The slope of the eigenenergies with the phononic occupation is therefore independent of temperature and we will remove the explicit dependence of the terms when considering the slope. Eq. ( 12) can thus be written

∆ε n (T ) = 3N m ∂ε n ∂n m n m (T ) + 1 2 . ( 15
)
Then, the zero-point motion contribution (ZPR) is simply the shift of the eigenenergies at T = 0, that is

∆ε n (T = 0) = 1 2 3N m ∂ε n ∂n m , (16)
highlighting that ε n (T = 0) = ε n [0]
. Additionally, from Eqs. ( 5) and ( 6) we can see that the phonon frequency goes as

M -1/2 κ . The U m,κ goes as M -1/2 κ
and the zU m,κ must have the dimension of length because it is an ionic displacement (see Eq. ( 12)). Therefore, the z must have the dimension of M 1/2 κ times length. All of this leads to the fact that the ZPR of Eq. ( 12) goes as (ω m M κ ) -1 , i.e. as M -1/2 κ . This "isotopic effect" allows for experimental measurements of the zero-point renormalization by substituting atoms with heavier isotopes, as explained in the review paper of Cardona and Thewalt 5 .

B. Eigenenergies as derivatives with respect to electronic occupation numbers fn

Following Janak, we extend to fractional occupations the first-principle Born-Oppenheimer energy E BO . Janak's theorem 61 then gives

ε n = ∂E BO ∂f n , (17)
where it has to be noted that the Janak theorem breaks down, within DFT, if the exact exchange-correlation functional is used. In that case, the total energy is not a continuous function of the electronic occupation anymore and eigenenergies must be defined as difference of total energies with integral occupations numbers.

We now complement the Born-Oppenheimer energy with phonon energy and entropy at the harmonic level. The energy becomes

E(T ) = E BO + E vib (T ) = E BO + 3N m ω m n m (T ) + 1 2 , ( 18
)
where E BO is the total energy without phonon and electron-phonon contributions. Taking into account the vibrational entropy gives a Helmholtz free energy 62

F (T ) E(T ) -T S vib (T ), (19)
where S vib is the vibrational entropy,

S vib (T ) = k B 3N m (1 + n m (T )) ln(1 + n m (T )) -n m (T ) ln(n m (T )) . ( 20
)
We now show that the temperature-dependent eigenenergies can be obtained from the extension of Janak's theorem to finite phonon temperature

ε n (T ) = ∂F (T ) ∂f n . (21)
For sake of simplicity, we neglect the dependence of the electronic occupations f n on electronic temperature. Actually, the explicit treatment of the electron system at finite temperature (e.g. using the Mermin functional 63 ), supposing (wrongly) the adiabatic approximation to be still valid, would not change the remaining of the paper. Using this definition of ε n (T ), Eq. ( 18), Eq. ( 19) and taking into account the dependence of the phonon frequencies on electronic occupation numbers as well as the dependence of phonon occupation numbers on electronic occupation numbers, the change of eigenenergies due to electron-phonon interaction becomes

∆ε n (T ) = 3N m ∂ω m ∂f n n m (T ) + 1 2 + ω m ∂n m (T ) ∂f n -T ∂S vib ∂n m (T ) ∂n m (T ) ∂f n . (22)
Substituting Eq. ( 20) for S vib into Eq. ( 22) gives

∆ε n (T ) = 3N m ∂ω m ∂f n n m (T ) + 1 2 + ω m ∂n m (T ) ∂f n -k B T ∂n m (T ) ∂f n ln 1 + n m (T ) n m (T ) . (23)
Using Eq. ( 13), the last 2 terms in the sum cancel out.

We thus obtain

∆ε n (T ) = 3N m ∂ω m ∂f n n m (T ) + 1 2 . ( 24
)
This is a first important result of the present paper. To the authors knowledge, it was never derived starting from the free energy Eq. ( 19) and the finite temperature extension of Janak's theorem Eq. (21). Identification of Eq. ( 24) with Eq. ( 15) obviously yields

∂ε n ∂n m = ∂ω m ∂f n . (25)
This link can be more rigorously established as follows. The derivative of the phonon frequency is retrieved from deriving Eq. ( 5) with respect to electronic occupation

2ω m ∂ω m ∂f n = κα κ ′ γ ∂C κα κ ′ γ ∂f n U * m,κ ′ γ U m,κα , (26)
where the derivative of the displacements U m,κα with respect to the occupations f n does not contribute thanks to the Hellmann-Feynman theorem 64,65 . The temperature dependence of the eigenenergies is obtained by substituting the preceding result for ∂ωm ∂fn into Eq. ( 24)

∆ε n (T ) = 3N m 1 2ω m κα κ ′ γ U * m,κ ′ γ ∂C κα κ ′ γ ∂f n U m,κα n m (T ) + 1 2 . (27)
Using Janak's theorem, we can reformulate the derivative of the IFC into derivatives of eigenenergies

∂C κα κ ′ γ ∂f n = ∂ ∂f n ∂ 2 E BO ∂R κα ∂R κ ′ γ = ∂ 2 ε n ∂R κα ∂R κ ′ γ . ( 28
)
Substituting the above expression into Eq. ( 27), we obtain

∆ε n (T ) = 3N m 1 2ω m κα κ ′ γ U * m,κ ′ γ ∂ 2 ε n ∂R κα ∂R κ ′ γ U m,κα n m (T ) + 1 2 . (29)
The double sum over atomic position displacements is actually the second-order derivative of the eigenenergy with respect to the normal mode,

κα κ ′ γ U * m,κ ′ γ ∂ 2 ε n ∂R κα ∂R κ ′ γ U m,κα = d 2 dz 2 εn [zU m,κ ] z=0 , (30)
so that we recover, from Eqs. ( 14), ( 25) and (30) the following relations

∂ε n ∂n m = 1 2ω m d 2 dz 2 εn [zU m,κ ] z=0 = 1 2ω m κα κ ′ γ U * m,κ ′ γ ∂ 2 ε n ∂R κα ∂R κ ′ γ U m,κα = ∂ω m ∂f n . (31)
We thus obtain a demonstration of the Brook's theorem in the first-principle context.

C. The Fan and Debye-Waller contributions

In order to establish the links with the Fan, Debye-Waller and AHC approaches, we now analyze Eq. ( 29) in more detail, focusing on the second-order derivative of the eigenenergies with respect to two atomic displacements. We can obtain it from perturbation theory. We start from

ε n = Ψ n | Ĥ |Ψ n , (32)
and differentiate it, also using the Hellmann-Feynman theorem. At the equilibrium geometry,

∂ε n ∂R κα = Ψ (0) n ∂ Ĥ R κα Ψ (0) n . (33)
Eq. ( 33) can be derived a second time, with respect to another atomic displacement. An equivalent result can be obtained by switching the two atomic displacements. Both can be combined and deliver an expression that is real and explicitly symmetric with respect to the indices κα and κ ′ γ

∂ 2 ε n ∂R κα ∂R κ ′ γ = Ψ (0) n ∂ 2 Ĥ ∂R κα R κ ′ γ Ψ (0) n + 1 2 ∂Ψ n ∂R κα ∂ Ĥ ∂R κ ′ γ Ψ (0) n + (κα) ↔ (κ ′ γ) + (c.c.) , (34)
where (κα) ↔ (κ ′ γ) stands for the previous term in which the indices κα and κ ′ γ have been exchanged, and where (c.c.) stands for the complex conjugate of the previous term.

The contribution from the second-order perturbation of the Hamiltonian,

∂ 2 Ĥ ∂RκαR κ ′ γ
, gives the Debye-Waller (DW) term 14,66 of the semi-empirical approach. The other bracketed term originates from the first-order modifications of the wavefunction and corresponds to the Fan term when considered in many-body perturbation theory. As mentioned in the introduction, the complementarity of the two terms for the description of the eigenenergy renormalization due to the electron-phonon interaction, although obvious in the present derivation, was first shown in 1974 by Baumann 21 .

III. FROM THE ADIABATIC APPROXIMATION TO THE ALLEN-HEINE-CARDONA THEORY

In this section, we relate Eq. ( 29) to the AHC theory, as formulated 80 by Giustino et al. 46 . We also carefully treat the RIA and unveil the terms it neglects. Furthermore, we discuss the translational invariance and its consequences. In order to focus on solids (in view of the application and testing for diamond and solid-state systems in general), we work with periodic systems and introduce the related definitions. More information in this respect is provided in appendix A 2 and A 3.

A. Wavevectors and translational invariance

We introduce a specific notation for the derivative of an arbitrary quantity X, that depends on the atomic coordinates, with respect to a collective displacement of atoms characterized by a wavevector q

∂X ∂R κα (q) = 1 N BvK l e iq•R l ∂X ∂R lκα , (35)
where R lκα is the coordinate along the α axis of the atom κ in the cell l, and N BvK is the number of primitive cells of the periodic system defined by the Born-von Karman boundary conditions [START_REF] Born | Dynamical Theory of Crystal Lattices[END_REF] .

The eigendisplacement vectors are solution of the dynamical equation

κα C κα κ ′ γ (q)U m,κα (q) = M κ ′ ω 2 mq U m,κ ′ γ (q), (36)
where the Fourier transform of the IFC appears

C κα κ ′ γ (q) = 1 N BvK ll ′ C lκα l ′ κ ′ γ e -iq•(R l -R l ′ ) = l ′ C 0κα l ′ κ ′ γ e iq•R l ′ =N BvK ∂ 2 E BO ∂R κ ′ α (-q)∂R κ ′′ β (q) , (37)
where R l is a translation vector of the Bravais lattice. The eigendisplacement vectors fulfill the following normalization condition

κα M κ U * m,κα (q)U m ′ ,κα (q) = δ mm ′ . ( 38
)
All properties of a crystal, including its eigenenergies and their derivatives, must be invariant upon a uniform translation δ. Therefore,

ε nk [{R lκ }] = ε nk [{R lκ + δ}] (39) ∂ε nk ∂R κ ′ α (Γ) [{R lκ }] = ∂ε nk ∂R κ ′ α (Γ) [{R lκ + δ}]. (40)
In Eq. 40, we have taken the derivative with respect to a collective displacement of atoms that does not break the translation symmetry (q = Γ) to avoid any problem with the Bloch theorem, and keep the Bloch notation nk.

By Taylor expanding the right hand side of Eq. ( 40), we obtain

∂ε nk ∂R κ ′ α (Γ) [{R lκ }] = ∂ε nk ∂R κ ′ α (Γ) [{R lκ }] + κ ′′ β δ β ∂ 2 ε nk ∂R κ ′ α (Γ)∂R κ ′′ β (Γ) [{R lκ }] + O(δ 2 ), (41)
where δ β is the β component of the vector δ. For the equality to hold for all δ, every term of order one and higher in the series must be identically zero

β δ β κ ′′ ∂ 2 ε nk ∂R κ ′ α (Γ)∂R κ ′′ β (Γ) [{R lκ }] =0 ∀ δ β ∈ R ⇒ κ ′′ ∂ 2 ε nk ∂R κ ′ α (Γ)∂R κ ′′ β (Γ) [{R l κ }] =0. ( 42
)
B. The temperature dependence in the adiabatic harmonic approximation for the solid periodic case

Eq. ( 29) can be generalized to the periodic case, with a discretized integral over q (N q is the number of wavevectors used to sample the Brillouin zone)

∆ε nk (T ) = 1 N q q 3N m 1 2ω mq κα κ ′ γ ll ′ ∂ 2 ε nk ∂R lκα ∂R l ′ κ ′ γ e -iq•(R l -R l ′ ) U * m,κ ′ γ (q)U m,κα (q) n mq (T ) + 1 2 . (43)
Similarly, the generalization of Eq. ( 15) leads to

∆ε nk (T ) = 1 N q q 3N m ∂ε nk ∂n mq n mq (T ) + 1 2 , (44)
where

∂ε nk ∂n mq = 1 2ω mq κα κ ′ γ ll ′ ∂ 2 ε nk ∂R lκα ∂R l ′ κ ′ γ e -iq•(R l -R l ′ ) U * m,κ ′ γ (q)U m,κα (q) (45)
and we will focus on the latter quantity, which represents the change of eigenenergy due to a specific phonon mode.

To split this expression in a Fan and a Debye-Waller contribution, we substitute the extension to periodic sys-tem of Eq. ( 34) in it and retrieve

∂ε nk ∂n mq = 1 2ω mq κα κ ′ γ U * m,κ ′ γ (q)U m,κα (q) u (0) nk ∂ 2 Ĥk,k ∂R κα (-q)∂R κ ′ γ (q) u (0) nk + 1 2 ∂u nk ∂R κα (q) ∂ Ĥk,k ∂R κ ′ γ (q) u (0) nk + (κα) ↔ (κ ′ γ) + (c.c.) , (46)
where Ĥk,k is defined through Eq. (A11) applied to the Hamiltonian. This allows us to introduce the following notation for the DW and Fan contributions related to band n and wavevector k (we skip the n and k indices, which should not be confusing in the present context)

D κα κ ′ γ (q) u (0) nk ∂ 2 Ĥk,k ∂R κα (-q)∂R κ ′ γ (q) u (0) nk (47)
F κα κ ′ γ (q) 1 2 ∂u nk ∂R κα (q) ∂ Ĥk,k ∂R κ ′ γ (q) u (0) nk +(κα) ↔ (κ ′ γ) + (c.c.) . (48)
The change of eigenenergy due to a specific phonon mode, Eq. ( 46), thus becomes

∂ε nk ∂n mq = 1 2ω mq κα κ ′ γ D κα κ ′ γ (q) + F κα κ ′ γ (q) U * m,κ ′ γ (q)U m,κα (q). ( 49
)
With the same notations, the translational invariance Eq. ( 42) reads

κ ′ D κα κ ′ γ (Γ) + F κα κ ′ γ (Γ) = 0. ∀ α, γ, κ (50)
In particular, one can multiply this expression by any expression independent of κ ′ , and still get zero. This gives us some freedom on the form of the added terms which are chosen in such a way that the D terms will cancel out in the rigid-ion approximation (see next subsection). Moreover, we can also sum over κ instead of κ ′ in Eq. ( 50) in such a way that the resulting expression is Hermitian

∂ε nk ∂n mq = 1 2ω mq κα κ ′ γ D κα κ ′ γ (q) + F κα κ ′ γ (q) U * m,κ ′ γ (q)U m,κα (q) -D κα κ ′ γ (Γ) + F κα κ ′ γ (Γ) 1 2 U * m,κγ (q)U m,κα (q) + U * m,κ ′ γ (q)U m κ ′ α (q) . ( 51
)
Beyond the adiabatic and harmonic approximations, we have not made any additional approximation until now. The Debye-Waller term that we have obtained, Eq. ( 47), invokes the second-order derivative of the firstprinciple Hamiltonian. In most DFPT procedures, calculations of the phonon band structure rely only on the evaluation of the first-order derivative of the selfconsistent DFT Hamiltonian and wavefunctions. Indeed, while the expression leading to such phonon band structure calculations include a second-order derivative with respect to the non-self-consistent electron-ion potential, the latter term does not depend on second-order derivative of the wavefunction and Hamiltonian. Thus, the Debye-Waller term is not a by-product of a phonon band structure calculation. We will now show how the rigid-ion approximation allows us to compute the Debye-Waller term without computing the second-order derivative of the self-consistent first-principle Hamiltonian.

C. The rigid-ion approximation

In the case of semi-empirical potentials, it is natural to suppose that the Hamiltonian depend on potentials created independently by each nucleus, screened by electrons attached to them. In this case, as pointed by Allen and Heine 22 , the numerical burden of computing a second-order derivative of the Hamiltonian can be completely avoided. A rigid-ion Hamiltonian has the following form

Ĥri = T + lκ V κ (r -R lκ ). ( 52
)
so that its mixed (off-site) second-order derivatives vanish

∂ 2 Ĥri ∂R lκα ∂R l ′ κ ′ γ = 0 if κ = κ ′ as well as if l = l ′ . ( 53
)
Within DFT, imposing such properties amounts to neglecting the effect of the density variation due to the displacement of one atom on the screening of the potential created by the displacement of another atom. Using the notation of Eq. ( 47), the rigid-ion approximation implies that D κα κ ′ γ (q) in RIA is

D RIA κα κ ′ γ (q) = D RIA κα κγ (q)δ κ,κ ′ = D RIA κα κγ (Γ)δ κ,κ ′ (54)
i.e. the non-site-diagonal Debye-Waller contributions vanish.

If we now apply Eq. ( 54) to Eq. ( 51), even if the latter is derived from a DFT Hamiltonian, only two terms remain, which we call the Fan term and the diagonal Debye-Waller term in the rigid-ion approximation

∂ε RIA nk ∂n mq = ∂ε FAN nk ∂n mq + ∂ε DDWRIA nk ∂n mq (55)
with

∂ε FAN nk ∂n mq = 1 2ω mq κα κ ′ γ F κα κ ′ γ (q)U * m,κ ′ γ (q)U m,κα (q) (56) ∂ε DDWRIA nk ∂n mq = -1 4ω mq κα κ ′ γ F κα κ ′ γ (Γ) U * m,κγ (q)U m,κα (q) + U * m,κ ′ γ (q)U m κ ′ α (q) , ( 57
)
where the DDW RIA actually originates from the secondorder contribution of displacing the same atom, although it assumes a Fan-like form. We can see that within RIA all DW-type terms disappear making the calculations easier to perform (only first-order derivative of the Hamiltonian have to be computed). All Fan-like contributions can be derived from DFPT, as explained in the appendix A 3. In practice, in the AHC theory as formulated by Giustino et al 46

, F κα κ ′ γ (q)
is obtained using Eq. ( 48) and

∂u nk ∂R κα (q) ∂ Ĥk,k ∂R κ ′ γ (q) u (0) nk = ∞ ′ n ′ =1 u (0) nk ∂ Ĥk,k ∂Rκα(-q) u (0) n ′ k+q u (0) n ′ k+q ∂ Ĥk,k ∂R κ ′ γ (q) u (0) nk ε (0) nk -ε (0) n ′ k+q , (58)
where the infinite sum over bands is truncated in numerical calculations, while the prime after the sum symbol indicates that the terms with a vanishing denominator (such situation always occurs at Γ when n = n ′ ) have to be excluded. In the resulting expression, one needs to sum over a large number of empty bands since the first-order wavefunction is expressed in the sum-over-states form. As shown by Sternheimer [START_REF] Sternheimer | [END_REF] , the summation over highly energetic bands can be replaced by the solution of a linear equation. This equation can then be solved iteratively with the same techniques as the ones of the DFPT approach used to calculate phonon eigenvectors and eigenenergies 47,48,69,70 . The resulting expression for the first-order wavefunction is detailed in appendix A 3 of this paper (Eq. (A25)) and leads to an alternative form of Eq. ( 58), proposed in Ref. 53

∂u nk ∂R κα (q) ∂ Ĥk,k ∂R κ ′ γ (q) u (0) nk = M ′ n ′ =1 u (0) nk ∂ Ĥk,k ∂Rκα(-q) u (0) n ′ k+q u (0) n ′ k+q ∂ Ĥk,k ∂R κ ′ γ (q) u (0) nk ε (0) nk -ε (0) n ′ k+q + P ck+q ∂u nk ∂R κα (q) ∂ Ĥk,k ∂R κ ′ γ (q) u (0) nk ( 59
)
This formulation is independent of the value of M, which is usually taken to be slightly larger than the number of bands for which the electron-phonon renormalization is sought. It removes the cumbersome sum over states, and results in a significant speed up of the calculations 53 as well as the elimination of the convergence study on the truncation of the sum. The complete expression for the change of eigenenergies due to electronphonon interactions, in the RIA, is obtained from the combination of Eqs. ( 44), ( 48), ( 55) and (59).

D. Beyond the rigid-ion approximation

One can actually analyze the full expression for the derivative of the eigenenergies with respect to phonon occupation numbers, Eq. ( 49), and split it into the sum of the two following contributions

∂ε nk ∂n mq = ∂ε FAN nk ∂n mq + ∂ε DW nk ∂n mq , ( 60
)
where the first term has already been identified and with

∂ε DW nk ∂n mq = ∂ε DDW nk ∂n mq + ∂ε NDDW nk ∂n mq , (61)
where the DW term has been divided into a DDW and a non-diagonal DW contribution (NDDW) term defined as simply the diagonal and non-diagonal in κ, κ ′ parts of the full DW term of Eq. ( 49)

∂ε DDW nk ∂n mq = 1 2ω mq κα κ ′ γ δ κ ′ κ D κα κ ′ γ (q)U * m,κ ′ γ (q)U m,κα (q) (62) ∂ε NDDW nk ∂n mq = 1 2ω mq κα κ ′ γ (1 -δ κκ ′ )D κα κ ′ γ (q)U * m,κ ′ γ (q)U m,κα (q), (63)
where the NDDW term is the mixed derivative of two different atoms in the same cell as well as its replicas.

Finally, the DDW term can be divided even further into a diagonal rigid-ion approximation contribution (DDW RIA ) already identified and a diagonal non-rigidion approximation contribution (DDW NRIA )

∂ε DDW nk ∂n mq = ∂ε DDWRIA nk ∂n mq + ∂ε DDWNRIA nk ∂n mq , (64)
where the DDW NRIA term can be obtained as the re-maining components from Eq. ( 51)

∂ε DDWNRIA nk ∂n mq = 1 2ω mq κα κ ′ γ δ κκ ′ D κα κ ′ γ (q) -D κα κ ′ γ (Γ) 1 2 U * m,κγ (q)U m,κα (q) + U * m,κ ′ γ (q)U m,κ ′ α (q) , (65)
where the DDW NRIA contains only the second-order derivative of the Hamiltonian with respect to two coordinates belonging to different periodic replicas of the same atom. Interestingly, although this term and the DDW RIA originate from the diagonal Debye-Waller term, its expression includes a non-diagonal contribution. We argue that this definition of the diagonal terms is mathematically more convenient as this definitions holds outside of the RIA. An alternative definition, previously used in the literature, defines the DDW as the diagonal term in κ, κ ′ at q = Γ which contain the same derivative taken twice with respect to the same atom. We may note that our FAN and DDW RIA terms are the same as Eq. 16 and Eq. 15 of Ref. 53, co-authored by some of us. Eq. 22 of the latter reference regroups our DDW NRIA and NDDW terms into a term that was called non-diagonal Debye-Waller. Incidentally, this paper studied isolated molecules, for which the only wavevector to be considered is q = Γ. Thus the DDW NRIA contribution always vanished in that case. We believe that the updated definitions are more general since our DDW term is really diagonal in κ, κ ′ , not only in the isolated molecule case, but also in the periodic case, while our NDDW term is purely off-diagonal, as implied by their names and in contrast to the definition of Ref. 53.

IV. TEMPERATURE DEPENDENCE FROM FINITE DIFFERENCES OVER ATOMIC DISPLACEMENTS

The temperature dependence can also be computed through a FP approach where, in a supercell, a set of self-consistent first-principles calculations are done with atoms displaced slightly from their equilibrium positions. The change of force due to atomic displacement allows one to construct the IFC, from which the phonon frequencies and eigenvectors can be deduced.

In a DFT approach (GW behaves similarly, with an electronic self-energy replacing the DFT exchangecorrelation potential), the Hamiltonian is the sum of the kinetic energy operator T and the Kohn-Sham potential VKS [ρ]. The latter can further be split in the sum of potentials generated by each ion V lκ (r -R lκ ) and the Hartree and exchange-correlation (Hxc) potential VHxc [ρ] generated by the electronic density of the system

Ĥ = T + VKS [ρ] = T + lκ V lκ (r -R lκ ) + VHxc [ρ]. (66)
Thus, the change of the Hamiltonian is related to the change of the density only through VHxc [ρ]. The change of the density due to the displacement of one atom being in general affected by the displacement of another atom, the second order derivative of the Hamiltonian with respect to the displacement of two different atoms will contain contributions from VHxc [ρ], but not from T and V lκ (r -R lκ ) since they have the form of Eq. ( 52). This can be seen as a consequence of the fact that, unlike the bare ionic potential, the Hxc potential is screened.

We can use the dielectric function ǫ to describe this effect more rigorously. The change of Kohn-Sham potential can then be expressed as the change of ionic (bare) potential, screened by the inverse dielectric function ǫ-1

∂ VKS ∂R lκ = ǫ-1 ∂ ∂R lκ l ′ κ ′ Vl ′ κ ′ . (67)
As the dielectric function depends on all the atomic positions, frozen-phonon DFT calculations do include nondiagonal terms. Since, in contrast, the AHC formalism (and thus DFPT) neglects non-diagonal contributions through the RIA, it becomes interesting to devise a scheme that allows one to obtain these contributions from FP calculations. One can then assess the validity of the AHC formalism and obtain the magnitude of the contributions neglected by the RIA. We start from Eq. ( 45) and compute the second order derivative of the eigenenergies through a FP approach based on second-order finite differences

∂ε nk ∂n mq = 1 2ω mq lκα l ′ κ ′ γ ∂ 2 ε nk ∂R lκα ∂R l ′ κ ′ γ e -iq•(R l -R l ′ ) U * m,κ ′ γ (q)U m,κα (q) = 1 2ω mq ∂ 2 ∂h 2 ε nk R lκα = R (0) lκα + hU m,κα (q)e -iq•R l h=0 , (68)
where h sets the amplitude of the FP displacement along the normal mode (m, q), R

lκα are the equilibrium position of the atoms, the eigenvalues ε nk are evaluated with the atoms displaced along the phonon mode (m, q) and the l, κ and α indices are iterated upon within the {}.

We also compute, from a FP approach, the first NRIA contribution to the eigenenergies renormalization, that is, the DDW NRIA contribution (Eq. ( 65))

∂ε DDWNRIA nk ∂n mq = 1 2ω mq κ ∂ 2 ∂h 2 κ u (0) nk Ĥk,k R lκ ′ α = R (0) lκ ′ α + h κ ′ U m,κ ′ α (q)e -iq•R l u (0) nk h=0 - κ ∂ 2 ∂h 2 κ u (0) nk Ĥk,k R lκ ′ α = R (0) lκ ′ α + h κ ′ U m,κ ′ α (q) u (0) nk h=0 - γ ∂ 2 ∂s∂t u (0) nk VHxc R lκ ′ α = R (0) lκ ′ α +sU * m,κ ′ γ (q)U m,κ ′ α (q) + tδ αγ u (0) nk h=0 + κ ∂ 2 ∂h 2 κ u (0) nk VHxc R lκ ′ α = R (0) lκ ′ α + h κ ′ U m,κ ′ α (q) u (0) nk h=0 , (69)
where {h κ }, s and t are scalars introduced for finite differences purposes. The first term computed through finite differences is a collective displacement of the atoms in a primitive cell. The second order derivative of this term is non-zero because the curvature of the displaced potential is evaluated at fixed equilibrium ion position through the unperturbed periodic part of the wavefunction u (0)

nk .

We finally proceed to calculate from a FP approach the second NRIA contribution to the renormalization, that is, the NDDW contribution

∂ε NDDW nk ∂n mq = 1 2ω mq κκ ′ (κ =κ ′ ) ∂ 2 ∂h κ ∂h κ ′ u (0) nk VHxc R lκ ′′ α = R (0) lκ ′′ α + h κ ′′ U m,κ ′′ α (q)e -iq•R l u (0) nk h=0
.

Antonius et al 39 computed the temperature dependence of diamond using finite differences (Eq. ( 68)) in the many-body GW framework, which led to an additional 200 meV ZPR of the diamond bandgap with respect to the AHC value, closer to the experimental bandgap. In this section, we compare AHC and and FP results for some chosen q-wavevector contributions to the ZPR of the diamond bandgap. We study the q-wavevector contributions to the ZPR instead of the full ZPR since the converged q-point integration required by the full ZPR is computationally out of reach for FP (it requires 70x70x70 q-point grids and associated supercells). For the full qpoint integration within AHC, please see Ref. 54.

The calculation of structural properties in this work are based on Density Functional Theory (DFT) [START_REF] Martin | Electronic Structure. Basic Theory and Practical Methods[END_REF] within the local density approximation (LDA) 71,72 . We use a planewave basis set to represent the wavefunctions and account for the core-valence interaction using normconserving pseudopotentials 73 . The valence electrons of carbon treated explicitly in our ab-initio calculations are 2s 2 2p 2 . All the calculations are done using the ABINIT software package 74 .

Convergences studies with a tolerance of 0.5mHa per atom on the total energy led to the use of a 6x6x6 Γ-centered Monkhorst-Pack sampling 75 of the Brillouin zone and an energy cut-off of 30 Hartree for the truncation of the planewave basis set. The lattice parameter of 6.675 Bohr was obtained by structural optimization.

The phonon frequencies were calculated using Eq. (5) for the DFPT method and a second-order derivative of the total energy with respect to atomic displacements converged with a Richardson interpolation of order 4 for the finite-difference method. Also, for the latter method and the q = L point, a 2x2x2 supercell with a 3x3x3 Γ-centered Monkhorst-Pack sampling were used to ensure the k-point sampling remained equivalent to the one used for the DFPT calculation on the primitive cell. The comparison between FP and DFPT for the phonon frequencies is given in Table I. The discrepancies between the two methods remain within 0.2% for the two q-points considered, which demonstrates both the equivalence of the two methods and the convergence of our calculations.

We now assess the accuracy of our AHC implementation and quantify the impact of the RIA. The contribution of the same two q-points to the ZPR at k = Γ and k = L is given in Table II. This table presents the renormalization of the different eigenenergies due to electronphonon coupling for the 4 first distinguishable bands of diamonds. The Γ 25 ′ valence band maximum is threefold degenerate as well as the Γ 15 conduction band. In diamond, the conduction band minimum is located between the k = Γ and k = X points. At the k = L, the valence L 1 and conduction L 3 ′ bands are only doubly degenerate. It can also be noted that our results show the contribution of some q-points to the ZPR with values in very close agreements (within 2 meV) with Figure 1 of Ref. 39 that show the electron-phonon coupling energies (there is therefore a conversion factor of 1/2).

The AHC results are split into Fan and DDW RIA contributions that are computed using Eq. ( 56) and Eq. ( 57), respectively.

The importance of the rigid-ion approximation (RIA) can be deduced by computing the DDW NRIA and NDDW terms through finite difference calculations using Eq. ( 69) and Eq. ( 70), respectively. It can be noticed that the DDW NRIA and the NDDW terms are equal at q = Γ, due to the T 2g optical mode of the diamond crystal for which U m,κα = -U m,κ ′ α when κ = κ ′ , as expected from Eqs. ( 65) and (63).

The impact of the RIA remains below 23% and is usually much smaller for wavevectors other than the zonecenter one. To assess completely the validity of the approximation, we should do a full q-point integration on the BZ. However, since the ZPR converges extremely slowly with the number of q-points 54 , this would require huge supercell calculations, which are currently computationally out of reach.

Finally, the sum of the AHC and NRIA contributions can be compared to the finite difference calculations done using Eq. (68). All the FP calculations are also converged with a Richardson interpolation of order 4. The discrepancies remain below 7 µeV in absolute value. Such discrepancies can be attributed to numerical noise or anharmonicity. Indeed, the finite displacements selected in the FP method ensure a good compromise between these two sources of error, so that both are present, but as small as allowed by our convergence criteria.

VI. CONCLUSIONS

The renormalization of eigenenergies due to electronphonon coupling can be computed by different methods. In this paper, we have reviewed three of them: the first-principle molecular dynamics method, the frozenphonon (FP) method and the Allen-Heine-Cardona (AHC) method based on density-functional perturbation theory. The two first methods are equivalent within the adiabatic harmonic approximation while the third is only equivalent when the rigid-ion approximation (RIA) is also performed.

The theory's key ingredient is the second-order derivative of the eigenenergies with respect to two atomic displacements. This derivative gives rise to a term stemming from the first-order modification of the wavefunction called the Fan term and a term corresponding to a second-order perturbation of the Hamiltonian called the Debye-Waller (DW) term. Although the two terms were discovered separately in the 50's, there was a lot of confusion in the litterature until Baumann realised in 1974 the complementarity of these two terms for the computation of the zero-point motion renormalization (ZPR). The present paper compared in detail two (AHC/DFPT and FP) of the three approaches still used today to calculate the ZPR. We considered the first for its efficiency in the computation of the ZPR at arbitrary q-points, crucial for periodic system, and the later, to go beyond the RIA and study its impact on the calculated ZPR.

Also, in this paper, we derived Brook's theorem in the first-principle context and obtained an expression for the eigenenergy renormalization from the finite temperature extension of Janak's theorem (see Eqs. (31) and ( 24)).

We also rederived how, within the RIA, the translational invariance (Eq. ( 42)) allows to express the DW contribution in terms of first-order derivatives of the Hamiltonian only.

A major contribution of this paper is the clarification of the terms appearing beyond the RIA made in the AHC theory. The DW term is divided into a diagonal rigid-ion approximation contribution (DDW RIA ), a diagonal nonrigid-ion approximation contribution (DDW NRIA ), and a non-diagonal DW contribution (NDDW). This allows all the term's definition to be coherent with their respective names (see Eqs. ( 56), ( 57), ( 65) and ( 63)).

Nevertheless, due to the computational limitation related to the evaluation of the second-order derivative of the Hamiltonian, the NRIA terms cannot be computed with DFPT. Therefore we derived the equations related to the FP approach, which allows us to numerically evaluate the NRIA terms in the case of bulk diamond (see Eqs. ( 68), ( 69) and ( 70)).

For the diamond phonon frequencies, the discrepancy between the FP and AHC approaches is below 0.2%. The differences of ZPR between the two methods are always below 7µeV in absolute value, which strengthens our confidence in the theory and numerical implementation within the ABINIT software. These small differences are attributed to the unavoidable numerical noise and anharmonicity in the FP calculations. The impact of the RIA is also evaluated for two q-points (Γ and L) and is found to be as large as 23% for the Γ zone center contribution to the ZPR but is usually much smaller for other q-points contributions.

al 78 ). 80 Although we agree with Eqs. (2-7) of Allen and Cardona 23 , we think that their Eq. ( 8) is erroneous, since the polarization vectors in the Debye-Waller term of this equation have different atomic indices κ and κ ′ , which is not coherent with their Eq. ( 5) or with Ref. 46, where the polarization vectors have the same atomic index κ. On the contrary, we perfectly agree with the formulation of Giustino et al. 46 .

TABLE I :

 I Comparison of phonon frequencies computed from FP and from DFPT.

	q-point	Mode	ω DFPT [Ha] ω FP [Ha] diff. [%]
	Γ	LO+TO	0.00606474	0.00605226	0.2058
	L	TA	0.00250256	0.00250111	0.0580
		LA	0.00494158	0.00494344	0.0376
		TO	0.00564786	0.00565244	0.0810
		LO	0.00577597	0.00577016	0.1007
	V. COMPARISON BETWEEN AHC/DFPT AND
		FINITE-DIFFERENCE RESULTS

TABLE II :

 II Comparison between the AHC/DFPT and FP results for the contributions of specific q-points to the ZPR of several bands of diamond. The column "diff." is the sum of the AHC and the NRIA contributions minus the FP result.

	q-point k-point		AHC [meV]		NRIA [meV]		sum all	FP	diff.
			Fan	DDWRIA sumAHC DDWNRIA NDDW % NRIA	[meV]	[meV]	[µeV]
	Γ	Γ1	-32.8174	20.2868 -12.5306 0.7250	0.7250 13.09	-11.0806 -11.0809	0.308
		Γ 25 ′	-332.4265 357.2564	24.8300 1.7991	1.7991 12.66	28.4282	28.4289 -0.638
		Γ15	-330.4363 316.2021 -14.2342 0.1924	0.1924	2.78	-13.8494 -13.8497	0.338
		Γ 2 ′	-63.3087	32.3760 -30.9327 0.1499	0.1499	0.98	-30.6328 -30.6335	0.748
		L 2 ′	-67.8598	46.8783 -20.9814 1.1410	1.1410 12.20	-18.6993 -18.6999	0.554
		L1	-146.0806 129.4769 -16.6037 0.5663	0.5663	7.32	-15.4710 -15.4714	0.367
		L 3 ′	-311.2306 321.3287	10.0981 1.4805	1.4805 22.67	13.0590	13.0592 -0.307
		L3	-473.0115 292.4656 -180.5458 0.0779	0.0779	0.09	-180.3900 -180.3937	6.977
	L	Γ1	-116.4278	62.6966 -53.7312 3.2318	-2.3250	1.72	-52.8256 -52.8245	1.104
		Γ 25 ′	-922.8240 1104.1052	181.2812 6.0491	-3.7542	1.25	183.5761	183.5771 -1.017
		Γ15 -1250.8082 977.2263 -273.5819 1.0840	-2.0900	0.37	-274.5878 -274.5881	0.244
		Γ 2 ′	-407.6022 100.0584 -307.5438 0.3437	-2.1920	0.60	-309.3921 -309.3973	5.131
		L 2 ′	-234.2353 144.8781 -89.3572 4.3856	-2.8996	1.69	-87.8712 -87.8728	1.542
		L1	-620.7070 400.1500 -220.5570 2.6651	-2.1582	0.23	-220.0501 -220.0525	2.359
		L 3 ′ -1018.9788 993.0698 -25.9090 5.1401	-3.3984	7.21	-24.1674 -24.1672 -0.210
		L3	-740.6821 903.8683	163.1862 0.7991	-1.8919	0.67	162.0934	162.0935 -0.142

VII. ACKNOWLEDGEMENTS

This work was supported by the FRS-FNRS through a FRIA grant (S.P.) and a FNRS grant (Y.G.). Moreover, A. M. would like to acknowledge financial support from the Futuro in Ricerca grant No. RBFR12SW0J of the Italian Ministry of Education, University and Research. The authors would like to thank Yann Pouillon and Jean-Michel Beuken for their valuable technical support and help with the test and build systems of ABINIT. Computational ressources have been provided by the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CECI) funded by the Fonds de la Recherche Scientifique de Belgique (FRS-FNRS) under Grant No. 2.5020.11.

Appendix A: Technicalities

Integrals of phonon wavefunctions and powers of the position operator

We evaluate the integrals present in Eq. (7), with the Taylor expansion of the eigenenergy.

In the harmonic approximation, the phonon wavefunctions are obtained by solving the Schrödinger equation for the harmonic oscillator. We obtain [START_REF] Griffiths | Introduction to Quantum Mechanics[END_REF] χ sm,m (z) = ω m π

where H sm (ξ) = (-1) sm e ξ 2 d sm dξ sm e -ξ 2 is the Hermite polynomial and ξ = √ ω m z is a dimensionless position variable. Hermite polynomials satisfy the following orthonormality condition

The phonon wavefunctions are thus normalized, with

The first-order integral cancels out, as an odd function integrates to zero.

Finally, the square of the dimensionless position operator in second quantization can be expressed as

Hence, the second-order in the Taylor expansion yields

2. Convention for the unperturbed periodic system

Following the same convention as Gonze 48 , the unperturbed wavefunction can be obtained as the product of a phase factor and a periodic function (Bloch's theorem)

where N BvK is the number of unit cells repeated in the Born-von Karman periodic box, Ω 0 volume of the unperturbed unit cell, n the band index and k label the wave vector of the wavefunction. The periodic part of the Bloch wavefunction, in Eq. (A8), is subject to the following orthonormalization condition

where the scalar product of periodic functions is defined as

For a generic operator we follow the following convention

(A11)

Perturbation theory

The perturbation theory for a periodic system is based on the idea that the solution of a reference system (usually the equilibrium ground state one-body Schrödinger equation) is known. In general, the perturbation can be incommensurate with the periodic system and characterized by a wave vector q. If the amplitude of the perturbation is characterized by a small scalar parameter λ, then any observable X(λ) can be expressed as a power serie X(λ) = X (0) + (λX (1) q + λ * X

(1)

-q,-q ) + • • • , (A12) were we use the superscript notation as a shorthand for derivatives:

The perturbed Schrödinger equation now depends explicitely on the parameter λ

where n is the band index, k is a wavevector in the Brillouin zone. The solutions |Ψ nk (λ) must fulfill the normalization condition Ψ nk (λ)|Ψ nk (λ) = 1. The Hamiltonian Ĥ(λ) depends parametrically on the atomic position R lκ of the atom κ in the cell l.

The translated first-order wave function becomes

where R a is a vector of the real space lattice. We can then factorize out the phase factor to map the incommensurate problem into a problem commensirate with periodicity of the unperturbed one. To this end, we introduce the periodic first-order wave functions

nk,q (r). (A15)

We define S Mk+q as the space of the M low-lying k + q states, M being larger or equal to the index of the highest band for which we aim to compute the temperaturedependent behaviour. We define P Mk+q as the projector on S Mk+q . We also define P ck+q as the projector on the subspace of k + q ground-state wavefunctions complementary to S Mk+q . The first-order derivatives of wavefunctions can be split in two contributions, one that is contained inside S Mk+q and one that belongs to its complimentary, u

nk,q = P Mk+q u

(1)

nk,q + P ck+q u

(1) nk,q .

(A16)

can be easily computed using standard perturbation theory,

nk,q , we want to avoid the summation over an infinite number of states. So, in the spirit of DFPT, we minimize

where E

(2) -q,q {u (0) , u (1) } is given in Eq. 42 of Ref. 48, under the constraint of parallel-transport gauge [START_REF] Gonze | [END_REF] ,

with the Lagrange parameters

The minimum of this expression with respect to variations of P ck+q u

(1) nk,q leads to the canonical Euler-Lagrange equation

We can then pre-multiply on each side by u

where n ′′ ∈ [1, M ] and get

This leads to the equation

We can then substitute Eq. (A23) inside Eq. (A21) and get

and finally

Supposing H (1) has already be determined by the usual DFPT self-consistency loop over occupied states only (in which M = N val ), the u (1) nk,q are found by solving Eq. (A25), combined with Eqs. (A16) and (A17). * Electronic address: samuel.ponce@uclouvain.be 1 W. G. Aulbur, L. Jönsson, and J. W.