Marco Fronzi

M Hussein

N Assadi

Dorian A H Hanaor

Theoretical insights into the hydrophobicity of low index CeO2 surfaces

Keywords: CeO2 surfaces, contact angle, density functional theory, GGA+U

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Understanding and controlling surface wettability in ceramic materials is paramount towards their implementation in a wide range of high-value novel biomedical and chemical applications. Hydrophobic coatings have been routinely utilised to control the reactivity and mechanical performance of surfaces while hydrophilicity can be used to facilitate microfluidics, self-cleaning behaviour and thermal coatings [1,[START_REF] Peng | Heat Mass Transfer, Experimental investigation on steam condensation heat transfer enhancement with vertically patterned hydrophobic-hydrophilic hybrid surfaces[END_REF]. Harnessing the intrinsic wettability of unmodified oxide surfaces rather than polymeric coatings or nanostructuring is a valuable pathway towards thermally and chemically robust surface functionality. Furthermore, the wetting characteristics of an oxide surface are indicative of more general surface features, including defect chemistry, energetics and lattice termination, with implications towards surface localised phenomena in general. The wettability of surfaces is generally parametrised in terms of the contact angle of water at the surface (Θ). Typically, materials are defined as super-hydrophilic if Θ is close to 0°, and super-hydrophobic if Θ is larger than 150°.

Recent experiments have identified robust intrinsic hydrophobicity of surfaces in rare-earth oxides (REOs), giving a prospect to multiple novel applications [START_REF] Azimi | [END_REF][4][5][START_REF] Reitz | Facile and general Synthesis of thermally stable ordered mesoporous rare-earth oxide ceramic thin films with uniform mid-size to large-size pores and strong crystalline texture[END_REF]. Moreover, the wetting characteristics of this class of material are of importance towards interpreting and improving their catalytic activity. In particular, hydrophobicity confers resistance to water deactivation at catalyst surfaces and enhances the adsorption of organic compounds. Consequently, hydrophobicity, or organophilicity, is frequently associated with higher performance and is often a desired trait in applications involving the oxidation of organic compounds and selective synthesis [START_REF] Nowak | A comparison of contact angle measurement techniques applied to highly porous catalyst supports[END_REF][START_REF] Chuang | catalytic oxidation of formaldehyde over hydrophobic catalysts[END_REF][START_REF] Wu | VOC deep oxidation over Pt catalysts using hydrophobic supports[END_REF][START_REF] Dacquin | Interdependent lateral interactions, hydrophobicity and acid strength and their influence on the catalytic activity of nanoporous sulfonic acid silicas[END_REF][START_REF] Francesco | New Simple Hydrophobic Proline Derivatives as Highly Active and stereoselective catalysts for the direct asymmetric Aldol reaction in aqueous medium[END_REF]. The lower levels of hydroxyl ion driven surface defects associated with hydrophobicity are further known to enhance luminescence [START_REF] Zheng | Formation and properties of hydrophobic CeO2 nanoparticles[END_REF]. In recent years, interest in hydrophobic surfaces has been further driven by heat transfer applications, as dropwise water condensation has been found to be a highly effective heat transfer mechanism. As intrinsically hydrophobic ceramic materials, REOs are particularly attractive for such applications, relative to surface functionalised materials, owing to their thermal and chemical robustness.

Because of its diverse functionality and the facile tunability of its bulk-and surface properties, cerium oxide, crystallising in a cubic fluorite structure, is the most widely studied REO [START_REF] Trovarelli | Catalytic properties of ceria and CeO2-containing materials[END_REF][START_REF] Choudhury | Annealing temperature and oxygen-vacancy-dependent variation of lattice strain, band gap and luminescence properties of CeO2 nanoparticles[END_REF]. CeO2 is often studied in nanoparticle, pristine and doped forms and has been extensively studied for applications in catalysis and biomedicine, in which the material's interaction with adsorbates, in general, and wettability in particular, plays a key role in performance. Studies have shown that adsorption behaviour of various molecular and ionic species diverges significantly between the different crystallographic planes of CeO2, highlighting the importance of application-targeted morphological design to enhance the expression of high-performance surfaces in ceria-based catalysts [START_REF] Pengren | durable anti-icing coatings based on CeO2 microparticles[END_REF][START_REF] Huang | Facetselective response of trigger molecule to CeO2 {110} for up-regulating oxidase-like activity[END_REF][START_REF] Tam | Non-wetting nickel-cerium oxide composite coatings with remarkable wear stability[END_REF].

Hydrophobicity and wetting angles of up to 116° have been reported for unmodified surfaces of CeO2 nanomaterials, which tend to be dominated by (111) plane [START_REF] Zheng | Formation and properties of hydrophobic CeO2 nanoparticles[END_REF]. Due to its intrinsic hydrophobicity, CeO2 has been proposed for applications in robust anti-icing and wear resistant coatings [START_REF] Agarwal | ChemSusChem, Exposed surfaces on shape-controlled ceria nanoparticles revealed through AC-TEM and water-gas shift reactivity[END_REF]. Although many experimental investigations have been conducted into the physical properties of this oxide, the relative interactions of its low index surfaces with water remain poorly understood. To gain fundamental mechanistic insights into the behaviour of applied CeO2 surfaces, here we use a theoreticalcomputational methodology to investigate the hydrophobicity of the principal surfaces of CeO2, namely the terminated (111), ( 110) and (100) planes, which are dominant in morphologies of octahedra, nanorods and cubes respectively [START_REF] Huang | Facetselective response of trigger molecule to CeO2 {110} for up-regulating oxidase-like activity[END_REF][START_REF] Yan | Template-free hydrothermal synthesis of CeO2 nano-octahedrons and nanorods: Investigation of the morphology evolution[END_REF]. In the present work we parametrise surface wettability to obtain a measurable quantity that can be compared with experimental results. From the adsorption energies of a water bilayer, we compute the contact angles Θ to compare the surfaces' hydrophobicity [START_REF] Zheng | Formation and properties of hydrophobic CeO2 nanoparticles[END_REF].

Methodology

The calculations presented in this work were performed using unrestricted density functional theory approach, within the generalised gradient approximation (GGA) based on the Perdew-Burke-Ernzerhof (PBE) parameterisation for the exchange-correlation functional [START_REF] Perdew | [END_REF]. We used periodic boundary conditions and a plane-wave basis set as implemented in VASP code [START_REF] Kresse | Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements[END_REF][START_REF] Kresse | From ultrasoft pseudopotentials to the projector augmented-wave method[END_REF][START_REF] Raybaud | Ab initio density functional studies of transition-metal sulphides: I[END_REF][24], where the projector augmented wave pseudopotential is used to replace the core electrons. The cut-off energy for the valence electrons was set at 500 eV. We included a Hubbard correction with a Ueff(Ce4f) of 4.2 eV to improve the electron localisation effects [25]. The inclusion of the Hubbard U term in the DFT+U approach corrects the resulting electronic structure of the reduced CeO2, favouring the insulating state rather than the incorrect metallic solution. However, this approach, in general, suffers from a strong linear dependence of the energetics on the value chosen for the parameter U and also on the localised projector functions chosen that define the U-dependent energy term. Here, to circumvent this problem, we verified that the energetics of the systems calculated using Ueff=4.2 eV agree with the previously calibrated results [START_REF] Fronzi | and reduced CeO2(111) surface: a first-principles investigation[END_REF], validating our choice of U. Since, in the CeO2 and H2O system, the effect of van der Waals correction is considerably smaller than the dominant localisation effect of f electrons [START_REF] Fernandez-Torre | Insight into the adsorption of water on the clean CeO2 (111) surface with van der Waals and hybrid density functionals[END_REF], the calculations were limited to the GGA+U framework.

We examined here the three principle low index surfaces of CeO2, corresponding to oxygen terminated (111), ( 110) and (100) planes. These surfaces were selected due to their significance in the surface-driven functionality of CeO2 across the various applications and growth morphologies of this material [START_REF] Yang | [END_REF][START_REF] Désaunay | Surface-dependent oxidation of H2 on CeO2 surfaces[END_REF]. We constructed these surfaces by cleaving the fluorite-type cubic CeO2 structure [START_REF] Loschen | study of cerium oxides: Dependence on the effective U parameter[END_REF]. The simulated supercells consisted of the cleaved surface slabs with a (2u×2v) expansion and a depth of three unit cells resulting in a 144 atom supercell for (100), a 96 atoms supercell for (110) and a 60 atom supercell for (111) surfaces. To find the ground state configuration of the adsorbed water layers on the CeO2 surface, we ran quenching ab initio molecular dynamics simulations with a target temperature of 20 K with steps of 0.1 fs. The molecular dynamics simulation was effective in finding the initial structures for geometry optimisation. Full geometry optimisation was then carried out on the equilibrated structures, with convergence criteria for the energy and forces of 10 -5 eV and 10 -2 eV/Å, respectively.

Results

Hydrophobicity

We calculated the adsorption energy (EAd) of a water double layer composed of eight molecules for (100) and (111) surfaces and sixteen molecules for the (110) surface by subtracting the total DFT energy of the supercell containing both the CeO2 slab and the water double layer from the summed total energies of the individual CeO2 surface slabs and the standalone water double layers. Since the (110) surface has a relatively larger surface cross-section (21.137 Å 2 ) than the (100) and (111) surfaces (14.947 Å 2 ), to fully cover the interface, a larger number of H2O molecules were required. Furthermore, (110) has three undercoordinated surface ions; two O ions and a Ce ion while (100) and (111) surfaces have only two undercoordinated surface ions; two O ions for the (100) surface, and an O ion and a Ce ion for the (111) surface. Consequently, (110) surface required twice as many water molecules to saturate all undercoordinated surface ions. The configurations of the water CeO2 surfaces and the double layers are presented in Figure 1. According to earlier works, the calculated adsorption energy of the water to the CeO2 surface does not significantly change if the thickness of the water slab is increased beyond two layers [START_REF] Carchini | Understanding and Tuning the Intrinsic Hydrophobicity of rare-earth oxides: A DFT+U Study[END_REF]. Based on this approximation we can use the optimised structures and their total energies calculated here to predict the contact angle between water and CeO2 surfaces. Young's Equation defines the contact angle between a given surface and a liquid, schematically shown in Figure 2(a), as follows [32]:

cos(𝛩) = 𝛾 𝑆𝑜𝑙𝑖𝑑-𝐺𝑎𝑠 -𝛾 𝑆𝑜𝑙𝑖𝑑-𝐿𝑖𝑞𝑢𝑖𝑑 𝛾 𝐿𝑖𝑞𝑢𝑖𝑑-𝐺𝑎𝑠 , ( 1
)
where γ is the interface energy. In our theoretical framework, these quantities are calculated in the absence of thermal excitation, and therefore we approximate here that the liquid at an atomistic level acts similarly to ice, an approximation justified experimentally by Van Oss et al. [START_REF] Van Oss | Determination of contact angles and pore sizes of porous media by column and thin layer wicking[END_REF]. Consequently, we consider the interfacial energies with ice rather than liquid water. As a second approximation, we use here vacuum conditions as a substitute for the gas phase, as simulation cells are cleaved by vacuum rather than water vapour. The magnitude of the error resulting from this approximation was shown to be minor by Carchini et al. [START_REF] Carchini | Understanding and Tuning the Intrinsic Hydrophobicity of rare-earth oxides: A DFT+U Study[END_REF]. By applying these two approximations; we obtain a new equation for the contact angle as cos(𝛩) ≃

𝛾 𝐶𝑒𝑂 2 -𝑣𝑎𝑐𝑢𝑢𝑚 -𝛾 𝐶𝑒𝑂 2 -𝑖𝑐𝑒 𝛾 𝑖𝑐𝑒-𝑣𝑎𝑐𝑢𝑢𝑚 . ( 2
)
The total energy of the supercell containing both the water double layer and ceria surface cell shown in Figure 2(b) can be decomposed as:

𝐸 𝐶𝑒𝑂 2 +𝑖𝑐𝑒 𝑇𝑜𝑡 = 𝐴 • 𝛾 𝐶𝑒𝑂 2 -𝑣𝑎𝑐𝑢𝑢𝑚 + 𝐴 • 𝛾 𝑖𝑐𝑒-𝑣𝑎𝑐𝑢𝑢𝑚 + 𝐴 • 𝛾 𝐶𝑒𝑂 2 -𝑖𝑐𝑒 + 𝑛𝐸 𝐼𝑐𝑒(𝐵𝑢𝑙𝑘) 𝑇𝑜𝑡 + 𝑚𝐸 𝐶𝑒𝑂 2 (𝐵𝑢𝑙𝑘) 𝑇𝑜𝑡 (3)
Here, 𝐸 𝐶𝑒𝑂 2 +𝑖𝑐𝑒 𝑇𝑜𝑡 is the total energy of the supercell containng the CeO2 and adsorbed water,

𝐸 𝑖𝑐𝑒(𝐵𝑢𝑙𝑘) 𝑇𝑜𝑡
is the total energy of bulk ice in stable hexagonal form, 𝐸 𝐶𝑒𝑂 2 (𝐵𝑢𝑙𝑘) 𝑇𝑜𝑡 is the total energy of bulk CeO2 and γ are interfacial energies indicated by their respective subscripts. The coeficients n and m correspond to the molecular ratios given by the stoichiometry of the supercell. A schematic representation of interface energies is given in Figure 2. 𝛾 𝐶𝑒𝑂 2 -𝑣𝑎𝑐𝑢𝑢𝑚 can be calculated from the difference of bulk and surface slabs as:

2𝐴 • 𝛾 𝐶𝑒𝑂 2 -𝑣𝑎𝑐𝑢𝑢𝑚 = 𝐸 𝐶𝑒𝑂 2 -𝑠𝑙𝑎𝑏 𝑇𝑜𝑡 -𝑚𝐸 𝐵𝑢𝑙𝑘-𝐶𝑒𝑂 2 𝑇𝑜𝑡 , (4)
where 𝐸 𝐶𝑒𝑂 2 -𝑠𝑙𝑎𝑏 𝑇𝑜𝑡 is the total energy of the CeO2 slab for a given surface. Similarly, 𝛾 𝑖𝑐𝑒-𝑣𝑎𝑐𝑢𝑢𝑚 , as illustrated in Figure 2(c), can be calculated as (5) On the basis of these interfacial energy values, we evaluate Equation ( 2) and obtain the adsorption energies and contact angles of water with the three surfaces as summarised in Table 1. The attraction between the water double layer and the CeO2 (110) surface is the strongest followed by the that of the (100) plane, with the (111) surface exhibiting the lowest adsorption energy. These results show that the (111) surface is strongly hydrophobic as its Θ value substantially exceeds 90°. The (110) surface can, on the other hand, be considered weakly hydrophilic as indicated by its calculated contact angle of Θ = 64.09° which is smaller than 90°. The (100) surface is slightly hydrophobic as Θ for this plane is marginally larger than 90°. The value obtained for the (111) surface is slightly larger than the earlier obtained Θ value of 99.9° calculated by Carchini et al. however it is closer to the experimentally obtained value (103°) reported in the same study [START_REF] Carchini | Understanding and Tuning the Intrinsic Hydrophobicity of rare-earth oxides: A DFT+U Study[END_REF]. To understand the origin of the different levels of hydrophobicity/hydrophilicity of CeO2 surfaces, we analysed the interaction between those hydrogen atoms closest to the CeO2 surface and the surface oxygen atoms by calculating the crystal orbital overlap populations (COOP). COOP is calculated by multiplying overlapping electronic populations by their corresponding density of states. Positive and negative values indicate bonding and antibonding interactions respectively. The COOP plots, presented in Figure 3, show strong bonding of H from the closest adsorbed water layer to the lattice oxygen at the (110) surface as indicated by the peak at ~-5.5 eV. For the (100) surface, the interaction of H with O results in more moderate bonding states at ~-2 eV and ~-7 eV. At the (111) surface, on the other hand, there are only minor bonding states close to the Fermi level with peaks at ~-1 eV with no strong bonding peaks anywhere deeper in the valence band. We thus observe here a correlation between the contact angle and the intensity of the bonding states, where more intense H-O bonding peaks deeper in the valence band indicate higher hydrophilic water ceria interaction. Indeed, the (110) surface has the most intense H-O bonding peak and the (111) the least intense, in line with the calculated contact angles that indicate a hydrophilic and hydrophobic behaviour for the (110) and the (111) surfaces, respectively. Furthermore, the O-H interaction at (100) and (110) surfaces results in antibonding states at ~-1 eV indicating that the surface To analyse the nature of bonding in the CeO2 outermost layer at the water/CeO2 interface, we calculated the COOP corresponding the overlap between neighbouring surface Ce-Ce and Ce-O pairs in the CeO2 slab as presented in Figure 4. Accordingly, the Ce-Ce COOP plots for all surfaces are dominated by antibonding states for energies between ~-5 eV and ~-2 eV. Additionally, there are two minor bonding Ce-Ce peaks for (111) at ∼-1.5 eV and ∼-4.3 eV. For the (110) surface, there is one minor bonding region at -5.0 eV. The Ce-Ce COOP for the (100) surface, however, shows no bonding states. Given that the Ce-Ce interaction is dominated by antibonding states, surface stability must, therefore, be driven by Ce-O or O-H interactions. The Ce-O interactions at the (110) and (111) surfaces generally show strong bonding character for deeper energies of the valence band (E <~-1.5 eV). The (100) surface, on the other hand, has two broad antibonding states at E = ~-1.5 eV and E = ~-3 eV that suggest the instability of this termination in vacuum conditions in the absence of O-H bonding states when fully covered with the H2O layer as shown in Figure 3.

Surface Thermodynamic Analysis

As Ce 4+ can be readily reduced to Ce 3+ in CeO2 alongside the liberation of O, the effects of such vacancies merit consideration and have been partially examined experimentally by high-resolution scanning force microscopy [34]. Here we calculated the formation energy of an oxygen vacancy at the (100), ( 110) and (111) surfaces according to the following equation:

where 𝐸 𝑆𝑢𝑟𝑓+𝑉 𝑂 𝑇𝑜𝑡 is the total energy of the surface following vacancy formation (reduced surface), 𝐸 𝑆𝑢𝑟𝑓 𝑇𝑜𝑡 is the total energy of the pristine (stoichiometric) surface and 𝐸 𝑂 2 𝑇𝑜𝑡 is the total energy of an oxygen molecule calculated in a vacuum box. The results are presented in the first row of Table 2. The formation energy of VO on the (100) surface has a negative value which indicates that bare O terminated (100) surface is not stable without water adsorption, as implied also by the COOP analysis of the preceding section. Specifically, given the two broad antibonding Ce-O COOP peaks in (100) shown in Figure 4(a), a negative formation energy for VO is not surprising at this surface. However, the formation energy of vacating the same O site at the (100) surface covered with the water double layer, marked with a circle in Figure 1(a), is 0.608 eV, leading to the conclusion that O terminated (100) surface is stablised in aqueous conditions. The (110) and ( 111) surfaces, on the other hand, are stable under any conditions of water adsorption as the VO formation energy has a positive values on these surfaces.

As O terminated (100) surface cannot be stabilised by partial H2O covering, we focus on the interaction of a single water molecule with (110) and ( 111) surfaces under the stoichiometric and reduced conditions to understand how these surfaces behave when partially covered with water. Water adsorption energies were calculated according to the following equation:

𝐸 𝐴𝑑𝑠 = 𝐸 𝑆𝑢𝑟𝑓+𝐻 2 𝑂 𝑇𝑜𝑡 -𝐸 𝑆𝑢𝑟𝑓 𝑇𝑜𝑡 -𝐸 𝐻 2 𝑂 𝑇𝑜𝑡 , (7)
where 𝐸 𝑆𝑢𝑟𝑓+𝐻 2 𝑂 𝑇𝑜𝑡 is the total energy of the surface upon water adsorption, 𝐸 𝑆𝑢𝑟𝑓 𝑇𝑜𝑡 is the total energy of the initial surface which can be either stoichiometric or reduced, and 𝐸 𝐻 2 𝑂 𝑇𝑜𝑡 is the total energy of a water molecule calculated in a vacuum box. Furthermore, for every scenario, two different configurations were considered; first molecular (or intact) adsorption of the water molecule in which an intact water molecule is adsorbed at a Ce site and second, a dissociated water molecule which consists of a HO group adsorbed on Ce site and a H atom adsorbed on a nearby O site. The relaxed structures of all configurations are presented in Figure 5.

As shown in Table 2, by comparing the molecular water adsorption on both surfaces, we observe that EAds(H2O) is lower at the (110) surface. Furthermore, the formation energy of an O vacancy was significantly larger at CeO2 (111) surface than that of (110) surface. Higher VO formation energy, in turn, results in lower molecular and dissociative water adsorption energies onto the reduced (111) surface. That is because bonding between the reduced (111) surface with the O atom of the adsorbed water molecule compensates for the chemical instability that was created by the removal of an oxygen atom at the (111) surface. The Results obtained here are in general agreement with the earlier GGA calculations [START_REF] Fronzi | and reduced CeO2(111) surface: a first-principles investigation[END_REF]. We then calculated surface energies, within the framework of the ab-initio thermodynamics, with the inclusion of the chemical potentials of water and oxygen [START_REF] Fronzi | and reduced CeO2(111) surface: a first-principles investigation[END_REF][START_REF] Fronzi | CHx adsorption (x = 1-4) and thermodynamic stability on the CeO2(111) surface: a first-principles investigation[END_REF]. We considered chemical potentials as independent variables so that a wider range of conditions could be extrapolated. Assuming that the gas reservoir exchanges particles with the surface without affecting its overall chemical potential, the surface energy can be calculated as follows:

𝛾(𝑝, 𝑇) = 1 𝐴 [𝐺 -∑ 𝑁 𝑖 𝜇 𝑖 (𝑝, 𝑇)], (8)
where A is the interface area, G is the Gibbs free energy of the crystal, ∑ 𝑁 𝑖 𝜇 𝑖 (𝑝, 𝑇) is the chemical potential of the atomic species in the system and the term Ni is the total number of atoms of the species i. We assume that the interface is in equilibrium with an environment formed by a mixture of two different chemical species in the gas-phase, namely O2 and H2O, which are considered to be in nonequilibrium with one another. We measure the chemical potentials with respect to the calculated energy of the isolated molecules, and thus define the relative chemical potentials to be 𝛥𝜇 𝑖 (𝑝, 𝑇) = 𝜇 𝑖 -𝐸 𝑖 . By expanding Equation 8to obtain a relationship between 𝛾(𝑝, 𝑇), 𝛾 𝑂 (𝑝, 𝑇) and 𝛾 𝐻 2 𝑂 (𝑝, 𝑇), where water and oxygen chemical potential are treated as independent variables, we obtain phase diagrmas of stable surface configurations as shown in Figure 6 in which each colored area indicates which surface termination is stable with respect to the other competing surfaces for any given values of chemical potentials. The borders of different areas indicate the transitions between the stability of two surface terminations, and give us information on the conditions for such transitions. For CeO2(100) we obtain that water termination is stable under any environmental conditions, confirming the conclusions reached in the previous section. We therefore do not include this phase diagram in Figure 6, since it would not bring any useful information. Furthermore, we observe that for the (111) surface there is no direct transition between the region corresponding the stability of water adsorbed at the stoichiometric surface and that corresponding to the bare reduced surface (denoted by 111-VO). For the (110) surface there is no direct transition between the stability field of the pristine surface and that of water terminated reduced surface. These phase diagrams thus demonstrate that for the (111) surface the presence of water stabilises surface oxygen vacancies, shifting the stability of VO towards higher values of the oxygen chemical potential. This is shown in the top panel of Figure 6 where H2O-VO, water adsorbed on a reduced surface, is stable for higher ΔμO compared to the water adsorbed on the bare reduced surface that contains VO. This point can be utilised experimentally as the stability of oxygen vacancy results in high chemical activity for ceria surface through the Mars-van Krevelen mechanism in which redox reactions are mediated by oxygen vacancies [START_REF] Doornkamp | The universal character of the Mars and Van Krevelen mechanism[END_REF]. Furthermore, on reduced ceria surfaces, superoxides could, in principle, be formed when an electron, trapped at a Ce 3+ site, is transferred to an adsorbed O2 molecule which is a critical mechanism in a wide variety of reaction pathways [37]. Unlike the (111) surface, for CeO2 (110), the presence of water molecules does not stabilise surface VO, and therefore no increase in surface chemical activity is anticipated in the presence of high H2O chemical potential. Finally, we revisited the contact angles taking into account the presence of VO at CeO2 and H2O interface, marked with circles in Figure 1. We repeated the calculations of Section 3.1 using 4u×4v extended supercells for which one VO was introduced to all interfaces. The VO concentration in the extended supercell was 3.125% for (100), (110) surfaces and 6.25% for the (111) surface. The contact angles for the O deficient surfaces was calculated to be 118.59° for (100) surface, 74.30° for (110) surface and 160.15° for (111) surface. Although we anticipate that Θ is sensitive to the VO concentration, these calculations, nonetheless, demonstrate the O vacancy on the low index CeO2 surfaces, in general, increases the hydrophobicity.

Discussion

Numerous rare earth oxides have been found to exhibit intrinsic hydrophobicity, i.e. without the use of chemical surface treatments. Of these REOs, one can consider CeO2 to be the most attractive, due to its relative availability and diverse functionality in different forms and microstructures. The origins of hydrophobicity in CeO2 have been the subject of discussions in recent literature [START_REF] Azimi | [END_REF]38]. By density functional theory (DFT) and crystal orbital overlap population (COOP), we have demonstrated here that the hydrophilicity or hydrophobicity of the different surfaces of rare earth oxides, such as CeO2, can be understood in light of their electronic and crystallographic structures and the consequent bonding states. The weaker H-O bonding at the CeO2(111) surface in the form of minor shallow bonding states indicated in COOP analysis, corresponds to a highly hydrophobic surface, consistent with the numerous experimental studies on this material and confirming the large contact angle calculated here [START_REF] Kabir | Structure and properties of hydrophobic CeO2-x coatings synthesized by reactive magnetron sputtering for biomedical applications[END_REF][START_REF] Kropp | the (111) and (100[END_REF]. The deeper bonding states found in (100) and strong bonding states for (110) correspond to the greater wettability of these planes.

The dissociative adsorption of water at the (111) and (100) surfaces has been observed experimentally, with analysis showing dissociation is relatively more favoured on (100) surfaces, with dissociation being enhanced for both cases in reduced materials [START_REF] Mullins | 2(100) and CeO2(111) thin films[END_REF][START_REF] Chen | Reactivity of hydroxyls and water on a CeO2(111) thin film surface: The role of oxygen vacancy[END_REF], in qualitative agreement with the results of the present study. Further studies, by atom force microscopy, have also shown the dissociation of water to be driven by the presence of oxygen vacancies at surface sites on the lower energy (111) planes, as hydroxide substitutes for surface oxygen [34].

Numerous experimental studies, mainly related to CO oxidation, have observed poorer catalytic activities at (111) planes in CeO2 relative to the higher energy (100) and (110) planes [START_REF] Tana | Morphology-dependent redox and catalytic properties of CeO2 nanostructures: Nanowires[END_REF]. In contrast, as a support for catalytic nanoparticles (Au, Cu, Pt) in water gas shift reactions, (111) plane of ceria is found to impart enhanced performance [START_REF] Rodriguez | Water gas shift reaction on Cu and Au nanoparticles supported on CeO2(111) and ZnO(000-1): Intrinsic activity and importance of support interactions[END_REF][START_REF] Baron | gold with cerium oxide supports: CeO2(111) thin films vs CeOx Nanoparticles[END_REF]. Whether wettability of a surface contributes or detracts from a heterogeneous catalyst activity depends on diverse mechanistic aspects, and therefore the design of CeO2 structures remains non-trivial. However, the findings here, indicating the stabilisation of surface oxygen vacancies in the presence of molecularly adsorbed water on (111) surfaces, suggesting that for catalytic pathways in aqueous systems involving vacancy mediated redox reactions, the increased prominence of this plane (such as in octahedral nanocrystals) may be advantageous.

According to our findings, the negative value of the O vacancy formation energy on the polar (100) surface demonstrates that this surface is intrinsically unstable. Experimentally, various stabilising surface reconstructions have been reported that are highly dependent on the sample preparation and treatment history. Capdevila-Cortada et al. showed that surface disorder created by partial removal of O could kinetically stabilise the bare (100) surface through fast surface diffusion process among nearly degenerate VO configurations. [46]. Judging from the positive value of the VO formation energy on fully water covered (100) surface, we infer that full water coverage can also provide a stabilisation mechanism for CeO2 (100) surface.

It has also been argued that a larger average Ce-Ce distance at the interface results in higher levels of hydrophilicity, implying that the Ce geometry at the interface is the main driving force behind the wettability of the system [START_REF] Carchini | Understanding and Tuning the Intrinsic Hydrophobicity of rare-earth oxides: A DFT+U Study[END_REF]. In contrast, as shown in Table 1, in the most hydrophobic surface i.e. the (111) surface, the Ce atoms have six Ce nearest neighbours with a Ce-Ce distance of 3.86 Å, whereas for the most hydrophilic (110) surface, Ce has only two nearest neighbours with Ce-Ce distance of 3.86 and two second nearest neighbours at 5.46 Å. Also, in (100) surface Ce has only four nearest neighbours with a distance of 3.86 Å and four second nearest neighbours with a distance of 5.47 Å. Consequently, our results suggest that large Ce-Ce distances are positively correlated with hydrophobicity rather than hydrophilicity. Therefore the previously stated correlation between Ce-Ce distances and hydrophobicity does not constitute a general trend, and instead, a positive correlation can be established only through the analysis of the specific electronic configuration and COOP [START_REF] Carchini | Understanding and Tuning the Intrinsic Hydrophobicity of rare-earth oxides: A DFT+U Study[END_REF].

Conclusions

CeO2 (111) surface was found to be the most hydrophobic with Θ = 112.53° followed by (100) with Θ = 93.91°. In contrast, the CeO2 (110) surface, was found the most hydrophilic with Θ = 64.09°. The contact angles increase for all three surfaces when O vacancy is present. Moreover, the O terminated (100) surface was found to be unstable in the absence of full coverage by molecular H2O. In the absence of water, CeO2 (100) surfaces will be, at least, partially Ce terminated. This finding implies that CeO2 (100) surfaces formed in an aqueous environment (through sol-gel synthesis for example) are likely to exhibit different surface properties from materials prepared using vacuum based methods (e.g. pulsed laser deposition). Finally, the thermodynamic analysis demonstrates that water adsorption stabilises oxygen vacancies on CeO2 (111) surfaces, shifting the stability of VO towards higher values of the oxygen chemical potential. The increased concentration of VO at higher oxygen chemical potentials is expected to improve the catalytic activity of (111) surface, while the opposite trend was found for CeO2 (110).

Figure 1 .

 1 Figure 1. (Color online) Schematic representation of a water double layer adsorption onto CeO2(100) (a), CeO2(110) (b) and CeO2(111) (c) surfaces. Yellow, red and blue spheres represent cerium, oxygen and hydrogen atoms, respectively.

  Fronzi, M., Assadi, M. H. N., & Hanaor, D. A. H (2019). Theoretical insights into the hydrophobicity of low index CeO2 surfaces. Applied Surface Science. 478, 68-74 https://doi.org/10.1016/j.apsusc.2019.01.208. 6 oxygens dissociate the water molecules and form separated H and H-O groups. Water dissociation within the first layer of water is evident from the relaxed structures presented in Figure 1(a) and (b).

Figure 2 .

 2 Figure 2. Schematic representation of A) water adsorption on a CeO2 surface B), C) and D) the models used in our calculations. Here we indicate the interfaces used in Equations 2-5, where each number indicates a specific interface as follow:1 the solid-gas, 2 the liquid-gas, 3 the solid-liquid, 4 the CeO2-ice, 5 the CeO2-vacuum, 6 ice-vacuum.

Figure 3 . 7 Figure 4 .

 374 Figure 3. (Color online) COOP of CeO2 (100) (yellow line), CeO2 (110) (red line) and CeO2 (111) (blue line) corresponding to the closest H from the water double layer to the CeO2 surface and the outermost O from CeO2 surface. The vertical line at 0 eV represents the Fermi level energy.

9 Figure 5 .Table 2 .

 952 Figure 5. (Color online) Schematic representation of molecular and dissociative water adsorption onto CeO2(110) (left) and (111) (right) surfaces. From top to bottom, molecular adsorption onto stoichiometric surfaces (a) and (b), molecular adsorption onto reduced surfaces (c) and (d), dissociative adsorption onto stoichiometric surfaces (e) and (f) and dissociative adsorption onto reduced surfaces (g) and (h). Yellow, red and blue spheres represent cerium, oxygen and hydrogen atoms, respectively. Table 2. Vacancy formation as well as water molecular (associative) and dissociative adsorption energies at CeO2 (111) and CeO2 (110) surfaces (calculated with to the water in the gas phase), schematically represented in Fig.5. Adsorption energies were not calculated for (100) surface due to the lack of stability of O termination in this surface when not fully covered with H2O.

Figure 6 .

 6 Figure 6. (Color online) Surface phase diagram of stable structures of CeO2(111) and CeO2(110) in equilibrium with an oxidising and "humid environment", as a function of ΔμO and 𝚫𝝁 𝑯 𝟐 𝑶 in the gas phase.

  Fronzi, M., Assadi, M. H. N., & Hanaor, D. A. H (2019). Theoretical insights into the hydrophobicity of low index CeO2 surfaces. Applied Surface Science. 478, 68-74 https://doi.org/10.1016/j.apsusc.2019.01.208.

	4

  Fronzi, M., Assadi, M. H. N., & Hanaor, D. A. H (2019). Theoretical insights into the hydrophobicity of low index CeO2 surfaces. Applied Surface Science. 478, 68-74 https://doi.org/10.1016/j.apsusc.2019.01.208.

	2𝛾 𝑖𝑐𝑒-𝑣𝑎𝑐𝑢𝑢𝑚 = 𝐸 𝑖𝑐𝑒-𝑠𝑙𝑎𝑏 𝑇𝑜𝑡	-𝑚𝐸 𝐼𝑐𝑒(𝑏𝑢𝑙𝑘) 𝑇𝑜𝑡	.
			5

Table 1 . Water double layer adsorption energies (EAd) at CeO2 (100), (110) and (111) surfaces calculated with respect to the isolated water double layer (DL), and respective CeO2 slabs are presented in the first two rows. Values are given per unit area and per adsorbed molecule. Furthermore, the calculated Θ and the distance of Ce ions from their neighbouring Ce ions are at the CeO2/water interface are also given.

 1

	Configuration	(100)	(110)	(111)
	EAd(H2O-DL) (eV/Å 2 )	-0.023	-0.038	-0.014
	EAd(H2O-DL) (eV/H2O)	-0.345	-0.404	-0.179
	Θ	93.91°	64.09°	112.53°
	Ce-Ce	3.86 Å	3.86, 5.46 Å	3.86 Å

  Fronzi, M., Assadi, M. H. N., & Hanaor, D. A. H (2019). Theoretical insights into the hydrophobicity of low index CeO2 surfaces. Applied Surface Science. 478, 68-74 https://doi.org/10.1016/j.apsusc.2019.01.208.

	𝐸 𝑉 𝑂 = 𝐸 𝑆𝑢𝑟𝑓+𝑉 𝑂 𝑇𝑜𝑡	-𝐸 𝑆𝑢𝑟𝑓 𝑇𝑜𝑡 +	1 2 𝐸 𝑂 2 𝑇𝑜𝑡 ,
			8

  Fronzi, M., Assadi, M. H. N., & Hanaor, D. A. H (2019). Theoretical insights into the hydrophobicity of low index CeO2 surfaces. Applied Surface Science. 478, 68-74 https://doi.org/10.1016/j.apsusc.2019.01.208.

	13

Acknowledgements

The financial support was provided by the National Natural Science Foundation of China (Grant No. 51323011). Computational resources were provided by the Center for Computational Sciences at the University of Tsukuba.


Thomas Milcent
email: thomas.milcent@u-bordeaux.fr

Emmanuel Maitre
email: emmanuel.maitre@imag.fr

EULERIAN MODEL OF IMMERSED ELASTIC SURFACES WITH FULL MEMBRANE ELASTICITY *

Keywords: Level set, fluid structure, Eulerian, immersed interface, Navier-Stokes, membrane elasticity AMS subject classifications. 74F10, 76D05, 65M06

We introduce an Eulerian model for the coupling of a fluid governed by the Navier-Stokes equations, with an immersed interface endowed with full membrane elasticity (i.e., including shear effects). We show numerical evidences of its ability to account for large displacements/shear in a relatively simple way, avoiding some drawbacks of Lagrangian representation.

Introduction

In biomechanics of immersed membranes, the ability for a numerical model and method to handle shear efficiently is a challenging issue. Human red blood cells (RBCs) are among the simplest animal cells, since they have no nucleus nor organelle. They are made, like vesicles, of a phospholipid bilayer, plus a protein network (of spectrin), the cytoskeleton. For simple vesicles, the involved energy is a curvature energy, which is minimized with an area and enclosed volume constraint. For that case, Eulerian and Lagrangian models were developed in the last decade and proved to successfully reproduce the dynamics of such objects in flow [START_REF] Beaucourt | Steady to unsteady dynamics of a vesicle in a flow[END_REF][START_REF] Thierry Biben | Three-dimensional vesicles under shear flow: Numerical study of dynamics and phase diagram[END_REF][START_REF] Thierry Biben | Phase-field approach to three-dimensional vesicle dynamics[END_REF][START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF][START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF][START_REF] Doyeux | Simulation of two-fluid flows using a finite element/level set method. Application to bubbles and vesicle dynamics[END_REF]. For the red blood cells, the cytoskeleton provides extra resistance to membrane shear. Capsules are usually constituted of a liquid drop protected by a thin elastic membrane. They are used for applications in cosmetic, food, and pharmaceutical industries. They are more rigid objects than vesicles, so the curvature energy is small compared to membrane elasticity, including shear. In the modeling of vesicles or capsules, taking into account shear efficiently is of paramount interest.

On one hand, the definition of local shear variation on a surface is apparently easy in Lagrangian coordinates and, in the context of capsules, was formalized decades ago [START_REF] Barthes-Biesel | Motion of a spherical microcapsule freely suspended in a linear shear flow[END_REF][START_REF] Barthes-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF][START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF]. However, when it comes to practical computations, the unavoidable singularity of Lagrangian parametrization of closed surfaces brings high complications [START_REF] Walter | Ellipsoidal capsules in simple shear flow: prolate versus oblate initial shapes[END_REF][START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF]. Moreover, large deformation and volume conservation in this Lagrangian setting also raise numerous difficulties, which were studied by several authors in the framework of the immersed-boundary method [START_REF] Lee | An immersed interface method for incompressible navier-stokes equations[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]. On the other hand, interface-capturing methods using level-set, while more intrinsic and dimension-independent regarding surface localization, are usually reported as unable to capture tangential motions (see, however, some connected studies on surface flows and multicomponent vesicles [START_REF] Knut | A diffuse-interface method for two-phase flows with soluble surfactants[END_REF][START_REF] John S Lowengrub | Phase-field modeling of the dynamics of multicomponent vesicles: Spinodal decomposition, coarsening, budding, and fission[END_REF][START_REF] Xu | A level-set method for interfacial flows with surfactant[END_REF]). This is due to the fact that the transport equation for the level-set function, namely ∂ t φ + u • ∇φ = 0, while recording some information on the surface area variations (enabling a complex fluid formulation of fluid-interface coupling [START_REF] Beale | Locally corrected semi-lagrangian methods for stokes flow with moving elastic interfaces[END_REF][START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF][START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF]), just ignore any tangential component of u, since ∇φ is normal to the interface. However, full 3D Eulerian elasticity and its coupling with fluid mechanics has already successfully been studied and implemented numerically by several teams both in the incompressible [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF][START_REF] Maitre | Applications of level set methods in computational biophysics[END_REF][START_REF] Richter | Finite elements for fluid-structure interaction in ALE and fully Eulerian coordinates[END_REF][START_REF] Richter | A fully Eulerian formulation for fluid-structure-interaction problems[END_REF][START_REF] Valkov | Eulerian method for fluid-structure interaction and submerged solid-solid contact problems[END_REF][START_REF] Zhang | An Eulerian numerical method for fluidsolid interaction[END_REF] and compressible [START_REF] Gorsse | A simple cartesian scheme for compressible multimaterials[END_REF][START_REF] Kamrin | Reference map technique for finitestrain elasticity and fluid-solid interaction[END_REF] cases. The method relies on backward characteristics (also called a reference map by some authors) which map back material points to their initial positions. This concept was also introduced early for point correspondence in level-set applied to image processing [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF] and more recently to improve numerical accuracy of advection equations (e.g. level-set motions) [START_REF] Kohno | A new method for the level set equation using a hierarchical-gradient truncation and remapping technique[END_REF][START_REF] Mercier | The characteristic mapping method for the linear advection of arbitrary sets[END_REF].

This article is the first attempt, to the best of our knowledge, to write a full Eulerian model of an immersed interface endowed with full membrane elasticity (i.e. including shear effects) and to provide numerical evidences of its ability to account for large displacements/shear in a relatively simple way.

In Section 2, we recall the basic definitions of interface capturing via level-set and backward characteristics. In Section 3, we provide simple examples illustrating the possible drawbacks of a Lagrangian approach to record shear variations. In Section 4, we build two invariants that will be used to record the full membrane elasticity. The first one records area change (in both compressible and incompressible settings), while the second measures shear. In Section 5, we consider the fluid-structure coupling problem: an immersed membrane into an incompressible fluid. The whole model is reduced, following [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF], to a Navier-Stokes equation with source terms expressed as a combination of partial derivatives of the backward characteristics, coupled with the vectorial transport equation giving these characteristics. We provide numerical illustrations of the ability of the numerical implementation of our model to simulate the dynamic relaxation of an initially sheared sphere. As the notion of shear is best understood through simple tests, we provide in Appendix A analytical examples of velocity fields under which the deformed continuous medium actually experiences a shear, and show that our invariants behave as expected. Appendix B deals with some technical lemmas.

Preliminaries

Forward and backward characteristics.

Let Ω 0 ⊂ R 3 be the reference configuration of a continuous medium and assume that this medium is deformed by a smooth map X : R 3 × R + -→ R 3 (the forward characteristics) to Ω t = X(Ω 0 ,t). A velocity field u : R 3 × R + -→ R 3 is naturally associated with X:

∂ t X(ξ,t) = u(X(ξ,t),t), X(ξ,0) = ξ, ξ ∈ Ω 0 . (2.1)
We introduce the backward characteristics Y : R 3 × R + -→ R 3 by the formula Y (X(ξ,t),t) = ξ (see Figure 2.1). The physical interpretation of Y (x,t) is the position at time 0 of a material particle lying in x at time t and moving at speed u. The derivative of this relation with respect to t and ξ in turn, gives with (2.1),

∂ t Y + (u • ∇ x )Y = 0, Y (x,0) = x, x ∈ Ω t , (2.2) and [∇ ξ X(ξ,t)] = [∇ x Y (x,t)] -1 for x = X(ξ,t). (2.3)
In the following, we use the notation ∇ ξ for the gradient with respect to ξ and we denote by ∇ the gradient with respect to x.

Eulerian representation of surfaces.

We consider a surface Γ t captured by a level set function φ :

R 3 × R + -→ R Γ t = {x ∈ R 3 / φ(x,t) = 0}.
The initial position of the interface Γ 0 is associated to a given function φ 0 (thus we set φ(•,0) = φ 0 ). We assume that the interface is advected by the velocity field u, such that φ verifies

∂ t φ + u • ∇φ = 0.
(2.4)

We have, with (2.2), that the solution of (2.4) is given by φ(x,t) = φ 0 (Y (x,t)).

(2.5)

Let the normal in the reference configuration be denoted by n 0 (ξ) for ξ ∈ Γ 0 and the normal in the deformed configuration by n(x,t) for x ∈ Γ t (see Figure 2.1). These normals are expressed in terms of the normalized gradient of the associated level set:

n(x,t) = ∇φ(x,t) |∇φ(x,t)| , n 0 (ξ) = ∇ ξ φ 0 (ξ) |∇ ξ φ 0 (ξ)| . (2.6)
The relation (2.5) gives ∇φ(x,t) = [∇Y (x,t)] T ∇φ 0 (Y (x,t)). Therefore, with the definitions (2.6),

n 0 (Y (x,t)) = [∇Y (x,t)] -T n(x,t) |[∇Y (x,t)] -T n(x,t)| . (2.7) With (2.
3), we get the Lagrangian equivalent of (2.7):

n(X(ξ,t),t) = [∇ ξ X(ξ,t)] -T n 0 (ξ) |[∇ ξ X(ξ,t)] -T n 0 (ξ)| . (2.8) Y (x,t) = ξ x = X(ξ,t) n 0 (ξ) n(x,t) Initial configuration Ω 0 Deformed configuration Ω t Γ 0 Γ t

Can we measure surface shear variation with a Lagrangian description of the surface?

We present in this section two intuitive methods to measure shear variation with a parametrization of the surface. The first is based on the angle between two normalized vectors, and the second uses the invariants of the metric tensor. We will show that, even in the simplest case of a plane deformed in a shear velocity field, the two methods fail because the shear variation depends on the choice of initial parametrization. For general closed surfaces, more difficulties arise from the singularity of Lagrangian parametrizations. The first method cannot be applied because it is impossible to construct the normalized vectors (hairy ball theorem). We will show that, for a sphere deformed in a shear velocity field, the invariant of the second method is discontinuous at the poles. In any case, the definition of shear depends completely on the choice of the parametrization, and hence there is no canonical reference state. The limitations of these methods advocate for an immersion of the surface in R 3 using a Eulerian description to capture shear.

Two Lagrangian methods to capture shear.

Let Γ 0 be the initial surface described by the parametrization Ψ 0 : U -→ R 3 , where U is an open subset of R 2 . A parametrization of the deformed surface Γ t is given by Ψ :

U × R + -→ R 3 , where Ψ(θ 1 ,θ 2 ,t) = X(Ψ 0 (θ 1 ,θ 2 ),t
) and the deformed vectors associated to the parametrization are given by a i (t) = ∂ θi Ψ(θ 1 ,θ 2 ,t). The first intuitive method to capture shear is to consider the angle between the two normalized deformed vectors

angle(t) = arccos a 1 (t) |a 1 (t)| • a 2 (t) |a 2 (t)| .
The second is to consider the metric tensor:

M(t) := a 1 (t) • a 1 (t) a 1 (t) • a 2 (t) a 1 (t) • a 2 (t) a 2 (t) • a 2 (t)
which is classically used to compute area, angles, and length of curves on the surface.

The trace of M(t) captures shear deformation but also contains information on local area. Dividing the trace by the square root of the determinant keeps only shear effects, as will be illustrated below. We thus introduce the following invariants:

Z 1 (t) = det(M(t)), Z 2 (t) = Tr(M(t)) 2 det(M(t)) , ( 3.1)
where Z 1 (t) correspond to the local area, and we will show that Z 2 (t) is a good Lagrangian candidate to capture the local shear. In the following, we are interested in the local variation of area and shear between time 0 and t. Therefore, we consider the quantity angle(t)/angle(0) and the invariants Z i (t)/Z i (0). In the following examples we consider the shear velocity field (see Figure A.10) with the associated forward characteristics given by

u(x,y,z,t) = ⎛ ⎝ 0 x 0 ⎞ ⎠ , X(ξ 1 ,ξ 2 ,ξ 3 ,t) = ⎛ ⎝ ξ 1 ξ 2 + tξ 1 ξ 3 ⎞ ⎠ .
3.2. Case of the plane: orthogonal parametrization. Let Γ 0 be the plane {(x,y,z) ∈ R 3 /z = 0} which is parametrized by Ψ 0 : R 2 -→ R 3 given by

Ψ 0 : (θ 1 ,θ 2 ) → ⎛ ⎝ cos(α)θ 1 -sin(α)θ 2 sin(α)θ 1 + cos(α)θ 2 0 ⎞ ⎠ .
where α ∈ [0,2π] is a parameter. The deformed surface Γ t is still geometrically Γ 0 in this case, and the deformed vectors are given by

a 1 (t) = ⎛ ⎝ cos(α) sin(α) + tcos(α) 0 ⎞ ⎠ , a 2 (t) = ⎛ ⎝ -sin(α) cos(α) -tsin(α) 0 ⎞ ⎠ .
At t = 0, these vectors are orthonormal and have an α angle with the canonical basis (see Figure 3.1). We can show that, for this example, we have angle(t) angle(0) = 2 π arccos t(2cos(2α)tsin(2α))

2 (1 + tsin(2α) + t 2 cos(α) 2 )(1 -tsin(2α) + t 2 sin(α) 2 ) . (3.2)
This quantity depends on the angle parameter α (even if all pair of initial vectors are orthonormal). Hence, this first method is not appropriate to measure the local shear variation. For the second method, the invariants (3.1) in the case of the plane are given by

Z 1 (t) Z 1 (0) = 1, Z 2 (t) Z 2 (0) = 1+ t 2 2 . (3.3)
Hence, there is no area variation and the shear is constant in space and increases with time, which correspond intuitively to what we can expect of a shear deformation. The second invariant is independent of α, so this second method seems appropriate to capture shear.

• α a 1 (0) a 2 (0) • β a 1 (0) a 2 (0)

Case of the plane: non orthogonal parametrization.

Let Γ 0 be the same plane {(x,y,z) ∈ R 3 /z = 0} described with another parametrization Ψ 0 : R 2 -→ R 3 given by

Ψ 0 : (θ 1 ,θ 2 ) → ⎛ ⎝ θ 1 + cos(β)θ 2 sin(β)θ 2 0 ⎞ ⎠ ,
where β ∈ [0,2π] is a parameter. The deformed vectors are given by

a 1 (t) = ⎛ ⎝ 1 t 0 ⎞ ⎠ , a 2 (t) = ⎛ ⎝ cos(β) sin(β) + tcos(β) 0 ⎞ ⎠ .
At t = 0, these vectors are not orthogonal and have a β angle (see Figure 3.1). For the second method, the invariants (3.1) are given by

Z 1 (t) Z 1 (0) = 1, Z 2 (t) Z 2 (0) = 1 2 1+ t 2 + cos(β) 2 + (tcos(β) + sin(β)) 2 . (3.4)
There is no area variation, and the shear depends on the parameter β. Therefore, the second method fails to capture local shear variation because it depends on the choice of the parametrization.

3.4.

Case of the sphere: orthogonal parametrization. Let Γ 0 be the sphere

{(x,y,z) ∈ R 3 /x 2 + y 2 + z 2 = 1} which is parametrized by Ψ 0 : [0,π] × [0,2π] -→ R 3 given by Ψ 0 : (θ 1 ,θ 2 ) → ⎛ ⎝ sin(θ 1 )cos(θ 2 ) sin(θ 1 )sin(θ 2 ) cos(θ 1 ) ⎞ ⎠ .
The deformed vectors are given by

a 1 (t) = ⎛ ⎝ cos(θ 1 )cos(θ 2 ) cos(θ 1 )sin(θ 2 ) + cos(θ 1 )cos(θ 2 )t -sin(θ 1 ) ⎞ ⎠ , a 2 (t) = ⎛ ⎝ -sin(θ 1 )sin(θ 2 ) sin(θ 1 )cos(θ 2 ) -sin(θ 1 )sin(θ 2 )t 0 ⎞ ⎠ .
At t = 0, these vectors are orthogonal (see Figure 3.2). The first method with angle between vectors also fails in this more complex case. Indeed, it is impossible to construct a continuous tangent vector field on a sphere (hairy ball theorem). For the second method, the invariants (3.1) in the case of the sphere are given by

Z 1 (t) Z 1 (0) = 1 -tsin(θ 1 ) 2 sin(2θ 2 ) + t 2 sin(θ 1 ) 2 sin(θ 2 ) 2 , (3.5) Z 2 (t) Z 2 (0) = 1 + sin(θ 1 ) 2 + tcos(2θ 1 )sin(2θ 2 ) + t 2 2 (1 + cos(2θ 1 )cos(2θ 2 )) (1 + sin(θ 1 ) 2 ) 1 -tsin(θ 1 ) 2 sin(2θ 2 ) + t 2 sin(θ 1 ) 2 sin(θ 2 ) 2 . ( 3.6)
The first invariant (3.5) is well defined and represents the local area variation. The limit at the poles (θ

1 = 0 or θ 1 = π) of the second invariant (3.6) is cos(θ 2 ) 2 + (sin(θ 2 ) + tcos(θ 2 )
) 2 and depends of θ 2 . Therefore, there is a discontinuity of the second invariant at the poles. We can reparametrize periodically the surface, but the algorithms are more complex. Moreover, the results will always depend on the choice of the parametrization, even in zones where there is no singularity.

To tackle this problem, we propose immersing the surface in R 3 and projecting on it the 3D deformations. In this case, we use the flat metric of R 3 instead of the metric of the surface and thus avoid the singularities of parametrizations. Also, the immersion in R 3 will allow us to get a reference state, which does not depend on the parametrization. Furthermore, we will introduce Z i , the Eulerian equivalent of the Lagrangian invariants Z i . We will show that Z 1 is the classical and well-known local-area variation. We will demonstrate with a variety of analytical illustrations in Appendix A that Z 2 is a good candidate to record the local shear variation. We want to measure the deformations on the surface Γ t . Following the Lagrangian description of [START_REF] Barthes-Biesel | Motion of a spherical microcapsule freely suspended in a linear shear flow[END_REF][START_REF] Barthes-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF][START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF], we introduce the tensor

M (X(ξ,t),t) := [∇ ξ X(ξ,t)][I -n 0 (ξ) ⊗ n 0 (ξ)]. ( 4.1)
Let v(ξ) be a vector defined at the point ξ ∈ Γ 0 . This vector is first projected with [In 0 (ξ) ⊗ n 0 (ξ)] to v τ (ξ) ∈ T ξ Γ 0 , the tangent plane of Γ 0 at ξ. Then the vector v τ (ξ) is deformed with X in the vector [∇ ξ X(ξ,t)]v τ (ξ) at X(ξ,t). A first property is that Mv is already in T X(ξ,t) Γ t , the tangent plane of Γ t at X(ξ,t). Indeed, using (2.8) and

v τ (ξ) • n 0 (ξ) = 0, we have (M (X(ξ,t),t)v(ξ)) • n(X(ξ,t),t) = ([∇ ξ X] T n(X(ξ,t),t)) • v τ (ξ)) = 0.
The tensor (4.1) is written in its Eulerian form with (2.3)

M (x,t) := [∇Y (x,t)] -1 [I -n 0 (Y (x,t)) ⊗ n 0 (Y (x,t))].
The associate Cauchy-Green tensor is defined as follows (In 0 ⊗ n 0 is a projector hence involutive):

A := MM T = [∇Y ] -1 (I -n 0 (Y ) ⊗ n 0 (Y ))[∇Y ] -T .
After introducing the right Cauchy-Green tensor

B = [∇ ξ X][∇ ξ X] T = [∇Y ] -1 [∇Y ] -T , (4.2) using (2.7), the relation A(v ⊗ v)A T = (Av) ⊗ (Av), and |[∇Y ] -T n| 2 = (Bn) • n, we get A = B - (Bn) ⊗ (Bn) (Bn) • n . (4.3)

Invariants of A.

According to (4.3), we have

An = 0. (4.4)
Thus, 0 is an eigenvalue and det(A) = 0. The other invariants are Tr(A) and Tr(Cof(A)) = 1 2 (Tr(A) 2 -Tr(A 2 )). As A is real and symmetric, there exists an orthonormal basis of eigenvectors. Moreover, A is positive since Ax • x = |M T x| 2 ≥ 0, so we denote its eigenvalues by 0, λ 2 1 and λ 2 2 . Therefore, Tr(A) = λ 2 1 + λ 2 2 and Tr(Cof(A)) = (λ 1 λ 2 ) 2 . We introduce the following invariants:

Z 1 = Tr(Cof(A)) = |λ 1 λ 2 |, ( 4.5
)

Z 2 = Tr(A) 2 Tr(Cof(A)) = 1 2
λ 1 λ 2 + λ 2 λ 1 . (4.6)
In the reference configuration (usually at initial time in our case), we have that A(0) = In 0 ⊗ n 0 so Tr(A(0)) = 2 and Tr(Cof(A(0))) = 1. Under the smoothness assumption on u, there exists a time T that for t < T , Tr(A) and Tr(Cof(A)) do not vanish. Therefore, the invariant Z 2 is well defined. We have that Z 1 = Z 2 = 1 at initial time and the inequalities Z 1 ≥ 0, Z 2 ≥ 1. Note that Z 1 and Z 2 are both Eulerian fields, and in 5.2 we will derive equations to compute them. We will prove in Appendix B that Z 1 is exactly equal to the local area variation of the surface (even if incompressibility is not assumed on u). We will demonstrate in Appendix A that Z 2 is a good candidate to record the local shear variation.

Immersed membrane in a fluid

In this section, we apply the above considerations to build an Eulerian fluidstructure model of an immersed elastic interface, which accounts for the full membrane elasticity through an energy involving invariants Z 1 and Z 2 . We derive Eulerian equations verified by the two invariants, and the forces induced by the energy variation during the structure motion. The final model is recast as a Navier-Stokes system with a source term, the latter depending on space derivatives of Z 1 and Z 2 , coupled with their Eulerian equations. This is a generalization to full membrane elasticity of a simple level-set formulation of the immersed boundary method that was introduced previously and that we recall now.

Mathematical model.

We consider an elastic membrane Γ t immersed into a fluid governed by, for example incompressible Navier-Stokes equations, lying in a bounded domain Ω of R 3 , for a time interval [0,T ]. The interface Γ t is captured by a level-set function φ : Ω× [0,T ] → R. Based on previous works on Eulerian description of immersed interface [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF][START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF], we write the coupling between the immersed surface and the surrounding fluid under the formulation

⎧ ⎪ ⎨ ⎪ ⎩ ∂ t u + u • ∇u -1 Re Δu + ∇p = f (φ) on Ω× [0,T ] divu = 0 on Ω × [0,T ] ∂ t φ + u • ∇φ = 0 on Ω × [0,T ]
where f (φ) accounts for the elastic force imparted on fluid by the immersed surface. Dependence on φ means here dependence on φ and its derivatives. We have shown that, in the incompressible case, |∇φ| captures the local area variation and that the force associated to the regularized energy

E = Ω E(|∇φ|) 1 ε ζ φ ε dx (5.1)
can be expressed by the formula

f (φ) = div E (|∇φ|)|∇φ| I - ∇φ ⊗ ∇φ |∇φ| 2 1 ε ζ φ ε (5.2)
where r → E(r) is the elastic constitutive law, ε is the width of the interface, and ζ is a cut-off function used to spread the interface near {φ = 0}. However, in order to capture the full membrane energy, the function φ is not sufficient. Indeed a velocity field u tangential to Γ t would verify u • ∇φ = 0 and therefore would not change the value of φ in its evolution equation. In order to capture some information on how points are moving on the surface, one way is to use the backward characteristics Y of the velocity field u.

We propose the new model:

⎧ ⎪ ⎨ ⎪ ⎩ ∂ t u + u • ∇u -1 Re Δu + ∇p = F 1 (φ,Y ) + F 2 (φ,Y ) on Ω× [0,T ] divu = 0 on Ω × [0,T ] ∂ t Y + u • ∇Y = 0 on Ω × [0,T ] (5.3)
where F i is the force associated with the invariant Z i , which expressions will be derived in the next section. In most cases, equations (5.3) will be supplemented with homogeneous initial and boundary conditions on the velocity, and with initial condition on Y given by the identity

u = 0 on ∂Ω × [0,T ], u= 0 on Ω × {0}, Y (0,x) = x on Ω.
Under such initial and boundary conditions, characteristics remain for all time inside Ω. Thus the initial and deformed whole continuous medium (fluid and structure) always occupies the same domain of R 3 (thus there is no need to introduce an initial and deformed configuration as in 2.1).

Moreover, with such initial and boundary conditions, there is no need to solve an equation on φ, since the knowledge of Y and φ 0 is sufficient. Indeed, the solution to the transport equation on φ is given by the following (see (2.5)):

φ(x,t) = φ 0 (Y (x,t)).
However if we consider an initially deformed interface (see (5.20) in the numerical tests), we could choose to set Y (0,•) to some maps not equal to the identity. In that case, we need to solve the transport equation on φ

∂ t φ + u • ∇φ = 0.
(5.4)

Remark 5.1. One could wonder about the optimality of our representation, regarding the number of fields to capture full membrane deformation. Note that, to capture the interface position itself, we need one field (a level-set function φ). In the general case, to record the area change, we need to compute Z 1 (or |∇φ| and J; see (B.1) in Appendix B), and therefore an extra field to get the shear. Thus, as we solve for Y , which has three components, and obtain φ, Z 1 and Z 2 from it, we are optimal. In the incompressible case, area change is recorded by |∇φ| alone. We are not able, for the moment, to express shear by only introducing another field that we could compute without involving Y .

Eulerian equations on Z i .

In order to build an Eulerian model, we need to find evolution equations associated with Eulerian invariants. Proposition 5.1. Under the smoothness assumption made on u, the invariants verify

∂ t Z 1 + u • ∇Z 1 = Z 1 [∇u] : C 1 , C 1 = I -n ⊗ n, (5.5) ∂ t Z 2 + u • ∇Z 2 = Z 2 [∇u] : C 2 , C 2 = 2A Tr(A) -(I -n ⊗ n).
(5.6)

Proof. We first focus on the equations related to A. Taking the gradient of (2.2), we get

∂ t [∇Y ] + u • ∇([∇Y ]) = -[∇Y ][∇u],
and its inverse verifies

∂ t ([∇Y ] -1 ) + u • ∇([∇Y ] -1 ) = [∇u][∇Y ] -1 .
(5.7)

The equation (4.2) on B is then given by

∂ t B + u • ∇B = [∇u]B + B[∇u] T . (5.8)
Using (2.2) we get that

∂ t (I -n 0 (Y ) ⊗ n 0 (Y )) + u • ∇(I -n 0 (Y ) ⊗ n 0 (Y )) = 0.
This relation together with (5.7) gives

∂ t M + u • ∇M = [∇u]M, ∂ t M T + u • ∇(M T ) = M T [∇u] T .
(5.9)

Using the relations (5.9), we get

∂ t A + u • ∇A = [∇u]A + A[∇u] T , ( 5.10)
where the initial condition is given by A(0) = In 0 ⊗ n 0 . Note that B verifies the same equation (5.8) with a different initial condition B(0) = I. Following (5.10), we get

∂ t Tr(A) + u • ∇Tr(A) = 2[∇u] : A (5.11)
and

∂ t Tr(Cof(A)) + u • ∇Tr(Cof(A)) = 2[ATr(A) -A 2 ] : [∇u].
(5.12) Using (5.12), we get

∂ t Z 1 + u • ∇Z 1 = Z 1 [ATr(A) -A 2 ] Tr(Cof(A)) : [∇u].
We show in Lemma B.1 of Appendix B that ATr(A) -A 2 = Tr(Cof(A))[In ⊗ n]. Using (5.11) and the equation on Z 1 , we get the equation for Z 2 . We have that C i n = 0 because An = 0 (see 4.4). We will prove in Proposition B.3 of Appendix B that Z 1 captures the local area variation. In Appendix A, we provide many illustrations to explain why Z 2 is a relevant measure of shear variation.

Computation of the elastic forces F i .

We introduce the regularized energy

E i = Ω E i (Z i ) 1 ε ζ φ ε dx. (5.13)
Here, E i is the constitutive law associated to the invariant Z i .

Proposition 5.2. The time variation of E i , using the principle of virtual power

∂ t E i = - Ω F i • u dx (5.14)
corresponds to the following force:

F i = ∇ E i (Z i ) 1 ε ζ φ ε + div E i (Z i )Z i C i 1 ε ζ φ ε .
(5.15)

Proof. The derivative with respect to t gives

∂ t E i = Ω E i (Z i )(Z i ) t 1 ε ζ φ ε dx + Ω E i (Z i ) 1 ε 2 ζ φ ε φ t dx.
Using the transport equation on φ and (5.5), (5.6) we get

∂ t E i = Ω E i (Z i )(-u • ∇Z i + [∇u] : Z i C i ) 1 ε ζ φ ε dx + Ω E i (Z i ) 1 ε 2 ζ φ ε (-u • ∇φ) dx.
Integrating the second term by parts gives the following (the integral on ∂Ω vanishes since ζ( φ ε ) = 0 on ∂Ω):

∂ t E i = - Ω u • ∇(E i (Z i )) 1 ε ζ φ ε + div E i (Z i )Z i C i 1 ε ζ φ ε • u + E i (Z i )u • ∇ 1 ε ζ φ ε dx.
Grouping the first and last term and using (5.14) leads to the expression (5.15).

In the incompressible case, the gradient can be forgotten up to a redefinition of pressure and Z 1 = |∇φ| thanks to (B.1). Therefore, we find the previous result (5.2). We have that

div E i (Z i )Z i C i 1 ε ζ φ ε = div(E i (Z i )Z i C i ) 1 ε ζ φ ε + E i (Z i )Z i C i n |∇φ| ε 2 ζ φ ε .
The last term vanishes because C i n = 0. With the identity div(Av) = div(A T ) • v + A T : ∇v and the symmetry of C i , we get

div(E i (Z i )Z i C i n) = 0 = div(E i (Z i )Z i C i ) • n + E i (Z i )Z i C i : [∇n].
With the definition of C 1 in (5.5) and with (5.6), we find the normal part of forces F i (without the pressure term):

F 1 • n = -Tr([∇n])E 1 (Z 1 )Z 1 1 ε ζ φ ε . (5.16)
With the relation 2Z 1 Z 2 = Tr(A) (see (4.6)), we get

F 2 • n = -(2Tr(A[∇n]) -Tr(A)Tr([∇n])) E 2 (Z 2 ) 2Z 1 1 ε ζ φ ε .
(5.17)

In the previous formulas we used the relation [

I -n ⊗ n] : [∇n] = Tr([∇n]) (where [∇n] T n = 0 since |n| 2 = 1
) which is by definition the surface mean curvature. Hence the normal force associated to the area variation (5.16) only depends on the geometry through the surface curvature. The normal force associated to the shear variation (5.17) not only depends on the geometry but also on the deformations on the surface through A. In the general case, the normal part of the shear force does not vanish.

Numerical method.

The fluid-structure equations (5.3) and (5.4) are discretized with finite difference schemes on a staggered grid (see Figure 5.1 for a 2D configuration) with the classical projection method. Let Δt be the time step and u n ,p n ,φ n ,Y n the time discretization of the variables at t n = nΔt. The semi-discretization in time is given by

• • • • • • • • • p,φ,Y × × u 1 × × u 2 Δx 2 Δx 1
Step 1:

u -u n Δt + u n • ∇u n - 1 R e Δu n = F 1 (φ n ,Y n ) + F 2 (φ n ,Y n ) Step 2: Δp n+1 = div(u ) Δt Step 3: u n+1 = u -Δt∇p n+1
Step 4:

φ n+1 -φ n Δt + u n+1 • ∇φ n = 0 Y n+1 -Y n Δt + u n+1 • ∇Y n = 0
In Step 1, a prediction of the velocity is computed with an explicit Euler scheme in time and classical central staggered velocity-pressure schemes of order two for the convection, the diffusion, and the source terms. An explicit scheme in time is used for the diffusion because the Reynolds number is large enough. In Step 2, the Poisson equation for the pressure is performed with the Fishpack library [START_REF] Adams | Fishpack -efficient FORTRAN subprograms for the solution of separable elliptic partial differential equations[END_REF] (we use homogeneous Neumann boundary conditions). In Step 3, the the velocity is corrected to enforce the incompressibility condition. In Step 4, the transport equations are discretiezd with an explicit Euler scheme in time and a WENO5 scheme in space [START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF].

We do not perform the redistancing. Instead, we use the renormalization φ |∇φ| to measure the distance to interface. Thus,

|∇φ| 1 ε ζ φ ε is replaced by 1 ε ζ φ |∇φ|ε .
This approach was proved in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF] to be efficient from the point of view of both volume conservation and interface force calculations. For the cut-off function, we considered the expression ζ(r) = 1 2 (1 + cos(πr)) on [-1,1] and ζ(r) = 0 elsewhere. We use in our simulations the linear elastic laws

E 1 (r) = λ(r -1) E 2 (r) = μ(r -1), (5.18)
where λ and μ are the elastic stretching and shear modulus (however the model is still strongly nonlinear due to the geometric nonlinearities and the coupling with Navier-Stokes equations).

Numerical tests.

In this section, we provide evidences of numerical convergence of our method on the illustrative test case of a sheared elastic sphere. The domain Q = [-1,1] 3 is discretized on a Cartesian mesh with 128 points in each direction. We choose in our simulation a Reynolds number R e = 100, an elastic stretching modulus λ = 1, and a shear modulus μ = 0.1. The parameter ε is equal to 3.5Δx in the simulations where Δx is the grid size. We take the time step Δt = 1.310 -3 s. We impose zero velocity for the initial and boundary conditions. The immersed initial surface is a sphere, so that

φ 0 (x,y,z) = x 2 + y 2 + z 2 -0.5, (5.19)
but which is not in its rest state, since material points have been moved, so that Y (x,y,z,0) = (xcos(t 0 z) + y sin(t 0 z),-xsin(t 0 z) + y cos(t 0 z),z). (5.20) This corresponds to a deformation of the sphere when a 3D circular shear (see the expression of Y for TC4 in Table A.1 and Figure A.8) is applied on opposite directions on the north and south poles (with respect to the z-axis) until t = t 0 . Here, we take t 0 = π. While this initial deformation was (artificially) imposed with no area variation (actually, the sheared surface is still geometrically a sphere), the area will change when it starts to relax, so the force F 1 is also involved. The motion is however driven by F 2 .

The numerical results at different time steps are presented in Figure 5.2. For plotting purposes, we represented on the deformed surface a checker board pattern which was tracked with markers to visualize the shear relaxation. In Figure 5.3, we plot the velocity magnitude at t = 0.5s while the relaxation takes place. Due to the large imposed shear, the surface undergoes a complex deformation involving some small folds. These kinds of ripples were also observed in [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF] in the simulation of a capsule in a simple shear flow. A nice feature of our numerical code is its ability to relax toward a stable solution without any curvature energy. Some grid effects are, however, present for large time (see the last image of Figure 5.2). The use of an unstructured mesh and finite-element solver to remove this effect is presently under study. To provide more quantitative results, we plotted in Figure 5.4 and Figure 5.5 pressure slices following the x-and z-axes during this relaxation for t = 0.1 and t = 1.2. The computations were performed with N points in each direction for three resolutions with N = 64,128,256. Accordingly, we chose ε = 1.75Δx,3.5Δx,7Δx, where Δx is the corresponding grid size, so that the numerical interface width remains constant. We observe numerical convergence. Likewise, Figure 5.6 depicts the variation of vertical radius through time.

As the flow is incompressible, the map Y (•,t) : Ω → Ω conserves volume; that is, det∇Y = 1 in the continuous case. However, after time and space discretization, and due to numerical errors introduce while solving for the transport equations on Y , such a nonlinear combination of derivatives will not fully respect this constraint. On Figure 5.7 and Figure 5.8, we depict the RMS (Root Mean Square) of det∇Y -1 as a function of time, both for the whole domain and on the interface. Thus, the plotted quantities are, respectively, Fig. 5.8: footnotesizeTime evolution, for N = 64,128,256 of t → RM S Γt (t) defined in (5.21).

RM S Ω (t) = 1 |Ω| Ω |det∇Y (x,t) -1| 2 dx 1 2
and RM S Γt (t) = 1 |Γ t | Γt |det∇Y (x,t) -1| 2 ds 1 2
.

(5.21)

While N = 64 is clearly under-resolved, the figures show numerical convergence of those RMS as N increases.

Conclusion

In this paper, we introduced a way to account for membrane shear in a completely Eulerian framework. The idea was to build two invariants from algebraic combinations of the space derivatives of the backward characteristics of the continuous medium velocity field. The first invariant accounts for changes of local area of the membrane, while the second records the local membrane shear. We showed how to compute those invariants in an Eulerian way and produced some evidences that the second invariant indeed captures the membrane shear. The fluid-structure coupling is therefore rephrased as a Navier-Stokes equation with source term involving the two invariants, coupled to a vectorial transport equation for the backward characteristics. The formulation is therefore very versatile since it could be implemented as an add-on to a preexisting fluid solver. In this first work, we presented outputs of a 3D numerical code using finite differences and a projection method on a staggered mesh, which simulate the relaxation of a spherical membrane which has been sheared from its poles. Our model and code are able to capture ripples due to the bucking of membrane in a stable way without any extra smoothing curvature energy. Future work will focus on the application of this method to model red blood cells and capsules in 3D flows.

Appendix A. Analytical illustrations for Z 2 .

We present some analytical illustrations to show that Z 2 is intuitively a good measure of the local shear variation of a surface. In all the test cases, we define an initial surface Γ 0 = {φ 0 = 0} and a velocity field u which will move material points of this surface. Then we compute the backward characteristics Y and the invariants Z 1 and Z 2 to see how these Eulerian quantities could record this motion. The test cases and the results are described in Table A.1. A.1. Deformations where Γ t = Γ 0 .

TC φ 0 (x,y,z) u Y (x,y,z,t) Z 1 Z 2 1 z ⎛ ⎝ x αy 0 ⎞ ⎠ ⎛ ⎝ xe -t ye -αt z ⎞ ⎠ e t(1+α) ch(t(1 -α)) 2 z ⎛ ⎝ βy x 0 ⎞ ⎠ ⎛ ⎝ xch(t √ β) -y √ βsh(t √ β) -x √ β sh(t √ β) + ych(t √ β) z ⎞ ⎠ 1 1 + (1+β) 2 2β sh 2 (t √ β) 3 x 2 + y 2 -1 ⎛ ⎝ -yz xz 0 ⎞ ⎠ ⎛ ⎝ xcos(tz) + ysin(tz) -xsin(tz) + ycos(tz) z ⎞ ⎠ 1 1 + t 2 2 (x 2 + y 2 ) 4 x 2 + y 2 + z 2 -1 ⎛ ⎝ -yz xz 0 ⎞ ⎠ ⎛ ⎝ xcos(tz) + ysin(tz) -xsin(tz) + ycos(tz) z ⎞ ⎠ 1 1 + t 2 (x 2 +y 2 ) 2 2(x 2 +y 2 +z 2 ) 5 x 2 + y 2 + z 2 -1 ⎛ ⎝ x y z ⎞ ⎠ ⎛ ⎝ e -t x e -t y e -t z ⎞ ⎠ e 2t 1 6 max(|x|,|y|,|z|) -1 ⎛ ⎝ 0 x 0 ⎞ ⎠ ⎛ ⎝ x y -tx z ⎞ ⎠ See (A.1)-(A.3) 7 x 2 + y 2 + z 2 -1 ⎛ ⎝ 0 x 0 ⎞ ⎠ ⎛ ⎝ x y -tx z ⎞ ⎠ See (A.4) See (A.5)
In this first series of test cases (TC1 to TC4), while material points could have moved, globally the initial and deformed surfaces are the same (Γ 0 = Γ t ). Moreover, the deformations are uniform in space in the sense that the invariants Z i on the surface do not depend on the spatial variables (except for TC4).

What we refer to as the α and β deformations (TC1 and TC2) are 2D, and the initial surface is the plane z = 0. The velocity field of each deformation is plotted in the figures below with the corresponding values of Z 1 and Z 2 . The β = -1 deformation is a rotation, and, as expected (slipping from hyperbolic to circular functions with √ β identified to i ∈ C), there is no area variation and shear variation (see The α = 1 deformation is a pure dilatation, and, as expected, there is only area variation (see Figure A.2). The β = 0, β = 1, and α = -1 deformations correspond to different shear transformations, and, as expected, there is only shear variation (see figures A.3, A.4, and A.5). Note that, for the β = 0 test case, we take the limit β -→ 0 the results for Y and Z i . Moreover, for the β = 0 deformation, we found the same invariants as in the Lagrangian case (see (3.3)). The α = 0 deformation is an uniaxial deformation and there is area and shear variation (see In each plane {z = α}, the velocity is a rotation of magnitude α. In these test cases, there is no area variation but pure shear variation. For TC3, Z 2 = 1+ t 2 2 is constant on the surface (x 2 + y 2 = 1 on the cylinder). This test case is clearly a 3D generalization on a cylinder of the 2D β = 0 deformation, and this is why we found the same invariants. For TC4, Z 2 = 1+ t 2 2 (1z 2 ) 2 on the surface and depends only on the height z which is intuitive (x 2 + y 2 + z 2 = 1 on the sphere). A.2. Deformations where Γ t = Γ 0 .

We present now three test cases (TC5-7) where the deformed surface Γ t is different from the initial shape Γ 0 . For the 3D dilatation deformation (TC5 see Figure A.9), the initial surface is a sphere. In this test case we have only area variation, as expected.

In the last test cases, the same shear velocity field (see Figure A.10) is applied on a cube and on a sphere. For the test case TC6, the deformed shape of the cube (see We have the relation det(B) = J 2 , so the expression of Z 1 (B.2) reduces to (B.1). This recovers the result of [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF], where the area variation captured by |∇φ| in the incompressible case (J = 1).

B.3.

Proof that Z 1 is the local area variation. Proposition B.3. Let Γ 0 a smooth surface, deformed with smooth forward characteristics in Γ t = X(Γ 0 ,t). We denote by u the associated smooth velocity field. Let a parametrization of the surface Γ t be given by (θ 1 ,θ 2 ) → Ψ(θ 1 ,θ 2 ,t), where Ψ : This result holds for all function g. Therefore, (B.3) holds, and Z 1 measures the local area variation (in the compressible and incompressible cases).

U × R + -→ R 3 is

Fig. 2 . 1 :

 21 Fig. 2.1: Initial and deformed configurations.

Fig. 3 . 1 :

 31 Fig. 3.1: Orthogonal (left) and non orthogonal (right) parametrization of the plane.

Fig. 3 . 2 :

 32 Fig. 3.2: Orthogonal parametrization of the sphere.

4 .

 4 Surface deformation in the Eulerian frame 4.1. The surface tensor A.

Fig. 5 . 1 :

 51 Fig. 5.1: Staggered grid with position of unknowns.

Fig. 5 . 2 :

 52 Fig. 5.2: Numerical simulation of the relaxation of a sphere due to shear force at times t = 0, 0.5s, 1s, 1.5s, 2s, 2.5s, 3s, 3.5, 4s, 4.5s, 5s, 9s. Lagrangian markers are used for plotting a checkered pattern coming back to an undeformed state when shear decreases.

Fig. 5 . 3 :

 53 Fig. 5.3: Velocity magnitude at t = 0.5s.

Fig. 5 . 4 :

 54 Fig. 5.4: Pressure along the x-axis at time t = 0.1 (Left); and t = 1.2 (Right) for N = 64,128,256.

Fig. 5 . 5 :

 55 Fig. 5.5: Pressure along the z-axis at time t = 0.1 (Left) and t = 1.2 (Right) for N = 64,128,256.

Fig. 5 . 6 :

 56 Fig. 5.6: Vertical radius up to t = 30 (Left); zoomed to t ∈ [0,5] for N = 64,128,256.(Right).

Fig. 5 . 7 :

 57 Fig. 5.7: Time evolution for N = 64,128,256 of t → RM S Ω (t) defined in (5.21).

  Figure A.1).

1 Fig. A. 1 :

 11 Fig. A.1: Velocity field (-y,x,0) of the β = -1 deformation (TC2).

5 Z1 = e 2t Z2 = 1 Fig. A. 2 :

 512 Fig. A.2: Velocity field (x,y,0) of the α = 1 deformation (TC1).

t 2 2 Fig. A. 3 :

 23 Fig. A.3: Velocity field (0,x,0) of the β = 0 deformation (TC2).

Fig. A. 4 :

 4 Fig. A.4: Velocity field (y,x,0) of the β = 1 deformation (TC2).

Figure A. 6 )

 6 . It can be quite surprising at first sight, but, when a surface is stretched in two directions with different magnitude (α = 1), we do have shear and therefore as expected Z 2 = ch(t(1α)) = 1. In the "3D circular shear" test cases the initial surfaces are a cylinder (TC3 see Figure A.7) and a sphere (TC4 see Figure A.8).

Fig. A. 5 :

 5 Fig. A.5: Velocity field (x,-y,0) of the α = -1 deformation (TC1).

Fig. A. 6 : 1 Z2 = 1+ t 2 2 (x 2 + y 2 )Fig. A. 7 : 2 2(x 2 +y 2 +z 2 )Fig. A. 8 :

 612227228 Fig. A.6: Velocity field (x,0,0) of the α = 0 deformation (TC1).

Figure A. 11 )A 2 .B. 2 . 2 and 2 . 2 =

 1122222 is given by the zero level set of φ(x,y,z,t) = max(|x|,|y -tx|,|z|) -1.The calculations are done independently on each plane that composes the cube. We have the following results:On the faces {x = ±1} Z 1 = 1, Z 2 = 1 (A.1)On the faces {ytx= ±1} Z 1 = 1+ t 2 Z 2 = 2+ t 2 2 √ 1+ t 2 (A.2)we have A i3 = 0. Since A is a symmetric tensor and B is an orthonormal basis, the matrix A is symmetric, and its structure is given by )2 + (A 12 ) 2 A 11 A 12 + A 12 A 22 0 A 11 A 12 + A 12 A 22 (A 12 ) 2 + (A 22 )In the basis B , we have In ⊗ n = and A have the same invariants soTr(A) = Tr(A ) = A 11 + A 22 Tr(Cof(A)) = Tr(Cof(A )) = A 11 A 22 -(A 12 )It is now easy to check that A Tr(A )-(A ) 2 = Tr(Cof(A ))[In ⊗ n],and the lemma is proved. Equivalent expression for Z 1 . In this subsection, we give an equivalent expression of Z 1 which links the present formulation with our previous work[START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF].Proposition B.2. Let J = det(∇ ξ X) = det(∇Y ) -1 be the volume ratio. ThenZ 1 = J |∇φ| |∇φ 0 (Y )| . (B.1)Proof. The invariant Z 1 is defined by (4.5)Z 1 = Tr(Cof(A)). (B.2)With the definition (4.3), we getTr(A) 2 = Tr(B) -(B 2 n) • n (Bn) • n Tr(A 2 ) = Tr(B 2 ) -2 (B 3 n) • n (Bn) • n + (B 2 n) • n (Bn) • nUsing the Cayley-Hamilton theorem B 3 -Tr(B)B 2 + Tr(Cof(B))Bdet(B)I = 0, we getTr(Cof(A)) = 1 (Bn) • n (B 3 n) • n -Tr(B)(B 2 n) • n + Tr(Cof(B))(Bn) • n = det(B) (Bn) • n .The gradient of (2.5) gives the relation ∇φ = [∇Y ] T ∇φ 0 (Y ), so (Bn) • n = [∇Y ] -T ∇φ |∇φ| |∇φ 0 (Y )| 2 |∇φ| 2 .

Table A . 1 :

 A1 Table of test cases.

  smooth on an open set U of R 2 . The local area variation verifies|Ψ ,1 (θ 1 ,θ 2 ,t) ∧ Ψ ,2 (θ 1 ,θ 2 ,t)| |Ψ ,1 (θ 1 ,θ 2 ,0) ∧ Ψ ,2 (θ 1 ,θ 2 ,0)| = Z 1 (Ψ(θ 1 ,θ 2 ,t),t) Z 1 (Ψ(θ 1 ,θ 2 ,0),0) = Z 1 (x,t) Z 1 (Y (x,t),0) (B.3) for x = Ψ(θ 1 ,θ 2 ,t) and with the notationΨ ,i (θ 1 ,θ 2 ,t) := ∂ θi Ψ(θ 1 ,θ 2 ,t). Let f : R 3 × R + -→ R be a smooth function. The Reynolds formula for surfaces reads Since, Ψ(θ 1 ,θ 2 ,0) = Y (Ψ(θ 1 ,θ 2 ,t),t) the equation (B.5) becomes ,θ 2 ,0)) Z 1 (Ψ(θ 1 ,θ 2 ,t),t) |Ψ ,1 (θ 1 ,θ 2 ,t) ∧ Ψ ,2 (θ 1 ,θ 2 ,t)| dθ 1 dθ 2 (Ψ(θ 1 ,θ 2 ,0),0) |Ψ ,1 (θ 1 ,θ 2 ,0) ∧ Ψ ,2 (θ 1 ,θ 2 ,0)| dθ 1 dθ 2 .

	Proof.
			d dt Γt	f ds =
					the
	previous integral becomes
	1 Z 1	(∂ (B.4)
	The first term vanishes with (2.2). The second term is also zero with (5.5). Therefore,
				d dt Γt	g(Y (x,t)) Z 1 (x,t)	ds = 0.	(B.5)
	=	U	g(Ψ(θ 1 ,θ 2 ,0)) Z 1

Γt ∂ t f + u • ∇f + f [∇u] : [In ⊗ n] ds. Let g : R 3 -→ R be smooth function and let f (x,t) = g(Y (x,t)) Z1(x,t) . The expression in t (g(Y )) + u • ∇(g(Y ))) -g(Y ) (Z 1 ) 2 (∂ t Z 1 + u • ∇Z 1 -Z 1 [∇u] : [In ⊗ n]). U g(Ψ(θ 1

On the faces {z = ±1}

On the faces corresponding to {x = ±1}, there is no area and shear variation (the faces are just translated). For the faces of equations {ytx = ±1}, the faces are stretched in only one direction and there is area and shear variation as in the 2D α = 0 deformation. For the faces {z = ±1},there is only shear variation as in the 2D β = 0 deformation. For the test case TC7, the deformed shape of the sphere Γ t is given by the zero level set of

We have that

In Figure A.12, iso-contours of Z 1 and Z 2 are plotted on the deformed surface Γ t . The results are intuitive: the area variation reaches its maximum on points that are orthogonal to the y = tx plane while the shear variation is larger along the z-axis. The deformation and the initial surface are the same as in the Lagrangian test case in Section 3.3. In the formula (A.5), there are no discontinuities at the poles (x 2 + (ytx) 2 + z 2 = 1 on the surface Γ t ), unlike in the Lagrangian formula (3.6).

Appendix B.

B.1. Technical lemma.

Lemma B.1. The following identity holds:

Proof.

It is equivalent to prove the lemma in any basis of R 3 . We consider the orthonormal basis of R 3 B = (e 1 ,e 2 ,e 3 ) = (τ 1 ,τ 2 ,n), where (τ 1 ,τ 2 ) is an orthonormal basis of the tangent plane orthogonal to n. Let A be the matrix of the tensor A in the basis B , and we denote by A ij its coefficients. We have with (4.4) that An = 0; hence