Dennis Trujillo

Ashok Gurung

Jiacheng Yu

Sanjeev K Nayak

S Pamir Alpay

Pierre-Eymeric Janolin
email: pierre-eymeric.janolin@centralesupelec.fr

Data-Driven Methods for Discovery of Next-Generation Electrostrictive Materials

Keywords: Electrostriction, data mining, first-principles calculation, halides

All dielectrics exhibit electrostriction, i.e., display a quadratic strain response to an electric field compared to the linear strain dependence of piezoelectrics. As such, there is significant interest in discovering new electrostrictors with enhanced electrostrictive coefficients, especially as electrostrictors can exhibit effective piezoelectricity when a bias electric field is applied. We present the results of a study combining data mining and first-principles computations that indicate that there exists a group of iodides, bromides, and chlorides that have electrostrictive coefficients exceeding 10 m 4 C -2 which are substantially higher than typical oxide electrostrictive ceramics and polymers. The corresponding effective piezoelectric voltage coefficients are three orders of magnitude larger than lead zirconate titanate. Trujillo et al., Discovery of Novel Electrostrictors, page 2

INTRODUCTION

The development of novel materials exhibiting exceptional strain response in the presence of an external electric field is imperative to the development of tunable electromechanical devices [START_REF] Fanson | Articulating fold mirror for the Wide-Field/Planetary Camera II[END_REF][START_REF] Uchino | Electrostrictive Actuators: Materials and Applications[END_REF][START_REF] Li | Electrostrictive effect in ferroelectrics: An alternative approach to improve piezoelectricity[END_REF][START_REF] Uchino | Electrostrictive effect in perovskites and its transducer applications[END_REF][START_REF] Newnham | Electrostriction: Nonlinear electromechanical coupling in solid dielectrics[END_REF][START_REF] Yuan | Giant Electrostriction of Soft Nanocomposites Based on Liquid Crystalline Graphene[END_REF][START_REF] Li | Giant thermally-enhanced electrostriction and polar surface phase in La2Mo2O9 oxygen ion conductors[END_REF][START_REF] Korobko | Giant Electrostriction in Gd-Doped Ceria[END_REF] . Promising applications involve next-generation transducers [START_REF] Chen | Electrostrictive transducers for medical ultrasonic applications[END_REF] , ultra-efficient actuators [START_REF] Fanson | Articulating fold mirror for the Wide-Field/Planetary Camera II[END_REF][START_REF] Uchino | Electrostrictive Actuators: Materials and Applications[END_REF][START_REF] Damjanovic | Electrostrictive and Piezoelectric Materials for Actuator Applications[END_REF] , sonars [START_REF] Pilgrim | Electrostrictive Ceramics for Sonar Projectors[END_REF] , and low power logic devices [START_REF] Anam | Proposal for an electrostrictive logic device with the epitaxial oxide heterostructure[END_REF] with application spaces in medicine, astronomy, and consumer electronics. Recent discoveries of 'giant' electrostrictors, for example, Gd-doped ceria [START_REF] Korobko | Giant Electrostriction in Gd-Doped Ceria[END_REF] , and (Nb, Y) doping in bismuth oxide [START_REF] Yavo | Large Nonclassical Electrostriction in (Y, Nb)-Stabilized δ-Bi2O3[END_REF] and lanthanum molybdenum oxide [START_REF] Li | Giant thermally-enhanced electrostriction and polar surface phase in La2Mo2O9 oxygen ion conductors[END_REF] , soft nano-composite based on carbon nanotubes (CNT) [START_REF] Luna | Giant Electrostrictive Response and Piezoresistivity of Emulsion Templated Nanocomposites[END_REF] and liquid crystalline graphene [START_REF] Yuan | Giant Electrostriction of Soft Nanocomposites Based on Liquid Crystalline Graphene[END_REF] have opened promising applications for superior electromechanical properties, and also as a replacement of lead-based Pb(Ti 1-x Zr x )O 3 (PZT) solid solutions that are widely used in current applications. Of particular interest is the possibility to induce an effective piezoelectric behavior under bias electric field, generating extraordinarily large piezoelectric coefficients (d 33

~

pC N -1 ) [START_REF] Park | Induced giant piezoelectricity in centrosymmetric oxides[END_REF] . These 'giant' electrostrictors have been reviewed recently [START_REF] Yu | Defining 'Giant' Electrostriction[END_REF] . There is a clear advantage in utilizing electrostrictive materials that exhibit a quadratic strain in response to an electric field. This may allow for a more significant strain response for a given applied external field, as well as take advantage of the anhysteretic nature of electrostriction for highprecision displacements.

This effect is mediated via an applied electric field ( ) and the dielectric polarization ( ) vectors as described by:

(

(

where are the strain components, and are piezoelectric tensors, and and are fourth-rank tensors describing the electrostrictive response. In Equation (1), the total strain (at constant stress) is expressed in terms of an applied electric field, whereas in Equation (2) it is given as a function of the polarization. The M tensor is of particular interest for actuator applications whereas Q is used as a design criterion in sensing applications.

The search for novel electrostrictors has historically appeared bounded by the universal empirical relation connecting the electrostrictive coefficient to the ratio of the elastic compressibility over the permittivity [START_REF] Newnham | Electrostriction: Nonlinear electromechanical coupling in solid dielectrics[END_REF] until the discovery of 'giant' electrostrictors that are several orders of magnitude greater than the empirically expected values [START_REF] Yu | Defining 'Giant' Electrostriction[END_REF] . These include [START_REF] Park | Induced giant piezoelectricity in centrosymmetric oxides[END_REF] , lead halide perovskites ( m 2 V -2 and 1266 -1417 m 4 C - 2 ) [START_REF] Chen | Large electrostrictive response in lead halide perovskites[END_REF] , and Ce 0.8 Gd 0.2 O 1.9 ( m 2 V -2 ) [START_REF] Yuan | Giant Electrostriction of Soft Nanocomposites Based on Liquid Crystalline Graphene[END_REF][START_REF] Li | Giant thermally-enhanced electrostriction and polar surface phase in La2Mo2O9 oxygen ion conductors[END_REF][START_REF] Korobko | Giant Electrostriction in Gd-Doped Ceria[END_REF] . The electrostriction coefficients of these materials are three orders larger than typical perovskites, including BaTiO 3 , PbTiO 3 , SrTiO 3 , and Pb(Ti 1-x Zr x )O 3 for which is in the range 0.02 -0.07 m 4 C - 2 5 .

materials such as La 2 Mo 2 O 9 ( m 2 V -2 and m 4 C -2 ), (Nb,Y)-stabilized Bi 2 O 3 ( - m 2 V -2 and 83-273 m 4 C -2 ) 13 , CeO 2-x ( 300 m 2 V -2 )
Using large open materials databases composed of accurately determined materials properties, the programmatic search for materials exhibiting specific properties becomes possible. In line with data-driven methods, advanced modeling approaches across length scales, and the successful insertion of the concept of materials genomics that allows for rapid development and deployment of novel materials [START_REF] Tanaka | Data-centric science for materials innovation[END_REF][START_REF] Suh | Evolving the Materials Genome: How Machine Learning Is Fueling the Next Generation of Materials Discovery[END_REF] , we employ here a combination of data mining and first-principles modeling in the discovery of materials exhibiting exceptional electrostrictive coefficients. Our results show several non-oxide inorganic compounds with electrostrictive properties of one or two orders of magnitude higher than existing oxide ceramics and with comparable or higher electrostriction coefficients of polymers.

RESULTS

Assessment of Electrostrictive Properties

The tensorial notation that has been employed in Equations ( 1) and ( 2) can be reduced to the following criteria to assess the magnitude of the electrostrictive properties of a wide variety of materials (see Supplementary Information), which is also defined as the Newnham's proxy , based on the approach presented in Refs. [START_REF] Newnham | Electrostriction: Nonlinear electromechanical coupling in solid dielectrics[END_REF][START_REF] Damjanovic | Electrostrictive and Piezoelectric Materials for Actuator Applications[END_REF] . can be derived from as below:

(3)

Here is the average eigenvalue of the dielectric tensor with , the relative permittivity of the material, and the permittivity of vacuum. in Equation ( 3) is the Reuss elastic compressibility given as [START_REF] Toonder | The relation between single crystal elasticity and the effective elastic behaviour of polycrystalline materials: theory, measurement and computation[END_REF][START_REF] De | Charting the complete elastic properties of inorganic crystalline compounds[END_REF] :

(

The Supplementary Information provides a more detailed treatment to define constitutive relations for the electrostrictive coefficients and . Equations ( 3) and (4) follow from specific electrical and mechanical boundary conditions corresponding to hydrostatic pressure/expansion and a single electric field component (uniaxial applied electric field in a short-circuit configuration) as discussed in the Supplementary Information. We note that Equation (3) describes the criterion used by Newnham et al. in their landmark paper that classified electrostriction for a range of materials over an extensive response space [START_REF] Newnham | Electrostriction: Nonlinear electromechanical coupling in solid dielectrics[END_REF] . We add here an additional criterion for M h [START_REF] Yu | Defining 'Giant' Electrostriction[END_REF][START_REF] Makagon | All-Solid-State Electro-Chemo-Mechanical Actuator Operating at Room Temperature[END_REF] . Although closely related, the quantities M h and Q h are not interchangeable as converse effects (Supplementary Information). Subtle differences in how these are defined give rise to slightly different criteria for the proxy derivation used to approximate M h and Q h . M h is readily accessible physical parameter in experiments.

In addition to the static empirical relations in Equation (3) and Equation ( 4), it is advantageous to use the differential forms of Q h and M h as described in the Supplementary Information, which is essentially the dielectric response with change in external hydrostatic pressure ( ). In the following, Q h and M h derived from our density functional theory (DFT) calculations are defined as and , respectively. Thus,

It must be noted that Equations ( 6) and ( 7) capture the mechanical response as differential in pressure ) and the equations are independent of Equation (3) and Equation (4).

DFT has been used to study electrostriction properties of compounds having high structural symmetry [START_REF] Jiang | Electrostriction coefficient of ferroelectric materials from ab initio computation[END_REF][START_REF] Kornev | Electrostriction Coefficients of GaN, AlN, MgO and ZnO in the Wurtzite Structure from First-Principles[END_REF] . Of particular interest in this regard is the study by Tanner et al. wherein the electrostrictive coefficients of cubic dielectric materials were computed using the differentials of the lattice parameters, which is similar to differentials of hydrostatic pressure for cubic crystal system [START_REF] Tanner | Optimized Methodology for the Calculation of Electrostriction from First-Principles[END_REF] . We have adopted the latter method using Equation ( 6) and Equation (7).

In our computations, the trace of external stress tensor, defined as pstress, is treated as the hydrostatic pressure with the conventional notation that positive and negative pstress values mean hydrostatic expansion and compression, respectively. This approach allows for uniform evaluation of electrostrictive coefficients over all seven crystal systems independent of their crystal structure.

Data Mining

To assign proxies for selecting materials exhibiting high electrostriction, we utilize a set of relationships between the dielectric response, the elastic coefficients, the electric field-induced strain, and the electrostrictive properties of materials as described in the previous section. A similar methodology developed by Newnham et al. has shown a semi-empirical connection between the elastic compliance, the dielectric response, and the hydrostatic electrostrictive

coefficient Q h 4,5 .
Existing databases of calculated functional properties provide easy access to data that can be used to shortlist materials with exceptional properties via a simple database query of related proxies [START_REF] Sahoo | Electronic and optical properties of zinc based hybrid organic-inorganic compounds[END_REF][START_REF] Nayak | Insight into point defects and impurities in titanium from first principles[END_REF] . The Materials Project (MP) serves as an ideal platform for novel materials discovery, given the large number of calculated material properties and optimized crystal structures contained in the database [START_REF] Jain | The Materials Project: A materials genome approach to accelerating materials innovation[END_REF] . The database contains 124,515 inorganic compounds with properties calculated via medium / high accuracy DFT methods. The MP team has calculated the elastic tensors for 13,751 and the dielectric tensors for 4,892 inorganic materials so far [START_REF] Jain | The Materials Project: A materials genome approach to accelerating materials innovation[END_REF] and made these publicly available [START_REF] Toonder | The relation between single crystal elasticity and the effective elastic behaviour of polycrystalline materials: theory, measurement and computation[END_REF] . In addition, there are currently 1,974 materials for which MP has calculated both the dielectric and elastic properties. This set of materials is ideal for data mining and training machine learning models to predict electrostrictive response. Starting from this set, we down-select using a set of simple property filters to find new promising candidates for electrostriction applications. We limit the materials by considering non-elemental (> 1 element) dielectrics with a bandgap > 0.5 eV and seek compounds that should exhibit relative stability as determined via the energy above the hull (< 0.1 eV). The schematic representation of the workflow is shown in Fig. (1).

These materials are then separated into centrosymmetric and non-centrosymmetric datasets based on their point group symmetries. This discretization serves the purpose of allowing us to consider materials lacking an inversion symmetry (non-centrosymmetric), which exhibit a spontaneous polarization via a piezoelectric response resulting from the dielectric and electrostrictive contributions separately from those that solely display an electrostrictive response (see Fig. 1). Taking these latter materials, it is then possible to consider the ratio of for each material and rank these according to the proxies described in Equations ( 3) and (4). To reduce the dielectric and elastic tensors to a scalar value for each data point, the inverse of the Reuss average bulk modulus was used in place of the elastic compressibility ( ) and the averaged dielectric tensor eigenvalue in place of , thus is considered a replacement for and a proxy for predicting . Alternatively, is the expression used for querying for high . These proxies are utilized throughout the rest of the study in combination with DFT calculated values for shortlisting materials which should exhibit exceptional electrostrictive strain response to an applied external electric field.

First Principles (DFT) Validation

The accuracy of DFT for materials properties that are determined at the electronic and atomic levels has been well tested; it provides a platform to bridge the gap between fundamental theory, high-throughput computational data science, and experiments [START_REF] Trujillo | Surface structure and energetics of low index facets of bismuth ferrite[END_REF][START_REF] Sahoo | Atomistic Insights into the Hydrogen Oxidation Reaction of Palladium-Ceria Bifunctional Catalysts for Anion-Exchange Membrane Fuel Cells[END_REF][START_REF] Adeagbo | Theoretical investigation of iron incorporation in hexagonal barium titanate[END_REF][START_REF] Sahoo | Segregation and ordering in binary transition metal clusters[END_REF][START_REF] Esquinazi | Defect-Induced Magnetism in Nonmagnetic Oxides: Basic Principles, Experimental Evidence, and Possible Devices with ZnO and TiO2[END_REF][START_REF] Sahoo | Surface phase diagrams of titanium in Oxygen, Nitrogen and Hydrogen environments: A first principles analysis[END_REF][START_REF] Nayak | Chromium point defects in hexagonal BaTiO3: A comparative study of firstprinciples calculations and experiments[END_REF][START_REF] Sahoo | Magnetic properties of small Pt-capped Fe, Co, and Ni clusters: A density functional theory study[END_REF][START_REF] Curry | Optical response of nickel-based superalloy Inconel-718 for applications in Trujillo et al., Discovery of Novel Electrostrictors, page 18 additive manufacturing[END_REF] . Here, we perform first-principles calculations with DFT using the plane-wave pseudopotential method [START_REF] Jones | The density functional formalism, its applications and prospects[END_REF] . The structural units of the compounds are taken from the MP database. We have used the PBEsol exchange-correlation functional [START_REF] Csonka | Assessing the performance of recent density functionals for bulk solids[END_REF] in place of the standard General Gradient Approximation (GGA) in the MP. PBEsol is a revised Perdew-Burke-Ernzerhof GGA chosen due to the enhanced prediction of equilibrium properties of densely packed solids, such as the dielectric and elastic tensors for insulating materials compared to GGA [START_REF] Tanner | Optimized Methodology for the Calculation of Electrostriction from First-Principles[END_REF] . Geometric optimization was carried out to obtain the equilibrium lattice parameters and corresponding atomic positions. The k-points grid size was set to 9 × 9 × 9, but models with atoms larger than 20 were treated with 7 × 7 × 7 grid size. The kinetic energy cutoff for the plane waves is 800 eV. The symmetry tag was switched on for the calculations so that the system did not deviate from the initialized structures.

Non-spin polarized calculations were performed as most compounds were not expected to have a magnetic ground state. Geometrical optimization for the models was carried out with the tolerance for total energy convergence set to 10 -5 eV and was performed in five cycles. The relaxed model from one cycle was then fed as the initial model for the next cycle to ensure a thorough relaxation. The calculations were performed with the Vienna ab initio simulation package (VASP) [START_REF] Kresse | Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF] .

The elastic compliance tensor was calculated by performing six finite distortions of the supercells and constructing the Hessian matrix from total energies obtained from density functional perturbation theory (DFPT) [START_REF] Baroni | Green's-function approach to linear response in solids[END_REF] . The static ion-clamped dielectric matrix was computed using DFPT via the linear Sternheimer equations [START_REF] Gajdoš | Linear optical properties in the projector-augmented wave methodology[END_REF][START_REF] Baroni | Ab initio calculation of the macroscopic dielectric constant in silicon[END_REF] . This method bypasses the complexity of computing the dielectric function through relaxation under finite electric fields. Equations ( 6) and ( 7) rely on the change in dielectric response as a function of hydrostatic pressure [START_REF] Francis | Finite basis set corrections to total energy pseudopotential calculations[END_REF] . The DFT computation of the dielectric function was done at 25 kbar pstress above and 25 kbar pstress below the equilibrium volume for which the pstress is zero. The mechanical perturbation of the order ± 25 kbar is well within the harmonic approximation and provides minimal deviation of the unit cell symmetry. We found slight variations in the equilibrium lattice parameters when compared to the MP. This is due to differences in the pseudopotentials that were utilized in this study in the relaxation procedures.

DISCUSSION

Applying the selection criteria developed in Section III has allowed us to shortlist a set of materials that may potentially exhibit a significant electrostrictive strain response. values estimated from the MP data using Equations ( 3) and ( 4) are compared with the experimental data taken from Refs. [START_REF] Li | Electrostrictive effect in ferroelectrics: An alternative approach to improve piezoelectricity[END_REF][START_REF] Uchino | Electrostrictive effect in perovskites and its transducer applications[END_REF][START_REF] Newnham | Electrostriction: Nonlinear electromechanical coupling in solid dielectrics[END_REF][START_REF] Uchino | Electrostriction and Its Interrelation with Other Anharmonic Properties of Materials[END_REF] in Fig. , respectively). We also plot a comparison of with and with in Fig. 3(a) and Fig. 3(b), respectively. There are significant deviations from the linear correlation that follows from the Newnham semi-empirical correlation between and . The result ascertains that Newnham's semi-empirical relation, is not a universal measure that is applicable across all materials.

The most striking aspect of the top compounds is that they are halides, chlorides, bromides, and iodides. We anticipated this due to recent observations of high electrostrictive properties in fluorides such as CaF 2 , BaF 2 , SrF 2 , LiF, and NaF with electrostrictive responses around 0.36 -0.46 m 4 C -2 3-5,47 . It is interesting to note that none of these compounds made our top ten obtained either via data mining of MP or from further refinement of MP results through high-precision DFT. The materials that are discovered through our approach are mainly binary non-cubic semi-metal halides. We note that three materials appear common to the top ten To explain these results and to identify the origin of the giant electrostrictive properties of materials displayed in Tables 1 and2, we carried out a detailed analysis of bonding, charge distribution, and hybridization for selected compounds. For example, we have examined GaBr

where is the breakdown voltage. This quantity essentially determines the maximum amount of strain that can be obtained through an applied electric field which is also true for electrostrictors.

For stress sensing, the piezoelectric charge coefficient (d) relates the applied stress ( ) to the dielectric displacement generated ( ) in short-circuit conditions whereas the piezoelectric voltage coefficient (g) relates the applied stress to the generated electric field ( ).

There is no direct way of computing . We used the breakdown voltage expression from Refs. [START_REF] Kim | Machine Learning Assisted Predictions of Intrinsic Dielectric Breakdown Strength of ABX3 Perovskites[END_REF][START_REF] Kim | From Organized High-Throughput Data to Phenomenological Theory using Machine Learning: The Example of Dielectric Breakdown[END_REF] :

(

where is the maximum phonon frequency at the -point that were calculated from the linear response approach and is the DFT obtained band gap. Although the band gap is underestimated in DFT, we do not expect it to affect the results of because Equation ( 10) is a machine-learned relationship that already accounts for the underestimated band gaps.

Table 3 lists the top ten materials for , , and . The magnitudes of for these materials are smaller than what generally is accepted as promising in the field of linear dielectrics, i.e., few hundreds of m V -1 . The values of , on the other hand, are three orders of magnitude larger than those commonly appear in linear dielectrics applications. This complements the extremely large effective piezoelectric coefficients obtained on ceria (d 33 ~ pC N -1 ), a known 'giant' electrostrictor, and underlines the potential of electrostrictors as an alternative to piezoelectrics. As for the top candidates, we find 70% of the materials common in the top ten list of and . The data signifies further that fluorides, chlorides, bromides, and iodides are promising material systems for both linear and quadratic dielectric response.

We have used a data mining backed by first principles calculations approach to down select materials for the discovery of novel electrostrictive materials. Our analysis shows that the halidesfluorides, chlorides, bromides, and iodides, are promising materials systems for high electrostriction response. However, unlike linear dielectrics and ferroelectric glasses, they do not follow Newnham's relation. We note that some of these materials may have significantly larger electrostrictive coefficients ( ) than even electroactive polymers. Discovery of electrostrictors with electrostrictive coefficients comparable to the 'giant' electrostrictors, where the physics is governed by active defects, presents the prospect of applications that are limited to lowfrequency response as is the case for 'giant' electrostrictors. The corresponding effective piezoelectric voltage coefficients obtained under bias electric field surpass classical piezoelectrics by three orders of magnitude and constitute an appealing alternative to piezoelectrics. We believe that this combination of computational and informatics-based approach will provide the necessary guidance for future experimental work, resulting in a new generation of electrostrictors with exceptionally high property coefficients.

METHODS

The Materials Project database was used for sourcing, compiling, and curating the crystal structure data and the existing elastic compliance and permittivity. The electrostriction coefficient was computed based on Newnham's proxy. It provided early-stage learning for shortlisting a list of compounds for carrying out more detailed first-principles calculations. A set of top candidates was selected for first-principles calculations. A set of randomly selected compounds was also included in the list. Detailed density functional theory (DFT) calculations using the Vienna Ab initio Simulation Package (VASP) were carried out with generalized gradient approximation (GGA) as the exchange-correlation functional and a high-density kpoints mesh. The electrostatic coefficient was calculated as a differential of pressure, where DFT calculations were applied with pressure constraint with the pstress tag. 7), Eq. ( 8), and Eq (9), respectively. The unit of , , and are MV/m, m/V and, Vm/N, respectively. The data for of top-ten materials are tabulated. All calculations were carried out at 0K and therefore on materials in their solid states.

Material

4 C 2 .

 42 3 , one of the top candidates with = 560.53 m 4 C -2 and = m 2 V -2 . Its oxide counterpart Ga 2 O 3 has significantly lower electrostriction response factor, = 0.02 m As such, Ga 2 O 3 did not make it to our lists even within the top 50 materials. The crystal structure and the atomic charge densities of GaBr 3 and Ga 2 O 3 are shown in Fig. 4 (a) and Fig. 4 (b), respectively. The atom resolved density of states (DOS) of the respective materials are plotted in Fig. 4 (c) and Fig. 4 (d). Despite each having 24 electrons per formula units in their valence band, the DOS is quite different. The bandgap observed from the DOS of the two compounds are 3.36 eV for GaBr 3 and 2.04 eV for Ga 2 O 3 , which is greater than 0.5 eV, a criterion used for our data mining. This corroborates with the charge separation observed between the cation and anion sites shown in Fig. 4 (a) and Fig. 4 (b) for GaBr 3 and Ga 2 O 3 , respectively. Owing to the large band gaps, the two materials have similar order of magnitude permittivities ( 2.86 for GaBr 3 and 9.28 for Ga 2 O 3 ). The overall energy spread of both the valence and conduction bands of GaBr 3 is smaller than for Ga 2 O 3 . In addition,

Figure 1 .

 1 Figure 1. Flow chart showing data curating and first-principles computation on datamined compounds. Data-mining Materials Project database and DFT computation workflow for search of high electrostriction materials.

Figure 3 .

 3 Figure 3. Comparison of electrostrictive coefficients of and of obtained from proxy and DFT computation. Comparison of (a) versus and (b) versus obtained from direct application of empirical formula to the Materials Project data and our DFT computations.

Figure 4 .

 4 Figure 4. Electronic structure of selected compounds. The crystal structure and the charge densities of GaBr 3 and Ga 2 O 3 are shown in (a) and (b). The atom projected electronic density of states for GaBr 3 and Ga 2 O 3 are shown in (c) and (d).

Figure 5 .

 5 Figure 5. Plot of permittivity versus compressibility obtained from DFT calculations. Electrostrictive coefficient as in Equation (3) falls along the quadrants I and III with larger values appearing is quadrant-III and small values in quadrant-I.

Figure 1 .

 1 Figure 1.

Figure 2 .

 2 Figure 2.

Figure 3 .

 3 Figure 3.

Figure 4 .

 4 Figure 4.

Figure 5 .

 5 Figure 5.

  ten high materials appear in both proxy and DFT data sets. These are: SiCl 4 , BCl 3 , PCl 3 , GaBr 3 , and WCl 6 . On the other hand, three out of the top ten materials appear in both the proxy and DFT data sets. These are SbSI, SbCl 2 , and WCl 6 . Materials in the top ten list of

	and		include MgCl 2 , MgI 2 , SiCl 4 , WCl 6 , and PCl 3 . Likewise, materials common to
	the top ten list of	and	are WCl 6 , GaBr 3 , and SbCl 3 . Overall, the data supports that the
	halide groups, especially the chlorides, the bromides, and the iodides are promising materials for
	electrostriction functionality.
	There are important differences in the output value of	that should be noted. For
	example, our	values for WCl 6 and PBr 3 are 22.83 m 4 C -2 and 20.99 m 4 C -2 , respectively,
	while	for these compounds is significantly larger: 560.53 m 4 C -2 for WCl 6 and 393.32 m 4
				2. It is observed that specialized polymers have generally larger
	than ceramics or glasses. Comparing	of materials from the MP dataset, one finds larger	in
	several ceramics and covalent salts that surpass polymers by almost 400% if the data from Ref. 3
	corresponding to polyurethane elastomer is considered an outlier. The top electrostrictive
	materials from MP are explored in more detail using more accurate DFT computations. The
	values based on Equation (6) display even more promising results with electrostrictive
	coefficients over 500% larger than typical electrostrictive polymers.
	The top ten	(	and	) and	(	and	) candidates from MP and

DFT calculations are listed in Tables

1 and 2

, respectively. It is encouraging to find that five out of our top C -2 for PBr 3 . Other materials such as GaCl 3 have similar values from the empirical and DFT relations (20.97 m 4 C -2 and 17.78 m 4 C -2

  GaBr 3 , WCl 6 , and SbCl 3 . These compounds have not yet been explored with respect to their electrostrictive properties. We certainly hope that our results will generate more experimental work to verify our predictions. We note that not all the compounds in Table1are viable materials, some are liquid at room temperature (e.g., SiCl 4 , PCl 3 , PBr 3 O), some are gases or volatile at room pressure (e.g., BCl 3 , SF 6 , or WCl 6 ). But others such as GaBr 3 with of 560.53 m 4 C -2 show significant promise.

Trujillo et al., Discovery of Novel Electrostrictors, page 11 and top ten datasets:

Table 1 .

 1 The list of top ten materials predicted from Newnham's proxy relation ( ) and our DFT calculations ( ). The unit of is m 4 /C 2 . All calculations were carried out at 0K and therefore on materials in their solid states.

	Material

Table 2 .

 2 The list of top ten materials predicted from and our DFT calculations ( ). The unit of is m 2 /V 2 . All calculations were carried out at 0K and therefore on materials in their solid states.

	Material

Trujillo et al., Discovery of Novel Electrostrictors, page 27

Table 3 .

 3 Values of dielectric breakdown, , piezoelectric coefficients , and computed from calculated from Eq. (

ACKNOWLEDGEMENTS

Authors AG, SKN, and SPA thank the Air Force Research Laboratory, Materials and Manufacturing Directorate (AFRL/RXMS) for support via Contract No. FA8650-18-C-5700. JY and PEJ have been supported by the ANR-19-ASTR-0024-01 and ANR-20-CE08-0012-1 grants.

DATA AVAILABILITY

The data supporting the results are available upon reasonable request to the authors.

within these bands, the DOS is composed of narrow peaks for GaBr 3 whereas Ga 2 O 3 exhibits a more continuous dispersion of its DOS.

The consequence of DOS differences is fundamentally connected with difference in elastic properties and thus on the compressibility , which for GaBr 3 is 22.67 GPa -1 . This value is three orders of magnitude larger than of Ga 2 O 3 (0.06 GPa -1 ). The plot of versus is shown to gain an insight into which physical factors play more important role in the derived properties. Fig. ( 5) is a four-dimensional plot with and as the absicissa and ordinate. The color coding differentiates the crystal system of the compounds, and the circle diameter scales according to the value of . Note that the values of and are independently computed using DFT. The figure is divided into four quadrants. Quadrant-II covers the region of large but small , and has large . Conversely, the quadrant-IV covers the region of small but large , and has smaller . These seem consistent with the semi-empirical relation of Equation ( 1).

However, we observe that quadrant-III is filled with data points showing that there are several materials with low and low that make up the chemical space. On the other hand, quadrant-I is almost empty. The result does not reflect the same trend as shown in Refs. [START_REF] Newnham | Electrostriction: Nonlinear electromechanical coupling in solid dielectrics[END_REF][START_REF] Damjanovic | Electrostrictive and Piezoelectric Materials for Actuator Applications[END_REF] , indicating that electrostrictive properties cannot be generalized using Newnham's proxy defined through Equations ( 3) and ( 4). The and derived from the Materials Project data are shown in Supplementary Figure 1 for all the datamined compounds and is provided in the Supplementary Information.

As a cross-check of the quality of data obtained from our study, we calculated the effective piezoelectric coefficients and from the derived and . The electrostriction coefficients are related to piezoelectric charge (d) and voltage (g) coefficients as:

AUTHOR CONTRIBUTION

The concepts presented herein were conceived by PEJ and SPA. DT performed the data mining process and provided the initial analysis. AG carried out DFT calculations under the supervision of SKN. All authors analyzed the data and prepared the manuscript for publication. DT and AG contributed equally for this work.

COMPETING INTERESTS

The Authors declare no Competing Financial or Non-Financial Interests


Strain-engineered divergent electrostriction in KTaO 3 Daniel S. P. Tanner, 1, 2, * Pierre-Eymeric Janolin, 2 and Eric Bousquet 1 1 Université de Liège, Q-MAT, CESAM, Institut de Physique 2 Université Paris-Saclay, CentraleSupélec, CNRS, Laboratoire SPMS, 91190 Gif-sur-Yvette, France (Dated: July 29, 2022)

We investigate the electrostrictive response across a ferroelectric phase transition from firstprinciples calculations and show that M , the field-induced electrostrictive tensor, controlling the amplitude of the electric-field induced strain, can be made arbitrarily large through strain engineering. We take as a case study the epitaxial strain-induced transition from para-to ferro-oelectricity of KTaO3. We show that the magnitude of the field-induced electrostriction diverges with the permitivity at the transition, hence exhibiting giant responses through a calculation of both the M and Q electrostrictive tensors. We explain the origin of this giant electrostrictive response in KTaO3 using a microscopic decomposition of the electrostriction coefficients, and use this understanding to propose design rules for the development of future giant electrostrictors for electromechanical applications. Finally, we introduce a further means to calculate electrostriction, specific to ferroelectrics, and not yet utilised in the literature.

Electrostriction is an electromechanical coupling present in all dielectrics which describes the quadratic induction of a strain or stress, due to an electric or polarisation field. Despite its ubiquity, it is often dwarfed in amplitude by its linear counterpart, piezoelectricity, which constitutes the primary response in most electromechanical systems. However, renewed interest in electrostriction has grown due to the recent discovery of materials which exhibit "giant" electrostrictive responses, [1][2][3] capable of competing with piezoelectrics in terms of responses. These electrostrictors exhibit a host of attractive advantages over piezoelectrics such as: (i) minimal hysteresis in the strain-field response, (ii) increased stability of induced deformations, (iii) smaller response time, (iv) enhanced temperature stability, and (v) materials which are lead-free and environmentally friendly. [4][5][6] In ferroelectric materials, there is a consensus that the large electromechanical response of ferroelectric materials is mostly due to their piezoelectric rather than electrostrictive response. Additionally, it is often accepted from experimental studies [5][6][7] and phenomenological treatment [START_REF] Strukov | Ferroelectric Phenomena in Crystals[END_REF] that the electrostriction remains relatively small and unchanged across ferroelectric phase transitions. However, these studies report the electrostrictive Q ijkl coefficient (relating strain to the polarisation squared) while a few reports [START_REF] Isupov | [END_REF]10] found that the electrostrictive M ijkl tensor (relating strain to the electric field squared) increases rapidly with temperature together with the permittivity near phase transitions. However, it is believed that the M ijkl response is nevertheless negligibly small with respect to the piezoelectric response, and therefore useless in terms of technological applications [START_REF] Shuvalov | Modern Crystallography IV Physical: Properties of Crystals[END_REF]. Hence, the question remains of whether inducing phase transitions where the permittivity diverges enables to increase the electrostrictive response significantly enough to make electrostriction in-teresting for electromechanical applications.

In this Letter, through a thorough examination of electrostriction at the strain-induced ferroelectric phase transition (FEPT) in KTaO 3 from first-principles calculations, we demonstrate a giant electrostrictive response can indeed be obtained in ferroelectric materials at the phase transition, confirming that, though the components of the polarisation tensor Q do not change much, the out-of-plane electric-field tensor components diverge and reaches values that make electrostriction interesting for electromechanical applications. In addition, through a microscopic decomposition of the electrostriction tensors, we show that it is the soft transverse optical mode that is responsible for the large electrostriction.

To proceed, we note that while one may consider electrostriction as the "strain/stress induced quadratically by an electric/polarisation field", one may also view it as the theromodynamically equivalent "linear change in susceptibility/inverse susceptibility induced by a stress or a strain" [START_REF] Devonshire | [END_REF]13] shown below:

1 0 ∂η ij ∂X kl = -2Q ijkl 1 0 ∂η ij 0 ∂x kl = 2q ijkl 0 ∂χ ij ∂X kl = 2M ijkl 0 ∂χ ij ∂x kl = -2m ijkl . (1)
See the supplemental material for details of the theromodynamical derivation of these equations.

[14] In our previous work [13], we outlined the advantages of this method of calculation, amongst which are greater computational efficiency and robustness. Viewed from the perspective of Eqs.(1), we see that at a strain-or stress-induced phase transition, involving the divergence of the susceptibility, some electrostrictive tensor components must become large. [START_REF] Jona | Ferroelectric Cristals[END_REF] To compute Eqs. (1) we have used the DFT package ABINIT (version 8.6.1) [START_REF] Gonze | [END_REF], with k-point grid densities of 8×8×8 and a plane wave cutoff energy of 75 Ha, to ensure convergence of electrostrictive coefficients of about 1%. The PseudoDojo [17] normconserving pseudopotentials were used, and the exchange-correlation functional was treated using the generalised gradient approximation of Purdue, Burke, and Ernzerhof, modified for solids (PBEsol). [18] We note that obtaining the electrostriction tensors in this manner using eq.1 and finite differences, will yield an "improper" electrostriction tensor containing contributions described by Nelson and Lax [19]; the "proper" tensor may be obtained from this by the addition of permitivity components to certain electrostriction components. [13,19] As a test case, we will analyse the electrostrictive response at the strain-induced FEPT present in KTaO 3 . KTaO 3 is a highly suitable material for such a study as it is an incipient ferroelectric which undergoes a ferroelectric phase transition due to relatively small axial stress. A FEPT is reported at 0.6 GPa by Uwe [7] and for a strain of 1% by Tyunina [20] at low temperature (4.2 K), which can be studied experimentally on single crystals or epitaxial thin films. In addition, KTaO 3 has a wide range of applications [20][21][22][23][24], including the recent discoveries of superconductivity [22,23] or dynamic multiferroism [24], making it a highly technologically relevant material.

We present in Table I our calculated ground state properties of KTaO 3 along with previous PBEsol calculations from the literature, as well as experimental measurements. We obtain a bandgap in good agreement with previous PBEsol theory, but it is underestimated with respect to the experimental values as expected from the PBEsol functional, though we can reproduce the indirect gap from R point to Γ-point. For the lattice and elastic constants good agreement is obtained with both experiment and previous theory. While no previous experimental or theoretical works have determined the M ijkl electrostrictive coefficients, for the Q ijkl coefficients, we obtain reasonable agreement with previous experiment [6,7], with the disagreement between our computed values and the experimental ones being equivalent to the variances amongst the experimental values.

To investigate electrostriction through the FEPT, we have calculated the properties of KTaO 3 over a range of biaxial compressive and tensile strains from -0.7% to 0.7%, as may be achieved through epitaxial growth on various substrates. This was done by imposing the strain components x 11 = x 22 = (a sa 0 )/a 0 where a 0 is the relaxed lattice constant of KTaO 3 and a s is the lattice constant of the substrate, and by relaxing the z direction of strain.

In the top part of Fig. 1 we show the impact of these various epitaxial strains on the permitivity tensor. Consistent with previous works on KTaO 3 , the epitaxial strain induces a FEPT [7,20], and we see it is in fact the tensile strains which drive the polar displacement. On the left, the middle panel reveals that the tensile strain induced in the z-direction by the compressive epitaxial xy

433 Q11 (m 4 C -2 ) 0.1256 - 0.087 a Q12 (m 4 C -2 ) -0.0313 - -0.023 a Q h (m 4 C -2 ) 0.063 - 0.041 a ,0.052 b M11 (m 2 V -2 ) 1.98x10 -18 - - M12 (m 2 V -2 ) -4.70x10 -19 - - M h (m 2 V -2 ) 1.05x10 -18
--TABLE I. Bulk paraelectric ground state KTaO3 electronic, structural, elastic, and electrostrictive properties obtained via calculation and measurement. a=Ref [7]; b=Ref [6]; c=Ref [25]; d=Ref [26;e=Ref[27]; f=Ref [28]; g=Ref [START_REF] Weber | Handbook of optical materials, 1 st edition[END_REF]; h=Ref [START_REF] Cabuk | [END_REF]; i=Ref [31] (soft mode at 300K, others at 10K); j=Ref [32] (RT). *In h the LO2 mode is mislabled as the TO3 mode due to the fact that LO2 is higher in frequency than the non polar TO3 and LO3 modes.

strain leads to the softening of the transverse optic (TO) mode polar in z, and a divergence of 33 , with the FEPT occurring at 0.4% epitaxial xy strain; on the right, the epitaxial tensile strain leads to the simultaneous softening of the polar TO modes in x and y, and a divergence of 11 and 22 , with the transition occurring at around 0.26% epitaxial strain. This earlier transition point is due to the fact that in this case the tensile strain is being directly imposed on the x and y coordinates, whereas with the compressive strain, the tensile strain in the zdirection occurs more slowly due to a Poisson effect. This occurrence of the FEPT at the rather small compressive and tensile strains of 0.4 and 0.26% is consistent with the incipient nature of the FE in KTaO 3 .

The bottom panel of Fig. 1 shows that the electrostrictive tensor components M 33 and M 13 also diverge at the FEPT, contrary to reports in the literature that electrostriction is constant at a phase transition [33]. To compute the electrostriction, at each fixed epitaxial strain, we apply a further strain in the z-direction, of magnitude of -0.005, to the already relaxed c lattice vector, and calculate again the permitivity using DFPT. The coefficient m 33 is obtained from the derivative of z with respect to the applied x 33 ; and the rate of change of x w.r.t. x 33 yields the coefficient m 13 . The technologically relevant M 33 and M 13 , are then derived by appropriate consideration of the epitaxial boundary con- x 33 (%)

ǫ 33 ǫ 11 ω z (cm -1
)

ω x (cm -1
)

ω z ω x M 33 (nm 2 V -2 ) -M 13 (nm 2 V -2
)

x 11 = x 22 (%)
where s ij are elastic susceptibilities, for which we have used the equilibrium values, noting that epitaxial strains of < ±1 are unlikely to induce significant elastic nonlinearities [34]. Focusing first on the compressive epitaxial strain induced transition (at -0.4%), we observe spikes in both M 33 and M 13 ; with M 33 reaching 1x10 -15 m 4 /C 2 , a value rarely observed except for some polymers, and M 13 reaching the smaller, though still considerable value of -1.5x10 -16 m 4 /C 2 . This means that in a compressively strained epitaxial layer, an electric field will produce a significant tensile strain in the z-direction, and also a compressive strain parallel to the plane of the layer. To put this in context, the peak M 33 is equivalent to an effective piezoelectric coefficient of d ef f 33 = 60, 000 pm/V, which compares very favourably with d 33 of up to 600 pm/V [35] values in classic piezoelectric material PZT. This large tensile z-strain could have applications in micro-fluidics, adjustable mirrors etc., where the leadfree nature of KTaO 3 makes it particularly attractive when compared with other lead-containing materials.

Looking at the tensile-strain-induced transition (at 0.26%), we find again a divergence of the coefficients, though this time both reaching much higher values, with M 13 reaching -5x10 -15 m 4 /C 2 , and M 33 reaching 2x10 -15 m 4 /C 2 . In this case an electric field in z will produce considerably more strain in-plane/xy than in the case of compressive epitaxial strain.

Figure 1 thus demonstrates that electrostrition in leadfree KTaO 3 , brought near the phase transition through epitaxial growth on an appropriate substrate, may be an ideal active material candidate for future electromechanical devices. In addition, we note that the values obtained of the order of 10 -15 m 4 /C 2 , compare very favorably with known giant electrostrictors, such as (Nb,Y)-doped Bi 2 O 3 (M≈10 -17 m 4 /C 2 ) [36], Gd-doped ceria (M≈10 -17 to 10 -16 m 4 /C 2 ) [37] and La 2 Mo 3 O 9 (M=1.5x10 -18 m 4 /C 2 ) [START_REF] Yu | Giant electrostriction[END_REF] without being restricted to low frequencies.

To better understand this giant electrostriction, we can profit from the particular manner in which the methods based on eqs.1, and introduced in Ref. 13, determine the electrostrictive coefficients. Through the decomposition of the permitivity tensor into its mode by mode contributions, we are also able to decompose the electrostriction tensor. This decomposition is illustrated for the x , z and the electrostriction tensor components m 33 and m 13 in Fig. 2. The figure shows that all four TO modes contribute to 11 and m 13 , whilst only three contribute to 33 and m 33 . However, amongst all of these contributions, the only that matter near the transitions are those of the soft TO1 modes, with all others being negligible and indistinguishable from zero in the plot. We note also that in Fig. 2, the behaviour of the m ij is much simpler than that of the corresponding M ij shown in the bottom panel of Fig. 1, with m 33 and m 13 each exhibiting only one peak; however, we see that it is the mixing of these coefficients, according to Eq.2, which leads to the double peaked M ij .

To clarify why some previous experimental observations reported that the Q electrostrictive response is constant across FEPT, [6] we have also examined the evolution of the Q 33 coefficient of our strained KTaO 3 . Figure 3 shows the development of η 33 and Q 33 at the compressive strain induced phase transition. We plot η 33 as the Q 33 is obtained by its derivative w.r.t. the compressive z-strain. Contrarily to the permitivity and M 33 , a regular, non-divergent behaviour for both η 33 and Q 33 . In fact, when approaching the FEPT at -0.4% epitaxial strain, the coefficient Q 33 evinces a constant behaviour, and whilst the coefficient decreases as the compressive strain increases beyond the transition, it is a very small change, i.e. < 3% over the considered strain range. x 11 = x 22 (%)

η 33 Q 33 FIG.
3. η33 and Q33 plotted against epitaxial strain. Q33 is obtained by applying a further z-strain of -0.05% at each epitaxial strain and computing the rate of change of 33 with respect to stress. See text.

It might surprising that the M 33 electrostrictive coefficient diverges at FEPT while the Q 33 is almost unaffected but this makes sense when we consider eqs. (1), as η = 1 /χ where we have that Q = -1 /χ 2∂χ /∂X. Hence, the diverging susceptibility (i.e. large ∂χ /∂X) is compensated by the 1 /χ 2 part such that they balance each other out to give a nearly constant Q coefficient at the FEPT, while the M coefficient is directly proportional to ∂χ /∂X and increases with it. Therefore, whether or not elec-trostriction is claimed to diverge and become huge at the FEPT [START_REF] Isupov | [END_REF] or not [6,7] depends on whether one analyses the M and m coefficients (which must diverge), or the Q and q coefficients (which do not) for continuous transitions.

Finally, we demonstrate that the ferroelectric phase offers a unique approach to the calculation of the electrostrictive coefficients, not yet implemented in DFT. This is based on the relation of the electrostriction tensor to the spontaneous polarisation and piezoelectricity tensor: [START_REF] Devonshire | [END_REF][39]

g ijm = 2Q ijmn P n , (3)
where g ijm is the piezoelectric tensor which describes the linear coupling of strains to polarisations. This means that in the range of strain where KTaO 3 is ferroelectric (below -0.4%), we can determine the Q coefficient through the spontaneous polarization (calculated through the Berry phase method) and the piezoelectric coefficient (directly accessible from DFPT). Just below the critical strain for FEPT (-0.4%), we obtained a spontaneous polarisation of P 3 = -1.57 µC/cm 2 and a piezoelectricity tensor components g 33 and g 31 of -4.015 and 1.035 10 -3 m 2 /C, respectively. Using eq. ( 3), this yields to Q 11 = 0.1278 m 4 C -2 and Q 12 = -0.03295 m 4 C -2 , in good agreement with the Q values of 0.1278 and -0.02986 m 4 C -2 , evaluated at the same point using eqs.1. This method to obtain the Q ijkl coefficients has advantages over both the methodology we have presented in Ref. 13, and the finite field methodologies found in the literature [START_REF] Jiang | [END_REF][41][42], the primary amongst which is that no finite differences or relaxations were required, just the ground-state based Berry phase and DFPT calculations. It also circumvents issues associated with applying an electric field to the paraelectric phase of a perovskite [START_REF] Jiang | [END_REF], and more general shortcomings associated with finite field methodologies [13].

In conclusion, we have performed a thorough investigation of electrostriction at the ferroelectric phase transition. We have outlined a route towards giant electrostriction through epitaxial strain engineering, demonstrating an electrostrictive M 33 coefficient equivalent to an effective d 33 coefficient of 60, 000 pm/V. The magnitude of this response, the ease with which the studied system may be grown, and the advantages of lead-free electrostrictive materials over piezoelectric materials more generally mark this as a significant scientific and technological advancement. In addition, through our treatment of both the Q ijkl and M ijkl tensors, we underline that whether one will observe a diverging and giant electrostriction, or a relatively constant and stable response, depends on whether one measures the Q or the M tensor. Finally, we have introduced a new method by which electrostriction may be computed in a ferroelectric using only the Berry phase technique and DFPT, without the need for relaxations or finite differences. Our results also question if the origin of giant electrostriction [43] e.g. in LAMOX, doped ceria, or Bi 2 O 3 , is coming from an intrinsic softening of a polar mode, like in KTaO 3 or if its due to intrinsic electronic effects or other extrinsic effects. This will require further experimental and theoretical researches to further understand giant electrostrictors.

13 FIG. 1 . 33 M

 13133 FIG.1. Variation of with epitaxial strain (x11 = x22; X33 = 0) of: permitivity tensor components 11 and 33 (top); lowest x and z polar phonon modes, ωx and ωz (mid); and electrostriction tensor components M33 and M13 (bottom). Compressive epitaxial strain induces a tensile strain induced phase transition in the z-direction, and tensile epitaxial strain induces a phase transition in the x=y direction.

13 FIG. 2 . 33 Q 33 (m 4 /C 2 )

 132333342 FIG.2. Decomposition of (top) permitivity tensor components 33 and 11, and (bottom) electrostrictive tensor components m33 and m13 into their mode by mode and electronic contributions. For both the permitivity and electrostriction, the primary contribution is from the softening TO mode.

ACKNOWLEDGEMENTS

Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (C ÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region and using the DECI resource BEM based in Poland at Wroc law with support from the PRACE OFFSPRING project. EB acknowledges FNRS for support and DT aknowledge ULiege Euraxess support. This work was also performed using HPC resources from the "Mésocentre" computing centre of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/). Financial support is acknowledged from public grants overseen by the French National Research Agency (ANR) in the ANR-20-CE08-0012-1 project and as part of the ASTRID program (ANR-19-AST-0024-02). This article is based upon work from COST Action OPERA -EurOPean NEtwork foR Innovative and Advanced Epitaxy, supported by COST (European Cooperation in Science and Technology.


Daniel Bennett

Daniel Tanner

Philippe Ghosez

Pierre-Eymeric Janolin

Eric Bousquet

Generalized relation between electromechanical responses at fixed voltage and fixed electric field

We present a general relation between the electromechanical couplings of infinitesimal strain and electric field to arbitrary order, measured at fixed voltage and at fixed electric field. We show that the improper response at fixed field can be written as the strain derivative of the n th order susceptibility tensor, and the proper response at fixed voltage drop can be written as the response at fixed field plus corrections for dilations and 90 • rotations induced by strain. Our theory correctly reproduces the proper piezoelectric response and we go beyond with the electrostrictive response. We present first-principles calculations of the improper electrostrictive response at fixed field, and illustrate how the correction is used to obtain the proper response at fixed voltage. This distinction is of high importance given the recent interest in giant electrostrictors exhibiting electromechanical responses as large as the piezoelectric ones.

INTRODUCTION

The calculation of the piezoelectric response from firstprinciples density functional theory (DFT) simulations is nowadays standard, either using finite differences [1] or density functional perturbation theory (DFPT) [2]. However, care should be taken, as two possible piezoelectric responses can be obtained: the so-called "proper" and "improper" responses [1][2][3][4]. The improper response is computed from DFT when calculations are performed at fixed electric field, whereas the proper response is typically measured experimentally, where the voltage and not the electric field is held fixed [1,5]. The distinction between fixed field and fixed voltage drop becomes necessary when strains are introduced.

This problem is sketched in Fig. 1: we have a dielectric material of thickness d sandwiched between two metal plates, across which a potential drop ∆V is applied. The voltage drop is related to the electric field in the material via

∆V = Ed . (1)
If a small strain η appears in the material, changing the thickness to (1 + η)d, Eq. (1) becomes:

∆V (η) = E(1 + η)d . (2)
If the field E is held fixed, which is typically the case in firstprinciples calculations, the potential difference has to change by η∆V . Thus when strain is introduced we need to compensate for this effect in order to hold the potential difference fixed. If the potential drop ∆V is fixed, then it is the reduced field E ′ = (1 + η)E rather than E which is held fixed in Eq. 2.

The voltage drop in each calculation will then be

∆V = (1 + η) -1 E ′ • (1 + η)d = E ′ d (3)
which is independent of η.

In addition to the strain dependence of the relationship between E and ∆V (Eq. 2), changes in the volume induced by * dbennett@uliege.be d

∆V (1 + η)d ∆V (1 + η)∆V FIG.
1. Sketch of a dielectric material of thickness d sandwiched between two metal plates. A potential difference ∆V is applied, resulting in an electric field of E = ∆V d inside the material. When a strain of η is applied to the material, the field changes to E = ∆V (1+η)d if the potential difference is held fixed (experiment), or the potential difference changes to (1 + η)∆V if the field is held fixed (first-principles).

strain can lead to different ways of measuring the energy per unit volume. Depending on whether calculations are done at fixed E or ∆V , and whether the energy per unit volume is measured with respect to the equilibrium reference volume (Ω 0 ) or the perturbed volume (Ω), there are four possible ways to define an electromechanical response, see Fig. 2: (i) the improper response, measured at fixed E and equilibrium volume Ω 0 , (ii) the proper response, measured at fixed ∆V and strained volume Ω, (iii) an intermediate situation with fixed E and Ω, and (iv) a second intermediate situation with fixed ∆V and Ω 0 are used. This means that great care should be taken when comparing first-principles calculations of any electromechanical response to experimental measurements.

Although this problem is well-known in the context of piezoelectricity, its generalization to other electromechanical responses measured at fixed voltage and fixed electric field has not been reported. This generalization is becoming relevant, with increasing interest in higher-order couplings such as the electro-optic response and electrostriction. Indeed, giant electrostrictors have recently been shown to give rise to electromechanical responses as large as piezoelectric ones [6,7], or even much larger: up to 2 × 10 5 pC/N [8]. However, no strict definition nor derivation of the proper versus improper responses have been given for electrostriction. In this paper, we derive an analogous relation for electrostriction and more generally any electromechanical response to an infinitesimal strain.

This paper is organized as follows: first, we review the relation between piezoelectric tensors measured at fixed electric field and fixed potential difference. The same methodology is then applied to obtain a corresponding relation for the electrostrictive tensor, and a general electromechanical coupling between an infinitesimal strain and electric field to arbitrary order. Detailed derivations are available in the Supplementary Material. We then use these relations to correct recent calculations of the electrostrictive response of rocksalt crystals in the literature [9], focusing on MgO, and re-examine the comparison with experimental measurements.

THEORY

For a crystal with mechanical degrees of freedom such as phonon displacements u and macroscopic strain deformations η, which is subject to an applied electric field E [10][11][12], the properties close to the equilibrium state can be described by Toupin's electric enthalpy density [13,14]:

F Ω 0 = F cell -E • P + 1 2 ε 0 E 2 , (4)
where Ω 0 is the volume of the reference state ( η = 0), F cell is the zero-field Kohn-Sham energy per unit volume of the system, and P is the total polarization of the system [15,16]. An homogeneous deformation results in a change in the positions r in the reference state: ∆r α = ηαβ r β . The symmetric part of η is simply the strain tensor, η. In Ref. [2] it was shown that the proper piezoelectric tensor is symmetric under the exchange of the strain indices. Thus, although we consider general deformations η in this paper, we use the terms strain and deformation interchangeably.

The various clamped responses (u = 0) are defined as the derivatives of Eq. ( 4) about the equilibrium state: { η, E} = 0. For example, the piezoelectric tensor e is defined as [2]:

e (i) αβ γ ≡ - ∂ 2 F Ω 0 ∂ E α ηβγ { η,E}=0 = ∂ P α ∂ ηβγ , (5)
where the superscript (i) denotes that this is the improper tensor, measured at fixed electric field. In order to obtain the proper tensor, at a fixed potential difference, two corrections must be made to Eq. ( 4). First, the electrostatic energy must be scaled by Ω Ω 0 , where Ω = Ω 0 det(I + η) is the volume after a strain deformation η. Secondly, in order to hold the potential difference fixed, we make the following change of variables to "reduced coordinates" for the electric field, polarization and displacement field D [17]:

E ′ = (I + η) E P ′ = Ω Ω 0 (I + η) -1 P D ′ = Ω Ω 0 (I + η) -1 D , (6)
which are equivalent to some of the reduced fields defined in Ref. [5], but here the units of each field are preserved [START_REF]E ′ is equivalent to ε = a T E ′ but preserving the units of electric field, where the columns of a are the lattice vectors (including strain). P ′ is equivalent to p = Ωb T P, where the columns of b are the reciprocal lattice vectors[END_REF]. Thus, the enthalpy Eq. ( 4) becomes, in reduced coordinates:

F ′ Ω = F cell -E ′ • P ′ + Ω Ω 0 1 2 ε 0 E ′T g -1 E ′ . ( 7
)
where g = (I + η) T (I + η). By differentiating we can see that

D ′ α = -∂ F ∂ E ′ α
, where we used D = ε 0 E + P, although note that

D ′ ̸ = ε 0 E ′ + P ′ .
The proper piezoelectric tensor, measured at fixed voltage drop, is defined as:

e αβ γ ≡ - ∂ 2 F ′ Ω ∂ E ′ α ∂ ηβγ { η,E ′ }=0 = ∂ P ′ α ∂ ηβγ . ( 8
)
The term linear in E is absent in Eq. ( 8) as the derivation is carried out about the equilibrium state ( η, E → 0). Summation convention is assumed, using Latin letters for dummy indices and Greek letters for free indices. Next, we explicitly evaluate the derivative in order to write Eq. ( 8) in terms of Eq. ( 5):

e αβ γ = ∂ det (I + η) (I + η) -1 α j P j ∂ ηβγ = ∂ (1 + ηii ) (δ α j -ηα j ) P j + O( η2 ) ∂ ηβγ = ∂ P α ∂ ηβγ + P α δ β γ -P γ δ αβ + ηii ∂ P α ∂ ηβγ -ηα j ∂ P j ∂ ηβγ . ( 9
)
Going from the first line to the second line, we used the expansion (I + η) -1 = I -η +O( η2 ), truncating to linear order, assuming the strain is infinitesimal. For the same reason, we use det(I + η) = 1 + tr ( η) + O( η2 ) to simplify the volume term. Note that it is important to do this before differentiating, otherwise contributions which are quadratic in strain will be retained. After differentiating, we set E, η → 0, because the derivatives are defined about the equilibrium state with zero strain and applied field. The relation between piezoelectric tensors measured at fixed voltage and at fixed field is then:

e αβ γ = e (i) αβ γ + P α δ β γ -P γ δ αβ . (10)
As discussed in Ref. 1, the first term is a correction for dilations of the polarization induced by strain. The second term is a correction for rotations of the polarization by 90 • , i.e. permutations of the indices. Additionally, a problem unique to the piezoelectric tensor is that Eq. ( 5) is sensitive to the branch on which the polarization is measured (quantum of polarization), but this is remedied by using Eq. (10), which is branch invariant [1].

Electrostriction

The improper strain electrostrictive tensor m, measured at fixed electric field, is defined as follows:

m (i) αβ γδ ≡ - 1 2 ∂ 3 F Ω 0 ∂ E α ∂ E β ∂ ηγδ { η,E}=0 = ∂ ε (i) αβ ∂ ηγδ , (11)
where

ε (i) αβ ≡ ∂ D α ∂ E β
is the improper permittivity tensor, measured at fixed electric field calculated about the equilibrium state. The proper electrostrictive tensor, measured at fixed voltage, is defined as follows:

m αβ γδ ≡ - 1 2 ∂ 3 F ′ Ω ∂ E ′ α ∂ E ′ β ∂ ηγδ { η,E ′ }=0 = ∂ ε αβ ∂ ηγδ , (12)
where

ε αβ ≡ ∂ D ′ α ∂ E ′ β
is the proper permittivity tensor, measured at fixed voltage. Immediately we see that Eqs. ( 11) and ( 12) are analogous to Eqs. ( 5) and ( 8), but with polarization replaced by dielectric permittivity. As before, we must write Eq. ( 12) in terms of Eq. (11). In order to do this, we must first relate the proper and improper permittivity tensors:

ε αβ = (1 + ηii )ε (i) αβ -ε (i) αi ηiβ -ηα j ε (i) jβ + O( η2 ) . (13)
Note that ε αβ = ε (i)

αβ at η = 0, which is expected. However, it is clear that their derivatives with respect to strain will not be equal at η = 0. Inserting Eq. ( 13) into Eq. ( 12):

m αβ γδ = ∂ ε αβ ∂ ηγδ = ∂ ∂ ηγδ (1 + ηii )ε (i) αβ -ε (i) αi ηiβ -ηα j ε (i) jβ + O( η2 ) = ∂ ε (i) αβ ∂ ηγδ + ε (i) αβ δ γδ -ε (i) αγ δ β δ -ε (i) δ β δ γα . (14)
Hence, we obtain an analogous relation between the proper and the improper response to Eq. ( 10) for the electrostriction:

m αβ γδ = m (i) αβ γδ + ε αβ δ γδ -ε αγ δ β δ -ε δ β δ γα , (15)
where we have dropped the (i) superscript on the permittivities because we have set η → 0 and hence ε (i) = ε. The corrections are similar to those appearing in the proper-improper relation for the piezoelectric tensors: the first term is a correction for changes in the permittivity induced by strain. The second term is a correction for permutations of the indices. Note that there are two of these terms, since ε is a rank 2 tensor.

Generalization to all orders of electromechanical coupling

The proper-improper relation can be generalized to electromechanical couplings which are linear in strain and to any order in electric field. The (improper) coupling between linear strain and n th order electric field is described by the following tensor:

A (i) E 1 ...E n η1 η2 ≡ - 1 n! ∂ n+1 F Ω 0 ∂ E E 1 . . . ∂ E E n ∂ η η1 η2 { η,E}=0 = ∂ χ (i) E 1 ...E n ∂ η η1 η2 , (16)
measured at fixed field, where n = 1 corresponds to the piezoelectric tensor and n = 2 corresponds to the electrostrictive tensor. We can see that in general A is the first strain derivative of the n th order susceptibility tensor:

χ (i) E 1 ...E n ≡              P E 1 , n = 1 ∂ D E 1 ∂ E E 2 , n = 2 ∂ n-1 P E 1 ∂ E E 2 ...∂ E En , n > 2 (17)
defined to be the polarization and permittivity for n = 1, 2, respectively. The proper tensor, measured at fixed voltage, is defined as:

A E 1 ...E n η1 η2 ≡ - 1 n! ∂ n F ′ Ω ∂ E ′ E 1 . . . ∂ E ′ E n ∂ η η1 η2 { η,E ′ }=0 = ∂ χ E 1 ...E n ∂ η η1 η2 , (18)
where

χ E 1 ...E n ≡                P ′ E 1 , n = 1 ∂ D ′ E 1 ∂ E ′ E 2 , n = 2 ∂ n-1 P ′ E 1 ∂ E ′ E 2 ...∂ E ′ En , n > 2 . ( 19
)
Again, the aim is to write Eq. ( 18) in terms of Eq. ( 16). In order to do this, we must first obtain the relation between the n th order susceptibility tensors:

χ E 1 ...E n = ∂ n-1 P ′ E 1 ∂ E ′ E 2 . . . ∂ E ′ E n = n ∏ i=2 (I + η) -1 k i E i ∂ n-1 P ′ E 1 ∂ E k 2 . . . ∂ E k n = n ∏ i=1 (I + η) -1 k i E i det(I + η)χ (i) k 1 ...k n (20)
Expanding det(I + η) and the products of (I + η) -1 , and truncating to linear order, we get the following expression:

χ E 1 ...E n = (1 + ηii )χ (i) E 1 ...E n - n ∑ i=1 χ (i) E 1 ...k i ...E n ηk i E i . (21)
Eq. ( 21) is the relation between n th order susceptibility tensors in reduced and unreduced coordinates at finite strain. It is reassuring to see that when η → 0 we have χ

E 1 ...E n = χ (i) E 1 .
..E n , as expected. Obtaining the proper-improper relation for A is now straightforward:

A E 1 ...E n η1 η2 = ∂ ∂ η η1 η2 (1 + ηii )χ (i) E 1 ...E n - n ∑ i=1 χ (i) E 1 ...k i ...E n ηk i E i = A (i) E 1 ...E n η1 η2 + δ η1 η2 χ E 1 ...E n - n ∑ i=1 χ E 1 ... η1 ...E n δ η2 E i ( 22
)
where η1 is in the i th position. There is one term for dilations induced by strain and n terms corresponding to permutations of the n indices of the susceptibility. For n = 1, 2, Eqs. (10) and ( 15) are reproduced, respectively.

First-principles calculations

We illustrate the proper-improper relation for electrostriction using first-principles calculations of the electrostrictive

ε ∞ ε (i) 11 ε (i) 22 ε 11 ε 22 ε 0 η 11 (%) FIG.
3. Improper (relative) dielectric permittivity of MgO as a function of applied strain η 11 (dashed lines) [9], and the (relative) proper permittivity obtained using Eq. ( 13) (solid lines), for (a): the electronic permittivity ε ∞ (clamped ions) and (b): the relaxed-ion permittivity ε 0 .

response of MgO, LiCl, LiF and NaCl. Typically, it is the electrostrictive response to stress which is measured experimentally, because it can be measured at fixed voltage [6]. The electrostrictive response to a stress is defined as

M αβ γδ = 1 2 ∂ ε αβ ∂ X γδ
, and is related to the electrostrictive response to strain used in the previous section by M αβ γδ = s αβ µν m µνγδ , where s is the compliance tensor, the inverse of the elastic tensor. In order to obtain the electrostrictive response to stress at fixed voltage, we first measure the electrostrictive response to strain m, apply the correction Eq. ( 15), and then contract with the compliance tensor to obtain M.

Fig. 3 shows the permittivity of MgO versus strain as obtained with the ABINIT code [22,23], following the methodology of Ref. [9]. Norm conserving pseudopotentials from PseudoDojo were used [24], and the PBEsol functional was used to treat exchange-correlation interactions [25]. A cutoff energy of 50 Ha was used to truncate the plane-wave basis, and a k-point grid of 8 × 8 × 8 was used to sample the Brillouin zone. The dielectric permittivity was then calculated using DFPT. The electronic or clamped-ion permittivity ε ∞ was first obtained, and then phonon calculations were performed [20] and c= Ref. [21].

in order to obtain the relaxed-ion permittivity ε 0 . The dashed lines show the components of the permittivity tensor obtained directly from ABINIT, which are in unreduced units and therefore improper. The components of the proper dielectric tensor, indicated by the solid lines, are obtained using Eq. ( 13). Note that for ε 22 , the proper and improper values are identical, which can be verified using Eq. ( 13). The improper electrostrictive response m (i) is obtained by taking the slope of ε (i) about η = 0, and the proper response m is obtained by taking the slope of ε about η = 0. The two intermediate values, obtained by either using reduced variables for the electric field or units of strained volume, were calculated using the corrections in Fig. 4 in order to illustrate their individual contributions to the total correction. The electrostrictive response to stress M was then obtained by contracting m with the compliance tensor in each case. In Table I we show a comparison of the electrostrictive responses for the proper, improper and intermediate cases, as well as experimental measurements.

DISCUSSION AND CONCLUSIONS

In this paper, we highlighted the difficulties in comparing electromechanical responses obtained from first-principles calculations to those obtained from experimental measurements. In the context of piezoelectricity, this problem is known as the 'proper-improper' relation, which may not be intuitive as 'proper' and 'improper' are ambiguous terms. We believe it is more physically insightful to reformulate the problem in terms of responses measured at fixed electric field and fixed voltage drop. Similar to the well-known relation for piezoelectricity, we have derived the relations for electrostrictive and general electromechanical responses to infinitesimal strain, measured at fixed field and at fixed voltage drop.

Depending on whether the field or voltage is held fixed, and whether the energy is given per unstrained or strained unit volume, there are four different ways to define an electromechanical response: the so-called 'improper' (fixed field) and 'proper' (fixed voltage) responses, and two intermediate ones.

We summarize these possibilities in Fig. 4 to help the reader in identifying each case. In order to correctly compare measurements from first-principles calculations with experimental measurements, it is important to identify which response from Fig. 4 has been measured in each case and, if necessary, make the appropriate correction.

Because the electrostrictive response is proportional to a third derivative of the free energy (Eq. ( 11)), it can be expressed in several different ways:

m αβ γδ = ∂ ε αβ ∂ η γδ = ∂ 2 X γδ ∂ E α ∂ E β = ∂ e β γδ ∂ E α . (23)
Therefore, the electrostrictive response can be calculated from first-principles using three different approaches:

(i) by calculating the linear evolution of the dielectric permittivity with respect to stress or strain.

(ii) by calculating the quadratic stress induced by an applied electric field.

(iii) by calculating the linear of evolution of the piezoelectric tensor in response to an applied electric field.

Second derivatives of the free energy can be calculated from first-principles calculations using second-order DFPT, and in

(P, E) → (P , E ) Ω 0 ↓ Ω Improper Intermediate Intermediate Proper F Ω 0 = F cell -E • P + 1 2 0 E 2 e αβγ = e (i) αβγ m αβγδ = m (i) αβγδ F Ω 0 = F cell - Ω 0 Ω E • P + 1 2 0 E T g -1 E e αβγ = e (i) αβγ -P γ δ αβ m αβγδ = m (i) αβγδ -αγ δ βδ -δβ δ αγ F Ω = F cell - Ω Ω 0 E • P + 1 2 0 E 2 e αβγ = e (i) αβγ + P α δ βγ m αβγδ = m (i) αβγδ + αβ δ γδ F Ω = F cell -E • P + Ω Ω 0 1 2 0 E T g -1 E
e αβγ = e (i) 4. Summary of the relations between the free energies, piezoelectric and electrostrictive responses for the four possible scenarios illustrated in Fig. 2.

αβγ + P α δ βγ -P γ δ αβ m αβγδ = m (i) αβγδ + αβ δ γδ -αγ δ βδ -δβ δ αγ FIG.
ABINIT, the proper piezoelectric response is automatically obtained because reduced units (E ′ , P ′ ) are used for the electric field and polarization. However, third-order DFPT for the electrostrictive response is not currently implemented in any widely available DFT code, and therefore in all three of the above methods, at least one of the derivatives must be evaluated using finite difference methods. Thus, a correction must be applied in order to obtain the proper response at fixed voltage.

With method (i), first introduced in Ref. [9] and used in this study, the electrostrictive response is obtained by measuring the derivative of the permittivity with respect to strain. The permittivity can easily be obtained from a single DFPT calculation, avoiding the need to perform a set of finite field calculations for each value of strain. Taking the derivative of the dielectric permittivity output from the code from finite differences will yield the improper response, i.e. Eq. ( 11) [9]. The proper response can be obtained by taking care to write the permittivity in terms of E ′ and P ′ , i.e. using Eq. ( 13),

m ≈ ε( η+ )-ε( η-) η+ -η-
, or by applying the correction (Eq. ( 15)) to the improper response. We propose that method (i) is the most straightforward method for calculating the electrostrictive response from first-principles.

With method (ii), the electrostrictive response M is obtained by measuring the quadratic strain response to an applied field, and measuring the curvature about zero field. This has previously been the most widely used method to measure electrostrictive responses from first-principles. However, if finite electric field calculations are performed [9,[START_REF] Kornev | AIP Conference Proceedings[END_REF][START_REF] Pedesseau | [END_REF], the improper response will be obtained and a correction will be needed. The proper response can be directly obtained by fixing the reduced electric field E ′ or the displacement field [9,28,29], both of which are possible in ABINIT, but not in most other widely available DFT codes. Additionally, a geometry relaxation calculation must be performed for each field value in order to find the induced strain, making this method more computationally expensive than method (i).

With method (iii), the electrostrictive response could be obtained by measuring the change in the piezoelectric tensor in response to an applied field [30]. To our knowledge, the electrostrictive response has not yet been calculated from firstprinciples in this way. However, if the proper piezoelectric response is obtained, and the reduced electric field E ′ or displacement field is held fixed, then the proper piezoelectric response should be obtained.

For the electronic permittivity, the sign of the slope changes when correcting from improper to proper. The reason for this is that the electronic permittivity is small and more sensitive to strain than the larger relaxed-ion permittivity (i.e. including phonon contributions), and the correction appears relatively larger.

In each case, the electrostrictive response is larger after correcting to the fixed voltage case. Agreement with experimental measurements is improved for the case of LiF, but not for MgO or NaCl. For the case of LiCl, the reported values in Ref. [21] are an order of magnitude larger than the typical values measured in other experimental studies. Experimental measurements of electrostrictive responses tend to vary significantly, in part due to the various measurement techniques. For example, the M 11 electrostrictive coefficient of MgO has been reported to be 2020 pm 2 /V 2 in Ref. [START_REF] Sundar | [END_REF] and about four times less (550 pm 2 /V 2 ) in Ref. [30]. For SrF 2 , the situation is even more dramatic as experimental reports disagree even on the sign of the longitudinal M coefficient: -1160 pm 2 /V 2 from Ref. [31] and +260 pm 2 /V 2 from Ref. [32]. As a consequence, the agreement with experimental values has to be considered with care. One typical source of error between first-principles calculations and experimental measurements is that first-principles calculations are typically done at zero temperature, whereas experimental measurements are performed at larger temperatures (in this case room temperature [START_REF] Sundar | [END_REF]). However, the evolution of the electrostrictive coefficients with respect to temperature is normally very small, provided a phase transition does not occur [7].

In spite of this, it is important to emphasise that the fixed field and fixed voltage responses are physically different quantities. When comparing first-principles calculations with experimental measurements, it is of primary importance to first ensure that the correct physical quantities are being compared.

  Illustration of the four possible ways to define electromechanical responses, depending on the quantity held fixed (E or ∆V ) and the units of volume used (Ω 0 or Ω). The improper response is obtained by fixing E and working in units of unstrained volume, Ω 0 , and the improper response is obtained by fixing ∆V and working in units of strained volume, Ω. Two additional 'intermediate' responses which are neither proper nor improper are possible by switching the quantity held fixed or the units of volume.

	Improper	Intermediate
	Fixed E Unstrained volume	Fixed ∆V Unstrained volume
	Intermediate	Proper
	Fixed E Strained volume	Fixed ∆V Strained volume
	FIG. 2.

TABLE I .

 I Measurements of the improper, proper intermediate electrostrictive responses M 11 ≡ M 1111 and M 12 ≡ M 2211 of MgO, LiCl, LiF and NaCl, compared to experimental measurements: a = Ref. [19], b = Ref.

ACKNOWLEDGEMENTS

The authors thank M. Stengel for helpful discussions. DB acknowledges funding from the University of Liége under special funds for research (IPD-STEMA fellowship programme). Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (C ÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region and using the DECI resource BEM based in Poland at Wrocław with support from the PRACE OFFSPRING project. EB acknowledges FNRS for support and DT aknowledge ULiége Euraxess support. This work was also performed using HPC resources from the "Mésocentre" computing centre of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/). Financial support is acknowledged from public grants overseen by the French National Research Agency (ANR) in the ANR-20-CE08-0012-1 project and as part of the ASTRID program (ANR-19-AST-0024-02).