Messaadi Mahassene

Kaid M'hamed

Ammam Abdelkader

REMOVAL OF CONGO RED BY A SYNTHESIZED LAYERED DOUBLE HYDROXIDE

Keywords: LDH, Congo red, Adsorption, Regeneration

The colored wastewaters from the textile industries may show toxic or carcinogenic eff ects on the organism, when discharged into the rivers and lakes, which are changing their biological life.

Therefore, in this study, the removal of Congo red by Zn-Al-SO 4 layered double hydroxide has been achieved. Zn-Al-SO 4 material was prepared by a facile co-precipitation method at constant pH, resulting in a suitable material for the adsorption of Congo red dye. The structure and morphology of the Zn-Al-SO 4 adsorbent were investigated using XRD, FT-IR, and BET and MEB techniques. Layered double hydroxides (LDHs) have been widely investigated in a wide range of applications in health, in the pharmaceutical industry and in the material of biotechnology industries. In this work the synthesis of Zn/Al double layered hydroxides by chemical co-precipitation method (molar ration 3). The samples were characterized and confirmed by X-ray diffraction (XRD), BET analysis, and infrared spectroscopy. A series of experiments was then carried out to study the influence on the adsorption capacity of certain parameters such as the amount of the adsorbent, the pH, the contact time, the initial dye concentration, the ionic strength effect and temperature. All the results obtained show that the adsorption kinetics of the dye on our material is well described by the second order model. The adsorption isotherms of the adsorbent / adsorbate systems studied are satisfactorily described by the Langmuir mathematical model. The intra-particle scattering model confirms the physisorption phenomenon Negative values of free energy prove that the phenomenon is physisorption. On the other hand, the thermodynamic study revealed that the adsorption is spontaneous and exothermic. The recovery of the material and reuse shows that the property of being able to regenerate is remarkable.

INTRODUCTION

Increasing environmental pollution caused by the rapid growth of technology is one of the serious concerns of global societies [START_REF] Boulaiche | Removal of heavy metals by chitin: equilibrium, kinetic and thermodynamic studies[END_REF]. The growth of civilization and industrialization has led to the destruction of water resources [START_REF] Järup | Hazards of heavy metal contamination[END_REF][START_REF] Niazi | Spectrophotometric determination of mercury in water samples after cloud point extraction using nonionic surfactant Triton X-114[END_REF][START_REF] Sajid | Removal of heavy metals and organic pollutants from water using dendritic polymers based adsorbents: a critical review[END_REF]. More than 10,000 dyes have been widely used in textile, rubber, paper, plastics, leather, cosmetic, pharmaceutical, and food industries, which generated huge volume of wastewater every year [START_REF] Tahmasebi | Application of mechanosynthesized azine-decorated zinc (II) metal-organic frameworks for highly efficient removal and extraction of some heavy-metal ions from aqueous samples: a comparative study[END_REF][START_REF] Luo | Removal of antimonite (Sb (III)) and antimonate (Sb (V)) from aqueous solution using carbon nanofibers that are decorated with zirconium oxide (ZrO2)[END_REF][START_REF] Qu | Effective and regenerable Ag/graphene adsorbent for Hg (II) removal from aqueous solution[END_REF][START_REF] Zhao | Treatment of lead contaminated water by a PVDF membrane that is modified by zirconium, phosphate and PVA[END_REF]. Annually the production of dyes is estimated at around 1.6 million tons which 10-15% of this amount is discharged to sewage. Exposure to these species causes skin irritation, respiratory problems, and cancer [START_REF] Carro | Interaction of heavy metals with Capretreated Sargassum muticum algal biomass: Characterization as a cation exchange process[END_REF]Wang et al.,2016). In addition, the produced colored wastewater from the dye manufacturing and textile industries e ff ects on biological life when is discharged into rivers and lakes without wastewater treatment. The long history of red coloring elements are figured out in two periods separated by a approximately recent year, 1858, marked by the invention of the first synthetic red dye. Previously, in the 19th century, foods were colored with natural red dyes, including cochineal red, beetroot red, madder, and metal salts as well. Yet, it had been merely following the invention of red dyes resulting from chemical synthesis that their use within the food industry spread from 1880. If the fields of application of red dyes remain very varied (the textile industry remains the main sector making use of these compounds), we note that since a few decades, the global food industry uses more and more significant amount of natural or artificial red dyes, in particular in confectionery, desserts, drinks, but also in cold cuts, fruits and vegetables and sugar. Otherwise, these compounds have a toxic character causing significant disturbances in the various natural mechanisms existing in the flora (selfpurification power of watercourses, inhibition of the growth of aquatic plants, etc.) and in the fauna (destruction of a category of fish, microorganisms, etc.) [START_REF] Alipour | Optimising the basic violet 16 adsorption from aqueous solutions by magnetic graphene oxide using the response surface model based on the Box-Behnken design[END_REF]) and can persist a long time in this environment. Besides to the toxicity and health risks of these pollutants, the wide kinds of red dyes usually have complex structure and non-biodegradable molecular due to their high molecular weight, so that is very hard to eliminate from environments [START_REF] Gong | Efficient removal of heavy metal ions from aqueous systems with the assembly of anisotropic layered double hydroxide nanocrystals@ carbon nanosphere[END_REF]. These are the reasons that the effective removal of these pollutants from water and their treatment before discharging into the water stream is a great subject of research nowadays [START_REF] Sajid | Removal of heavy metals and organic pollutants from water using dendritic polymers based adsorbents: a critical review[END_REF]. Diff erent methods are used to remove organic pollutants, such as dyes from wastewater [START_REF] Agarwal | Heavy metal removal from wastewater using various adsorbents: a review[END_REF], like ozonation [START_REF] Joshi | Synthesis and characterizations of polythiophene-Al2O3 based nanosorbent and its applications in the removal of Pb2+, Cd2+ and Zn2+ ions[END_REF], coagulation [START_REF] Chen | Selective removal of Cu (II) ions by using cation-exchange resin-supported polyethyleneimine (PEI) nanoclusters[END_REF], electrocoagulation [START_REF] Wang | Removal of cadmium ( Ⅱ) from aqueous solution by granular activated carbon supported magnesium hydroxide[END_REF] photocatalytic removal [START_REF] Ferrandon | Influence d'alcools aliphatiques sur l'adsorption d'acides aminés sur charbon actif[END_REF], membrane filter [START_REF] Xu | A new combined process for efficient removal of Cu (II) organic complexes from wastewater: Fe (III) displacement/UV degradation/alkaline precipitation[END_REF], and adsorption (Mnasri-Ghnimi et al., 2015; [START_REF] Naseem | Adsorptive removal of heavy metal ions using polystyrene-poly (N-isopropylmethacrylamide-acrylic acid) core/shell gel particles: adsorption isotherms and kinetic study[END_REF]. Among various methods, and using a suitable adsorbent, adsorption is widely used as a reasonable technique due to its simplicity, ease of use, low cost, high efficiency and the availability of the materials [START_REF] Lahreche | Application of Activated Carbon Adsorbents Prepared from Prickly Pear Fruit Seeds and a Conductive Polymer Matrix to Remove Congo Red from Aqueous Solutions[END_REF]. The objective of this work is to prepare and characterize (ZnAl-SO 4 ) by co precipitation method. This green adsorbent has been used to study the effects of different parameters on Congo red adsorption, such as pH of the aqueous solution, quantity of hybrid material, initial concentration of heavy metal ion, adsorption equilibrium time, stirring speed, salts effects and temperature.

MATERIAL AND METHOD Experimal Materials

The reagents used in this study were Zn Al-SO 4 synthesized in the laboratory (Kadari et al., 2017;[START_REF] Mahassene | The intercalation of Zn/Al HDL by diaminododecylphosphonic acid: synthesis and properties of absorption of cadmium[END_REF], the raw materials ZnSO 4 (Riedel De Haen), Al 2 (SO4) 3 (Riedel De Haen), nitric acid (64% Riedel De Haen), Congo red (98%, Sigma Aldrich), NaOH (99%, Riedel De Haen), Na 2 SO 4 , KNO 3 (99%, Riedel De Haen), UV/Vis spectra were obtained using a UV/Visible double beam type'' OPTEZEN 3220'' spectrometer to analyze the concentration of congo red in aqueous solution, a digital pH meter Consort C863 type.

Methods

Synthesis of the ZnAlSO 4

Zn-Al-SO 4 was prepared using a co-precipitation method similar to that described by Miyata [START_REF] Miyata | The Syntheses of Hydrotalcite-Like Compounds and Their Structures and Physico-Chemical Properties-I: The Systems Mg2+-Al3+-NO 3-, Mg2+-Al3+-Cl-, Mg2+-Al3+-ClO 4-, Ni2+-Al3+-Cl-and Zn2+-Al3+-Cl[END_REF] with a Zn/Al ratio of 3 a. In a reactor, 150 mL of an aqueous solution of metallic sulphates (0.75 M in ZnSO 4 and 0.25 M in Al 2 (SO 4 ) 3 was added drop wise, at 25°C and under magnetic stirring. 1.5 M NaOH was added to maintain the pH of the reaction mixture constant at 10.5 ± 0.1. When the reaction was achieved, the resultant gel was made under reflux at 60-70°C during 15 hours to permit crystal growth. After ripening, at room temperature for 24 h in mother liquors, the mixture was filtered, washed several times then dried, a white precipitate was obtained.

Adsorption experiment

The LDH of type ZnAl-SO 4 with a molar ratio of 3 has been prepared by the method of co-precipitation described by Miyata. In this work, we have chosen the method of spectrophotometry molecular absorption as colorful indicator for the determination of Congo red at λ max = 500 nm [START_REF] Guiza | Kinetics study of Congo red dye adsorption from aqueous solutions onto bentonite[END_REF] Several mixtures of solid phase (ZnAl-SO 4 ) and aqueous phase of the dye are subjected to magnetic stirring. After separating the two phases by centrifugation, the amount of Congo red remaining in the aqueous phase has been determined by spectrophotometry visible UV/V (double beam type'' OPTEZEN 3220'' spectrometer), The amounts of Congo red adsorbed are calculated from the following equation ( 1):

(1)

Where C0 and Ce are the initial and equilibrium concentrations of Congo red (mg L -1 ) respectively; m the amount of adsorbent (g) and V the volume of solution (L). The efficiency percentage (Y) of dye removal was calculated as follows (equation ( 2):

(2)

Characterization techniques XRD Analyses

Powder X-ray diffraction (XRD) data were collected on a PANAnalytical X'Pert Pro diffractometer in reflection mode at 40 kV and 40 mA using Cu Kα radiation (1.5405980 Å). Scans were recorded from 5°≤2θ≤100° with varying scan speeds and slit sizes. Samples were mounted on stainless steel sample holders.

FTIR Analyses

Fourier transform infrared (FTIR) spectra were recorded using a Perkin Elmer 16 PC-FTIR equipped with a thermostat to maintain the temperature of the sample at 25.0 ± 01 °C on KBr pellets in the range of 4000-400 cm -1 .

BET analyses

Specific surface areas and pore size were analyzed using the Brunauer-Emmett-Teller (BET) method. The samples were measured from the N2 adsorption and desorption isotherms at 77 K collected from a Quantachrome Autosorb-6 surface area and pore size analyzer.

UV/V analyses

A UV/Visible double beam type'' OPTEZEN 3220'' spectrometer was used for UV-vis spectral measurements. The accuracy and the precision of the method on the instrument were determined by measuring the absorbance at 500 nm against the prepared standards in the concentration range of (10 -3 to 10 -6 mol/l) [START_REF] Chen | Selective removal of Cu (II) ions by using cation-exchange resin-supported polyethyleneimine (PEI) nanoclusters[END_REF]. A standard curve for UV-vis measurement demonstrated high degree of accuracy with a coefficient of regression, R 2 =99.7 %, was used in the calculation of unknown UV-vis concentrations from absorbance readings.

RESULTS AND DISCUSSION

X-ray diffraction

The powder XRD pattern of Zn-Al-SO 4 is presented in figure (Figure 1). The diffraction peaks corresponding to the LDH phase are observed at the 2θ position 12.4, 20, 27.5, 37. Morphology SEM images of ZnAl-SO 4 -LDH are presented in Figure 4. The images clearly show that the sample existing lamellar particles with hexagon layer structure, which is the typical structure of the hydrotalcite-like material.

Effect of various parameters on Congo red dye removal efficiency

Different parameters have been investigated to study the retention yield of Congo red dye

Effect of stirring speed

We observe that the adsorption rate depends on the stirring speed. On the other hand, it is fast for t < 10 min, then it remains constant (Figure 5). The agitation speed affects the resistance to the transfer of material outside the clay particle, this with accordance with literature (Guiza et al., 2019).

Eff ect of pH on the efficiency of Congo red removal

The removal of Congo red has been investigated in the pH range from 3.5 to 12 by adding nitric acid HNO 3 in different solutions with an adsorbent dose of 0.025 g. According to the results shown in Figure 8, we noted that the removal of dye by the ZnAl-SO 4 is much more important at pH=8 with best adsorption capacity of 23 mg/g. That is to say the more acidic (pH = 3.5) is the solution, the less dye are extracted. This is probably due either to competition between the ionic organic pollutant dye in solution and the ions released by the nitric acid or to the non-stability of the LDH in acidic solution. (

) (4) 3
From these two equations, the following expression was obtained:

(5)

The calculation of certain thermodynamic parameters is essential in determining the nature of the retention process. The equilibrium constant Kd, can be calculated from the following relation: [START_REF] Chen | Selective removal of Cu (II) ions by using cation-exchange resin-supported polyethyleneimine (PEI) nanoclusters[END_REF] The capacity of adsorption (q) of the studied heavy metals by the LDH is determined by the following relation The thermodynamics parameters are listed in Table 2, where ΔG 0 values were obtained by plotting lnK against 1/T for the Langmuir data. The results of the thermodynamic study show that in (Table 2):

The negative ∆H°value indicates the exothermic nature of adsorption of the Congo red and was less than 10 kcal/mol indicating physical adsorption. The ∆G °values obtained were negative meaning that the adsorption was spontaneous and positive value of ∆S° suggested a slight increase in the randomness at the solid/solution interface during adsorption of dye on to the adsorbent.

Adsorption isotherm

The study of the adsorption isotherm is fundamental and plays an important role in determining the maximal capacity of adsorption. In order to adapt for the considered system, an adequate model that can reproduce the experimental results obtained, Langmuir, Freundlich and Dubinin-Radushkevich models have been considered.

In Figure 10, the quantities adsorbed according to the concentrations of the aqueous solution in equilibrium shows that the isotherm is of L-type which is obtained when the adsorption of the solvent is low and when the molecules are not oriented vertically but rather flat. In this case, the adsorption of the solute on the solid takes place in a monolayer.

Langmuir isotherm model

Mathematically, Linearization of Langmuir isotherm model is defined as:

(

) 8
Where qm is the monalayer capacity of the adsorbent, and KL is the Langmuir adsorption constant. qm and KL can be determined from the slope and intercept, respectively from the plot of 1/qe versus 1/Ce (figure 11).

Freundlich isotherm model

This model describes systems where the adsorption is done on heterogeneous surfaces with interactions between the adsorbed molecules; the linear form is given by: log q e =log Kf + 1/n log C e

The constant (Kf), due to the bond energy, and the heterogeneity factor (1/n) which measures the deviation from the linear part are determined from the plot log qe versus log C e . The experimental results obtained at ambient temperature are presented in Figure 12. The Freundlich equation ( 8) is applied for the adsorption of the Congo red dye on our prepared LDH.

Dubinin -Radushkevich Isotherm

Plotting ln q e = f (ε 2 ) gives us the line of the slope β, and qs is the y-intercept. The results obtained are gathered in Table 4 and Figure 13 shows a linear adjustment.

Kinetic study

The kinetic adsorption of Congo red dye by ZnAl-SO4 is studied at different contact times in the batch system and has been tested with pseudo first order pseudo second order and Intra particles diffusion model.

Pseudo-first-order model

The fitting results based on the pseudo-first order kinetic model are shown in Figure 14. The values for the kinetic models are shown in Table 5 [START_REF] Guiza | External mass transport process during the adsorption of fluoride from aqueous solution by activated clay[END_REF] Where qe and qt refer to the amount of mercury adsorbed per unit weight of Zn-Al-SO4 respectively at equilibrium and at any time:

Pseudo-second-order

Using the pseudo-second-order equation also, the adsorption process may be described. The differential equation is the following:

The results obtained (Figure 15) show that the pseudo-second-order model is more suitable to describe the kinetics of the retention of Congo red toxic dye by ZnAl-SO 4 .

Figure 15. Kinetics of pseudo second order for Congo red retention

Intra particles diffusion model

Solute transfer is generally characterized by either the external mass transfer step or intra particle diffusion or both. To study the existence of intraparticle diffusion during adsorption, the most widely used equation ( 12) is that given by (Abdel-Ghani et al ., 2016). From figure 16, we distinguish a double linearity reflecting the existence of two stages: firstly, the high adsorption rate, resulting from the diffusion of the outer film and through the boundary layer of the outer surface of the adsorbent. Second, the adsorption process is controlled by intra-particle diffusion, which is characterized by slowing the rate of adsorption, and it is known as the rate-limiting step. These two phases are involved simultaneously during adsorption.

Material Regeneration

The economic and environmental aspect of the use of adsorbent materials makes it important to reuse anionic clays, given their low cost and their ability to regenerate. After an adsorption with a yield R of almost 100%, desorption has been carried out using a solution of NaOH of known concentration C = 5×10 -3 M which is chosen after many trials using different NaOH concentrations.

The material used has been used for seven cycles showing best performance at each cycle.

CONCLUSION

Our prepared material (ZnAl-SO4) was found very effective adsorbent for the removal of Congo red dye from aqueous solutions. The synthetic method used to synthesize this adsorbent by chemical co-precipitation method is highly efficient and very low cost. This compound adsorbent showed that elimination can reach 100% under the optimal conditions and can be used at least for seven cycles. This can allow us to provide this adsorbent (Zn-AlSO4) for the waste water treatments.

  2, 42. and 47.1 having the respective d-values 10.94, 5.47, 3.66, 2.66 and 2,.18 and the respective miller indices are 003, 006, 015, 018, 110 and 113 characteristic of an LDH phase with a basal spacing of 7.98 A° and an inter metallic distance of 1.52 A°. The XRD patterns corresponding to the rhombohedral symmetry (space group, R-3m) are indexed in an hexagonal lattice. The (a) and (c) lattice parameters are estimated using (110) and (003) reflections, respectively.

Figure 1 .

 1 Figure 1. DRX pattern of ZnAl-SO4

Figure 2 .Figure 3 .

 23 Figure 2. FTIR spectrum of Zn-Al-S0 4

Figure 4 .

 4 Figure 4. SEM image of ZnAl-SO4

Figure 5 .

 5 Figure 5. Effect of stirring speed on the retention of Congo red

Figure 6 .

 6 Figure 6. Adsorbent amount effect on the Congo red retention

Figure 7 .

 7 Figure 7. Concentration effect on the retention of Congo red. T = 25C°, m = 0,025 g

Figure 8 .

 8 Figure 8. Effect of initial solution pH on the dye retention. T=25C°, m=0,025 [dye] = 55,74 mg/L

( 7 )

 7 qe: the capacity of the adsorption at equilibrium C0: initial concentration of Cong red in mg/L Ce: the concentration of the Cong red at equilibrium in mol/lL V: the volume of the solution of the treated Cong red (10 mL) M: the molar mass of Congo red dye m: the mass of ZnAl-SO4 (0.025 g) R: the perfect gas constant (R= 8.314 J mol -1 K -1 ) Kd: the distribution coefficient between the aqueous and solid phases

Figure 9 .

 9 Figure 9. lnKads fonction of 1/T for the congo red retention

Figure 10 .

 10 Figure 10. Graphical representation of the isotherm of the Congo red aqueous solution

Figure 11 .

 11 Figure 11. Langmuir linearization

Figure 12 .

 12 Figure 12. Freundlich linearization

Figure 13 .

 13 Figure 13. Dubinin-Radushkevich model for dye retention

Figure 14 .

 14 Figure 14. Kinetics of pseudo first-order for congo red retention

Figure 16 .

 16 Figure 16. Diffusion model of Congo red dye particles

Figure 17 .

 17 Figure 17. Recycling and reuse of ZnAl-SO4 towards adsorption of Congo red

Table 1

 1

	Zn/Al-LDH	62.5	0.29	18.6

. Pore properties and specification surface area of ZnAl-SO4 BETsurface area(m 2 g -1 )

Pore volume(cm 3 g -1 ) Average pore size(nm)

Table 2 .

 2 Adsorption Thermodynamic parameters at 303K

	Thermodynamic Parameters	-1 ΔG° (kcal. mol )	-1 ΔH°, (kcal.mol )	ΔS° (cal.m	-1 .kcal	-1 )
	Values	-1.55	-0.088		4.83

Table 3 .

 3 Langmuir and Freundlich constants for the elimination of Congo red on ZnAl-SO4

	Langmuir model	Freundlich model
	Qm (mg/g)	KL (l/mg)	1/n	KF
	22.22	0.52	-0.023	21.75

Acknowledgements. This work was supported by Dr Moulay Tahar University, Saida, Algeria and "PunchOrga" Network (Pole Universitaire Normand de Chimie Organique), the "Ministère de la Recherche et des Nouvelles Technologies," CNRS (Centre National de la Recherche Scientifique).


Yarkın Aybars Çetin

Benjamí Martorell

Francesc Serratosa

Neus Aguilera-Porta

Monica Calatayud

Analyzing the TiO 2 surface reactivity based on oxygen vacancies computed by DFT and DFTB methods

Keywords: nanoparticles, oxygen vacancy, titania, DFT+U, VASP, SCC-DFTB

Titanium dioxide is a key material in many fields, including technological, industrial and biomedical applications. Many of these applications are related to the surface reactivity of TiO2 and involve its reducibility properties. Recently titania has been banned as a food additive due to its (nano)toxicity, and the release of reactive oxygen species plays a crucial role in many toxicological mechanisms. Determining chemical descriptors that account for the extension of reduction is necessary to understand such processes and necessary for predicting the reactivity of an unknown system. In the present work, we compute a set of chemical descriptors for selected surfaces of anatase and rutile TiO2. The aim is twofold: we want to provide chemically meaningful information on the surface reactivity, and benchmark the descriptors for two ab initio schemes. To do so, we compute the oxygen vacancy formation energy, and the corresponding electronic structure, in four slab models with two different computational schemes (DFT+U and DFTB). In this way, we characterize the robustness of the dataset, with the purpose of scaling up to more realistic model systems such as nanoparticles or explicit solvent, which are too computationally demanding for state-of-the-art density functional theory approaches.

Introduction

TiO2 is a versatile material with applications in a variety of technological fields. [START_REF] Burello | A Theoretical Framework for Predicting the Oxidative Stress Potential of Oxide Nanoparticles[END_REF][START_REF] Liu | Development of Structure-Activity Relationship for Metal Oxide Nanoparticles[END_REF][START_REF] Xia | Comparison of the Mechanism of Toxicity of Zinc Oxide and Cerium Oxide Nanoparticles Based on Dissolution and Oxidative Stress Properties[END_REF][START_REF] Ropers | Titanium Dioxide as Food Additive[END_REF] From the construction to the food additive industry, TiO2 has been used widely with a decreasing trend of use after recent regulations. [START_REF] Blaznik | Use of Food Additive Titanium Dioxide (E171) before the Introduction of Regulatory Restrictions Due to Concern for Genotoxicity[END_REF] TiO2 is also an excipient in medicinal products and a pigment and thickener component in personal care products, [START_REF] Weir | Titanium Dioxide Nanoparticles in Food and Personal Care Products[END_REF][START_REF] Peters | Characterization of Titanium Dioxide Nanoparticles in Food Products: Analytical Methods To Define Nanoparticles[END_REF] and is found in mineral sunscreen products as UV filters. 8,9 Because of its wide range of applications and presence in everyday market products, TiO2 brings a reasonable further demand for knowledge related to the toxicology perspective. [START_REF] Elstner | Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties[END_REF] Among many different toxicological mechanisms, the oxidative potential of metal oxides is a key parameter to enable the prediction of Reactive Oxygen Species (ROS) formation, which, in turn, is related to the cytotoxicity of materials. In particular, statistical studies between the relationship of physicochemical properties of TiO2 (as well as other oxides or metal systems) with the toxic effect of the materials [START_REF] Rybińska-Fryca | Structure-Activity Prediction Networks (SAPNets): A Step beyond Nano-QSAR for Effective Implementation of the Safe-by-Design Concept[END_REF][START_REF] Mikolajczyk | Nano-QSAR Modeling for Ecosafe Design of Heterogeneous TiO2-Based Nano-Photocatalysts[END_REF][START_REF] Beyer | Unimolecular Reactions of Dihydrated Alkaline Earth Metal Dications M2+(H2O)2, M = Be, Mg, Ca, Sr, and Ba: Salt-Bridge Mechanism in the Proton-Transfer Reaction M2+(H2O)2 → MOH+ + H3O+[END_REF][START_REF] Boukhvalov | Development of Theoretical Descriptors for Cytotoxicity Evaluation of Metallic Nanoparticles[END_REF][15] have been extensively performed using different approximations, demonstrating that correlation. The relevance between ROS formation and the oxidative potential, which is directly related to the bandgap value (defined as the difference between the conduction and valence band energy levels), has been shown in previous works. Since the formation of oxygen vacancies in TiO2 has been connected to the reducibility, such a process has been extensively used to explain the properties of bulk and surfaces. 16- 22 Also related to reducibility, TiO2 has been reported in the literature to be phototoxic and dependent on several factors such as morphology [START_REF] Tong | Effects of Material Morphology on the Phototoxicity of Nano-TiO2 to Bacteria[END_REF] , size [START_REF] Xiong | Size of TiO2 Nanoparticles Influences Their Phototoxicity: An in Vitro Investigation[END_REF] or corona proteins [START_REF] Garvas | Protein Corona Prevents TiO2 Phototoxicity[END_REF] to cite a few. Understanding the mechanisms behind the light induced toxicity would require sophisticated theoretical methods and will not be addressed in the present work. Instead, the focus will be put on the characterization of oxygendeficient TiO2 systems to investigate the potential use of computational features as toxicity descriptors in realistic models.

Theoretical nanoparticles (NP) related studies are available for high-level electronic correlation and electronic structure methods. Some works report the comparison of different computational schemes such as density functional theory (DFT) based tight binding DFTB and Hybrid-DFT [START_REF] Selli | Modelling Realistic TiO2 Nanospheres: A Benchmark Study of SCC-DFTB against Hybrid DFT[END_REF][START_REF] Spiegelman | Density-Functional Tight-Binding: Basic Concepts and Applications to Molecules and Clusters[END_REF] . These studies enable a well-detailed NP property analysis, explaining theoretical insights of surface physics and chemistry by covering the different levels of solution methodologies, focusing on geometry, electronic structure and surface reactivity of non-reduced particles [START_REF] Rousseau | Theoretical Insights into the Surface Physics and Chemistry of Redox-Active Oxides[END_REF] . To the best of our knowledge, the reducibility and in particular the oxygen vacancy formation energy calculation, have not been addressed yet as regards the use of DFTB vs DFT schemes.

The state-of-the-art methodology for studying accurately titania materials is based on the electronic structure. The electronic structure is obtained from implementing quantum methods such as DFT, combined with a set of plane waves basis and pseudopotentials in periodic boundary conditions. The simulation package VASP [START_REF] Kresse | Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set[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 | Ab Initio Molecular-Dynamics Simulation of the Liquid-Metal-Amorphous-Semiconductor Transition in Germanium[END_REF][START_REF] Kresse | Ab Initio Molecular Dynamics for Liquid Metals[END_REF] is widely used among the different available codes in the materials science community. Using the pure DFT approximation allows the fast evaluation of energy through Kohn-Sham orbitals for systems up to approximately 500 atoms. To tackle nanotoxicity phenomena, it would be desirable to calculate more realistic systems such as nanoparticles of few nanometers, or to explicitly include solvent and coating. Note that in practice, the study of surfaces is done by creating slabs repeated in 3D and containing a vacuum space to prevent interaction between successive slabs. This computational approach is computationally challenging for the DFT method since the required atoms are much larger than 500. The implicit solvent [START_REF] Mathew | Implicit Solvation Model for Density-Functional Study of Nanocrystal Surfaces and Reaction Pathways[END_REF] , although not ideal, is an improvement to evaluate solvation effects.

On the other side, DFTB theory allows efficient processing of a system containing up to 10.000 atoms, which makes it useful to study TiO2 properties and toxicity. In this case, the theoretical model size becomes comparable to the experimental samples. Currently, DFTB can be solved by the DFTB+ 10,33 code package. DFTB+ numerical solver is a quantum mechanics based, open-source code. In order to assess the scalability between DFT and DFTB methods, and the reliability of the electronic structure provided by DFTB, a set of representative descriptors must be compared.

In the present work, we compute the oxygen vacancy formation energy and electronic structure for selected TiO2 stoichiometric surfaces using the two methods: DFT and DFTB. The first one is based on the plane-wave DFT, and the second is based on the self-consistent charge DFTB (SCC-DFTB). We compare the formation of oxygen vacancy's energetic and electronic features as a descriptor of titania surface reducibility. Four different surfaces are considered, namely anatase (101), anatase (100), anatase (001) and rutile (110). Our final goal is to compare the results obtained by the two approaches so as to assess the reliability and robustness of the results. Two frequently used computational settings are chosen. The aim of this comparison is ultimately to answer the following questions: What may serve as a foundation for describing nanotoxicity from a modeler point of view? What are the similarities and differences in the results obtained by the two methodologies? Note that the focus in the present work is put on the stability of the reduced, defective TiO2 models, therefore the nature and stability of ROS is not addressed.

Computational Methods

DFTB in DFTB+

The calculations used in this study to determine the oxygen formation vacancy process rely on the self-consistent charge (SCC) DFTB technique and the DFTB+ software code [START_REF] Elstner | Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties[END_REF][START_REF] Aradi | a Sparse Matrix-Based Implementation of the DFTB Method †[END_REF] . The key aspects of this concept originated from traditional Density Functional Theory (DFT), and the fundamentals of the methodology can be found in the literature [START_REF] Elstner | Self-Consistent-Charge Density-Functional Tight-Binding Method for Simulations of Complex Materials Properties[END_REF][START_REF] Koskinen | Density-Functional Tight-Binding for Beginners[END_REF][START_REF] Seifert | Tight-Binding Density Functional Theory: An Approximate Kohn-Sham DFT Scheme †[END_REF] . The Slater-Koster files, which contain all parameterized Vrep of the different atoms in the simulations and symmetry operation transformations for SCC-DFTB simulation experiments, are those of the TIORG set for Ti interactions with O, which can be found on the www.dftb.org web page with all details. TIORG set of interactions was adjusted using hybrid DFT functional B3LYP with the Pople-style double-ζ 6-31G(d). For the sake of comparison, calculations have also been performed on the MATSCI set of interactions (also available on www.dftb.org), which was adjusted from PW91 GGA DFT functional; results are displayed in Section 1 and 2 in the Supporting Information.

The total electronic energies were converged until a criterion of 10 -5 eV in the system's total energy (SCC Tolerance) was reached through SCC cycles. The included script component, hamiltonian, relies on the electronic optimization (spin-polarized) procedure of the desired system by recognising Slater-Koster libraries. Driver in the script is responsible for changing the geometry of the input structure during the calculation. Among the currently available other methods, the conjugated gradient option was employed for geometry optimisation with a convergence criterion of 10 -3 in the maximum force component. All simulations were performed in the framework of 3D periodic boundary conditions (pbc).

DFT in VASP

VASP evaluates an approximate solution to the many-body Schrödinger equation, either using DFT to solve the Kohn-Sham equations or using the Hartree-Fock (HF) approximation to solve the Roothaan equations. Hybrid functionals, which combine the HF method with density functional theory, are also implemented. Central variables like one-electron orbitals, electronic charge density, and local potential are represented in plane-wave basis sets in VASP. Electron-ion interactions are modelled using normconserving or ultrasoft pseudopotentials or the projector-augmented-wave technique.

VASP employs efficient iterative matrix diagonalisation approaches to find the electronic ground state, such as the residual minimisation method with blocked Davidson algorithms or direct inversion of the iterative subspace (RMM-DIIS).

In this work, the PBE Perdew-Berke Ernzerhof functional, spin-polarized, is employed [START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] . A planewave basis set of cutoff 400 eV described the valence electrons (H, Ti and O described by 1, 4 and 6 valence electrons, respectively), while the core electrons are represented by projector Augmented Wave (PAW) pseudopotentials [START_REF] Kresse | Ab Initio Molecular Dynamics for Liquid Metals[END_REF][START_REF] Blöchl | Projector Augmented-Wave Method[END_REF] . The 3d electrons of Ti were treated with the effective Hubbard onsite Coulomb interaction parameter (U' = U-J = 4 eV ) [START_REF] Dudarev | Electron-Energy-Loss Spectra and the Structural Stability of Nickel Oxide: An LSDA U Study[END_REF] . The Brillouin zone was sampled with a ~0.02 Å -1 spaced k-point grid. This computational setting has been successfully used to describe TiO2 rutile and anatase systems [START_REF] Wei | Understanding the Role of Rutile TiO 2 Surface Orientation on Molecular Hydrogen Activation[END_REF][START_REF] Wei | The Subsurface Diffusion of Hydrogen on Rutile TiO2 Surfaces: A Periodic DFT Study[END_REF] . The conjugate gradient algorithm was used for the geometry optimisation, with energy convergence criteria of 10 -3 eV for the ionic and 10 -4 eV for the electronic loop. All the atoms were allowed to relax. The same slabs, used in previous studies [START_REF] Wei | The Subsurface Diffusion of Hydrogen on Rutile TiO2 Surfaces: A Periodic DFT Study[END_REF][START_REF] Wei | Hydrogen Activation on Anatase TiO2: Effect of Surface Termination[END_REF][START_REF] Tielens | Characterization of Amorphous Silica Based Catalysts Using DFT Computational Methods[END_REF] , were used for VASP and DFTB calculations. They are available as Supporting information in the references cited.

Oxygen Vacancy Formation Energy

The energy cost associated to remove an oxygen from the stoichiometric systems EVO is calculated as:

2 1 2 2 2 2 1 2 1 2 1 2 n n n n VO MO O MO MO MO O E E E E - - = + = + - (1)
where MO2n-1 and MO2n correspond to oxygen deficient and stoichiometric system respectively, and O2 the oxygen gas-phase molecule calculated in a box. The lower the EVO value, the easier the reduction. The O-deficient systems are open-shell because the removal of a neutral oxygen leaves two electrons trapped in titanium sites. Thus, two titanium Ti 4+ sites become formally Ti 3+ . The structures were calculated for singlet (Nα-Nβ=0) and triplet (Nα-Nβ=2) configurations, the results are presented in Table 1 and S1 and discussed below.

Slab Models

Periodic 3D slabs cut out from the bulk rutile and anatase have been used. A vacuum of at least 12 Å was used to prevent interaction between slabs. The bare slabs are available in the supporting information of references [START_REF] Wei | The Subsurface Diffusion of Hydrogen on Rutile TiO2 Surfaces: A Periodic DFT Study[END_REF][START_REF] Wei | Hydrogen Activation on Anatase TiO2: Effect of Surface Termination[END_REF][START_REF] Tielens | Characterization of Amorphous Silica Based Catalysts Using DFT Computational Methods[END_REF] . Figure 1 displays the slabs models (supercells) used for the anatase (001), ( 100) and (101) surfaces as well as the rutile (110) slab. It also includes the labeling for each generated vacancy. The simulation cell of the anatase (001) plane had lattice parameters of a=b=7.74 Å and 12.73 Å in thickness. Energies were evaluated on a 3x3x1 reciprocal space grid. The systems consisted of 72 (24Ti, 48O) atoms. For anatase (100), parameters were a=7.74 Å, b=9.56 Å, and 13.55 Å in slab thickness, with a k-points grid of 2x3x1 and consisted of 96 (32Ti, 64O) atoms. Lastly, in the anatase group, anatase (101) surface had a=10.32 Å, b=7.74 Å, and thickness 13.16 Å. It consisted of the same k-points grid and atoms as anatase (100). We named those structures A-001, A-100 and A-101, respectively. The removal of an oxygen atom on anatase (001) gave oxygen deficiency of 1 in 48 atoms (2.08% oxygen deficiency); on anatase (100) and (101), it was 1 in 64 atoms (1.56% oxygen deficiency).

Finally, the rutile slab model of this study (R-110 32Ti, 64O) has 12.47 Å slab thickness, a 4x2x1 kpoints grid and a=5.92 Å and b=13.18 Å. When removing an oxygen atom, oxygen deficiency was 1 in 64 atoms (1.56% oxygen deficiency).

Results

In this section, we report the results obtained for oxygen vacancy formation energies and the electronic structure calculations (density of states and charges) for both methodologies. For the sake of comparison on methodologies, we present the relative stability of the studied structure using MatSci (derived from PW91 DFT) vs Tiorg implementations on DFTB in Table S1. That comparison shows that, because differences in energies are small (<0.2 eV) for some structures (electronic and geometric), the relative stability may easily change depending on the method used for calculations. Rutile shows in all studied cases a higher EVO formation energy than anatase, which is in agreement with the larger reducibility of anatase respect to rutile for DFT simulations. Nevertheless, for DFTB, rutile shows a competitive EVO value for reduction respect to anatase (only 0.15 eV lower than Osub for A-100 and A-101).

Energetics

In On the other hand, as a general rule, the most reducible site corresponds to the most exposed 2-fold coordinated oxygen atom is on the surface. For DFTB, the O2c energy is 5.62 for R-110 and 6.30 eV for A-101, and 5.72 for A-001 and 5.81 for A-100. However, for A-100 and A-101 the most reducible site for DFTB is the Osub, with only 5.77 and 5.78 eV; for VASP in A-(100), energy for Osub is only 0.04 eV higher than the O2c, the reason is that a strong relaxation moves the O2c in the surface to replace the vacant oxygen in the subsurface, leading to a geometry similar to O2c.

Comparison with literature results is not straightforward due to the different computational setups and the models used. In ref [START_REF] Ha | Oxygen Vacancies of Anatase(101): Extreme Sensitivity to the Density Functional Theory Method[END_REF] , which uses Quantum Espresso PBE+U among others, it is shown that the values and relative ordering of O-vacant A-101 models is sensitive to the computational setting, we obtain similar values and trends: 3.61 eV vs our 3.48 eV for O2C, 4.71 eV vs our 4.88 eV for O3C. In reference [START_REF] Bouzoubaa | Comparison of the Reduction of Metal Oxide Surfaces: TiO2-Anatase, TiO2-Rutile and SnO2-Rutile[END_REF] , with the PW91 functional and narrower slabs, half surface coverage and relaxation limited to neighbouring atoms, EVO values are found larger than our VASP values, especially for anatase slabs: A-001O2C is 4.60 eV vs our 2.80 eV; A-100 O2C is 4.63 eV vs our 3.32 eV, R-110 O2C is 4.21 eV vs our 3.63 eV. It can be concluded that a numerical comparison has to be done carefully despite the fact that the trends seem to be captured.

Relaxation

Removal of a lattice oxygen atom induces moderate relaxation of the atoms neighbouring the vacancy. The distance between 2 Ti across the vacancy increases, as shown in Figure 2, for both DFT and DFTB methods. For instance, the A-001 Ti-Ti distance goes from 3.870 Å in the bare slab to 4.615 Å and 4.325 Å (DFT and DFTB, respectively) in the presence of an O2c vacancy. The Ti-O bonds exhibit an increase of ~5% at the DFT level, where they express a decrease of the same order in the case of DFTB, concerning the unrelaxed distance. In the case of DFT, this is a consequence of the electronic redistribution, with two excess electrons that form Ti 3+ sites bigger in size than Ti 4+ . In the case of A-100 Osub and A-101 Osub2, the conjugate gradient algorithm used in VASP finds a minimum that resembles the O2c system: the O2c atom moves and fills the vacancy. In the case of DFTB, the situation is similar; however, Ti 3+ atoms tend to insert more into the lattice, leading to shorter Ti-O distances; these differences in trends will be explained by differences in charge distributions in the next section.

Charge Analysis

Table 2. Summary of charges on the two adjacent titanium atoms (three atoms for Osub structures) to oxygen vacancies on bare slab and on the most reducible structures for each method (in parenthesis, the most reducible structure). DFTB and VASP employ Mulliken-point and Bader [START_REF] Henkelman | A Fast and Robust Algorithm for Bader Decomposition of Charge Density[END_REF] Table 2 summarizes the charges on titanium atom for the most reducible structures respect to the bare slab found for each method. Although the formal charge on Ti atoms for TiO2 is +4, this material has a large covalent contribution, and charges on Ti are much lower than expected. In addition, the way charges are evaluated differs between the codes employed. For the VASP results, the Bader scheme was employed [START_REF] Henkelman | A Fast and Robust Algorithm for Bader Decomposition of Charge Density[END_REF] , whereas DFTB makes use the Mulliken-point charge scheme.

Thus, DFTB calculations have previously depicted values much lower, in the range of +0.76 to +0.86 |e|, and for oxygen atoms the charge range is only from -0.35 to -0.45|e| [START_REF] Fox | Bulk and Surface Properties of Rutile TiO2 from Self-Consistent-Charge Density Functional Tight Binding[END_REF] . Our results for bare slabs in DFTB show a similar trend, varying charges from +0.83 to +1.00 |e|, increasing the charge value from the outer-most Ti atoms to the most internal ones of slabs. When creating the oxygen vacancy, the charge on O atoms is reduced. The less affected charge is the one on Ti atom on rutile structure, which being symmetric both Ti atoms just reduce from +0.83 to +0.70 |e|. The reduction is more important on the A-001 structure, with the same two-fold coordinated oxygen vacancy as for rutile. In this case there is a clear reduction from +0.86 to +0.50 |e| on Ti atoms. This trend is even larger in the tricoordinated sublayer structures in A-100 and A-101 , where the most reduced Ti atoms decrease the charge to +0.26 and +0.33|e| respectively on the surface sites, being the other Ti atoms of the tricoordinated site less reduced (+0.70 to +0.81 |e|).

For the VASP results, the absolute values are much higher than for DFTB. Typical values for Ti sites are +2.66-+2.70 |e| for stoichiometric systems, and decrease to +2.20/+2.42 |e| for reduced systems.

Whereas the anatase reduced systems exhibit lower values of +2.20/+2.30 for Ti charges, the rutile slab shows slightly larger values, +2.42/+2.41. The decrease observed in the charge on Ti sites correponds to the reduction of the sites by the excess electrons. This is confirmed by the integration of the spin density which gives values close to 1 on each of the Ti sites with lower charge (those presented in Table 2). O values range between -1.30/-1.37 |e| for the surface stoichiometric sites, and increase to -1.40/-1.42 |e| for the defective slab models, showing a redictribution of the electron density due to the excess electrons.

A remark concerns the localization of the excess electrons in neighboring Ti sites. Although in our study the electronic optimization leads to the reduction of two Ti sites close to the oxygen vacancy, such solution is not always the least energy solution and often quasi-degenerate states appear [START_REF] Franchini | Polarons in Materials[END_REF][START_REF] Di Valentin | Reduced and N-Type Doped TiO2: Nature of Ti3+ Species[END_REF] . A systematic study where all the possible distributions are considered would be needed, and it is out of the scope of the present work that relies on a semi-quantitative comparison between DFTB and VASP.

Finally, it is important to note that the spin states consider either singlet (Nα-Nβ=0) or triplet (Nα-Nβ=2) solutions. Whereas the two methods agree for the O2C states being a triplet (Table 1), discrepancies are observed for the other vacancy types. Both DFTB settings, based on B3LYP (Tiorg) and PW91 (Matsci), seem to exhibit a preference for singlet states, although the difference in energy with triplet states may be small, such as in R-110 O2C where DFTB Tiorg shows 0.08 eV less stable the singlet, whereas Matsci stabilizes the singlet by 0.10 eV with respect to the triplet (Table S1). On the contrary, VASP seems to favor triplet states, with almost isoenergetic values for R-110 Osub1, A-001 Osub1 and O3C, A-100 Osub1, A-101 Osub2. In conclusion, the cases where the singlet-triplet spin states are close in energy, seem to be captured by the two methods, although a finer electronic state analysis is required to provide more accurate description.

Density of States

The electronic structure of titania systems changes drastically in the presence of oxygen vacancies due to the excess electrons transferred to titanium sites. As a result, the systems may exhibit open-shell spin configurations that should be reflected in the corresponding density of states (DOS), see for instance [START_REF] Helali | Scaling Reducibility of Metal Oxides[END_REF][START_REF] Di Valentin | Electronic Structure of Defect States in Hydroxylated and Reduced Rutile TiO2(110) Surfaces[END_REF][START_REF] Ha | Oxygen Vacancies of Anatase(101): Extreme Sensitivity to the Density Functional Theory Method[END_REF] . The main feature is the presence of states in the bandgap, which using DFT+U (or hybrid functionals) corresponds to localized electrons on two distinct Ti sites. For the sake of comparison of methods and how they exhibit electronic features, we include in Figure S2 of the Supporting Information the DOS of bulk 2x2x2 TiO2 at both DFT and DFTB levels, with GGA(+U) and hybrids implementations. It can be observed that the excess electrons occupy a Ti 3d level in the bottom of the conduction band (GGA Matsci DFTB), whereas it lies as an isolated state in the band gap for the methods using hybrid functionals (Tiorg DFTB and B3LYP VASP) as well as DFT+U VASP.

The VASP calculated DOS is in very good agreement with the literature as regards the open-shell nature of the reduced systems, and the presence of localized states in the gap. Figure 3 displays the total DOS for stoichiometric and defective A-001 slabs. It can be observed that the stoichiometric slab displays a bandgap of 1.6 eV in the region -0.3 to 1.3 eV. In the same region, all the oxygen-deficient slabs present localized states: O2c as a peak centred at 0.20 eV, O3c as a peak centered at -0.22, Osub as a peak centered at 0 eV and Obulk as a more delocalized state between -0.4 and 0.3 eV. The other slabs present similar features, in all cases the peaks in the gap corresponding to Ti 3d states carrying unpaired electrons in agreement with the Bader spin density analysis (not shown). Figure 4 shows the total DOS for DFTB for the A-001 structure (analogue to Figure 3). This time, the bare slab presents a wider bandgap of 2.8 eV, from -2.8 to -0.0 eV. Again, all the reduced systems show localized states in the gap, this time much closer to the conduction band than on Figure 3: O2c as a peak centred at -0.6 eV, O3c as a peak centered at -0.4 eV, Osub as a peak centered at -0.7 eV and Obulk a state at -0.5 eV. This trend is similar to DFT. A deeper analysis of the DOS and their projection on atomic species orbitals (PDOS) for DFTB can be found in the Figures S2 to S5 of Section 3 in Supporting Information; in general, the effect of removing oxygen atoms from the slabs is the reduction of DOS on the top site of the valence band and the appearance of the peak in the bandgap.

Discussion

Figure 5 displays the VASP and DFTB EVO values for selected O-deficient systems. The EVO values show that the surface termination is selective: A-001 presents the lower values followed by A-100, A-101 and R-110. This is consistent with the literature where anatase is found to be more reducible than rutile. However, to have a quantitative comparison, many effects need to be considered: computational setting, a structural model used, coverage, symmetry, etc.

The coordination of the O atom removed seems directly related to the energy cost: the lower coordinated O2c needs lower EVO values. However, subsurface positions are in some cases found more stable for DFTB and even competitive with O2c as in the case of A-100 where O2c is 3.32 eV and Osub 3.36 eV (VASP), 5.81 and 5.77 eV (DFTB). In all cases, O3c and Obulk values display the higher EVO values.

The comparison between VASP and DFTB values shows a systematic increase of EVO for DFTB compared with VASP. This could be due to the different setup and/or value for O2 molecule that is known to be badly described by DFT. Also, the VASP values are more dispersed than the DFTB ones (~1.5 eV range vs < 1 eV respectively). Interestingly, the trend of EVO obtained with VASP and DFTB is similar, as can be seen in Figure 5. The O2c and Osub positions show the lowest EVO values in all the cases. In the case of DFTB, the trend O2c-Osub might be inversed. Charge analysis shows that the formation of oxygen vacancies decreases the absolute charge on the atoms neighbouring the defect, irrespective of the method used. However, it can be observed that the VASP values are higher than the DFTB ones due to different schemes of partitioning the charge, as previously shown at DFT, hybrid and post-HF schemes [START_REF] Koch | On the Charge State of Titanium in Titanium Dioxide[END_REF] . However, in both cases, one clearly sees the reduction of adjacent Ti atoms to the oxygen vacancy.

As for the DOS, the two methods display the appearance of peaks in the bandgap upon reduction. However, the position of the peak is closer to the valence band for VASP, and closer to the conduction band for DFTB. This is associated with the functional used, PBE+U for the former vs the DFTB implementation using B3LYP for the latter. The states are localized in the Ti sites that become reduced. This effect is observed in the charge decrease, and the integrated Bader spin density for all the slabs studied.

DFTB and DFT+U calculations presented in this paper are based on B3LYP and GGA+U, respectively. Note bandgap states can be scrutinized through both of these computational methodologies. Besides, Evo trends and bandgap graphics obtained from VASP and DFTB seem to have similar values and be meaningful measures.

Charge delocalization analysis is valuable for its pivotal role to describe the relation among Ovac, reducibility and ROS production. The study supports DFT+U by showing that the U value can reproduce the experimental value to describe TiO2 nanoparticles´ electronic structure. [START_REF] Morales-García | On the Use of DFT+U to Describe the Electronic Structure of TiO2 Nanoparticles: (TiO2)35 as a Case Study[END_REF] We conclude that, regardless of the methodology, our study shows an apparent adjacent Ti atoms reduction to Ovac. This allows us to mention that both DFTB and DFT+U would consistently serve the purpose of being used for NP toxicology studies, in a semi-quantitative manner.

Finally, we want to make some considerations on how the oxygen vacancies in TiO2 may increase the toxicity of the material. On the one hand, the formation of Ovac may lead to the release of Reactive Oxidative Species due to the redox potential of Ti 3+ . Besides, the presence of oxygen vacancies introduces in the reduced material new sites that may be a source for reactivity. On the other hand, the inclusion of oxygen vacancies creates the states in the bandgap, close to the conduction band for DFTB; and, at the same time, the conduction band decreases its energy and gets much closer to the Fermi level. This fact, together with the lower Z/r ratio of Ti 3+ respect to Ti 4+ , may explain higher toxicity by the rule found by Liu et al. [START_REF] Liu | Development of Structure-Activity Relationship for Metal Oxide Nanoparticles[END_REF] , that shows that metal oxides showing particular conduction band energies and small Z/r ratios are much more toxic than other metal oxides. The relationship between oxygen vacancies formation energies and other reducibility descriptors (band gap, electron affinity…) has also been reported for other potential nanotoxic metal oxides such as ZnO, CeO2 or SiO2 [START_REF] Helali | Scaling Reducibility of Metal Oxides[END_REF][START_REF] Deml | Intrinsic Material Properties Dictating Oxygen Vacancy Formation Energetics in Metal Oxides[END_REF] .

Conclusions

In the present work, the comparison between two widely used computational methods, VASP (PBE+U, plane waves) and DFTB (SK files adjusted on a base of B3LYP, Pople-style double-ζ 6-31G(d) calculations ), has been conducted to analyze the formation of oxygen vacancies in TiO2 slab models. Both schemes provide qualitative agreement as regards the energetics of formation of vacancies and the sites involved: anatase are more reducible than rutile slabs, and O2c sites are the easiest to remove. However, some discrepancy occurs for A-100 and A-101 where the two methods do not predict the same site. The charge analysis provides similar results regarding the decrease of charge upon reduction, with different absolute charges though; this is likely due to the different ways of computing charges. Finally, the electronic structure is analyzed by obtaining DOS that, in the two methods, show similar behaviour: localized peaks in the gap of the reduced material. However, the relative position of the peak in the gap is significantly different when using DFT+U or DFTB schemes used in this work. In conclusion, the two computational schemes converge to the same qualitative results and trends, but some caution is needed to conduct a quantitative comparison, in particular as regards the accuracy for using of oxygen vacancy as a descriptor for toxicologic studies.

Figure 1 :

 1 Figure 1: Side view of the cells for (001), (100), (101) surfaces of the anatase (A) and (110) surface of rutile (R). Ti atoms in blue, O atoms in red. Oxygen sites that are removed for each vacancy are here labeled.

Figure 2 :

 2 Figure 2: Selected distances (in Å) before and after relaxation for the slab models of O2c vacancies used in the present work. The cross indicates the atom removed and the dashed circle the position of the vacancy. For the sake of simplicity, only the VASP values are shown, the DFTB are described in the text.

Figure 3 :

 3 Figure 3: Total density of states obtained for the A-001 slab model, stoichiometric (slab) and defective for DFT. The reduced systems show a peak in the middle of the gap, indicating localized excess electrons in the Ti sites (indicated with an arrow). The top of the valence band is aligned for all the systems for the sake of comparison.

Figure 4 :

 4 Figure 4: Total density of states obtained for the A-001 slab model, stoichiometric (slab) and defective for DFTB. The reduced systems show a peak in the middle of the gap, indicating localized excess electrons in the Ti sites (indicated with an arrow). The top of the valence band is aligned for all the systems for the sake of comparison.

Figure 5 :

 5 Figure 5: Schematic behaviour of the EVO for selected stable vacancies found for VASP and DFTB.

Table 1 .

 1 Summary the settings used. VASP calculations performed in this work are based on the Generalised Gradient Approximation (GGA+U), whereas DFTB implementation for the Slater-Koster files Tiorg used in this work are based on the hybrid B3LYP functional.

	Surface	Method	Obulk	Osub1	Osub2	O3c	O3c_2	O2c
	A-001	DFTB	6.44 (s)	6.01 (s)	/	6.03 (s)	/	5.72 (t)
		VASP	4.35(t)	3.88(s)	/	3.89(s)	/	2.80(t)
		Literature				5.46(s) a		4.60(s) a
	A-100	DFTB	6.38 (s)	5.77 (s)	6.15(s)	6.38 (s)	/	5.81 (t)
		VASP	4.51(t)	to O2c	4.49(t)	4.65(t)	/	3.32(t)
		Literature				5.56(s) a		4.63(t) a
	A-101	DFTB	6.47 (s)	6.26 (s)	5.78 (s)	6.50 (s)	6.31 (s)	6.30 (t)
		VASP	4.47(t)	3.94(s)	to O2c	4.88(s)	4.47(t)	3.48(t)
		Literature		4.50(t) b	4.43(st) b	4.71(t) b	4.62(t) b	3.61(t) b
	R-110	DFTB	7.02 (s)	6.49 (s)	5.97 (s)	6.59 (t)	/	5.62 (t)
		VASP	4.45(t)	4.91(st)	4.24(t)	4.41(t)	/	3.63(t)
		Literature				5.27 (t) a		4.21(t) a
	a Data from 43

of the EVO values obtained computed following equation 1 in eV. Labelling of Ovac is displayed in Figure

1

. In parenthesis, the most stable electronic state for the structure (s: singlet; t:triplet).

b Data from

[START_REF] Ha | Oxygen Vacancies of Anatase(101): Extreme Sensitivity to the Density Functional Theory Method[END_REF]

, spin polarized are labelled (t) in the table

Table 1

shows the EVO cost formation for all studied structures, for both at DFTB and DFT levels, for the most stable electronic structure in each case (singlet or triplet). Absolute values for DFTB are always larger (between 5.62 and 7.02 eV) than DFT (2.80 to 3.63 eV). This difference in absolute values comes from

  general, the Obulk vacancies show the largest energy for all structures and methods. For DFTB, this moves in a narrow range of 6.38 to 6.47 eV in anatase and 7.02 for rutile; for VASP, these values are 4.35 to 4.51 for anatase and 4.45 eV for rutile. Nevertheless, for A-100 and A-101, DFTB displays almost an equal value (6.38 and 6.47 eV, respectively for each surface) for the bulk oxygen and the tricoordinated (O3c) oxygen on the surface (6.38 and 6.50 eV), and VASP displays also differences for the Obulk (4.51 and 4.47 eV) and O3c (4.65 and 4.88 eV). Rutile (110) follows the general trend for DFTB, but VASP presents similar results to A-100 and A-101, being the O3c vacancy (4.41 eV) close to the bulk value (4.45 eV).

Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 814426 NanoInformaTIX (Development and Implementation of a Sustainable Modelling Platform for NanoInformatics). Dr. B. Diawara is acknowledged for the Modelview program. This work was performed using HPC resources from GENCI-CINES/IDRIS (Grant 2020-x2020082131 and 2021-x2021082131).

ASSOCIATED CONTENT


Clément Mayet
email: clement.mayet@lecnam.net

Denis Labrousse
email: denis.labrousse@satie.ens-cachan.fr

Rihab Bkekri
email: rbkekri@satie.ens-cachan.fr

Francis Roy
email: francis.roy@stellantis.com

Gaël Pongnot
email: gael.pongnot@ens-paris-saclay.fr

Energetic Macroscopic Representation and Inversion-Based Control of a Multi-Level Inverter with Integrated Battery for Electric Vehicles

Keywords: electric vehicles, integrated battery system, multilevel inverter, simulation and control, battery management system, Energetic Macroscopic Representation

This paper deals with the Energetic Macroscopic Representation (EMR) of a Multi-Level Inverter with Integrated Battery (MLI-IB) for electric vehicles. The MLI-IB consists of modular cascaded modules connected in series. Each module is composed by an individual battery and an H-bridge converter. An inversion-based control is deduced from the EMR to control the output voltage. In addition, the MLI-IB is managed to balance the state-of-charge of the batteries between the various modules.

I. INTRODUCTION

Energy saving, greenhouse gases emissions and fuels depletion are critical issues for the coming decades. The researchers and development activities related to transportation have therefore emphasized the development of high efficiency, clean, and reliable transportation systems [START_REF] Eshani | Modern Electric, Hybrid Electric, and Fuel Cell Vehicles: Fundamentals, Theory, and Design[END_REF]. Electric Vehicles (EVs), Hybrid Electric Vehicles (HEVs), and Fuel Cell Electric Vehicles (FCEVs) are expected solutions to increase future mobility while limiting environmental impacts [START_REF] Chan | The State of the Art of Electric, Hybrid, and Fuel Cell Vehicles[END_REF]. Nevertheless, the development of EVs faces different challenges such as high initial cost, short driving range (range anxiety) and long charging (refueling) time compared to conventional Internal Combustion Engine Vehicles (ICEVs), which justifies new research and development on these different topics [START_REF] Eshani | Modern Electric, Hybrid Electric, and Fuel Cell Vehicles: Fundamentals, Theory, and Design[END_REF], [START_REF] Chan | The State of the Art of Electric, Hybrid, and Fuel Cell Vehicles[END_REF].

Focusing on EVs, the traction system is achieved using one or more Electric Machines (EMs), which are supplied by power electronics converters and energy sources. For conventional EVs, a high voltage battery (typically 400 V) is usually considered as main energy source and is a combination of individual battery cells (for example, 3.6 V Li-ion cells) in series and in parallel. However, the reliability of the entire battery depends on the good operating conditions of each cell and requires a Battery Management System (BMS) [START_REF] Gandoman | Status and future perspectives of reliability assessment for electric vehiles[END_REF]. The role of the BMS is to passively balance the State-of-Charge (SoC) between the different cells, to manage the thermal aspect of the different cells, and eventually to shut down the system in case of failure [START_REF] Altaf | Load Management of Modular Battery Using Model Predictive Control: Thermal and State-fo-Charge Balancing[END_REF], [START_REF] Hannan | State-of-the-Art and Energy Management System of Lithium-Ion Batteries in Electric Vehicle Applications: Issues and Recommendations[END_REF]. One of the disadvantages of such a topology is that the entire battery pack is limited by the weakest cell in terms of capacity, which can cause the battery to shut down even if the other cells still have energy remaining. Another aspect is related to the necessity to use the technique of Pulse Width Modulation (PWM) to control the inverter and adapt the battery voltage to the supply voltage requested by the EM. This induces switching losses, which mainly depend on the switching frequency, the battery voltage, and the current [START_REF] Eskandari | Accurate Analytical Switching Loss Model for High Voltage SiC MOSFETs Includes Parasitics and Body Diode Rerverse Recovery Effects[END_REF].

The Multi-Level Inverter [START_REF] Rodriguez | Multilevel Converters: An Enabling Technology for High-Power Applications[END_REF]- [START_REF] Vasu | Asymmetric Cascaded H-Bridge Multilevel Inverter With Single DC Source per Phase[END_REF] with Integrated Battery (MLI-IB) is a potential candidate to solve the aforementioned issues. The MLI-IB integrates battery cells into cascaded H-bridge modules. It aims to fulfill the roles of phase voltage control, BMS, active SoC balancing, and can manage each cell individually [START_REF] Mayet | Simulation and Control of a New Integrated Battery System for Automotive Applications[END_REF]. If a cell fails, it can be disconnected from the system, ensuring a continuous operation (albeit in a degraded mode), and improving the reliability of the entire system. The MLI-IB allows removing the passive balancing system of conventional BMS by adding this function directly in the control. In addition, a specific control without PWM technique, which, added to the fact that the voltage of the individual battery used in the MLI-IB is lower than the high voltage of the battery pack used in conventional EVs, leads to drastically reduce the switching losses of power electronics devices. Another advantage of this topology is its modularity. Indeed, this offers opportunities for significant reconfiguration, but also for the recovery and reuse of used battery cells, which is a sustainable way of using them given the current global challenges of energy storage. However, MLI-IB is a more complex system than conventional battery system and must be controlled and assessed for future uses in automotive applications. A first study was carried out to propose a model and control of the MLI-IB, integrating the SoC balancing management [START_REF] Mayet | Simulation and Control of a New Integrated Battery System for Automotive Applications[END_REF].

The objective of this paper is to propose a new organization and control of the MLI-IB based on the Energetic Macroscopic Representation (EMR) methodology. The EMR is a graphical description tool that highlights the energy properties of complex systems while respecting physical causality [START_REF] Bouscayrol | Chapter 3: Graphic Formalism for the Control of Multi-Physical Energetic Systems: COG and EMR[END_REF], [START_REF] Iwasaki | Causality and model abstraction[END_REF]. Inversion rules have been defined to deduce the Inversion-Based Control (IBC) of the studied system. This methodology makes it possible to better highlight the different functions and levels of control of this new topology. This paper focuses on the MLI-IB. The traction system of the EV is outside the scope of this paper, but its model, EMR, and control can be found in [START_REF] Chan | Electric, Hybrid, and Fuel-Cell Vehicles: Architectures and Modeling[END_REF] as an example. In Section II, the architecture of the MLI-IB is presented. Then, its model, EMR and IBC are developed in Section III. The IBC of the MLI-IB requires defining energy management strategies, which will be explained in Section IV. Finally, the simulation results will be analyzed and discussed in Section V.

II. ARCHITECTURE OF THE MLI-IB

The MLI-IB is composed of 3 phases to supply the EM of the EV (Fig. 1). The objective of the MLI-IB is to impose the phase voltage requested by the control of the EM. Each phase j consists of q modules connected in series. The different modules are controlled to impose their output voltages according to a specific control and strategies which will be presented in the following sections. Each module is composed of an independent battery and two Unit Conversion Elements (UCEs). The battery can be composed by one or several cells in series and/or parallel. The two UCEs form an H-bridge converter which enables bidirectional power flow from each battery. Finally, each UCE is composed by two MOSFETs. The two UCEs are controlled to impose the reference output voltage of the module. This architecture is modular and can be easily reconfigured. It also offers high redundancy, which increases the reliability of the system and gives additional degrees of freedom that can be used in the control specially to add the function of management of the SoC of the battery cells. In the next section, the model, EMR and IBC of one phase of the MLI-IB are presented.

III. EMR AND IBC OF THE MLI-IB

This section describes the model of the MLI-IB, which is organized using the EMR. A control structure is then deduced.

A. EMR of the MLI-IB

The EMR is a graphical description tool, which highlight the energy properties of a system [START_REF] Bouscayrol | Chapter 3: Graphic Formalism for the Control of Multi-Physical Energetic Systems: COG and EMR[END_REF]. It organizes the system into interconnected basic elements: source of energy (green oval), accumulation of energy (orange crossed rectangle), mono (orange square) or multi (orange circle) domain conversion of energy, and distribution/coupling of energy (double orange square). All elements are connected according to the interaction principle: the product of the action and reaction variables leads to the power exchanged. Furthermore, all the components are described with respect to the physical causality [START_REF] Iwasaki | Causality and model abstraction[END_REF].

The EMR of the MLI-IB (Fig. 2) is composed of different elements that convert energy between the batteries (BAT) of the different modules and the EM. This EMR is given for one phase j and can be duplicated to represent the 3-phase system. The EM gear gear

ip-1 ip-2 ip-3 vp-3 vp-3-ref 3-phase MLI-IB Trac. sub. MLI-IB Phase 3 vp-2-ref MLI-IB Phase 2 vp-2 vp-1-ref MLI-IB Phase 1 vp-1 vp-j vm-j-q vm-j-q-ref vm-j-3 vm-j-3-ref vm-j-2 vm-j-2-ref vm-j-1 vm-j-1-ref ip-j Module q Module 3 Module 2 Module 1 Phase to module control vb-j-k vo-j-k-1 ii-j-k-1 vm-j-k ip-j ib-j-k io-j-k-1 UCE 1 UCE 2 ii-j-k-2 vo-j-k-2 io-j-k-2 vo-j-k-1-ref vo-j-k-2-ref One phase of the MLI-IB One module of the MLI-IB K1-j-k-p K2-j-k-p ii-j-k-p uc-j-k-p vo-j-k-p io-j-k-p vp-j-ref Module to UCE control vm-j-k-ref UCE control vo-j-k-p-ref vb-j-k One UCE Fig. 1. Traction subsystem of an EV supplied by a 3-phase MLI-IB vb-j-k BAT ib-j-k uc-j-k-2 vb-j-k ii-j-k-2 ii-j-k-1 vb-j-k uc-j-k-1 vo-j-k-2 io-j-k-2 vo-j-k-1 io-j-k-1 vm-j-k EM ip-j vp-j ip-j vp-j-ref vm-j-k-ref vo-j-k-2-ref vo-j-k-1-ref uc-j-k-1 Strategy Module Strategy Phase vo-j-k-1-ref vm-j-k'-ref UCEs Modules Phase SoCb-j-k (1) (2) (3) (3) (4)
(5, 6)

(7) (8) (9) (9)
(Fig. 3) (Fig. 3 source EM represents one phase of the EM and imposes the phase current ip-j to the MLI-IB. A coupling element is used to represent the serial connection of the different modules k. A vectorial representation is chosen. The variables vm-j-k and ip-j are vectors containing respectively the output voltages and currents of the different modules. The phase voltage vp-j is the sum of all the output voltages vm-j-k and the currents circulating in the different modules are the common phase current ip-j [START_REF] Eshani | Modern Electric, Hybrid Electric, and Fuel Cell Vehicles: Fundamentals, Theory, and Design[END_REF].

⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎧ = = = ⎣ ⎢ ⎢ ⎢ ⎡ ⋮ ⎦ ⎥ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎢ ⎡ ⋮ ⎦ ⎥ ⎥ ⎥ ⎤ = ⎣ ⎢ ⎢ ⎢ ⎡ ⋮ ⎦ ⎥ ⎥ ⎥ ⎤ (1)
In the following expressions, all variables are 1-by-q vectors, where each component corresponds to a module. Another element is used to represent the coupling between the two UCEs that exist in each module. The output voltages vm-j-k of the modules are the differences between vo-j-k-1 and vo-j-k-2, which are respectively the output voltages of the UCEs 1 and the UCEs 2 (2). The output currents of each UCEs are also given in [START_REF] Chan | The State of the Art of Electric, Hybrid, and Fuel Cell Vehicles[END_REF].

= - = = - (2)
All the UCEs 1 are represented by a single conversion element. An identical conversion element is used to represent all the UCEs 2. They are described by [START_REF] Gandoman | Status and future perspectives of reliability assessment for electric vehiles[END_REF], where p is 1 for the UCEs 1 and 2 for the UCEs 2. The output voltages vo-j-k-p depend on the voltages vb-j-k of the batteries and the control variables uc-j-k-p of the UCEs. Equivalent resistances Req-j-k-p are considered to take into account the losses of the UCEs. The input currents iij-k-p are also expressed in (3).

= × ! -" # × $ = × (3)
A last coupling element is used to represent the parallel connection on the battery of the two UCEs of all modules. The voltages vb-j-k of the batteries are thus common inputs for the UCEs 1 and the UCEs 2 (4). The currents ib-j-k of the batteries are then the sum of the input currents of the UCEs 1 and of the UCEs 2, as expressed in [START_REF] Altaf | Load Management of Modular Battery Using Model Predictive Control: Thermal and State-fo-Charge Balancing[END_REF].

% ! common ! = $ + $ (4)
Finally, the batteries of each module are represented by a single source, which imposes the voltages vb-j-k of the different batteries. A simple RE type battery model is used and is described by [START_REF] Hannan | State-of-the-Art and Energy Management System of Lithium-Ion Batteries in Electric Vehicle Applications: Issues and Recommendations[END_REF], where eb-j-k are the open circuit voltages and rb-j-k are the internal resistances of the batteries. Both parameters depend on the SoCs SoCb-j-k of the batteries. The SoCs are determined using a Coulomb counting method, which is expressed by [START_REF] Eskandari | Accurate Analytical Switching Loss Model for High Voltage SiC MOSFETs Includes Parasitics and Body Diode Rerverse Recovery Effects[END_REF] and where SoCb0-j-k is the initial SoCs and Cbj-k is the effective capacity of the batteries (in Ah).

! = + ! ,-./ ! 0 -1 ! ,-./ ! 0 × ! (5) -./ ! = -./ !2 - 1 3600. / ! 8 ! . 9: (6)

B. IBC of the MLI-IB

A local control, which is called IBC, is deduced from the EMR using inversion rules (see the blue part of Fig. 2) [START_REF] Bouscayrol | Chapter 3: Graphic Formalism for the Control of Multi-Physical Energetic Systems: COG and EMR[END_REF]. The objective of this control is to impose the output reference phase voltage vp-j-ref.

A tuning path is defined to link the tuning variables uc-j-k-1 and uc-j-k-2 to the objective of the control (vp-j). The control structure is obtained by inverting this tuning path step-by-step. Conversion elements are directly inverted and coupling elements require distribution of energy.

The inversion of the coupling element (1) requires sharing the energy between the different modules. Considering it is a coupling of dimension q, it is possible to directly impose the reference output voltages of q -1 modules in the vector vm-j-k'-ref.

The reference output voltage of the q th module can be deduced by [START_REF] Rodriguez | Multilevel Converters: An Enabling Technology for High-Power Applications[END_REF]. The vector containing the reference output voltages vm-j-kref of all the modules is thus obtained by merging the vector vm-jk'-ref and the voltage vm-j-q-ref.

The vector vm-j-k'-ref will be defined in Section IV.B based on a specific voltage generation strategy coupled with a SoC balancing management strategy.

;#< = ;#< - = ;#< (7)
The inversion of the second coupling element (2) requires sharing the energy between the UCE 1 and UCE 2 of each module. Considering it is a coupling of dimension 2, it is chosen to directly impose the reference output voltages of the UCEs 1 in the vector vo-j-k-1-ref. The reference output voltages of the UCEs 2 are thus defined by [START_REF] Malinowski | A survey on Cascaded Multilevel Inverters[END_REF]. The vector vo-j-k-1-ref will be defined in Section IV.A based on the operation strategy of the H-bridge converter, at the module level.

> ;#< = ;#< ;#< = ;#< - ;#< (8)
Finally, the conversion elements (3) representing the UCEs are directly inverted. Neglecting the losses, the duty cycles αj-k-p-ref of each UCE are defined by [START_REF] Babaei | A Single-Phase Cascaded Multilevel Inverter Based on a New Basic Unit With Reduced Number of Power Switches[END_REF]. They can take any values between 0 and 1 and can be connected to the control variables ucj-k-p in several ways according to the chosen model and control of the UCEs (Fig. 3).

kα-swi > αj-k-p-ref uc-j-k-p > αj-k-p-ref uc-j-k-p αj-k-p-ref uc-j-k-p

In case of an average model, the control variables are identical to the duty cycles. However, on the real system, the control variables must be binary signals that can be only 0 or 1. In this case, two different controls are proposed (see Fig. 3). The first is based on PWM technique whereas the second is a direct control where the duty cycles are compared to a constant ka-swi. ka-swi is a constant between 0 and 1 and corresponds to a threshold that determines the switching condition of the UCEs.

IV. ENERGY MANAGEMENT STRATEGIES

Based on the defined IBC, two strategies must be defined. First, the operation strategy of the H-bridge converter of each module is defined to distribute the energy between the two UCEs. Then, a specific voltage generation strategy of phase j, coupled with a SoC balancing management, is defined to share the energy between the different modules.

A. Operation strategy of an H-bridge converter at the module level

According to the operation strategy of the two UCEs that compose the H-bridge converter of a module (Table I), the reference output voltage vo-j-k-1-ref of the UCE 1 is used to impose a positive output voltage vm-j-k of a module while the reference output voltage vo-j-k-2-ref of the UCE 2 is kept to zero. On the contrary, the reference output voltage vo-j-k-2-ref of the UCE 2 is used to impose a negative output voltage vm-j-k of a module while the reference output voltage vo-j-k-1-ref of the UCE 1 is kept to zero. In addition, there are two possibilities for imposing the output voltage vm-j-k of a module to zero. In this case, this paper considers that both UCEs impose reference output voltages equal to zero. This operation strategy is described by directly imposing the reference output voltages of the UCEs 1 as defined in [START_REF] Vasu | Asymmetric Cascaded H-Bridge Multilevel Inverter With Single DC Source per Phase[END_REF]. As explained previously, in such a case, the reference output voltages of the UCEs 2 are deduced in the IBC by [START_REF] Malinowski | A survey on Cascaded Multilevel Inverters[END_REF].

%

;#< = ;#< when ;#< > 0 ;#< = 0 when ;#< ≤ 0 (10)

B. Voltage generation strategy at the phase level

As summarized in equation ( 1), the voltage vp-j of phase j is set by turning on or off some of the modules connected in series. The strategy must activate/deactivate a certain number of modules according to the reference phase voltage vp-j-ref and the actual voltages available within each module (i.e. voltages vb-j-k) that can be measured or estimated. There are two possible modes for the voltage generation strategy, namely unipolar and bipolar [START_REF] Altaf | Load Management of Modular Battery Using Model Predictive Control: Thermal and State-fo-Charge Balancing[END_REF], [START_REF] Rodriguez | Multilevel Converters: An Enabling Technology for High-Power Applications[END_REF]- [START_REF] Vasu | Asymmetric Cascaded H-Bridge Multilevel Inverter With Single DC Source per Phase[END_REF]. The unipolar mode does not allow polarity inversion of any module (i.e. batteries), which means that at any time instant either all batteries of the modules charging or all are discharging. The bipolar mode allows polarity inversion of some modules, which means that at any time instant it is possible to charge some cells while discharging the others. Only the unipolar mode is considered in this paper. Naturally, this architecture and this strategy involve unbalancing the SoCs of the batteries between the modules. A specific ranking strategy has thus been defined to adapt the order of activation of the different modules accordingly to their actual SoCs [START_REF] Mayet | Simulation and Control of a New Integrated Battery System for Automotive Applications[END_REF]. This ranking strategy takes on the role of SoC balancing management and is expressed through an adaptation matrix T(SoCb-j-k) [START_REF] Mayet | Simulation and Control of a New Integrated Battery System for Automotive Applications[END_REF]. This matrix is updated according to a specific period, which must be defined beforehand, or as soon as the maximum delta allowed between the SoCs is reached. The voltage generation strategy uses vectors ranked according to the SoCs of the batteries (variables named xxr-j-n) instead of a natural rank which depends on the topology (variables named xx-j-k) (see Fig. 1).

E F G; H = IJ,-./ ! 0K × F G F G = IJ,-./ ! 0K L × F G; H (11)
The first step of this strategy is thus to rank the voltages of the batteries according to their SoCs [START_REF] Bouscayrol | Chapter 3: Graphic Formalism for the Control of Multi-Physical Energetic Systems: COG and EMR[END_REF].

!; H = IJ,-./ ! 0K × ! (12)
Then, a simple comparison between the reference phase voltage vp-j-ref and a vector, which considers the cumulative sum of the voltages vbr-j-n of the batteries via the function cumsum [START_REF] Iwasaki | Causality and model abstraction[END_REF], leads to obtain the equivalent control variables at the level of the modules [START_REF] Mayet | Simulation and Control of a New Integrated Battery System for Automotive Applications[END_REF]. As the output voltages of the modules can be positive or negative, the sign of the reference output voltage must also be taken into account [START_REF] Chan | Electric, Hybrid, and Fuel-Cell Vehicles: Architectures and Modeling[END_REF].

; $HL H ;#< = ;#< > cumsum, !; H 0 (13)
; H ;#< = ; $HL H ;#< × sign, ;#< 0 (14) Then, the reference output voltages of the different modules are obtained by an Hadamard product between vectors ucmr-j-n-ref and vbr-j-n (15). However, in this way, the last module n that could be potentially connected is never considered because the modules are only activated when the summation of their voltages is just below the reference phase voltage as indicated in [START_REF] Iwasaki | Causality and model abstraction[END_REF]. This leads to always requesting the full battery voltage of a module and never an intermediary voltage. The intermediate reference output voltage of the latter module n can be defined as expressed in (16). At this point, the strategy has defined all the reference output voltages for all modules ranked according to their SoCs. Equation (17) reorganizes the vector according to the natural order (topology) excluding the q th module as explained in (7). Simulations are performed to prove the effectiveness of the MLI-IB in controlling the requested phase voltage and balancing the SoC of each battery. The simulated topology consists of 24 modules in series per phase (q = 24). Each module is made up of 4 battery cells (3.6 V Li-ion type) in series and an H-bridge converter. The EM requires a sinusoidal voltage supply with a frequency between 0 and 1 kHz, a voltage amplitude between 0 and 320 V, and an active power up to 58 kW for the 3-phase system. As an example, a complete model of an electric traction system can be found in [START_REF] Chan | Electric, Hybrid, and Fuel-Cell Vehicles: Architectures and Modeling[END_REF]. However, in this paper the MLI-IB is tested under steady-state conditions and not over a full driving cycle. Therefore, a sinusoidal phase voltage of 200 Hz with an amplitude of 80 V is simulated. The EM imposes a current of amplitude 150 A with a power factor or 0.90 (phase angle of -0.45 rad), which corresponds to absorbed active and reactive powers of 5.4 kW and 2.6 kVAR, respectively, on each phase.

A. Results 1: Ability to generate the requested phase voltage

This subsection presents the ability of the topology and its control to generate an output phase voltage as close as possible to the requested reference voltage. The ideal output voltage is obtained using the average model and control (see Fig. 3). In such a case, the output voltage is strictly equal to its reference as shown on Fig. 4. In this case, seven modules are needed and are sequentially turned on and off to generate the reference output phase voltage. With the average model and control, their output voltages vm-j-k are identical to their references vm-j-k-ref. However, the average output voltages of the modules shown on Fig. 4 cannot be directly imposed by the real modules and the control variables must be converted into binary signals. Two different techniques are considered in this paper (see Fig. 3).

The PWM technique allows generating a voltage with high frequency harmonics, which can be filtered if needed. This leads to a voltage waveform for the EM very similar to its reference (Fig. 5). However, this technique also induces high switching losses in the MOSFETs due to the high frequency switching. Another technique using direct control is tested (Fig. 6). In this case the PWM technique is not used, which highly reduces the number of switchings, and thus the losses. The module is activated when its reference voltage reaches a threshold, defined by the coefficient ka-swi (see Fig. 3). In this case, and depending on its activation, a module can only impose zero or ± the battery voltage as output voltage, which can cause its actual average output voltage to be more or less far from its reference. It means that, depending on the amplitude of the reference phase voltage and the voltage level of the batteries, the quality of the harmonics can highly vary. High phase voltage amplitude with low battery voltages will lead to low harmonics, while low phase voltage amplitude with high battery voltages will lead to high harmonics. Hybrid control, switching between these two techniques, could help reduce switching losses while keeping low harmonics. This aspect needs to be studied further.

B. Results 2: Ability to balance the SoCs of the different cells

This subsection presents the ability of the topology and its control to balance the SoCs of the different modules. As explained in Section IV.B, this topology naturally unbalances the SoCs of the different modules. A specific ranking strategy is thus implemented to adapt the order of use of the different modules, and therefore takes the role of battery balancing management strategy [START_REF] Mayet | Simulation and Control of a New Integrated Battery System for Automotive Applications[END_REF]. The situation described in subsection V.A is simulated over 1 s with an adaptation of the ranks every 40 ms (Fig. 7). However, this period could be adapted if necessary. The results show that all batteries start with a SoC of 80 %. They are sequentially discharged according to the activation/deactivation of the different modules. The proposed ranking strategy is able to choose appropriate ranks to maintain an acceptable balance between the different SoCs. In such conditions, all 24 modules are successively used to generate the requested phase voltage although the maximum number of modules activated at the same time is seven as mentioned in the previous subsection. Such a ranking strategy can be adapted according to different conditions. For example, the ranks can be defined according to the sign (charge/discharge) of the instantaneous or of the active power [START_REF] Mayet | Simulation and Control of a New Integrated Battery System for Automotive Applications[END_REF]. In addition, the condition for adapting the ranks can be a specific period of time (as in Fig. 7) or a maximum acceptable difference between the SoCs of the modules not to be exceeded. The choice of the ranking conditions has an impact on the losses within the batteries but also on the switching losses of the converters. Future studies will be carried out on this topic.

VI. CONCLUSION

Fig. 2 .

 2 Fig. 2. EMR and IBC of one phase of the MLI-IB.

Fig. 3 .

 3 Fig. 3. Inversion of the UCEs.

Fig. 4 .

 4 Fig. 4. Simulation results with the ideal average (reference) control.

Fig. 5 .

 5 Fig. 5. Simulation results with the PWM control.

1 Fig. 6 .

 16 Fig. 6. Simulation results with the direct control.

Fig. 7 .

 7 Fig. 7. Simulation results with the ranking strategy (SoCs balancing).

TABLE I

 I

	.	LOGICAL OPERATION OF AN H-BRIDGE CONVERTER
	Control variables of	Output voltages of	Output voltage of
	both UCEs	both UCEs	the module
	uc-j-k-1	uc-j-k-2	vo-j-k-1	vo-j-k-2	vm-j-k
	0	0	zero	zero	zero
	0	1	zero	positive	negative
	1	0	positive	zero	positive
	1	1	positive	positive	zero

This paper has been achieved within IBIS project, which is funded by ADEME (French Agency for the Ecological Transition) through the Program "Investissement d'Avenir (PIA)".

This paper investigates a Multi-Level Inverter with Integrated-Battery (MLI-IB) for automotive applications. The MLI-IB consists of modular cascaded modules connected in series. Each module is composed of an individual battery and an H-bridge converter. The objective of the MLI-IB is to supply the traction drive of electric vehicles by allowing better management of battery cells compared to a conventional battery system. The MLI-IB is modeled and represented using the Energetic Macroscopic Representation (EMR). Then, an Inversion-Based Control (IBC) is deduced to control and manage the MLI-IB. Preliminary simulation show promising results for future automotive applications. Further work needs to be done to study cell ageing, efficiency, and harmonic generation within such a topology.