Weiguang Zhou

Yanping Zhang

Mathieu Salanne
email: mathieu.salanne@sorbonne-universite.fr

Effects of fluoride salt addition to the physico-chemical properties of the MgCl2-NaCl-KCl heat transfer fluid : a molecular dynamics study

Keywords: Molecular dynamic simulation, Molten ternary chloride, Fluoride, Local structure, Thermodynamic properties

Concentrated solar plants are promising solutions for electricity production. In these plants, the heat transfer fluid plays an important role, and finding systems with good thermal properties is very important. In this regard, molten salts, and more particularly molten chlorides, are currently investigated. Experimental studies of these melts are difficult and expensive, so complementing them with simulations would allow to test a wider range of compositions. In this work, we show that classical molecular dynamics simulations are suitable for predicting the properties of a ternary salt composed of MgCl2, KCl and NaCl by extensive comparisons with experimental data (and previous simulations) on the density, heat capacity, viscosity and thermal conductivity. We then study the effect of adding fluoride ions in the melt on these properties in order to investigate the suitability of mixed chlorides-fluorides for future heat transfer fluids studies.

Introduction

Concentrated solar power (CSP) plant is a promising high-temperature solar power technology and is becoming one of the most massively planned solar power plant technology in the next few decades.

An electricity generation target of 630 GW for CSP technology in 2050 is set by International Energy Agency (IEA) [START_REF] Guillot | Corrosion effects between molten salts and thermal storage material for concentrated solar power plants[END_REF]. According to the US Department of Energy's sunshine initiative, increasing the competitiveness of CSP plants requires [START_REF] Guillot | Corrosion effects between molten salts and thermal storage material for concentrated solar power plants[END_REF]

H ij V B k , i j q q , , x y z L L L ij r r ij B ij h p C ij C ij D µ , , , ij ij ij ij dd dq dd dq f f d d l , , , ij ij ij n n n f c b n J t ( ) M t N i r i M D V A N , , , i j a b
, , , i j a b T , , x y z , , x y z P improve theoretical efficiency, and (2) effective thermal energy storage (TES) strategies to ensure dispatchability [START_REF] Ding | Progress in Research and Development of Molten Chloride Salt Technology for Next Generation Concentrated Solar Power Plants[END_REF]. Currently, most CSP plants are equipped with TES system, so they can supply dispatchable and low-cost electricity with abundant but intermittent solar energy. In a TES system, the heat transfer fluid (HTF) is particularly critical for overall performance and efficiency. The properties of an ideal HTF include: low melting point, high boiling point, thermal stability, low vapor pressure at high temperatures, low corrosion, low viscosity, high thermal conductivity, high specific heat capacity and low cost [START_REF] Vignarooban | Heat transfer fluids for concentrating solar power systems -A review[END_REF][START_REF] Pacio | Assessment of liquid metal technology status and research paths for their use as efficient heat transfer fluids in solar central receiver systems[END_REF][START_REF] Cordaro | Multicomponent molten salt mixtures based on nitrate/nitrite anions[END_REF]. Molten salts are the most widely studied HTFs due to their good properties [START_REF] Peng | The preparation and properties of multi-component molten salts[END_REF]. Molten nitrate salts are used in many large-scale current CSP plants, but they are limited to low temperatures, generally below 600℃ [START_REF] Herrmann | Survey of Thermal Energy Storage for Parabolic Trough Power Plants[END_REF][START_REF] Stern | High Temperature Properties and Decomposition of Inorganic Salts Part 3, Nitrates and Nitrites[END_REF]. Nowadays, hightemperature CSP plants are developing, even exceeding 950℃ [START_REF] Hoffschmidt | Performance Evaluation of the 200-kWth HiTRec-II Open Volumetric Air Receiver[END_REF]. Myers et al. [START_REF] Myers | Thermal energy storage using chloride salts and their eutectics[END_REF] compared the data and concluded that chloride salts and their eutectics, such as the eutectic formed by the ternary system MgCl2-NaCl-KCl (MNK), will be increasingly used in high-temperature TES system in CSP plants. It is therefore necessary to study the properties of molten chloride salts and their eutectics in the applications of HTF.

Molecular dynamics (MD) simulations have the advantages of low cost and high accuracy in acquiring the properties of molten salt at high temperature, while traditional experiments are costly and demanding. The main ingredient of these simulations is the interaction potential between the atoms, which can be estimated using an analytical force field (classical MD) or using first-principles electronic structure calculations, generally performed at the density functional theory (DFT) level (first-principles MD, FPMD). The latter being more accurate, it has been extensively been used to study molten chloride HTFs. Using this method, Liang et al. [START_REF] Liang | First-principles molecular dynamics simulations on the local structure and thermo-kinetic properties of molten magnesium chloride[END_REF] studied the local structure and thermo-kinetic properties of pure molten magnesium chloride. Li et al. [START_REF] Li | Unrevealing the thermophysical properties and microstructural evolution of MgCl2-NaCl-KCl eutectic: FPMD simulations and experimental measurements[END_REF] unrevealed the thermophysical properties and microstructural evolution of MNK eutectic by comparing experimental measurements and FPMD simulations. Xu et al. [START_REF] Xu | Powerful predictability of FPMD simulations for the phase transition behavior of NaCl-MgCl2 eutectic salt[END_REF] predicted the occurrence of structural transitions in the NaCl-MgCl2 eutectic salt, which also provide insights for the applications of molten salts in TES system. Li et al. [START_REF] Li | Insight into dynamic interaction of molten MgCl2-NaCl-KCl with impurity water via FPMD simulations[END_REF] studied the dynamic interaction of molten MNK with water impurity, which sheds insight for the development of salt purification technology.

However, FPMD has a large computational cost, which could be avoided using classical MD.

Simple molten chlorides, such as LiCl, NaCl, KCl and their mixtures have long been studied using the Born-Mayer-Huggins (BHM) potential [START_REF] Ciccotti | Transport properties of molten alkali halides[END_REF], using the parameters proposed by Fumi and Tosi [START_REF] Tosi | Ionic sizes and born repulsive parameters in the NaCl-type alkali halides-II: The generalized Huggins-Mayer form[END_REF]. Based on the same force field, Xie et al. [START_REF] Xie | Heat and mass transportation properties of binary chloride salt as a high-temperature heat storage and transfer media[END_REF] calculated the heat and mass transportation properties of LiCl-KCl and obtained results consistent with experimental literature data. However, the predictions may be improved by the use of polarizable force fields, such as the polarizable ion model (PIM) proposed by Madden and co-workers [START_REF] Wilson | Polarization effects in ionic systems from first principles[END_REF]. In the PIM, induced dipoles are introduced as additional degrees of freedom, which yields a better representation of complex ionic arrangements. Such force fields can be parameterized using DFT, ensuring a better transferability from one system to another [START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF]. Indeed, Wu et al. [START_REF] Wu | Molecular dynamics simulation on local structure and thermodynamic properties of molten ternary chlorides systems for thermal energy storage[END_REF] investigated structural and thermodynamic properties of LiCl-NaCl-KCl system with different KCl concentrations by classical MD using either the BHM or the PIM force field and compared the results with reference FPMD calculations.

They concluded that polarization effect increases the accuracy of the results. In parallel to the case of molten chlorides, PIM force fields were also developed for a large variety of molten fluorides [START_REF] Salanne | Including many-body effects in models for ionic liquids[END_REF].

Recently, new parameters were obtained for the MgCl2 system and its mixtures with alkali chlorides [START_REF] Wu | Elucidating Ionic Correlations Beyond Simple Charge Alternation in Molten MgCl2-KCl Mixtures[END_REF], which was validated against X-ray scattering experiments [START_REF] Wu | Elucidating Ionic Correlations Beyond Simple Charge Alternation in Molten MgCl2-KCl Mixtures[END_REF][START_REF] Wu | Temperature Dependence of Short and Intermediate Range Order in Molten MgCl2 and Its Mixture with KCl[END_REF]. This opens the door for the efficient simulation of HTFs based on this compound.

In this work, we first study the local structure and thermophysical properties of the MNK eutectic using PIM-based classical MD simulations. We then investigate the effect of additional fluoride ions, which are expected to deeply change the solvation shell of cations. The evaluated properties include radial distribution functions, density, specific heat capacity, self-diffusion coefficients, shear viscosity and thermal conductivity. All the simulations are performed under a wide range of temperatures. We have examined the precision of our MD simulation setup by comparing the simulated results for the ternary MNK with the available experimental data. After that, a small part of chloride ions in molten salt are replaced with fluoride ions, in order to compare and analyze the impact of fluoride ions on the local structure and thermophysical properties of molten salts.

Methodology

Force field

In this study, we use a pairwise additive force field, in which all the terms depend on the distance between the ions 𝑟 "# . It is given by:

(1)

The first term on the right side of Eq. ( 1) denotes the electrostatic interaction between ions. Then the second to fourth terms are Fumi-Tosi potential [START_REF] Tosi | Ionic sizes and born repulsive parameters in the NaCl-type alkali halides-II: The generalized Huggins-Mayer form[END_REF], in which the second term represents shortrange repulsion force, while the third and fourth term describe dipole-dipole and dipole-quadrupole dispersion interactions respectively. Tang-Toennies functions [START_REF] Tang | An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients[END_REF] are used to damp the dispersion terms (third and fourth term) at short range:

(2)

(3)
All the parameters used in the Fumi-Tosi potential are listed in Table 1.

6 8 3 exp i j ij ij j i ij ij ij ij ij ij ij ij ij dd ij dq ij n ij ij ij ij ij q q C D q r V r B r f r f r f r r r r r µ h × = + - - - + ( ) ( ) ( ) 6 0 1 exp ! k ij dd ij ij ij dd ij dd ij k d r f r d r k = = - - å ( ) ( ) ( ) 8 0 1 exp ! k ij dq ij ij ij dq ij dq ij k d r f r d r k = = - - å
The last term of Eq. ( 1) accounts for the polarization effects. Polarizable atoms carry an induced dipole 𝜇 " interacting with all the point charges and all the other dipoles [START_REF] Wilson | Polarization effects in ionic systems from first principles[END_REF]. The influence of polarization effect on the structure and characteristics of molten salt should not be ignored, in particular in the presence of divalent species such as the Mg 2+ cation [START_REF] Salanne | Polarization effects in ionic solids and melts[END_REF]. As mentioned in the introduction, studies have shown that the results of molecular simulation considering the polarizable ion model are more robust than with the simpler BHM model [START_REF] Lu | Thermal and Transport Properties of Molten Chloride Salts with Polarization Effect on Microstructure[END_REF]. The charge-dipole interactions are also damped at short-range with Tang-Toennies functions [START_REF] Tang | An improved simple model for the van der Waals potential based on universal damping functions for the dispersion coefficients[END_REF] defined as:

All the parameters of polarization term are listed in Table 2. [START_REF] Ishii | Transport coefficients and the Stokes-Einstein relation in molten alkali halides with polarisable ion model[END_REF]. The dipoles are calculated selfconsistently at each timestep of the simulation using the conjugate gradient method, so that the polarization terms have a many-body character despite the pairwise additive form of the potential in Eq (1).

( ) ( ) ( ) 0 1 exp ! k ij n n ij ij ij ij n ij n n ij k b r f r c b r k = = - - å 2.

Simulation details

All simulations in this study are performed using the MetalWalls (MW) software [START_REF] Marin-Laflèche | MetalWalls: A classical molecular dynamics software dedicated to the simulation of electrochemical systems[END_REF], which is a MD code dedicated to the modelling of electrochemical systems. The simulation is divided into two groups, a first one labeled "base salt", which consists in the ternary MgCl2-NaCl-KCl eutectic salts with a molar ratio of 45.37:32.96:21.67 mol%. The other salts include fluoride ions, that is they are made from the base salt in which 1% ,5%, 10% and 20% of chloride ions are replaced with fluoride ions. In the following they are labeled "MgNaKClF1%" , "MgNaKClF5%", "MgNaKClF10%" and "MgNaKClF20%" respectively. The corresponding atom numbers are provided in Table 3. Also in each simulation, the periodic boundary conditions method is employed to keep constant the number of particles and eliminate boundary effects in the simulation box. The Ewald summation [START_REF] Ewald | Die Berechnung optischer und elektrostatischer Gitterpotentiale[END_REF] method is used for long-range Coulomb potential as well as the charge-dipole and the dipole-dipole interactions. The velocity-Verlet algorithm is used to integrate the equations of motion in MetalWalls [START_REF] Swope | A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters[END_REF]. The experimental melting point of the ternary eutectic salts is about 656 K, so the preset temperatures are 800 K, 900 K, 1000 K, 1100 K and 1200 K for each system. Timestep for integrating the equations of motion is 0.5 fs. For all the simulations, the system are first simulated at ambient pressure in the NPT ensemble with Nosé-Hoover thermostat [START_REF] Nosé | A unified formulation of the constant temperature molecular dynamics methods[END_REF][START_REF] Hoover | Canonical dynamics: Equilibrium phase-space distributions[END_REF] and barostat [START_REF] Martyna | Constant pressure molecular dynamics algorithms[END_REF] at a constant pressure of 0.1 MPa for 500000 steps to reach equilibrium. After that, we keep the equilibrated system volume constant and simulate it in the NVT ensemble for at least 2000000 steps

(1 ns). The relaxation times are 0.5 ps for the barostat and 0.1 ps for the thermostat.

Evaluated properties 2.3.1 Radial distribution functions

The local structures of molten salts is described by the partial radial distribution functions (RDF), which are defined as [START_REF] Cordaro | Multicomponent molten salt mixtures based on nitrate/nitrite anions[END_REF] where 𝜌 D is the number density of species 𝛽, 𝑁 GD (𝑟) is the number of species 𝛽 locating in a sphere with 𝑟 as the radius and a specie 𝛼 as the center.

In this work, all RDF data was calculated using TRAVIS [START_REF] Brehm | TRAVIS -a free analyzer and visualizer for Monte Carlo and molecular dynamics trajectories[END_REF] which is a program package for analyzing and visualizing molecular dynamics trajectories.

Density

Density is one of the critical properties of molten salt. It is defined as [START_REF] Peng | The preparation and properties of multi-component molten salts[END_REF] Where 𝑁 is the number of particles, 𝑀 the molar mass and 𝑉 the equilibrium volume of the system at given temperature and ambient pressure in the NPT ensemble simulations, and 𝑁 M is Avogadro's constant.

Specific heat capacity

Specific heat capacity is calculated as [START_REF] Herrmann | Survey of Thermal Energy Storage for Parabolic Trough Power Plants[END_REF] Where 𝑈 is the total energy including kinetic energy and potential energy of the system. The variation of H is linear with respect to the temperature, so that 𝐶 O can be calculated from the slope of temperature enthalpy diagram of NPT simulations.

Shear viscosity

Shear viscosity (𝜇) is a critical property related to diffusion coefficients and structure, which partially determines the types of pipes and equipment used in the transportation and circulation of ( ) ( )

2 1 g 4 dN r r r dr ab ab b pr é ù = ê ú ë û A NM VN r = ( ) p p U PV H H C T T T D + ¶ D ae ö ae ö = » = ç ÷ ç ÷ ¶ D D è ø è ø
HTFs in CSP. The shear viscosities are calculated in equilibrium MD using the Green-Kubo relationship. Specifically, they are obtained by time integration of the autocorrelation function of the stress tensor obtained by simulation [START_REF] Allen | Computer simulatons of liquids[END_REF] is given as follows,

Where 𝑘 Q is Boltzmann constant and 𝑃 ST are the off-diagonal components of the stress tensor.

In this study, there are five independent anisotropic components of the stress tensor, including 2𝑧𝑧 -𝑥𝑥 -𝑦𝑦, 𝑥𝑥 -𝑦𝑦, 𝑥𝑦, 𝑥𝑧 and 𝑦𝑧. Fig. 1 shows the autocorrelation function of these five stress tensor components for the MgNaKCl system at 1000 K as an example. It can be seen that each autocorrelation functions first decreases with time, and then fluctuates around 0, indicating that its time correlation has decayed. An average over the five components is made in order to improve the statistics. It is found that the integral of the integral of this averaged function reaches a plateau value after about 2 ps, even if the individual autocorrelation functions remained quite noisy.

The viscosity is calculated from the value of this plateau. Note that 2 ps represents an estimate for the time scale of the structural relaxation of the system. ( ) ( )

0 0 xy xy B V P P t dt k T µ ¥ = ò 2.

Thermal conductivity

The thermal conductivity is calculated from Fourier's law using non-equilibrium MD (NEMD) method. [START_REF] Hoffschmidt | Performance Evaluation of the 200-kWth HiTRec-II Open Volumetric Air Receiver[END_REF] Where 𝐽 is the heat flux and ∇𝑇 is the temperature gradient. The NEMD method calculates thermal conductivity by imposing a heat flux to the simulation box. As illustrated on Fig. 2, the simulation box is divided into 20 slices evenly along 𝑧-axis. The 1 st slice is the hot region and 11 th slice is the cold one. Because of the periodic boundary conditions, there are two resulting temperature gradients across the cell. Technically, this is achieved by modifying the particle velocity to increase or decrease a certain amount of energy in the two regions [START_REF] Ikeshoji | Non-equilibrium molecular dynamics calculation of heat conduction in liquid and through liquid-gas interface[END_REF]. The temperature gradient ∇𝑇 is calculated from the average temperature of each slice. All the simulations start in the NVT ensemble until the system stabilizes at the set temperature. In a second step, the heat flux is imposed in NVE simulations for a total of 10 5 time steps to get the linear temperature gradient. Thermal conductivity is then calculated from Eq. ( 9).

Self-diffusion coefficient

The ion self-diffusion coefficient is calculated to study the transport properties of molten salt. In this work, they were calculated from time-dependent mean squared displacement (MSD) according to Einstein equation [START_REF] Allen | Computer simulatons of liquids[END_REF].

Where the angle brackets <…> denotes ensemble average and 𝑟 " is the position of ion 𝑖.

Results and discussions

Base salt

In order to study the local structure of the chloride at different temperatures above the melting point, a series of partial radial distribution functions were calculated. Fig. 3 shows the RDF curves of the cation-anion pair and anion-anion pair in the MgNaKCl simulation system at 800 K. It can be seen that all RDF curves start from 0: because of the short-range repulsion, there is no other ions within the ionic radius of the ions. Then, all curves display a well-defined first peak which characterizes the short-range order. It is followed by a series of peaks which intensity gradually decreases with the distance; the RDFs progressively flatten and tend towards 1, which corresponds to the typical structure of molten salts [START_REF] Lantelme | Application of the molecular dynamics method to a liquid system with long range forces (Molten NaCl)[END_REF]. ( ) ( ) ( )

2 0 i i M t r t r = - ( ) 1 6 dM t D dt =
More specifically, Fig. 3 shows the RDFs of Na + -Cl -, K + -Cl -, Mg 2+ -Cl -and Cl --Cl -ion pairs at 800 K. It can be seen the first peak of Mg 2+ -Cl -is significantly higher than that of Na + -Cl -and K + -Cl -, because of the stronger interaction between Mg 2+ -Cl -ion pairs. The first peak positions of each ion pair from FPMD simulation by Li et al. [START_REF] Li | Unrevealing the thermophysical properties and microstructural evolution of MgCl2-NaCl-KCl eutectic: FPMD simulations and experimental measurements[END_REF] are also marked in Fig. 3. In all cases the agreement is good for both the positions and the heights of the peaks, showing the ability of the PIM-based classical simulations to predict the correct structure of the melt. Fig. 4 shows RDFs of Na + -Cl -, Mg 2+ -Cl -, K + -Cl -and Cl --Cl -ion pairs at different temperatures. It can be seen from the partial enlarged view of the peak that as the temperature increases, the height of the peak decreases, which is due to the decrease of the density. Moreover, it can be found that the positions of the first peaks of these three cation-anion pairs are slightly shifted leftward as the temperature increases, which means that the distance between anion and cation becomes shorter at high temperatures. However, there is a slight shift to the right for the anion-anion pairs, which indicates that the distance between the ion clusters increases at high temperatures. Obviously, density and temperature are negatively correlated linearly. The densities decrease from 1.92 g/cm a to 1.64 g/cm a with an increase in temperature from 800 K to 1200 K. A linear fit of our results yield the following relationship. [START_REF] Li | Unrevealing the thermophysical properties and microstructural evolution of MgCl2-NaCl-KCl eutectic: FPMD simulations and experimental measurements[END_REF] The simulation results are slightly larger than the experimental results, and the overall deviation is 8.14%, which may be due to the fact that the force field parameters were derived from DFT calculations, and are not empirical [START_REF] Wu | Temperature Dependence of Short and Intermediate Range Order in Molten MgCl2 and Its Mixture with KCl[END_REF][START_REF] Ishii | Transport coefficients and the Stokes-Einstein relation in molten alkali halides with polarisable ion model[END_REF]. Indeed, despite the intrinsically more accurate character of FPMD, Li et al. obtained similar level of agreement with the experimental density in their previous study [START_REF] Li | Unrevealing the thermophysical properties and microstructural evolution of MgCl2-NaCl-KCl eutectic: FPMD simulations and experimental measurements[END_REF]. Fig. 6 shows the temperature variation enthalpy of the MgNaKCl system, obtained from the NPT simulations. The enthalpy varies linearly, which shows that the specific heat capacity is nearly constant among the whole temperature range based on equation [START_REF] Herrmann | Survey of Thermal Energy Storage for Parabolic Trough Power Plants[END_REF]. It is worth noting that some discrepancies are observed among experimental results about such a trend. In recent studies (on different salt compositions), Wang et al. reported changes of the heat capacity with temperature [START_REF] Wang | Thermophysical Properties Experimentally Tested for NaCl-KCl-MgCl2 Eutectic Molten Salt as a Next-Generation High-Temperature Heat Transfer Fluids in Concentrated Solar Power Systems[END_REF], while Xu et al. did not observe such a variation [START_REF] Xu | Basic properties of eutectic chloride salts NaCl-KCl-ZnCl2 and NaCl-KCl-MgCl2as HTFs and thermal storage media measured using simultaneous DSC-TGA[END_REF]. In their review of the literature, Villada et al.

showed that the latter behavior is more common, in agreement with our results [START_REF] Villada | Engineering molten MgCl2-KCl-NaCl salt for hightemperature thermal energy storage: Review on salt properties and corrosion control strategies[END_REF]. The calculated specific heat capacity of MgNaKCl is 1.09 J/(g • K), is in good agreement compared with the experimental result of 1.14 J/(g • K) by Li et al. [START_REF] Li | Unrevealing the thermophysical properties and microstructural evolution of MgCl2-NaCl-KCl eutectic: FPMD simulations and experimental measurements[END_REF] with a deviation of 4.4%. As shown in Fig. 7, the viscosity of the system decreases sharply as the temperature rises. Although no available experimental data on the viscosity of our studied composition in the literature, we compare the simulation results to the ones on MgCl2-NaCl-KCl (50-30-20 mol.%) by Li et al. [START_REF] Li | Survey and evaluation of equations for thermophysical properties of binary/ternary eutectic salts from NaCl, KCl, MgCl2, CaCl2, ZnCl2 for heat transfer and thermal storage fluids in CSP[END_REF] and MgCl2-NaCl-KCl (47.1-30.2-22.7 mol.%) by Villada et al. [START_REF] Villada | Engineering molten MgCl2-KCl-NaCl salt for hightemperature thermal energy storage: Review on salt properties and corrosion control strategies[END_REF]. Considering the different compositions of the molten salt, the agreement is good, except for the lowest studied temperature at which the predicted viscosity is too large. For the NEMD calculation of the thermal conductivity, parallelepipedal simulation cells were generated by replicating several times the initial cubic cell along the z-axis. We performed several tests in which the simulation cell was replicated between 2 and 6 times. They yielded similar results, so that the cell with Lz = 2 Lx was further chosen in order to keep a low computational cost. The temperature profiles obtained for MgNaKCl with an applied heat flux of 2.32 × 10 ghi W/Å 2 at each temperature are plotted in Fig. 8. A linear variation is observed between the hot region and the cold region. The slopes of the temperature variation along the z dimension were extracted in order to compute the thermal conductivity of the system according to Eq. ( 9). As shown in Fig. 9, the thermal conductivity of MgNaKCl decreases with temperature, from 0.44 W/K • m at 800K to 0.36 W/K • m at 1200K. As for the viscosity, experimental data on this specific MgNaKCl composition is lacking. We thus compare the simulation results with the measurements of Wang et al. [START_REF] Wang | Thermophysical Properties Experimentally Tested for NaCl-KCl-MgCl2 Eutectic Molten Salt as a Next-Generation High-Temperature Heat Transfer Fluids in Concentrated Solar Power Systems[END_REF] for a (45.98%-15.11%-38.91% mol) composition; a good overall agreement is obtained.

Villada et al. [START_REF] Villada | Engineering molten MgCl2-KCl-NaCl salt for hightemperature thermal energy storage: Review on salt properties and corrosion control strategies[END_REF] also tried to estimate the conductivity of the (47.1-30.2-22.7 mol.%) composition using a linear mixing rule and the values provided in the literature for the pure salts (NaCl and KCl by Nagasaka et al. [START_REF] Nagasaka | Experimental determination of the thermal diffusivity of molten alkali halides by the forced Rayleigh scattering method. I. Molten LiCl, NaCl, KCl, RbCl, and CsCl[END_REF], MgCl2 by Lu et al. [START_REF] Lu | Thermal and Transport Properties of Molten Chloride Salts with Polarization Effect on Microstructure[END_REF]). Their results are slightly overestimated, which may be due to the use of mixing rule. As was demonstrated from simulations of binary molten fluorides [START_REF] Ishii | Thermal Conductivity of Molten Alkali Metal Fluorides (LiF, NaF, KF) and Their Mixtures[END_REF], the mixtures systematically display lower values for the thermal conductivity than pure salts. Fig. 10 shows the linear relationship between the MSD of each ion and time at 1000 K. According to Eq. ( 11), the diffusion coefficient of each ion at each temperature is calculated by the slope. It can be found that the order of the diffusion coefficient of each ion is D mn o > D q o > D rs t > D uv wo in Fig. 11, which is affected by the mass, ionic radius, the number of charges and ion interaction. As far as we know, there are few studies on the self-diffusion coefficient of ternary MgCl2-NaCl-KCl molten salts, but by comparing with the experimental data of NaCl by Janz [START_REF] Janz | II.C -METAL OXIDE-MOLTEN SALT SYSTEMS[END_REF] , the experimental data of KCl by Bockris et al. [START_REF] Bockris | Self-Diffusion and Structure in Molten Group II Chlorides1[END_REF] and the FPMD data of MgCl2 by Liang et al. [START_REF] Liang | First-principles molecular dynamics simulations on the local structure and thermo-kinetic properties of molten magnesium chloride[END_REF] , the differences are small and the trend is consistent. 11 shows that the self-diffusion coefficient increases with increasing temperature, which means that the movement of ions will be easier. This agrees with previous result: as the temperature rises, the structure of the molten salt becomes loose and diffusion becomes easier. The self-diffusion coefficient of divalent cation Mg 2+ is smaller than of monovalent species due to the stronger Mg 2+ -Cl -interaction. are compared with previous experimental (Janz [START_REF] Janz | II.C -METAL OXIDE-MOLTEN SALT SYSTEMS[END_REF], Bockris [START_REF] Bockris | Self-Diffusion and Structure in Molten Group II Chlorides1[END_REF]) and FPMD (Liang [START_REF] Liang | First-principles molecular dynamics simulations on the local structure and thermo-kinetic properties of molten magnesium chloride[END_REF]) studies on the pure salts.

Effect of Fluoride salt addition

Since the computed properties on the ternary molten chloride are in agreement with experimental data, it is possible to make predictions on new salt compositions. In this study, four systems were simulated, in which 1%, 5%, 10% and 20% of the chloride ions of the base salt were replaced by fluoride ions to study their effect on the local structure and thermodynamic properties.

For the local structure, the RDFs of cation-anion and anion-anion pairs are shown in Fig. 12. At first sight, the interactions seem to be stronger for the Na + -Cl -pair than for the Na + -F -one. A similar effect is observed for K + -Cl -and K + -F -. On the contrary, there is stronger interaction between Mg 2+ -F -than Mg 2+ -Cl -. These observations can therefore be rationalized by the fact that the Mg ion attracts more fluoride ions in its vicinity due to its divalent character, which explains the depletion of the latter in the solvation shell of Na + and K + .

For anion-anion pairs, the Cl --F -nearest-neighbor distance is shorter than the Cl --Cl -one, which is due to the smaller radius of F -. This indicates that the distance between the ion clusters becomes smaller upon addition of F -, and thus the system volume will be smaller, which can impact the density. The smaller distance between the ion clusters may also lead to larger shear viscosity.

As the proportion of fluoride ions increases, the height of the peak between Na + -Cl -and K + -Cl - gradually increases, while the one of the Mg 2+ -Cl -RDF gradually decreases. In contrast, as the proportion of F -increases, the interaction between Na + -F -and K + -F -gradually decreases, while the Mg 2+ -F -one gradually increases. Again, the main reason for this difference is that the interaction between Mg 2+ and F -is very strong, so as the proportion of F -increases, more Mg 2+ will surround F -, and the extra vacancies around Cl -will naturally be occupied by Na + and K + cations. We now analyze variations in the density. Fig. 13 shows the density of each system at 1000 K. We observe that replacing part of the chloride ions in MgNaKCl with fluoride ions leads to an increase of the density. The greater the proportion of fluoride, the greater the increase in density; the 20% fluoride system leads to an increase of the density from 1.78 g/cm a to 1.84 g/cm a with respect to the base salt. However, comparing the molar mass of chlorine and fluorine, Cl > F, so the overall mass of the system decreases as the replacement ratio increases. The reason for the increase in density must be related to changes in its volume and structure. From the previous RDFs analysis, the relatively strong interactions of Mg 2+ -F -pair makes the main contribution to volume reduction and structural changes of the system. For specific heat capacity (𝐶 x ), the results are 1.0976 J/(g • K) , 1.1119 J/(g • K) , 1.1213 J/(g • K) and 1.1509 J/(g • K) for MgNaKClF1%, MgNaKClF5%, MgNaKClF10% and MgNaKClF20%, respectively. Fig. 13 shows the specific heat capacity of each simulated systems.

Obviously, the adding of fluoride ions leads to an increase oof the 𝐶 x of the salt, and the larger the proportion of fluoride ion, the more 𝐶 x increases. Specific heat capacity increased by 0.69%, 2.00%, 2.86% and 5.58% respectively. As shown on Fig. 13, this change in the specific density is well correlated with the corresponding density increase with composition.

For viscosity, Fig. 14 shows that the influence of the concentration in fluoride ion is not monotonous; it appears that the effect is important at large concentration only. In addition, at high temperature, the fluoride ion addition has a smaller effect on the viscosity. Note however that the uncertainty on the calculation of the viscosity is relatively large, which may explain the absence of a clear trend with concentration. For thermal conductivity, Fig. 15 shows that the thermal conductivities of the three systems all decrease with increasing temperature. As the proportion of fluoride ions increases, the thermal conductivity of the system increases. The main reason is that the volume of the system becomes smaller, so that the heat transfer between the ions becomes easier. Compared with MgNaKCl, the thermal conductivity of MgNaKClF1%, MgNaKClF5%, MgNaKClF10% and MgNaKClF20% each have an average increase of 1.21% and 5.27%, 5.79% and 13.36%. Finally, Fig. 16 shows the diffusion coefficient of each particles at each temperature for the five systems. Consistent with the previous discussion, the self-diffusion coefficient of each particle in each system increases as the temperature rises. Fig. 17 shows the self-diffusion coefficients of the particles in the 5 different systems at 1000 K. It can be seen that as the proportion of fluoride ions increases, the self-diffusion coefficient of each particle generally decreases, but it is not monotonous.

The decrease is more pronounced at high temperatures. The impact of fluoride ion concentration on the various ions is also different. At 1200 K, the self-diffusion coefficient of Mg 2+ , Na + and K + decreased by 37.99%, 21.95% and 10.97% in the MgNaKClF20% system. For the base salt system, by comparing our simulations with the experimental results in the literature, we showed that the selected force field accounts properly for the interactions of the ternary chloride salt. The temperature dependence of each properties were determined. In a second step, this study discusses the effect of adding fluoride ions on the evaluated properties. In terms of local structure, changes are mainly due to the strong interaction between Mg 2+ -F -, which leads to more F -in the solvation shell of Mg 2+ instead of Cl -. From the perspective of HTF application, high specific heat capacity, low shear viscosity and high thermal conductivity are more important in the selected properties. The addition of fluoride ion is very beneficial to increase the specific heat capacity and thermal conductivity of the molten salt. 20% fluoride addition leads to a specific heat increase of 5.58% and an average thermal conductivity increase of 13.36%. However, the addition of fluoride salt also leads to a higher viscosity. Using this strategy for enhancing HTFs would therefore need to select the best compromise between the various target properties. It should be noted that the phase diagram of the mixed chloride-fluoride salts are not known, which would be a necessary step before a practical application of such salts in HTFs.

Figure 1 .

 1 Figure 1. Stress tensor autocorrelation function of MgNaKCl simulation system at 1000 K. The viscosity is

Figure 2 .

 2 Figure 2. Schematic representation of the simulation box in the NEMD simulations. The box limits are

Figure 3 .

 3 Figure 3. RDFs of cation-anion and anion-anion pairs at 800 K. The results are compared with FPMD

Figure 4 .

 4 Figure 4. RDFs of cation-anion and anion-anion pairs at each temperature

Figure 5 .

 5 Figure 5. Temperature dependence of density for MgNaKCl. The simulation results are compared with

Figure 6 .

 6 Figure 6. Variation of the calculated enthalpy with temperature for the base salt.

Figure 7 .

 7 Figure 7. Viscosity of MgNaKCl system of MgCl2-NaCl-KCl. The simulation results are compared with the

Figure 8 .

 8 Figure 8. Mean temperature profile along the simulation cell at each imposed temperature

Figure 9 .

 9 Figure 9. Thermal conductivities of MgNaKCl system in different temperature. Results are compared with

Figure 10 .

 10 Figure 10. MSD of each ion in MgNaKCl system at 1000 K

Figure 11 .

 11 Figure 11. Self-diffusion coefficient of each ions in MgNaKCl systems in different temperature. The results

Figure 12 .

 12 Figure 12.RDFs of cation-anion and anion-anion pairs for the fluoride ion containing systems

Figure 13 .

 13 Figure 13. Density (at 1000K) and specific heat capacity of each simulated systems

Figure 14 .

 14 Figure 14. viscosity of each simulated systems at each temperatures

Figure 15 .Figure 16 .

 1516 Figure 15. Thermal conductivity of each simulated systems on the whole range of temperatures

Figure 17 .

 17 Figure 17. Self-diffusion coefficient of each particles at 1000K in different systems

  increasing the operating temperature to

	Nomenclature		enthalpy
		Interaction potential	Boltzmann constant
		charge of ion	Cell length dimension in the x, y and z direction
		distance between particle i and j	density
	,	BMH repulsive parameters	Specific heat capacity
	,	Van der Waals parameters	shear viscosity
		Tang-Toennies parameters for Fumi-Tosi potential	thermal conductivity
		Tang-Toennies parameters for polarization effects	heat flux
		simulation time	mean squared displacement
		particle number	position of ion
		molar mass	self-diffusion coefficient
		equilibrium volume	Subscripts
		Avogadro's constant	particle
		temperature	dimension
		pressure

Table 1 . Parameters of the Fumi-Tosi potential (all in atomic units) [23,25].

 1

	Pairs	𝜂 "#	𝐵 "#	𝐶 "#	𝐷 "#	𝑑 )) "#	𝑑 )* "#
	Cl -Cl	1.797	275.1	140.0	280.0	1.7	1.7
	Cl -Na	1.726	67.5	47.4	187.3	1.7	1.7
	Na -Na	1.000	0.0	11.7	51.8	1.7	1.7
	Cl -K	1.833	208.6	29.4	128.3	1.7	1.7
	K -K	5.00	1.0	1.0	10.0	1.7	1.7
	K -Mg	5.00	1.0	0.0	0.0	1.7	1.7
	Mg -Mg	5.00	1.0	0.0	0.0	1.7	1.7
	Cl -Mg	1.732	68.6	10.0	20.0	1.7	1.7
	K -Na	5.00	1.0	0.0	0.0	1.7	1.7
	Mg -Na	5.00	1.0	0.0	0.0	1.7	1.7
	Cl -F	2.096	278.8	45.8	204.94	1.7	1.7
	K -F	2.04	138.8	3.87	38.7	1.9	1.9
	Na -F	1.974	52.83	13.25	88.15	1.9	1.9
	Mg-F	1.938	45.58	13.25	88.15	1.9	1.9
	F -F	2.44	282.3	15	150	1.9	1.9

Table 2 . Parameters for polarization term (all in atomic units). [c

 2

	charge	dipole	𝑏 -"#	𝑛	𝑐 -"#	species	polarizability
	Na	Cl	1.760	4	3.000
						F	7.894
	K	Cl	1.632	4	3.000
	Mg	Cl	1.873	4	2.875
						Cl	20.000
	Cl	Na	1.760	4	0.697
	F	Na	1.873	4	-0.190
						Na	0.890
	K	Na	1.760	4	0.000
	Mg	Na	1.760	4	0.000
						Mg	2.000
	Cl	K	1.632	4	0.917
	F	K	1.745	4	-0.310
						K	4.717
	Na	K	1.760	4	0.000
	Mg	K	1.760	4	0.000
	Na	F	1.837	4	2.543
	K	F	1.745	4	2.500
	Mg	F	1.918	4	1.858

Table 3 . Number of atoms for each simulated system

 3

	System		Na	K	Mg	Cl	F	Total
	Base salt	MgNaKCl	329	217	454	1454	-	2454
	Fluoride addition	MgNaKClF1%	329	217	454	1439	15	2454
		MgNaKClF5%	329	217	454	1379	75	2454
		MgNaKClF10%	329	217	454	1304	150	2454
		MgNaKClF20%	329	217	454	1154	300	2454

Acknowledgments

The authors acknowledge HPC resources granted by the HPCaVe Centre at Sorbonne Université.

Data availability

The data generated within this study are available on Zenodo at the following address: http://dx.doi.org/10.5281/zenodo.5045186.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

CRediT authorship contribution statement


Hassan Denawi

Mathieu Abel

Olivier Siri

Roland Hayn

Electronic and magnetic properties of metalorganic polymers with 4d and 5d-transition metal ions

published or not. The documents may come

Electronic and magnetic properties of metal-organic

polymers with 4d and 5d-transition metal ions

INTRODUCTION

The study of metal-ligand molecules or polymers has a long history and is the subject of organometallic chemistry [START_REF] Kealy | A New Type of Organo-Iron Compound[END_REF][START_REF] Stone | Perspectives in organometallic chemistry[END_REF][START_REF] Braga | Arene Clusters[END_REF][START_REF] Wadepohl | Benzene and Its Derivatives as Bridging Ligands in Transition-Metal Complexes[END_REF]. They are of highest importance for biological and technological processes and their investigation is on the upturn at present. The choice of the ligand is here a crucial issue and the monocyclic 2,5-oxoamino-1,4-benzoquinonemonoimines (ZQ, Figure 1), featuring a stable zwitterionic ground state, appeared in this context extremely appealing [START_REF] Siri | Unprecedented zwitterion in quinonoid chemistry[END_REF]. The charge separation in ZQ is due to an intramolecular proton transfer that generates two distinct 6π-electron subunits, i.e., a positively charged trimethine cyanine and a negatively charged trimethine oxonol, linked through two single bonds (see ZQ, Scheme 1).

According to the coupling principle of polymethine dyes [START_REF] Dähne | Coupling Principles in Organic Dyes[END_REF], such a ground-state electronic structure is highly stable, generating a large intrinsic dipole (ca. 10 D), that could be exploited for the elaboration of molecular films following adsorption [START_REF] Xiao | The interface bonding and orientation of a quinonoid zwitterion[END_REF][START_REF] Routaboul | Altering the static dipole on surfaces through chemistry: Molecular films of zwitterionic quinonoids[END_REF][START_REF] Simpson | Coverage-dependent interactions at the organics-metal interface: Quinonoid zwitterions on Au(111)[END_REF][START_REF] Simpson | Modulating Bond Lengths via Backdonation: A First-Principles Investigation of a Quinonoid Zwitterion Adsorbed to Coinage Metal Surfaces[END_REF][START_REF] Routaboul | Changes in molecular film metallicity with minor modifications of the constitutive quinonoid zwitterions[END_REF] or polymerization on metallic surfaces [START_REF] Kunkel | Dipole driven bonding schemes of quinonoid zwitterions on surfaces[END_REF][START_REF] Koudia | On-surface synthesis of covalent coordination polymers on micrometer scale[END_REF][START_REF] Denawi | On-Surface Synthesis of Spin Crossover Polymeric Chains[END_REF] along with interfaces [START_REF] Fang | Charge-assisted hydrogen bond-directed self-assembly of an amphiphilic zwitterionic quinonemonoimine at the liquid-solid interface[END_REF][START_REF] Behyan | Structural organization and phase behaviour of: Meta -substituted dioctadecylaminobenzoquinones at the air/water interface[END_REF] and, more recently, as interlayers in graphene transistors [START_REF] Mahmood | Tuning graphene transistors through: Ad hoc electrostatics induced by a nanometer-thick molecular underlayer[END_REF].

Importantly, the use of ZC as a ligand in coordination chemistry has been extensively investigated and clearly demonstrated its efficiency and versatility in solution (Scheme 1). In presence of 1 equiv. of metal precursors (TM), the reaction proceeds by two successive deprotonation/ metalation steps from ZC, affording mononuclear complexes 1-4 with TM = Pd [START_REF] Taquet | Stepwise synthesis, structures, and reactivity of mono-, di-, and trimetallic metal complexes with a 6π + 6π quinonoid zwitterion[END_REF], Cu [START_REF] Braunstein | Synthesis of bis(phosphinoferrocenyl) copper complexes from zwitterionic quinonoid ligands and their structural and redox properties †[END_REF], Ru [START_REF] Das | First structurally characterized mono-and dinuclear ruthenium complexes derived from zwitterionic quinonoid ligands[END_REF] and Pt [START_REF] Kar | Colour tuning by the stepwise synthesis of mononuclear and homo-and hetero-dinuclear platinum(II) complexes using a zwitterionic quinonoid ligand[END_REF]. Interestingly, mononuclear complexes 5-8 with TM/L ratio of 1/2 could be also obtained when 2 equiv. of ZC was used (M = Pd [START_REF] Taquet | Stepwise synthesis, structures, and reactivity of mono-, di-, and trimetallic metal complexes with a 6π + 6π quinonoid zwitterion[END_REF], Ni [START_REF] Taquet | Stepwise synthesis, structures, and reactivity of mono-, di-, and trimetallic metal complexes with a 6π + 6π quinonoid zwitterion[END_REF], Zn [START_REF] Braunstein | Synthesis and properties of copper quinonoid complexes for optical recording application[END_REF] and Cu [START_REF] Braunstein | Synthesis and properties of copper quinonoid complexes for optical recording application[END_REF]). In all cases, each quinoid unit (ZC) is linked to the metal center through two covalent and two coordinative bonds through N and O atoms. As a result, the oxidation state of the metal centers is always +II. Furthermore, the zwitterionic character of the ligand is lost, but its π electrons are now localized. In contrast, the condensation of 2 equiv. of TM precursors on ZC (instead of one) afforded the dinuclear species 9-11 (M = Pd [START_REF] Ghisolfi | Reversible switching of the coordination modes of a pyridine-functionalized quinonoid zwitterion; Its Di-and tetranuclear palladium complexes[END_REF], Ru [START_REF] Das | First structurally characterized mono-and dinuclear ruthenium complexes derived from zwitterionic quinonoid ligands[END_REF],

and Pt [START_REF] Kar | Colour tuning by the stepwise synthesis of mononuclear and homo-and hetero-dinuclear platinum(II) complexes using a zwitterionic quinonoid ligand[END_REF]) in which the π electrons are still delocalized, a key feature to form conductive materials.

Scheme 1

This metal (TM) -ligand (L) assembly was recently investigated by some of us on surface to obtain polymers with full delocalization of the electrons. The on-surface polymerization of Fe and zwitterionic quinone (ZQ) to long polymer chains was indeed successfully reported and attracted considerable attention [START_REF] Koudia | On-surface synthesis of covalent coordination polymers on micrometer scale[END_REF]. Remarkably, this polymerization reaction could not be realized in the gas phase or under wet conditions but only on the surface of noble metals like Ag or Au [START_REF] Denawi | On-Surface Synthesis of Spin Crossover Polymeric Chains[END_REF][START_REF] Denawi | Magnetic Polymer Chains of Iron and Zwitterionic Quinoidal Ligands on the Ag(111) Surface[END_REF]. The polymer chains mentioned above are interesting for spintronics applications in the prospect of spin-polarized one-dimensional (1D) electron transport. With the rapid miniaturization of the semiconductor industry, molecular-sized components will be used in the foreseeable future in hybrid electronics or spin electronics (spintronics) devices.

The Fe-ZQ polymers are also remarkable, since they can be adapted to very different substrates in connection with a predicted spin-state transition and they build different 1D or 2D arrangements [START_REF] Denawi | On-Surface Synthesis of Spin Crossover Polymeric Chains[END_REF][START_REF] Denawi | Magnetic Polymer Chains of Iron and Zwitterionic Quinoidal Ligands on the Ag(111) Surface[END_REF]. There is certainly no experimental difficulty to replace Fe by other transition metal (TM) atoms and the whole class of 3d TM-ZQ polymer chains was studied theoretically [START_REF] Denawi | Magnetic Polymer Chains of Transition Metal Atoms and Zwitterionic Quinone[END_REF] predicting most of them to be antiferromagnetic semiconductors with the exception of metallic behavior for Mn-PC.

But to fully use the capacity of magnetic organometallic polymer chains for spintronics applications, a high value of the magnetic anisotropy energy (MAE) would be crucial. It may stabilize a potential magnetic ordering and opens the possibility of molecular magnetic information storage. Magnetic monolayers of organometallic polymers could be interesting for two different classes of applications. A high MAE of single magnetic units of molecular size qualifies them as candidates for magnetic molecular memory where the best examples are found at present in rare earth ions packed into sandwich complexes [START_REF] Münzfeld | Synthesis, structures and magnetic properties of [(η9-C9H9)Ln(η8-C8H8)] super sandwich complexes[END_REF] or C 80 molecules [START_REF] Liu | Air-stable redox-active nanomagnets with lanthanide spins radical-bridged by a metal-metal bond[END_REF]. Alternatively, if it would be possible to obtain monolayers or thin films which are ferromagnetically ordered, they also could serve in magnetic tunnel junctions where a large magnetic anisotropy barrier is useful for spintronic devices [START_REF] Jiang | Origin of the large interfacial perpendicular magnetic anisotropy in MgO/Co2FeAl[END_REF].

For that reason, 4d and 5d TM materials become interesting, in the form of atomic chains [START_REF] Mokrousov | Giant magnetocrystalline anisotropies of 4d transition-metal monowires[END_REF][START_REF] García-Suárez | Anisotropic magnetoresistance in atomic chains of iridium and platinum from first principles[END_REF] or as Rh-Phthalocyanine nanojunctions [START_REF] Li | Rheniumphthalocyanine molecular nanojunction with high magnetic anisotropy and high spin filtering efficiency Rhenium-phthalocyanine molecular nanojunction with high magnetic anisotropy and high spin filtering efficiency[END_REF]. Since the amount of spin-orbit coupling increases with the atomic number, these 4d and 5d TM materials have usually higher MAE (up to the order of several tens of meV) than the corresponding 3d series. For that reason, we investigate here selected 4d and 5d TM-ZQ polymers and can indeed confirm high MAE's for some of them. Especially remarkable are the significant MAE's for Rh-ZQ and Re-ZQ.

Furthermore, we study in detail the structural, electronic, and magnetic properties of 4d and 5d TM-ZQ polymers.

METHOD

We used the Vienna Ab-initio Simulation Package (VASP) [START_REF] Truhlar | Ab-initio simulations of materials using VASP: Density-functional theory and beyond[END_REF] at the level of the spinpolarized generalized-gradient approximation in the form of the Perdew and Wang PW91 functional [START_REF] Perdew | Erratum: Atoms, molecules, solids, and surfaces: Applications of the generalized gradient approximation for exchange and correlation[END_REF][START_REF] Perdew | Generalized gradient approximation made simple[END_REF]. It is well known that the standard Density Functional Theory (DFT) has difficulties to describe correctly the electronic structure of transition metal compounds, especially for 3d but to a lesser amount also for 4d and 5d transition metal ions. For that reason, we apply here the Spin polarized Generalized-Gradient Approximation with Hubbard term U (SGGA+U). The SGGA+U corrections were introduced by Liechtenstein et al. [START_REF] Liechtenstein | Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators[END_REF] Unfortunately, the value of U is not very well known. We take in the following U = 5 eV and an exchange energy of J = 0.90 eV as it was justified for Fe-ZQ. That has to be understood as an upper bound for the U value since that parameter is certainly reduced in 4d and 5d systems with respect to 3d compounds. The necessity of the SGGA + U method for metal-organic compounds with transition metal ions is proven by many examples such as TM-TCNB [START_REF] Mabrouk | Ligand Influence on Local Magnetic Moments in Fe-Based Metal-Organic Networks[END_REF][START_REF] Mabrouk | Possibility of a ferromagnetic and conducting metal-organic network[END_REF] or TM-TCNQ [START_REF] Ma | Novel two-dimensional tetragonal monolayer: Metal-TCNQ networks[END_REF], TM-ZQ [START_REF] Denawi | Magnetic Polymer Chains of Transition Metal Atoms and Zwitterionic Quinone[END_REF] where TM is a transition metal. We also perform calculations without the Hubbard U correction, i.e. we apply the SGGA method, to compensate the uncertain knowledge of the U parameter. In such a way we can justify the calculated MAE values and estimate their uncertainty.

The interaction between the valence electrons and ionic cores was described within the framework of the projector augmented wave (PAW) method [START_REF] Blöchl | Projector augmented-wave method[END_REF][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 | From ultrasoft pseudopotentials to the projector augmentedwave method[END_REF][START_REF] Kresse | Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[END_REF]. The electronic wave functions were expanded in plane waves with a kinetic energy cutoff of 480 eV, the polymers were relaxed until the residual forces were below eV/Å and the convergence criteria for the energy deviations was eV. The relaxation of all atomic positions guaranties the structural stability of the polymers. The Gaussian smearing method was used in these calculations and a width of 0.01 eV was adopted. The Brillouin zone was determined by a set of 8×1×1 k-points in the unit cell using the Monkhorst-Pack points [START_REF] Pack | special points for Brillouin-zone integrations[END_REF]. The geometry optimization was performed by relaxing all atomic coordinates for a fixed lattice constant (see Figure 1). The geometry optimization was performed without spin-orbit coupling. But a relativistic calculation with spin-orbit coupling was necessary to calculate the orbital moments and the MAE's.

RESULTS AND DISCUSSION

A. Electronic and magnetic properties of TM-ZQ polymers

We first report the results of the SGGA+U calculation without spin-orbit (SO) coupling and summarize the results in Table 1. The valence of the TM can be deduced from the local occupation numbers and was found to be +2 for all investigated polymers. Also, all TM-ZQ compounds were magnetic with the exception of the two TM's with d 8 configuration, Pd-ZQ and Pt-ZQ, which have a nonmagnetic ground-state. Two possible spin states, namely a highspin solution and an intermediate spin solution, were found for Mo-ZQ by the variation of the lattice constant. That is visible in two different curves of total energy vs.

lattice constant for Mo-ZQ which are given in the supplemental material (Fig. S1) for all investigated materials. The energetically lowest state for Mo-ZQ is the high-spin state but that is an exception since all other investigated 4d and 5d TM-ZQ polymers realize the intermediate or no spin state in difference to the 3d series of TM-ZQ [START_REF] Denawi | Magnetic Polymer Chains of Transition Metal Atoms and Zwitterionic Quinone[END_REF]. The optimized lattice constants for TM-ZQ range from 15.74 to 16.30 Å, corresponding to the trend of the atomic radii of the metal which generally decrease with the atomic number in the 4d or 5d series. The same trend is also visible in the interatomic distances d TM-N and d TM as

shown in Table 1. S1 in the Supplementary material). That is especially true for all structural parameters (not all are shown). Also, the exchange energies E ex have the same sign and the same order of magnitude without or with U for all investigated compounds. As it can be expected, the modulus of E ex is larger in SGGA than in SGGA+U with the only exception Rh-ZQ. Different are the band gaps which are reduced in SGGA for the Ru-, Pd-, Re-and Pt-ZQ polymers and even absent (metallic behaviour) for the Mo-, Rh-and Ir-ZQ polymers while half-metallic for the Ag-ZQ polymer (see Table S1).

B. Density of states

In order to understand in more detail the magnetic properties of TM-ZQ polymers, we plot the density of states (DOS) and the atomic site projected density of state (PDOS) in Figures 2 and3. The density of states is an important physical quantity that is very often calculated in solid state physics. We now analyze the electronic structure, i.e., the total densities of states (DOS), of the ground state solution for the TM-ZQ polymer network more in detail to understand the behaviors of TM-ZQ polymers. The total DOS for these polymers are shown in Figure 2

C. Magnetic anisotropy energy

The electronic structures of TM-ZQ are completed by SOC calculations. The initial directions of magnetization are along the x, y, z orientations, where the chains run along the x-direction and the z-axis is perpendicular to the plan of the polymers. For the TM-ZQ polymer, three different directions of magnetization are considered in the SOC calculations. The energy difference for the magnetizations lying in the x-and z-directions defines the magnetic anisotropy energy (MAE) which for the magnetic chains with Mo, Ru, Rh, Re and Ir takes rather important values of 5.37, 4.80, 143.25, 93.45 and 64.56 meV/atom (see Table 2). We calculated the FM and AFM chains at zero temperature. It is known that the isotropic Heisenberg model allows no long-range order at any finite temperature in one or two dimensions. However, a large MAE can stabilize long-range magnetic order in low dimensions. Furthermore, a large MAE would lead to interesting magnetic properties even above a hypothetical long-range ordering temperature. The local magnetic centers react then independently of each other to a magnetic field. But for such large MAE's as reported here, there might exist a magnetic blocking temperature and magnetic hysteresis due to the large energy barriers between different directions of the local magnetic moments. As summarized in Table 2, the easy axis of TM-ZQ follows the direction of the lowest energy. In addition to the local spin moments M d , Table 2 shows also the atomic orbital magnetic moments M L . The SOC generally reduces the local spin moments by a few percent in 4d and 5d but induces a large orbital magnetic moment in the polymers Ru-, Rh-, and Ir-ZQ; see Taking similar distances as in Fe-ZQ on Au(110) [START_REF] Simpson | Coverage-dependent interactions at the organics-metal interface: Quinonoid zwitterions on Au(111)[END_REF] one can translate a MAE of 65 meV (for Ir-ZQ) per magnetic ion into a surface (or interface) anisotropy K i of about 25 mJ/m 2 which would be an extremely high value.

It is instructive to compare with the MAE's of the 3d series TM-ZQ [START_REF] Denawi | Magnetic Polymer Chains of Transition Metal Atoms and Zwitterionic Quinone[END_REF]. There, indeed, all MAE's are smaller than 6 meV/atom with the only exception V-ZQ which reaches an easy plane anisotropy of about -22 meV/atom comparable to the values reported above for the 4d/5d series. But the reason is completely different: in V-ZQ there are two 3d orbitals (3d xz and 3d yz ) which are nearly degenerate and partially filled, leading to a large orbital moment and a large MAE. Therefore, a high MAE is possible despite a small SOC coupling constant for V. A similar degeneration among the d orbitals was not observed in the 4d/5d series of TM-ZQ.

CONCLUSIONS

We investigated systematically the structural, electronic, and magnetic properties of organometallic polymers which are built of zwitterion quinone (ZQ) molecules and 4d/5d transition metals. We used fully relativistic SGGA and SGGA+U calculations and found several compounds with surprisingly high MAE. Especially remarkable is the easy axis behaviour of Rh-ZQ and Re-ZQ with MAE's of 143.25 and 93.45 meV/atom in SGGA+U+SOC, respectively (41.28 and 108.1 meV/atom without U correction). Such high values can lead to high blocking temperatures below which one can expect a magnetic hysteresis. It would be very interesting to verify our prediction experimentally. And it could be the basis to use these materials in future devices of magnetic molecular memory.

There is great variation of magnetic properties in the investigated material family: (Re, Ir)-ZQ polymers are FMs, (Mo, Ru, Rh, Ag)-ZQ polymers are AFM, and non-magnetic are Pt-ZQ and Pd-ZQ polymers. We also analyse the electronic structure of all compounds, and find magnetic compounds, which combine local magnetic moments with wide bands characteristic to conjugated polymer compounds. There is a strong variation of MAE within the 4d and 5d series. Most variations are good understood but not all. They demand further investigations.

There is a certain influence of the U parameter on the electronic and magnetic properties, the gaps are reduced, and some insulating compounds become metals. But the MAE's do not strongly depend on U (with the exception of Ir-ZQ). So, the prediction of a significant MAE for Rh-ZQ and Re-ZQ makes them highly interesting candidate materials for practical spintronics applications.

Review B. 13 (1976)

Figure 1 :

 1 Figure 1: Geometrical structure of a free-standing polymer chain of TM atoms and zwitterionic Quinone (C: grey, O: red, N: blue, H: white, and TM: light blue).

Figure 2 :Figure 3 )Figure 3 :

 233 Figure 2: Total DOS of the 1D transition-metal (TM) atoms zwitterionic quinoidal polymer with the PAW-SGGA+U method, when TM is (a) Mo (b) Ru, (c) Rh, (d) Pd, (e) Ag, (f) Re, (g) Ir, (h) Pt. To analyse further the magnetic properties, we plot the projected density of states (PDOS) of the d orbitals of the TM atom. For the TM +2 configurations of TM-ZQ, the d xy orbital is totally

Table 1 :

 1 Distance between the central TM Atoms and the N Atoms (d TM-N , in Å), Distance between the central TM Atoms and the O Atoms (d TM-O , in Å), Lattice constant (a, in Ǻ), Exchange Energies (E ex =E FM -E AFM per TM atom, in meV), Total Magnetic Moments (M per TM atom, in ), Local Magnetic Moments of the d orbital at the TM atoms (M d per TM atom, in ), Local Magnetic Moments at the TM atoms (Mm per TM atom, in ), Energy BandGaps (spin up (E a ) and spin down (E b ), in eV), Total Energy gaps (Eg,, in eV), for 1D TMquinone polymer in the SGGA+U method for the TM-ZQ polymer.

			Mo	Ru	Rh		Pd	Ag	Re	Ir	Pt
	d TM-N	2.18	2.06	2.01	2.00	2.10	2.06	1.97	1.99
	d TM-O	2.15	2.06	2.03	2.03	2.20	2.03	1.99	2.02
	a (Ǻ)	16.48	16.03	15.797	15.74	16.3	16.11	15.75	15.74
	E ex		48.950	4.663	19.229	-	14.546	-3.463	-3.468	-
	M		4	2	1		0	1	3	1	0
	M d		3.357	1.725	0.863	0	0.473	3.084	1.125	0
	M m		3.493	1.816	0.873	0	0.471	3.250	1.143	0
	E a		0.35	0.63	0.85	0.76	1.2	0.35	0.41	0.41
	E b		1.2	0.69	0.94	0.76	0.35	1.07	0.31	0.41
	E g (FM)	0.43	0.63	0.85	0.76	0.35	0.08	0.05	0.41
	E g (AFM)	0.88	0.69	1		-	0.56	0.37	0.00	-
	the	exchange energy				which is the difference between the
	antiferromagnetic and ferromagnetic states. The two polymers (Re, Ir)-ZQ have ferromagnetic
	ground states (	is negative), with an energy difference between the FM and AFM states of
	-3.463 and -3.468 meV/TM respectively. Further, the polymers (Mo, Ru, Rh, Ag)-ZQ have an
	antiferromagnetic ground states (	is positive) which is more stable than the FM state by an
	energy difference of 48.950, 4.663, 19.229 and 14.546 meV per TM atom, respectively, as
	summarized in

Based on the most stable structures of TM-ZQ, we systematically study their magnetic and electronic properties. The corresponding magnetic moments are listed in Table

1

, as well as

Table 1 .

 1 We also give in Table1the band gap values of all compounds which

	are all direct band gap semiconductors within the SGGA+U density functional. The band gaps
	of Re-ZQ and Ir-ZQ polymers are very small (0.08 and 0.05 eV respectively). The physical
	parameters which were calculated by SGGA+U are confirmed by the SGGA method (without
	U, see Table

Table 2 .

 2 We verified the large MAE's calculated with SGGA+U+SO by fully relativistic SGGA calculations, or SGGA+SO (see TableS2in the supplementary material). In general, the MAE's without the U correction confirm the results summarized in Table2. First of all, they confirm the easy axis behaviour of most compounds. The only exception is Ag-ZQ for which the MAE changes sign between SGGA+SO and SGGA+U+SO, but the magnitude of MAE is very small in that case since we have just one hole in the d shell. Interesting are the large MAE's for Re-ZQ (108.1 meV/atom) and Rh-ZQ (41.8 meV/atom) confirming the large values of Table2where the U correction is included. The MAE of Ir-ZQ reduces from 64.6 meV/atom (with U) to 11.3 meV/atom (without U) but even 11.3 meV/atom is a rather large MAE for d ions. For Mo-ZQ and Ru-ZQ, the SGGA+SO values are even larger than in SGGA+U+SO.The two organo-metallic polymers Re-ZQ and Ir-ZQ have a ferromagnetic nearest neighbour exchange, a very small gap and a large MAE value. Like Fe-ZQ they could probably be deposited on metallic surfaces and build well-ordered monolayers. If they also would show ferromagnetic long range order they could eventually serve in magnetic tunnel junctions.

Table 2 :

 2 Magnetic anisotropy energy (MAE, per TM atom in meV), total magnetic moments M i , local magnetic moments of the d orbitals at the TM site M d (i) and orbital moments (M L (i)per TM atom in ) for selected 1D TM-ZQ polymers of the 4d and 5d series. The data result from SGGA+U calculations including SO coupling with the magnetic moments directed towards the i={x,y,z} direction.

		Mo	Ru	Rh	Pd	Ag	Re	Ir	Pt
	MAE (E x -E z )	5.37	4.80	143.25	-	-0.47	93.45	64.56	-
	M x	3.99	1.63	0.56	-	0.99	2.74	0.78	-
	M y	3.99	1.92	0.32	-	0.99	2.70	0.47	-
	M z	3.99	1.73	0.34	-	1	2.63	0.47	-
	M d (x)	3.348	1.411	0.470	-	0.470	2.883	0.976	-
	M d (y)	3.344	1.651	0.264	-	0.470	2.809	0.636	-
	M d (z)	3.343	1.496	0.285	-	0.473	2.779	0.410	-
	M L (x)	-0.056	0.425	0.632	-	0.036	0.047	1.144	-
	M L (y)	-0.056	0.600	0.217	-	0.032	0.024	0.340	-
	M L (z)	-0.082	0.099	0.143	-	0.114	-0.052	0.093	-

  5188-5192. doi:10.1103/PhysRevB.16.1748.

ACKNOWLEDGEMENTS

This work was supported by the Computer resources of the Centre Informatique National de l'Enseignement Supérieur (CINES), Project No. A0020906873 and the High Performance Computing (HPC) resources of Aix-Marseille University financed by the project Equip@Meso (ANR-10-EQPX-29-01). We thank M. Koudia for helpful discussions.

AUTHOR INFORMATION

Corresponding Author *E-mail: hassan.denawi@gmail.com (H.D.)

Conflict of interest

The authors declare that they have no conflict of interest.