Ali El Boutaybi
email: ali.el-boutaybi@c2n.upsaclay.fr

Panagiotis Karamanis

Thomas Maroutian

Sylvia Matzen

Laurent Vivien

Philippe Lecoeur

Michel Rérat

Electro-optic properties of ZrO 2 , HfO 2 and LiNbO 3 ferroelectric phases: A comprehensive and comparative study with density functional theory

We report the Pockels electro-optic properties of ZrO 2 and HfO 2 orthorhombic Pbc2 1 and rhombohedral R3m ferroelectric phases, and we compare them to the well-known LiNbO 3 Pockels material from density functional theory calculations using the CRYSTAL suite of program. Specifically, three essential processes are explicitly investigated: The electronic, the ionic (or vibrational), and the piezoelectric contributions. Our calculations reveal that the ionic part coming from the lowfrequency phonon modes contributes the most to the electro-optic coefficients of rhombohedral LiNbO 3 and of orthorhombic ZrO 2 and HfO 2 . While these low-frequency modes show zero contribution to Pockels coefficients for the rhombohedral phase of the latter compounds.

Introduction

Silicon photonics platforms are becoming indispensable for developing low-cost integrated photonic circuits. Electro-optic (EO) effect is one of the main properties required in plethora applications including optical communications [START_REF] Mohapatra | A giant electrooptic effect using polarizable dark states[END_REF][START_REF] Melikyan | High-speed plasmonic phase modulators[END_REF][START_REF] Rueda | Resonant electro-optic frequency comb[END_REF][START_REF] Abel | Large Pockels effect in micro-and nanostructured barium titanate integrated on silicon[END_REF], quantum computing [START_REF] O'brien | Photonic quantum technologies[END_REF][START_REF] Wang | Integrated photonic quantum technologies[END_REF], and neuromorphic applications [START_REF] George | Neuromorphic photonics with electro-absorption modulators[END_REF][START_REF] Shen | Deep learning with coherent nanophotonic circuits[END_REF][START_REF] Offrein | Prospects for photonic implementations of neuromorphic devices and systems[END_REF]. The principal physical process behind these applications is the so-called Pockels or linear EO effect, which is related to the modulation of the refractive index of a material 1 arXiv:2209.12536v1 [cond-mat.mtrl-sci] 26 Sep 2022 under an applied external electric field [START_REF] Lines | Principles and Applications of Ferroelectrics and Related Materials[END_REF]. Nevertheless, due to its centrosymmetric atomistic crystalline structure, pristine silicon features vanishing second-order optical nonlinearities; thus, it is unsuitable for EO applications. Bearing in mind that second-order nonlinear optical (NLO) responses strictly manifest in non-centrosymmetric materials, ferroelectric (FE) crystalline solids have emerged as the main materials of choice in EO applications. Such materials exhibit an intrinsic electric polarization that, in addition, can be controlled and tuned in a reversible manner by the application of an external electric field.

In the realm of EO applications, the most studied FE materials are currently ABO 3 perovskitetype crystalline solids. A representative member of this family is lithium niobate (LiNbO 3 , LNO), a ferroelectric material featuring a strong Pockels crystalline bulk coefficient (r 33 ) of about 30 pm/V [START_REF] Turner | High-frequency electro-optic coefficients of lithium niobate[END_REF]. LNO is widely applied in the telecommunications industry [START_REF] Wooten | A review of lithium niobate modulators for fiber-optic communications systems[END_REF] and specifically in the fabrication of optical modulators. Despite LNO's complexity of integration on silicon substrates, several studies showed its compatibility with silicon photonics [START_REF] Wang | Integrated lithium niobate electro-optic modulators operating at CMOScompatible voltages[END_REF][START_REF] Payam | Heterogeneous lithium niobate photonics on silicon substrates[END_REF]. An additional appealing perovskitetype crystalline solid for EO applications is tetragonal ferroelectric barium titanate (BaTiO 3 , BTO), featuring one of the highest Pockels bulk coefficients (r 42 ≈ 1300 pm/V) [START_REF] Zgonik | Dielectric, elastic, piezoelectric, electro-optic, and elasto-optic tensors of BaTiO 3 crystals[END_REF]. However, the growth of high-quality BTO thin films on silicon without affecting its strong Pockels coefficient is still challenging and of high cost, although high-quality BTO films can be fabricated via molecular beam epitaxy (MBE), [START_REF] Abel | Large Pockels effect in micro-and nanostructured barium titanate integrated on silicon[END_REF][START_REF] Messner | Plasmonic Ferroelectric Modulators[END_REF][START_REF] Xiong | Active Silicon Integrated Nanophotonics: Ferroelectric BaTiO 3 Devices[END_REF]. In addition to tetragonal BTO, its rhombohedral R3m phase also gained attention for cryogenic technologies as, for instance, for quantum computing [START_REF] Paoletta | Pockels effect in low-temperature rhombohedral BaTiO 3[END_REF][START_REF] Eltes | An integrated optical modulator operating at cryogenic temperatures[END_REF].

Recently, it has been demonstrated that among HfO 2 -ZrO 2 -based thin films, orthorhombic Pbc2 1 [START_REF] Böscke | Ferroelectricity in hafnium oxide thin films[END_REF] and rhombohedral R3m [START_REF] Wei | A rhombohedral ferroelectric phase in epitaxially strained Hf 0.5 Zr 0.5 O 2 thin films[END_REF] phases exhibit ferroelectric behavior. This is a significant outcome in terms of practical applications because both HfO 2 and ZrO 2 have an excellent silicon compatibility [START_REF] Ali | High Endurance Ferroelectric Hafnium Oxide-Based FeFET Memory Without Retention Penalty[END_REF][START_REF] Müller | Ferroelectric Zr 0.5 Hf 0.5 O 2 thin films for nonvolatile memory applications[END_REF][START_REF] Park | Thin Hf x Zr 1-x O 2 Films: A New Lead-Free System for Electrostatic Supercapacitors with Large Energy Storage Density and Robust Thermal Stability[END_REF][START_REF] Silva | Advances in dielectric thin films for energy storage applications, revealing the promise of group iv binary oxides[END_REF]. Combined with their ferroelectric nature, this might lead to significant breakthroughs in the realm of integrated photonic circuits that could significantly reduce the fabrication cost. Nevertheless, there is still much work to be done since the film quality and the ferroelectricity of ZrO 2 and HfO 2 are not well optimized at high film thicknesses. This is critical for their optimal application due to the inverse proportionality between the optical losses and the film thickness [START_REF] Fork | Optical losses in ferroelectric oxide thin films: Is there light at the end of the tunnel[END_REF][START_REF] Sando | Epitaxial ferroelectric oxide thin films for optical applications[END_REF]. In the case of the orthorhombic Pbc2 1 (o-phase), which is usually reported in polycrystalline thin films, ferroelectricity has been demonstrated in a large range of thickness from 1 nm to 1µm [START_REF] Suraj S Cheema | Enhanced ferroelectricity in ultrathin films grown directly on silicon[END_REF][START_REF] Polakowski | Ferroelectricity in undoped hafnium oxide[END_REF][START_REF] Starschich | Ferroelectric and piezoelectric properties of Hf 1-x Zr x O 2 and pure ZrO 2 films[END_REF][START_REF] Shimura | Preparation of 1 µm thick Y-doped HfO 2 ferroelectric films on (111) Pt/TiO x /SiO 2 /(001)Si substrates by a sputtering method and their ferroelectric and piezoelectric properties[END_REF]. This o-phase has also been reported in epitaxial thin films with a thickness thinner than 20 nm [START_REF] Lyu | Robust ferroelectricity in epitaxial Hf 1/2 Zr 1/2 O 2 thin films[END_REF][START_REF] Jiao | Ferroelectric Hf 0.5 Zr 0.5 O 2 thin films deposited epitaxially on (110)-oriented SrTiO 3[END_REF][START_REF] Yun | Intrinsic ferroelectricity in Y-doped HfO 2 thin films[END_REF]. Regarding the rhombohedral R3m (r-phase), the reported thicknesses are ranged from 5 to 40 nm [START_REF] Wei | A rhombohedral ferroelectric phase in epitaxially strained Hf 0.5 Zr 0.5 O 2 thin films[END_REF][START_REF] Silva | Wake-up Free Ferroelectric Rhombohedral Phase in Epitaxially Strained ZrO 2 Thin Films[END_REF][START_REF] Boutaybi | Stabilization of the epitaxial rhombohedral ferroelectric phase in ZrO 2 by surface energy[END_REF]. Indeed, for applications, ultrathin films are often advantageous; for example, in the case of nonvolatile memories and ferroelectric field-effect transistors [START_REF] Ali | High Endurance Ferroelectric Hafnium Oxide-Based FeFET Memory Without Retention Penalty[END_REF][START_REF] Müller | Ferroelectric Zr 0.5 Hf 0.5 O 2 thin films for nonvolatile memory applications[END_REF], energy storage [START_REF] Park | Thin Hf x Zr 1-x O 2 Films: A New Lead-Free System for Electrostatic Supercapacitors with Large Energy Storage Density and Robust Thermal Stability[END_REF][START_REF] Silva | Advances in dielectric thin films for energy storage applications, revealing the promise of group iv binary oxides[END_REF], negative capacitance [START_REF] Hoffmann | Direct Observation of Negative Capacitance in Polycrystalline Ferroelectric HfO 2[END_REF][START_REF] Hoffmann | Negative capacitance in hfo2-and zro2-based ferroelectrics[END_REF], and tunnel junctions [START_REF] Wei | Magnetic Tunnel Junctions Based on Ferroelectric Hf 0.5 Zr 0.5 O 2 Tunnel Barriers[END_REF][START_REF] Goh | The effect of the bottom electrode on ferroelectric tunnel junctions based on CMOS-compatible HfO 2[END_REF].

In terms of NLO properties, the studies on ZrO 2 -HfO 2 materials are scarce, leaving open the question of the potential of these compounds for optical applications. Indeed, to our knowledge, only two experimental studies can be found in the literature [START_REF] Kondo | Linear electro-optic effect in ferroelectric HfO 2 -based epitaxial thin films[END_REF][START_REF] Kondo | Influence of orientation on the electro-optic effect in epitaxial Y-doped HfO 2 ferroelectric thin films[END_REF]. They investigate the EO effect in Y 2 O 3 -doped HfO 2 , reporting an EO coefficient of about 0.67 pm/V for (111)-oriented thin-films, a modest value compared to the ones in LNO and BTO materials. However, the exact EO coefficients of the two experimentally reported ferroelectric phases in HfO 2 and ZrO 2 are not known yet. Therefore, this study's primary objective is to report and investigate the EO coefficients of HfO 2 and ZrO 2 materials.

In this work, we study the second-order nonlinear optical (NLO) properties of the orthorhombic Pbc2 1 and rhombohedral R3m phases of ZrO 2 and HfO 2 , relying on all-electron density functional theory (DFT) computations to report reliable EO coefficients. The obtained results are compared to the ones computed for the rhombohedral R3c LiNbO 3 (LNO) phase, for which experimental data are available allowing us to use LNO as a reference material. All reported EO coefficients have been determined from the microscopic second-order NLO responses of these systems that, in turn, they have been obtained via a coupled perturbed Kohn-Sham (CPKS) analytical approach developed and implemented in the CRYSTAL17 suite of quantum chemical programs by one of the authors of this study [START_REF] Dovesi | CRYS-TAL14: A program for the ab initio investigation of crystalline solids[END_REF][START_REF] Dovesi | Quantum-mechanical condensed matter simulations with CRYSTAL[END_REF]. It has already been established that this perturbative approach delivers high accuracy molecular and bulk electric properties, free from numerical errors [START_REF] Maschio | Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis[END_REF][START_REF] Maschio | Ab initio analytical infrared intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method[END_REF][START_REF] Maschio | Comment on "Ab initio analytical infrared intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method[END_REF][START_REF] Maschio | Calculation of the dynamic first electronic hyperpolarizability β(-ω σ ; ω 1 , ω 2 ) of periodic systems. Theory, validation, and application to multi-layer MoS 2[END_REF][START_REF] Rérat | Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab Initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme[END_REF]. To achieve the maximum possible relevance with the available or future experimental data, the reported coefficients comprise electronic, vibrational effects, and piezoelectric contributions, which have also been obtained analytically. Especially for comparison between theory and experiment, such contributions are of essential importance for a reliable prediction of the overall NLO response of a given crystalline or molecular material and the complete understanding of its origin. Our article is organized as follows: First, we give the theoretical background of our study, insisting on the close relationship between the macroscopic NLO susceptibilities and the microscopic hyper-polarizabilities of the unit shells. Then, we study the EO properties of LNO using CRYSTAL code, which will serve as the material of reference. Specifically, we report reliable values of its microscopic second-order responses and analyze the electronic, vibrational, and piezoelectric contributions. Finally, we report and discuss the Pockels EO coefficients of o-and r-phase ZrO 2 and HfO 2 , comparing them to those calculated for the R3c rhombohedral phase of LNO.

Theoretical framework and computational details

The microscopic linear and non-linear optical properties of a finite zero-dimensional system are defined in terms of a Taylor power series expansion with respect to an external applied electric field E j kl as

µ i = µ 0 + α tot i j E j + 1 2 β tot i j k E j E k + 1 6 γ tot i j kl E j E k E l + ... (1)
Where, i, j, k and l stand for Cartesian coordinates (x, y, z) , µ i and µ 0 are the induced and the permanent dipole moment, α tot i j is the total dipole electric polarizability, while β tot i j k and γ tot i j kl are the total first and second dipole hyperpolarizabilities, respectively. By the term total (tot) it is implied that both electronic (ele) and vibrational (vib) contributions [START_REF] Maschio | Calculation of the dynamic first electronic hyperpolarizability β(-ω σ ; ω 1 , ω 2 ) of periodic systems. Theory, validation, and application to multi-layer MoS 2[END_REF][START_REF] Rérat | Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab Initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme[END_REF][START_REF] Bishop | Molecular vibrational and rotational motion in static and dynamic electric fields[END_REF] are taken into account, Hence:

α tot i j = α el e i j + α vi b i j (2)
β tot i j k = β el e i j k + β vi b i j k (3)
For three-dimensional periodic crystalline lattices, equation (1) can be written as

P i = P 0 + χ (1) i j E j + χ (2) i j k E j E k + χ (3) i j kl E j E k E l + ... ( 4
)
Where P i is the unit-cell induced polarization while χ (1) i j , χ (2) i j k , and χ (3) i j kl stand for the total first order (linear), and the second-and the third-order (nonlinear) optical macroscopic bulk susceptibilities, respectively. The first and second coefficients of equation ( 4) are deduced from the total microscopic (hyper) polarizabilities per unit-cell volume (V ) as follows (in atomic units):

χ (1) i j = 4π V α tot i j
(5)

χ (2) i j k = 2π V β t ot i j k (6)
The total first-order susceptibility χ (1) i j can be used to deduce the relative dielectric tensor i j , which in turn gives access to bulk refractive indices n via the following expressions: i j = δ i j + χ (1) i j [START_REF] George | Neuromorphic photonics with electro-absorption modulators[END_REF]

n i i = 1/2 i i (8)
In the above equations, δ i j represents the elements of the identity matrix (δ i j = 1 for i = j , and 0 for i = j ). Note that we consider in this work a diagonal i j tensor so that the subscript i of i i corresponds to one of the three principal axes n i i of the ellipsoidal indicatrix. Now, to obtain the EO Pockels coefficients (r i j k ), as a first step, we consider the constant stress (or "clamped") conditions, wherein χ (2) i j k and the bulk refractive indices n i i given in equation ( 8) can be used to determine one out of the two terms contributing to the total r i j k coefficients. The clamped EO coefficient is denoted here as r S i j k , comprises both vibrational and electronic contributions, and is defined as:

r S i j k = -2 χ (2) i j k n 2 i i n 2 j j (9)
In the second step, we should model and consider the relevant acoustic contribution stemming from the piezoelectric effect. This effect depends on the frequency of the applied ac electric field, as it originates from the inverse piezoelectric effect, according to which the application of an ac oscillating field induces mechanical deformations on the crystalline lattice. The resulting contribution to the r i j k coefficients, denoted here as r p i j k , is defined from the photoelastic tensor P i j µν and converse piezoelectric tensor d kµν [START_REF] Erba | Piezo-optic tensor of crystals from quantum-mechanical calculations[END_REF][START_REF] Jazbinšek | Material tensor parameters of LiNbO 3 relevant for electro-and elasto-optics[END_REF] as: [START_REF] Lines | Principles and Applications of Ferroelectrics and Related Materials[END_REF] The calculation of the tensorial components of r p i j k is achieved via the computation of the photoelastic constants ∂η i j ∂S µν (P i j µν ) defined as the derivative of the impermeability tensor η i j with respect to strain S µν , which are the elements of the fourth rank photoelastic tensor, and of the converse piezoelectric tensor ∂S µν ∂E k (d kµν ) [START_REF] Erba | Photoelasticity of crystals from theoretical simulations[END_REF]. Finally, the total EO coefficients that we report in this work are calculated as r i j k = r S i j k + r p i j k . More detailed information about photoelastic calculations can be found elsewhere [START_REF] Erba | Piezo-optic tensor of crystals from quantum-mechanical calculations[END_REF][START_REF] Jazbinšek | Material tensor parameters of LiNbO 3 relevant for electro-and elasto-optics[END_REF][START_REF] Erba | Photoelasticity of crystals from theoretical simulations[END_REF].

r p i j k = 3 µ,ν=1 P i j µν d kµν
All computations of frequency-dependent properties have been performed with the CRYSTAL suite of programs [START_REF] Dovesi | Quantum-mechanical condensed matter simulations with CRYSTAL[END_REF]. In order to determine the dynamic vibrational and electronic (clamped nuclei) hyper-polarizabilities, a fully analytical approach relying on the coupled perturbed Hartree-Fock and Kohn-Sham method (CPHF / CPKS) has been applied [START_REF] Maschio | Comment on "Ab initio analytical infrared intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method[END_REF][START_REF] Maschio | Calculation of the dynamic first electronic hyperpolarizability β(-ω σ ; ω 1 , ω 2 ) of periodic systems. Theory, validation, and application to multi-layer MoS 2[END_REF][START_REF] Rérat | Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab Initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme[END_REF][START_REF] Pascale | The calculation of the vibrational frequencies of crystalline compounds and its implementation in the crystal code[END_REF][START_REF] Zicovich-Wilson | Calculation of the vibration frequencies of α-quartz: The effect of Hamiltonian and basis set[END_REF]. In the case of LNO, we considered three different DFT functionals to investigate the dependence on the DFT method of the properties of interest and to compare the results with available experimental data. That includes the Becke, three-parameter, Lee-Yang-Parr exchange-correlation functional (B3LYP) [START_REF] Becke | Density-functional thermochemistry. III. The role of exact exchange[END_REF][START_REF] Lee | Development of the colle-salvetti correlation-energy formula into a functional of the electron density[END_REF], the pure Perdew-Burke-Ernzerhof one (PBE) [START_REF] Perdew | Generalized gradient approximation made simple[END_REF], and with 25% Hartree-Fock (HF) exchange (PBE0) [START_REF] Adamo | Toward reliable density functional methods without adjustable parameters: The PBE0 model[END_REF]. In the case of ZrO 2 and HfO 2 (orthorhombic Pbc2 1 and rhombohedral R3m phases), only the B3LYP functional has been applied as it gives good results for LNO. The basis sets of LNO are taken from Refs. [START_REF] Dall'olio | Spontaneous polarization as a berry phase of the hartreefock wave function: The case of KNbO 3[END_REF], [START_REF] Vilela | BSSE-correction scheme for consistent gaussian basis sets of double-and triple-zeta valence with polarization quality for solidstate calculations[END_REF], and [START_REF] Mccall | Anisotropic optical properties of an oriented-emeraldine-base polymer and an emeraldine-hydrochloride-salt polymer[END_REF] for Nb, Li, and O, respectively. The basis sets for ZrO 2 are from Refs. [START_REF] Vilela | BSSE-correction scheme for consistent gaussian basis sets of double-and triple-zeta valence with polarization quality for solidstate calculations[END_REF] and [START_REF] Valenzano | Disclosing the Complex Structure of UiO-66 Metal-Organic Framework: A Synergic Combination of Experiment and Theory[END_REF] for O and Zr, respectively, while those for HfO 2 are taken from Refs. [START_REF] Bredow | Electronic and magnetic structure of ScMnO 3[END_REF] and [START_REF] Munoz Ramo | Spectroscopic properties of oxygen vacancies in monoclinic HfO 2 calculated with periodic and embedded cluster density functional theory[END_REF] for O and Hf, respectively.

Results and discussion

Nb

LNO rhombohedral R3c

We started our calculations with LNO, which features a ferroelectric R3c ground state at low temperatures. The crystal structure of the bulk system, built from ten atoms per unit cell, is shown in Figure 1a (lattice constants and atomic coordinates are given in Table S.1 in Supplemental Material (SM) [66]). Computed band gaps, dynamic refractive indices n i i and i i values of LNO determined with three functionals of different types, namely B3LYP, PBE0, and PBE at an optical frequency corresponding to 633 nm wavelength, are listed in Table 1. Apart from electronic contributions, which can be trivially computed with the CPKS approach, the reported values comprise as well the respective vibrational contributions computed at the same level of theory. The wavelength (633 nm) was chosen to compare our results to available data reported in the literature; both experimental and theoretical data are reported in Table 1. Note that the EO properties are also investigated up to the near-infrared wavelength range in section 3.2.

Starting from the computed bandgap of LNO, it is seen that out of the three methods used, B3LYP and PBE0 yield notably larger values as compared to the reported experimental values. On the other hand, the pure PBE functional follows the opposite trend. The frequently cited direct band gap value of LNO, mostly concluded from optical experiments [START_REF] Dhar | Optical properties of reduced lithium niobate single crystals[END_REF], is about 3.78 eV. However, this value should be considered with care due to electron-hole attraction effects, which might induce important underestimations in the final outcome of optical measurement, as already discussed in Refs. [START_REF] Thierfelder | Do we know the band gap of lithium niobate?[END_REF][START_REF] Nahm | First-principles study of microscopic properties of the Nb antisite in LiNbO 3 : Comparison to phenomenological polaron theory[END_REF]. Therefore, for LNO crystal, a band gap of 4.7 eV [START_REF] Thierfelder | Do we know the band gap of lithium niobate?[END_REF] should be considered more suitable for comparisons between theory and experiment. Bearing in mind the latter correction, the B3LYP functional and the associated basis set reproduce the band gap of LNO better than PBE.

Turning our attention to the refractive indices computed in this study, we also notice that both B3LYP and PBE0 functionals yield values that fall very close to the available experimental measurements. On the other hand, pure PBE overshoots the experimental reference values, yielding values close to the earlier Local-Density approximations (LDA) reported results. This behavior is not surprising since pure DFT functionals generally underestimate the gap returning overestimated cell polarizabilities α i j [START_REF] Bishop | Molecular vibrational and rotational motion in static and dynamic electric fields[END_REF] and dielectric constants. Hence, it is expected to affect the refractive indices (equation 8). Nonetheless, despite the apparent deviations in functional performance, all methods considered predict relatively stable birefringence values (denoted δn in Table 1), in good agreement with the available experimental measurements conducted on stoichiometric LNO crystals (cationic ratio Li/Nb ∼ 1) [START_REF] Schlarb | Refractive indices of lithium niobate as a function of wavelength and composition[END_REF][START_REF] Bergman | Curie temperature, birefringence, and phase-matching temperature variations in LiNbO 3 as function of melt stoichiometry[END_REF].

Table 1: Band gap, refractive indices n i i , the birefringence δn, and dielectric constants i i of LNO calculated at 633 nm wavelength using different functionals and compared to experimental results. Values in parentheses correspond to the static values (at "infinite" wavelength).

B3LYP

PBE0 PBE LDA [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF] exp [START_REF] Schlarb | Refractive indices of lithium niobate as a function of wavelength and composition[END_REF] exp [START_REF] Nelson | Refractive indices of congruently melting lithium niobate[END_REF] Gap (eV) The clamped EO coefficients r S i j k of LNO, together with the corresponding electronic and vibrational contributions computed with the B3LYP functional, are listed in Table 2. We also list previous published experimental measurements and theoretical data in the same table. All properties have been deducted from computed χ (2) i j k (-ω σ ; ω, 0) (for Pockels coefficient calculation, ω stands for a laser wavelength, ω = 0 for a static electric field, and ω σ is the sum of two wavelengths, see SM for more details [66]) obtained at the same level of theory, considering only the four EO independent tensorial components of LNO r 33 , r 13 , r 11 , and r 51 [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF]. The vibrational contribution is a result of the transverse optical phonon modes; they can be represented as 4A 1 + 5A 2 + 9E , where the A 1 and E modes are Raman and infrared (IR) active; in contrast, A 2 modes are silent. The A 1 modes are coupled to r 33 and r 13 , while the E modes are linked to r 11 and r 51 . The corresponding computed phonon frequencies, showing a good agreement with available experimental data, together with the second-order nonlinear susceptibilities χ (2) i j k (-ω σ ; ω, 0) at a wavelength of 633 nm, are given in SM [66]. Looking at the computed clamped EO coefficients presented in Table 2, the results obtained with B3LYP functional are in very good agreement with experimental results, notably for r 33 and r 13 coefficients. On the other hand, an overestimation of r 11 and an underestimation of r 51 , respectively, are revealed. Nevertheless, our calculations expose that the most dominant contribution out of the two considered for r S i j k is the vibrational contribution to the microscopic second hyperpolarizability of the unit cell. This result, relying on a hybrid DFT treatment of the wavefunction, confirms previous computations within the LDA approximation [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF], producing at the same time a very good estimation of the EO tensorial components with respect to the experimental data.

Table 2: r i j k (-ω σ ; ω, 0) EO tensor (pm/V) of LNO at 633 nm wavelength: electronic and ionic (clamped), and piezoelectric (unclamped) contributions, together with results from previous calculations and experiments.

B3LYP

LDA [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF] exp [START_REF] Turner | High-frequency electro-optic coefficients of lithium niobate[END_REF] We now consider the total unclamped EO coefficients r i j k of LNO, including the piezoelectric contribution, which is computed at the same level of theory as discussed above. We note that the photoelastic constants P i j µν , here computed at 633 nm, are found relatively independent on wavelength above 633 nm (see Figure S.2 in SM [66]). The same trend has been observed for MgO and CaWO 4 crystals from calculations at the PBE level [START_REF] Erba | Piezo-optic tensor of crystals from quantum-mechanical calculations[END_REF][START_REF] Erba | Photoelasticity of crystals from theoretical simulations[END_REF]. The piezoelectric contribution to the EO coefficients is listed in Table 2 as r P i j k . Looking at the total r i j k = r S i j k + r p i j k , an excellent agreement with experimental results is obtained. For the three coefficients r 33 , r 13 , and r 11 , the vibrational contribution is still dominant, while r 51 coefficient has the highest piezoelectric contribution (Table 2). Three functionals, B3LYP (Table 2), PBE0, and PBE (Table S.5), give a piezoelectric contribution to r 51 of around 11 pm/V, that represents around 33 % of r 51 unclamped value, lower than the 48 % found in Ref. [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF]. Nevertheless, the r 51 EO coefficient is found with the highest piezoelectric contribution from the four independent EO coefficients of LNO material, in line with the previous theoretical results [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF]. Experimentally, only a piezoelectric contribution of about 4.2 pm/V has been reported [START_REF] Turner | High-frequency electro-optic coefficients of lithium niobate[END_REF]. This disagreement might be due to the experimental method, which can not suppress the total piezoelectric contribution during the measurement of the clamped coefficient, as previously suggested in Ref. [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF].

It is worth noting that the PBE0 functional showed an overestimation of the piezoelectric contribution to r 33 (∼6 pm/V see Table S.5), which is higher than the value obtained with B3LYP and PBE functionals (∼1 pm/V close to the 1.8 pm/V found experimentally [START_REF] Turner | High-frequency electro-optic coefficients of lithium niobate[END_REF]). Moreover, pure PBE functional delivers a reasonable estimation of the EO coefficients (Table S.5), even though it underestimates the gap, thus, giving a high refractive index (Table 1) as discussed above. That is because underestimated gaps lead to overestimated susceptibility χ (2) i j k ; in return, this will compensate for the high values of the refractive indices when computing the EO coefficients r S i j k (see equation 5). The polarization was also computed for o-phase ZrO 2 and HfO 2 via Berry phase calculation using Quantum Espresso (QE) [START_REF] Giannozzi | QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials[END_REF], in order to verify that the optimized polar phases are consistent with the ones reported in literature [START_REF] Böscke | Ferroelectricity in hafnium oxide thin films[END_REF][START_REF] Materlik | The origin of ferroelectricity in Hf 1-x Zr x O 2 : A computational investigation and a surface energy model[END_REF]. Table 3 displays the band gap, refractive indices, and remanent polarization of o-phase ZrO 2 and HfO 2 ; Note that regarding optical properties, no experimental data are available for comparison. However, we did obtain polarization values similar to previous calculations [START_REF] Materlik | The origin of ferroelectricity in Hf 1-x Zr x O 2 : A computational investigation and a surface energy model[END_REF].

ZrO

It is well known that ZrO 2 and HfO 2 compounds have a wide band gap, as they are used for high-k applications. From Table 3 it can be seen that the band gap of o-HfO 2 is higher by about 1 eV compared to o-ZrO 2 . Another feature is that the refractive indices at 633 nm of o-ZrO 2 are higher than the ones of o-HfO 2 . This is interesting because a high refractive index is an essential element for silicon photonics applications in order to ensure a sufficient contrast with the refractive index of SiO 2 [START_REF] Messner | Plasmonic Ferroelectric Modulators[END_REF]. We also have to consider the vibrational modes of the Pbc2 1 phase, as these modes are directly linked and contribute to the EO coefficients. The primitive unit cell of Pbc2 1 phase with 12 atoms (4 formula units) results in 36 vibrational modes, 3 for translation modes, and 33 optical modes, having the following irreducible representation at the zone center Γ:

Γ = 8A 1 + 9A 2 + 8B 1 + 8B 2 (11)
The The orthorhombic Pbc2 1 with a 2mm point group has five independent EO coefficients [START_REF] Yariv | Photonics: Optical electronics in modern communications[END_REF]: r 13 , r 23 , r 33 , r 42 , and r 51 . Moreover, as given above, only A 1 , B 1 , and B 2 contribute to the EO coefficients in the Pbc2 1 phase. The A 1 modes are linked to r 13 , r 23 , and r 33 , the B 1 and B 2 modes are coupled to r 42 and r 51 , respectively. These five EO coefficients are computed at the same level of theory as discussed for LNO, and they are listed in Table 5 showing the three different contributions, electronic r el e i j k , vibrational r vi b i j k , and piezoelectric r p i j k . The first significant result is that the EO coefficients in o-ZrO 2 are higher and almost twice the ones in o-HfO 2 . The r 13 and r 33 coefficients have the highest values of the five independent EO elements. For both ZrO 2 and HfO 2 , these two coefficients have an opposite sign and are quite close in absolute values. Comparing the calculated EO results of o-ZrO 2 to those of the reference LNO material (Table 2), the r 13 of o-ZrO 2 is close to the r 13 of LNO, while the r 33 of o-ZrO 2 is around one-third of LNO. Unfortunately, no experiment results are available for pure o-ZrO 2 and o-HfO 2 . However, for Y 2 O 3 -doped HfO 2 o-phase, EO coefficients have been reported with effective values of about 0.46 pm/V and 0.67 pm/V for (100) and (111) oriented thin films, respectively [START_REF] Kondo | Linear electro-optic effect in ferroelectric HfO 2 -based epitaxial thin films[END_REF][START_REF] Kondo | Influence of orientation on the electro-optic effect in epitaxial Y-doped HfO 2 ferroelectric thin films[END_REF].

From Table 5, we also deduce that ionic is the dominant contribution to the EO coefficients in o-phase ZrO 2 and HfO 2 , as in LNO. The approach used here for EO calculations also provides further insights into this ionic contribution computed by evaluating the Raman and IR intensities [START_REF] Rérat | Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab Initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme[END_REF], as described in refs. [START_REF] Maschio | Comment on "Ab initio analytical infrared intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method[END_REF] and [START_REF] Maschio | Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis[END_REF], respectively. In the case of LNO crystal, we found that the two A 1 and E polar modes of lowest frequency contribute the most to r 33 and r 51 , respectively, amounting to 43% and 36% of the value of these coefficients. This observation is in good agreement with previous calculations [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF], and with the fact that for LNO, unstable phonon modes in the paraelectric phase transform into lower frequency and highly polar modes in the low-temperature ferroelectric phase [START_REF] Veithen | First-principles study of the electro-optic effect in ferroelectric oxides[END_REF]. Nevertheless, ZrO 2 and HfO 2 materials have numerous polymorphs [START_REF] Kisi | Crystal structures of zirconia phases and their inter-relation[END_REF] and other mechanisms that can stabilize the ferroelectric phase [START_REF] Materlik | The origin of ferroelectricity in Hf 1-x Zr x O 2 : A computational investigation and a surface energy model[END_REF][START_REF] Kisi | Influence of hydrostatic pressure on the t→o transformation in mg-psz studied by in situ neutron diffraction[END_REF], away from the classical paraelectricferroelectric phase transition picture. However, looking at the lowest phonon modes in o-ZrO 2 , the same trend is observed as in LNO material. For example, the lowest frequency A 1 TO (transverse optical) mode at 102 cm -1 (Table 4) is responsible for about 40% and more than 60% of the values of r 33 and r 13 in o-ZrO 2 , respectively. Table 5: The five independent elements of r i j k (-ω σ ; ω, 0) tensor for Pbc2 1 space group, together with their respective electronic, vibrational and piezoelectric contributions, computed at a wavelength of 633 nm. All values are given in pm/V. o-ZrO 2 o-HfO 2 r el e i j k

r vi b i j k r p i j k r i j k r el e i j k r vi b i j k r p i j k r i

Wavelength dependence of EO coefficients:

The EO values reported above are computed at 633 nm, while, for example, wavelengths for telecom applications are in the near-infrared region, typically around 1500 nm. Therefore, with a focus on the influence of laser wavelength on Pockels coefficients, the same calculations as performed at 633 nm were repeated at different wavelengths up to near-infrared. The results are displayed in Figure 2 for the unclamped r 33 and r 13 coefficients of LNO and o-ZrO 2 . ZrO 2 is here chosen instead of HfO 2 because of its higher EO coefficients compared to HfO 2 (Table 5). First, bearing in mind that the three contributions depend on wavelength, the electronic part is more sensitive to the electric field than the vibrational and piezoelectric contributions. Also, the electronic part depends strongly on the gap of the materials. Therefore, approaching the IR wavelengths, the electronic contribution decreases and converges to a static value [START_REF] Maschio | Calculation of the dynamic first electronic hyperpolarizability β(-ω σ ; ω 1 , ω 2 ) of periodic systems. Theory, validation, and application to multi-layer MoS 2[END_REF][START_REF] Rérat | Computation of Second Harmonic Generation for Crystalline Urea and KDP. An ab Initio Approach through the Coupled Perturbed Hartree-Fock/Kohn-Sham Scheme[END_REF]. This decrease is associated with a decrease in the refractive indices: At a wavelength of 1550 nm, the refractive index n zz is found to be 2.08 (2.19) for LNO (for ZrO 2 ), which is lower than the value obtained at 633 nm (Tables 1 and3). We found that the vibrational contribution remains almost constant from 633 nm to 1550 nm. The last contribution is piezoelectric, where the dependency on wavelength comes from the photoelastic coefficients, which remain relatively constant at high wavelengths (see Figure S1 in SM [66]). The general relative increase of the total EO coefficients observed in Figure 2 when increasing the wavelength is thus due to the decrease of the refractive indices (see equation 9). The EO coefficients of o-ZrO 2 are found lower than the ones of LNO at all investigated wavelengths. Nevertheless, in terms of applications, other parameters should be taken into consideration, such as the refractive index and the total dielectric constant, to define a figure of merit (n 3 r i j k / [START_REF] Sando | Epitaxial ferroelectric oxide thin films for optical applications[END_REF], with a total dielectric constant comprising both the electronic and the ionic contributions) that could bring o-phase ZrO 2 /HfO 2 as a competitive material for EO applications. Figure 3 displays the comparison of this figure of merit between LNO, HfO 2 and ZrO 2 materials. For r 33 coefficient, LNO has the highest figure of merit; however, the ones for ZrO 2 and HfO 2 are not negligible. Considering the compatibility of ZrO 2 and HfO 2 with silicon, these results are significant enough for further investigating different paths to increase the EO coefficients in these materials.

ZrO 2 -HfO 2 rhombohedral R3m

In this section, we focus on the EO properties of the rhombohedral R3m (space group No. 160) phase of ZrO 2 -HfO 2 . This phase is less studied than the orthorhombic phase and has only been reported in compressively strained epitaxial thin films [START_REF] Wei | A rhombohedral ferroelectric phase in epitaxially strained Hf 0.5 Zr 0.5 O 2 thin films[END_REF][START_REF] Boutaybi | Stabilization of the epitaxial rhombohedral ferroelectric phase in ZrO 2 by surface energy[END_REF]. For the EO investigation, we also used a unit cell of 12 atoms (4 formula units). One should note that in the case of the r-phase, an instability (with negative/imaginary frequency) of TO A 1 -TO 1 phonon mode was observed when the phonon was computed at the optimized unit cell volume for r-ZrO 2 (141 Å 3 ) and r-HfO 2 (136 Å 3 ), the optimized unit cells being given in SM [66]. For the calculation of EO coefficients and other properties to be reliable, one should have only real frequencies. To achieve such a condition for ZrO 2 , an experimental volume (∼132 Å 3 ) was imposed together with a rhombohedral angle equal to 89.56° [START_REF] Boutaybi | Stabilization of the epitaxial rhombohedral ferroelectric phase in ZrO 2 by surface energy[END_REF], forcing the imaginary phonon mode that is calculated at 82i cm -1 to move to 76 cm -1 at the experimental volume. While for HfO 2 , no r-phase unit cell volume was experimentally reported, the volume was reduced until achieving a positive phonon frequency, and the volume used here is 128 Å 3 . The computed band gap, refractive indices, and polarization are displayed in Table 6. The band gaps obtained for ZrO 2 and HfO 2 in the r-phase are significantly higher than in the o-phase. The remanent polarization for the r-phase was also computed via Berry Phase using QE. The value obtained for the polarization is lower than in the case of the o-phase, as expected from previous calculations [START_REF] Wei | A rhombohedral ferroelectric phase in epitaxially strained Hf 0.5 Zr 0.5 O 2 thin films[END_REF][START_REF] Silva | Wake-up Free Ferroelectric Rhombohedral Phase in Epitaxially Strained ZrO 2 Thin Films[END_REF]. Note that this polarization depends strongly on the compressive strain, which is directly related to the rhombohedral angle [START_REF] Wei | A rhombohedral ferroelectric phase in epitaxially strained Hf 0.5 Zr 0.5 O 2 thin films[END_REF][START_REF] Boutaybi | Stabilization of the epitaxial rhombohedral ferroelectric phase in ZrO 2 by surface energy[END_REF]. Here, the polarization is computed at a fixed angle using the rhombohedral primitive unit cell. In contrast, in Refs. [START_REF] Wei | A rhombohedral ferroelectric phase in epitaxially strained Hf 0.5 Zr 0.5 O 2 thin films[END_REF][START_REF] Silva | Wake-up Free Ferroelectric Rhombohedral Phase in Epitaxially Strained ZrO 2 Thin Films[END_REF] the polarization is computed with hexagonal unit cells and at different compressive strains.

Table 7: The four independent elements of r i j k (-ω σ ; ω, 0) tensor for R3m space group, together with their respective electronic, vibrational and piezoelectric contributions, computed at a wavelength of 633 nm. All values are given in pm/V. r-ZrO 2 r-HfO 2 r el e i j k As seen for o-phase (section 3.2), a unit cell with 12 atoms results in 36 vibrational modes, 3 for translation modes, and 33 for optical phonon modes. In the case of R3m space group, there are three different modes A 1 , A 2 , and E ; the irreducible representation at Γ zone center is:

r vi b i j k r p i j k r i j k r el e i j k r vi b i j k r p i j k r i j k r 33 -0.
Γ = 8A 1 + 3A 2 + 11E (12)
The A 1 and (degenerate) E modes are simultaneously Raman and IR active, while the A 2 modes are Raman active and IR inactive. The calculated A 1 mode frequencies are shown in Table 4 for ZrO 2 and HfO 2 , experimentally no data are available for ZrO 2 -HfO 2 in r-phase. The E and A 2 phonon frequencies are given in SM (Table S.8 in SM [66]) Following the same method as described in sections 3.1 and 3.2, the EO coefficients obtained for ZrO 2 and HfO 2 r-phase are given in Table 7. Comparing these EO coefficients to the ones of LNO and o-phase, r-phase shows very modest EO values. This trend is also observed in terms of remanent polarization obtained for r-phase compared to orthorhombic Pbc2 1 (Tables 3 and6). Strain dependence of EO coefficients: Because the ZrO 2 -HfO 2 R3m phase has a strong dependence on the compressive strain, r 33 EO coefficient of r-ZrO 2 is computed at different rhombohedral angles. The results are displayed in Figure 4, showing a linear dependence of r 33 with respect to the rhombohedral angle. Even at a high compressive strain (rhombohedral angle of 88°), the EO coefficient is still lower than the one calculated for o-ZrO 2 . Two main reasons cause the low EO coefficients of r-ZrO 2 compared to o-ZrO 2 . First is the relatively high band gap in the r-phase, which reduces the electronic contribution. The second reason is that the A 1 low-frequency phonon modes have negligible contributions to EO coefficients in the r-phase, while similar modes highly contribute to the ones of the o-phase (see section 3.2). In more details, for the r-phase, from 8 A 1 active modes (Table 4), only 3 modes at 220, 352, and 449 cm -1 contribute significantly to r 33 and r [START_REF] Wang | Integrated lithium niobate electro-optic modulators operating at CMOScompatible voltages[END_REF] . The modes at 220 and 352 cm -1 give almost the same value of the hyper-polarizabilities β 33 (equation 6), but with opposite signs, hence they cancel each other, therefore solely the mode at 449 cm -1 is contributing to r 33 . This is the most probable reason why the r-phase shows a lower EO coefficient than the o-phase. Finally, this gives r-phase a very low figure of merit compared to o-phase and to LNO as illustrated in Figure 3.

Conclusion

In conclusion, we investigated the electro-optic Pockels coefficients of the ZrO 2 and HfO 2 compounds in two noncentrosymmetric phases: rhombohedral R3m and orthorhombic Pbc2 1 . We determined that the o-phase has considerably higher EO coefficients than the r-phase. This remains true for a highly strained rhombohedral phase; the optical nonlinearity is weak compared to the o-phase. Moreover, r 33 and r 13 in the orthorhombic phase are higher in ZrO 2 than in HfO 2 , with values of about 12 and 7 pm/V for ZrO 2 and HfO 2 , respectively. Comparing the EO coefficients of o-ZrO 2 to rhombohedral LiNbO 3 reference material, the r 13 coefficient shows a comparable value to the r 13 of LiNbO 3 . In contrast, the r 33 coefficient of o-ZrO 2 is three times lower than the r 33 of LiNbO 3 . We observed that the polar and low-frequency phonon modes contribute most to the EO in the orthorhombic phase; the same trend is also observed in LiNbO 3 . Our numerical results suggest that the ferroelectric orthorhombic phase of ZrO 2 and HfO 2 -based thin films can be relevant for EO applications when considering their compatibility with silicon. Finally, considering the dominant contribution by the ionic part as identified in this work, we believe that doping ZrO 2 /HfO 2 with a suitable element to stabilize the orthorhombic phase but also increase the ionicity of these materials will lead to high EO coefficients. Thus, ZrO 2 /HfO 2 -based compounds may be promising candidates for EO properties and their use in photonic devices.

SUPPLEMENTAL MATERIAL EO coefficients of LiNbO 3 and ZrO 2 -HfO 2 materials

The electro-optic (EO) coefficients for 3m space group (SG) [START_REF] Yariv | Photonics: Optical electronics in modern communications[END_REF], thus for ZrO 2 -HfO 2 rhombohedral R3m (SG N°160) and LiNbO 3 R3c (SG N°161), have 4 independent EO coefficients with two possible representation:

m ⊥ x 1

Symmetry of β (SHG) and β (EO / dc-Pockels)

In the general case, the second order non-linearity β i j k is written as

β i j k = β i j k (-ω σ ; ω 1 , ω 2 ) (3)

Second harmonic generation (SHG)

For SHG, β i j k (ω σ = 2 ω and ω 1 = ω 2 = ω ) can be written as

β i j k = β i j k (-2ω; ω, ω) (4)
So in the case of second harmonic generation, for example with β y x y equal to β y y x , in general the symmetry given in (4) can simply be written as β i j k (-2ω; ω, ω) = β i k j (-2ω; ω, ω).

EO / dc-Pockels effect

For dc-Pockels effect, β i j k (ω σ = ω , ω 1 = ω, and ω 2 = 0) is written as

β i j k = β i j k (-ω; ω, 0) = β i k j (-ω; 0, ω) (5)
The same permutation as given above for SHG then leads to β xxz (-ω; ω,0) = β xzx (-ω; 0, ω). In other words, β xxz (-ω; ω,0) tensor element, where the static field is along z, is equal to β xzx (-ω; 0, ω), so that the static electric field is always applied along z direction. But β xxz (-ω; ω,0) is not equal to β xzx (-ω; ω, 0) because in the latter tensor element the static electric field is along x. Finally, this will ensure that r i j k coefficients respect the permutation r i j k = r j i k .

Results

LiNbO The refractive indices of LiNbO3 calculated using B3LYP functional, together with the electronic and vibrational contributions to the polarizability α i j , are reported in Table 2. The ionic contribution at 633 nm is almost null. At the Γ point, the optical phonons of LiNbO 3 can be classified according to the irreducible representation of the space group R3c into Γ = 4A 1 + 5A 2 + 9E . The A 1 and E modes are Raman and infrared active, while A 2 modes are inactive.

The calculated frequencies of the A 1 and E modes are shown in Table 3 and compared with previous theoretical and experimental data. There is a good agreement between our calculations and available experimental values. In particular, B3LYP functional calculated frequencies are close to the experimental ones, with only a relative underestimate of the first A 1 mode. This might be due to the anharmonicity of this mode: The experimental frequency is substantially higher than the one calculated in the harmonic approximation; Considering the one-dimensional non-interacting oscillator approximation, a good approximation to the experimental value was obtained by Inbar and Cohen [START_REF] Inbar | Comparison of the electronic structures and energetics of ferroelectric LiNbO 3 and LiTaO 3[END_REF]. Increasing HF exchange (PBE0) gives relatively higher frequency values, while decreasing exchange (PBE(10%) and PBE) brings lower values compared to the experimental results. In order to obtain the EO coefficients of LiNbO 3 using CRYSTAL17, we computed the electronic and vibrational contributions to the second hyperpolarizability β i j k (-ω σ ; ω, 0) at a wavelength of 633 nm as described in the main paper, then deduced the χ (2) i j k (-ω σ ; ω, 0) values that are given in Table 4 (see the main paper for the relation between the susceptibility χ (2) i j k and the hyperpolarizability β i j k ). For symmetry reasons, only 11 out of 27 elements of the rank three χ (2) i j k tensor are not null and only four elements are independent. These four elements are given in Table 4. We note that pure PBE functional gives the highest values of χ (2) i j k , that is because pure DFT-PBE underestimates the gap, as discussed in the main paper. Table 4: Non-zero elements of χ (2) i j k (-ω σ ; ω, 0) (pm/V) tensor in LiNbO 3 computed at a wavelength of 633 nm using different functionals, with their vibrational and electronic contributions.

B3LYP

PBE0 PBE(10%) PBE χ (2) i Table 5 gives the clamped and unclamped EO coefficients of LiNbO 3 , considering only the four independent elements of r i j k , and compares them with previous calculations and experiments. All functionals used in the present study agree well with experimental results, notably for the r 33 and r 13 coefficients. The clamped EO coefficients, taking in consideration the piezoelectric contribution, are also given in Table 4. As the experimental values of EO coefficients are usually measured at finite electric field wavelength, in Figure 1 we show the dependence of photoelastic coefficients on the electric field wavelength (nm). From 633 nm up to the near-infrared wavelength range, the two coefficients remain almost constant, while below 633 nm a decrease is observed. The orthorhombic Pbc2 1 phase with a unit of 12 atoms (4 formula units) results in 36 vibrational modes, 3 for translation modes, and 33 optical phonon modes having the following irreducible representation at the zone center Γ:

Γ = 8A 1 + 9A 2 + 8B 1 + 8B 2
The computed phonon frequencies are listed in Table 1 .

For the rhombohedral phase R3m, also with a unit of 12 atoms (4 formula units) resulting in 36 vibrational modes, 3 for translation modes and 33 optical phonon modes, the irreducible representation at the zone center Γ is:

Γ = 8A 1 + 3A 2 + 11E
The calculated frequencies are given in Table 1. Note that the computed E modes here are degenerate. The lowest E-TO1 (transverse optical) mode frequency depends strongly on the unit cell volume and can also be affected by other factors, like strain, oxygen vacancies, etc.

Figure 1 :

 1 Figure 1: The crystallographic structures of a) LiNbO 3 (space group R3c), b) rhombohedral ZrO 2 -HfO 2 (space group R3m), and c) orthorhombic ZrO 2 -HfO 2 (space group Pbc2 1 ).

2 -HfO 2 orthorhombic Pbc2 1

 21 In this section, we present electro-optics properties of the o-phase of ZrO 2 and HfO 2 compounds, with space group Pbc2 1 (No. 29). The optimization of Pbc2 1 ZrO 2 and HfO 2 is performed using 12 atoms in a unit cell. The optimized crystalline structure and atomic positions are given in SM[66].

2 Figure 2 :

 22 Figure 2: EO coefficients (unclamped) r 33 and r 13 of LNO and o-ZrO 2 at different wavelengths using B3LYP functional.

Figure 3 :

 3 Figure of merit (pm/V)

Figure 4 :

 4 Figure 4: r 33 EO coefficient of r-ZrO 2 computed at different angles.

  the case of Pbc2 1 (SG N°29 and 2mm point group), there are five independent electro-optic coefficients

Figure 1 :

 1 Figure 1: Independent photoelastic p 52 and p 21 coefficients of the LiNbO 3 crystal computed at B3LYP level as a function of the electric field wavelength (nm).

ZrO 2 -

 2 HfO 2

Table 3 :

 3 Band gap and refractive indices of orthorhombic Pbc2 1 ZrO 2 and HfO 2 at 633 nm. The polarization P b (along the b axis) from Berry phase calculation is also shown.

		o-ZrO 2	o-HfO 2
	gap (eV)	5.40	6.45
	n xx	2.204	2.059
	n y y	2.301	2.122
	n zz	2.236	2.081
	P b (µC/cm 2 ) 54 (58, [77] ) 49.5 (50, [77] )

Table 4 :

 4 A 1 , B 1 , and B 2 modes are Raman and IR active. Hence, they will contribute to the electro-optic coefficients, while A 2 modes are only Raman active. The calculated A 1 phonon frequencies are given in Table4, and the A 2 , B 1 , and B 2 phonon frequencies are given in SM[66]. Only a few experimental results concerning the phonon frequencies of Pbc2 1 o-phase are available for pure ZrO 2 and HfO 2 . They are displayed for o-ZrO 2 in Table4, showing a good agreement with the calculated phonon modes. Calculated frequencies (cm -1 ) of Raman and IR active A 1 phonon modes of ferroelectric ZrO 2 and HfO 2 using B3LYP functional.

			o-ZrO 2	o-HfO 2	r-ZrO 2	r-HfO 2
			exp [78]
	A 1	102	/	113	171	118
		186	200	153	220	171
		291	/	253	267	267
		315	320	295	352	284
		351	340	329	449	423
		387	/	369	534	537
		445	/	451	674	692
		566	580	589	755	770

Table 6 :

 6 Band gap and refractive indices of rhombohedral R3m ZrO 2 and HfO 2 at 633 nm. The rhombohedral polarization P r from Berry phase calculation is also shown.

		r-ZrO 2 r-HfO 2
	gap (eV)	5.79	6.77
	n xx	2.282	2.114
	n zz	2.271	2.107
	δn	0.011	0.007
	P r (µC/cm 2 )	6.8	5.1

Table 1 :

 1 The optimized geometry for LiNbO 3 , on the right the LiNbO 3 R3c structure.

	3
	LiNbO 3	a (Å)	b (Å)	c (Å)
		5.216	5.216	14.14
	Atom	x/a	y/b	z/c
	Li	0	0	0.2614
	Nb	0	0	-0.0188
	O	-0.0473 0.3435	0.0441
	V (Å 3 )	111

Table 2 :

 2 Refractive indices n i i of LiNbO 3 at 633 nm wavelength using B3LYP functional, showing the electronic and vibrational contributions to the polarizability α i j .

		electronic vibrational	n i i
	α xx	244.77	-0.60	2.257
	α y y	244.77	-0.60	2.257
	α zz	216.36	-0.48	2.149

Table 3 :

 3 Raman and IR active phonon mode frequencies (cm -1 ) of ferroelectric LiNbO 3 calculated with several functionals, compared with previous theoretical and experimental results.

		B3LYB PBE0 PBE(10%) PBE LDA [83] exp [84] exp [85]
	A1	232	251	242	234	208	256	/
		275	291	279	271	279	275	/
		340	346	332	322	344	332	/
		636	650	628	611	583	637	/
	E	143	152	144	140	151	/	155
		233	241	229	221	/	/	238
		246	259	250	245	236	/	265
		304	315	302	294	307	/	325
		344	367	252	343	334	/	/
		361	380	365	355	352	/	/
		427	434	414	402	432	/	431
		576	589	572	560	526	/	582
		655	674	654	641	617	/	668

Table 5 :

 5 Total r i j k (-ω σ ; ω, 0) EO tensor (pm/V) in LiNbO 3 computed at a wavelength of 633 nm using different functionals, together with previous calculations and experiments.

		B3LYB PBE0 PBE(10%)	PBE	LDA [72] exp [11] exp [74]
	clamped	r 33 30.64 28.02	30.45	30.87	26.9	30.8	34
		r 13 10.68	9.14	11.39	12.84	9.7	8.6	10.9
		r 11 -7.86 -6.06	-7.79	-9.59	4.6	3.4	/
		r 51 22.51 19.41	21.99	24.51	14.9	28	/
	unclamped							exp [75]
		r 33 31.59 34.32	30.90	31.20	27	32.6	/
		r 13 12.10 10.92	12.61	13.81	10.5	10	/
		r 11 -7.00 -5.00	-6.84	-10.20	7.5	6	9.9
		r 51 33.28 27.09	31.94	34.04	28.6	32.2	/

Table 6 :

 6 The optimized geometry for ZrO 2 orthorhombic Pbc2 1 .

	o-ZrO 2	a (Å)	b (Å)	c (Å)
		5.34	5.14	5.16
	Atom	x/a	y/b	z/c
	Zr	0.0274 0.2644 -0.2545
	O1	0.3602 0.0619 -0.1021
	O2	0.2262 -0.4688 -0.0064
	V (Å 3 ) 141.62

Table 7 :

 7 The optimized geometry for HfO 2 orthorhombic Pbc2 1 .

	o-ZrO 2	a (Å)	b (Å)	c (Å)
		5.28	5.07	5.10
	Atom	x/a	y/b	z/c
	Zr	0.0312 0.2651	-0.2551
	O1	0.3618 0.0638	-0.1051
	O2	0.2267 -0.4699 -0.00278
	V (Å 3 ) 136.52

Table 8 :

 8 Raman and IR phonon mode frequencies (cm -1 ) of ferroelectric ZrO 2 and HfO 2 calculated using B3LYP functional. Experimental data are also reported when available.

		o-ZrO 2		o-HfO 2		r-ZrO 2	r-HfO 2
		B3LYP exp [78]	B3LYP		B3LYP	B3LYP
	A 1	102	/	113	A 1	171	118
		186	200	153		220	171
		291	/	253		267	267
		315	320	295		352	284
		351	340	329		449	423
		387	/	369		534	537
		445	/	451		674	692
		566	580	589		755	770
	A 2	172	/	122	A 2	169	117
		180	/	134		332	312
		192	200	148		586	595
		295	/	303
		400	/	414	E	76	20
		482	/	496		165	114
		561	/	581		174	123
		634	/	654		233	187
		659	/	674		325	271
	B 1	188	/	141		358	304
		224	/	218		456	436
		298	/	238		536	540
		329	/	311		586	594
		402	/	395		683	702
		514	/	533		721	726
		545	550	554
		726	/	736
	B 2	165	/	121
		233	/	226
		289	/	269
		345	/	312
		399	/	369
		461	/	484
		622	620	637
		683	/	697

Acknowledgements: This work has received support from the Agence nationale de la recherche (ANR) under project FOIST (N°ANR-18-CE24-0030), and from the French national network RENAT-ECH for nanofabrication.