Shuang Tan

Justine Jean

Valerie Pralong

Simona Moldovan

Hailing Guo

Svetlana Mintova
email: mintova@ensicaen.fr.

Effect of particle size on the sodium ions utilization efficiency of zeolite-templated carbon as the anode in sodium ion battery

Keywords: Zeolite-templated carbon, nanosized zeolite, high sodium storage, carbon anode material

A series of zeolite-templated carbon (ZTC) with the particl z f 1μ , 70 , 10 h been prepared. The external surface area and diffusion path of Na + inside the ZTC was varied by controlling the size of the zeolite crystals used as templates. The three ZTC samples have been used to store sodium cations (Na + ) for the purpose of sodium-ion batteries (SIB). It was found that the ZTC with a particle size of 10 nm exhibited the highest first discharge capacity of 2116mAh•g -1 at a current density of 0.1mA•g -1 due to the highest external surface area and more complete filling of the internal pore volume of the ZTC. While the ZTC with 70nm particle size showed better reversible capacity of 180mAh•g -1 at a current density of 0.1mA•g -1 , and the 10nm size ZTC caused more electrolyte decomposition reaction and more solid electrolyte interphase (SEI) layer on the external surface of the crystals, resulting in a low first coulombic efficiency and poor Na + extraction.

Introduction

Due to the abundance of sodium in the earth's crust and the similarity of their energy storage systems to Lithium-ion batteries (LIBs), sodium-ion batteries (SIBs) have drawn a lot of attention as the next-generation energy storage technology [START_REF] Delmas | Sodium and Sodium-Ion Batteries: 50 Years of Research[END_REF] . However, owing to the large ionic diameter of Na + (1.02Å) compare to Li + (0.76Å), the massive volume expansion for accommodating Na + has posed serious challenges to the stability of the anode material [START_REF] Hwang | Sodium-ion batteries: present and future[END_REF] . Additionally, the electrochemical kinetics of SIBs are also slowed by the large ion diameter of Na + , which results in poor cycle performance and low capacity. Thus shortening the sodium ion diffusion path to increase its storage rate, along with constructing a stable and sufficient space to accommodate the sodium ion expansion are the main directions for designing high-performance SIB anode materials [START_REF] Slater | Sodium-Ion Batteries[END_REF] .

Due to their electrical conductivity, chemical stability, accessibility, and sustainability, carbon materials have been exploited as the anode in SIBs. Graphite, the most stable allotrope of carbon, is capable to provide a theoretical capacity of 372 mAh•g -1 for lithium ions due to the large layer spacing of 3.35 Å. However, the same space can only provide a theoretical capacity of 35 mAh•g -1 for sodium ions [START_REF] Wen | Expanded graphite as superior anode for sodium-ion batteries[END_REF] . Also, the rate performance is low caused by the larger ionic size and slower kinetics of sodium ions [START_REF] Wang | Nano-Structured Phosphorus Composite as High-Capacity Anode Materials for Lithium Batteries[END_REF] . In order to maximize their Na + storage capability and efficiency, researchers have looked into the storage of sodium in carbon materials with various pore configurations to facilitate the accessibility of Na + to void structures such as core-shell structure [START_REF] Sun | A Passionfruit-Like Carbon-Confined Cu2ZnSnS4 Anode for Ultralong-Life Sodium Storage[END_REF] , nanosheets [START_REF] Xia | Hard Carbon Nanosheets with Uniform Ultramicropores and Accessible Functional Groups Showing High Realistic Capacity and Superior Rate Performance for Sodium-Ion Storage[END_REF] , microtube [START_REF] Li | Hard Carbon Microtubes Made from Renewable Cotton as High-Performance Anode Material for Sodium-Ion Batteries[END_REF] , nanofiber [START_REF] Ren | High-performance sodium-ion storage: multi-channel carbon nanofiber freestanding anode contrived via ingenious solvent-induced phase separation[END_REF] , and porous carbon framework [START_REF] Hou | Carbon Quantum Dots and Their Derivative 3D Porous Carbon Frameworks for Sodium-Ion Batteries with Ultralong Cycle Life[END_REF][START_REF] Wang | A nitrogen-doped three-dimensional carbon framework for high performance sodium ion batteries[END_REF] . These examples demonstrate that constructing cavities to accommodate the volume expansion of sodium ions and increasing the external surface area of the material to facilitate the contact of sodium ions with the carbon material is a better way to store Na + . Recently, the storage behavior of Na + in nanoporous attracted considerable interest and it is believed that this could result in hard carbon anodes with prolonged low-potential plateaus and high capacity [START_REF] Li | Regulating Pore Structure of Hierarchical Porous Waste Cork-Derived Hard Carbon Anode for Enhanced Na Storage Performance[END_REF][START_REF] Matei Ghimbeu | Insights on the Na+ ion storage mechanism in hard carbon: Discrimination between the porosity, surface functional groups and defects[END_REF] .

Zeolite templated carbon (ZTC) is obtained by calcining carbon-containing precursors in the zeolite channels, the carbon replica of the zeolite channel structure is obtained under high temperature treatment conditions. It has been reported that the abundant and uniform porous of ZTC are crucial for storing small molecules such as CO 2 14 , H 2 15 , even Li +16 and electron [START_REF] Mostazo-López | Ultraporous nitrogen-doped zeolite-templated carbon for high power density aqueous-based supercapacitors[END_REF] . In our previous work [START_REF] Tan | Ordered sodium zeolite-templated carbon with high first discharge capacity for sodium battery application[END_REF] , the ZTC was used as anode materials of SIB and offered a first discharge capacity of up to 1714 mAh g -1 at a current density of 0.1 A•g -1 . The results suggested that sodium ions can be successfully stored inside the ZTC. Unfortunately, some Na + were trapped in the ZTC and unable to be released, leading to a significant capacity loss.

To investigate the storage behavior of Na + in ZTC and to further enhance the utilization of Na + , we reduced the size of the ZTC in order to facilitate the diffusion of Na + into the interior space of ZTC. Three ZTC samples with diameters of 1µm, 70 nm, and 10 nm were synthesized and evaluated as SiB anodes at various current densities. The contribution of different place in ZTC to Na + storage capacity was evaluated, and the impact of their particle size on the capacity was assessed. The results demonstrated that the nanoporous played a major role in the initial storage of Na + , and by reducing the particle size of the ZTC, a higher amount of Na + was stored but this also caused a significant capacity loss by enhancing the formation of the SEI layer. Experimental 1.1 Preparation of FAU zeolite crystals with particle size of 1 μm, 70 nm, and 10 nm.

The FAU zeolites with different particle sizes were prepared following the procedure published by Mintova [START_REF] Li | Regulating Pore Structure of Hierarchical Porous Waste Cork-Derived Hard Carbon Anode for Enhanced Na Storage Performance[END_REF] . In a typical synthesis, a clear suspension with the molar ratio of 8 Na 2 O: 0. O; the suspension was aged at ambient conditions for 24h and then subjected to crystallization at 100°C for 12h. The powder obtained after crystallization was labeled as X-FA (X=1μm, 70nm, and 10nm corresponding to the particle size of FAU zeolite crystals). The crystalline suspensions were filtered to collect the FAU zeolite crystals and then washed with deionized water and dried before use.

1.2 Preparation of FAU-ZTC with particle sizes of 1 μm, 70 nm, and 10 nm.

A vapor impregnation method described elsewhere [START_REF] Tan | Ordered sodium zeolite-templated carbon with high first discharge capacity for sodium battery application[END_REF] was used to prepare the FAU-ZTC with various particle sizes. In a typical synthesis, 5 g zeolite (the template) was placed in a vertical glass tube lined with quartz wool, and a stream of furfuryl alcohol (FA) vapor was produced under vacuum and flowed through the zeolite layer. After 1h, the zeolite template was saturated with the FA, and then the atmospheric pressure was restored and the zeolite was left to continue absorbing FA at ambient temperature for 24 hours. The zeolite samples loaded with FA were then washed with mesitylene to remove the excess of FA adsorbed on the external surface. The excess mesitylene was then removed by heating at 80°C for 12 hours. Then the FA-loaded FAU zeolite samples were heated at 150°C for 8 hours in nitrogen to pre-polymerize the FA. The pre-polymerized samples were then heated at 900°C for 4h with a heating rate of 5°C/min under a nitrogen atmosphere. Finally, the zeolite template was removed from the ZTC products by acid etching with a solution of 0.4M HF/0.2M HCl. Finally, the ZTC samples were washed with distilled water and dried at 100°C overnight. The obtained samples were labeled as X-FAU-ZTC (X=1μ , 70 , 10 ).

Characterization

The morphology and crystal size of the samples were characterized using a JEOL SEM-7900F scanning electron microscope (SEM) with an acceleration voltage of 5 kV. The Scanning Transmission Electron Microscopy (STEM) experiments were performed on a double corrected Analytical JEOL ARM200CF TEM operated at 200kV for imaging and 120kV for mapping. The TEM is equipped with a 100mm 2 Centurio Energy Dispersive X Rays (EDX) detector for the chemical analysis of the constituents. The samples were prepared under ultrasound treatment of powders in ethanol obtained by mechanical scratching of the cathode materials after the different treatments. The solution was drop cast on a holey Carbon membrane deposited on a 300 mesh Au grid adapted for the TEM experiments. Elemental maps (256 x 256px) were acquired using the Analysis Station (Jeol) software with a scanning speed of 0.02 ms/px in order to minimize the damaging effect of the electron beam and a spatial tracking algorithm was applied every 60 seconds. The same software was employed for the qualitative and quantitative analysis of the experimental data. X-ray powder diffraction (XRD) patterns were collected using a Bruker D8 Advance. The N 2 adsorption-desorption experiments were performed at 77 K using a Quantachrome Autosorb iQ3 gas-sorption analyzer. The micropore volume was calculated from the v-t plot. The pore size distribution was determined using the Quenched solid density functional theory (QSDFT) method.

Electrochemical test

The electrochemical performance of the ZTC samples was evaluated using half-cells assembled in a glove box. The powder samples were mixed with super P and poly vinylidene fluoride in a weight ratio of 8:1:1 in N-methyl-2-pyrrolidinone to form a homogeneous slurry, which were then deposited on a Cu foil and dried in a vacuum oven at 80°C for 12 h. Sodium metal foils were used as counter and reference electrodes. A solution of 1 mol•L -1 NaClO 4 in ethylene carbonate/diethyl carbonate (1/1; v/v) with 5 wt.% fluoroethylene carbonate was used as the electrolyte and Celgard 2400 was used as the separator. The active materials mass loaded on each round of Cu foil was about 0.7 mg•cm -2 . Cyclic voltammetric measurements with a scan rate of 0.1 mV•s -1 were performed on a CHI760D electrochemical workstation. Galvanostatic charge/discharge measurements were performed using a multichannel battery testing system (LAND CT2001A) in a voltage range of 0.01-3 V for anodes versus Na/Na + at different constant current densities.

The amount of Na + accommodated in supercages of the ZTC was calculated using the following equations:

(1) (2) where C is Capacity, F is Faraday constant, N a is Avogadro's constant, V micro is Micropores volume of ZTC, V pore is single supercage volume of ZTC, and R is pore diameter of ZTC. The structure of FAU-type zeolite (Figure S1A), and ZTC (Figure 1A) were determined by X-ray diffraction. All three samples of zeolites with various particle sizes exhibit characteristic Bragg peaks corresponding to the FAU-framework structure. As the particle size decreases, the intensity of the peaks decreased and the half-height width of the peak expands, indicating that the structure is less ordered. The non-negligible decrease in crystallinity led to a decrease in the specific surface area and micropore volume of the zeolite templates as confirmed by nitrogen physisorption analysis (Figure S2A, Table S1). In comparison to sample 1µm-FAU, the specific surface areas of 70nm-FAU and 10 nm-FAU decreased by 26% and 50%, respectively. Similar to other ZTC materials synthesized with furfuryl alcohol [START_REF] Zhang | Boosting Zn-Ion Energy Storage Capability of Hierarchically Porous Carbon by Promoting Chemical Adsorption[END_REF] , h pl 1μ -FAU-ZTC exhibits a weak peak at 6.3° 2θ h XRD p , wh ch pl c f 111 c y l pl f h FA z l [START_REF] Antoniou | Methane storage in zeolite-like carbon materials[END_REF][START_REF] Kwon | N-doped zeolite-templated carbon as a metal-free electrocatalyst for oxygen reduction[END_REF] . This c h h 1μ -FAU-ZTC sample has faithfully replicated the FAU zeolite structure. Only weak peaks or no peaks at the same position in the XRD patterns of samples 70 nm-FAU-ZTC and 10 nm-FAU-ZTC were observed as a consequence of the small particle size and partial collapse of nanosized zeolite template at high temperature [START_REF] Parmentier | Effect of the zeolite crystal size on the structure and properties of carbon replicas made by a nanocasting process[END_REF] . At the same time, no broad peaks representing the disordered graphitic structure of the amorphous c w 25° 2θ f ll h samples [START_REF] Shi | Rapidly reversible adsorption of methane with a high storage capacity on the zeolite templated carbons with glucose as carbon precursors[END_REF] , suggesting that the ZTC samples have formed a self-standing porous carbon skeleton by replicating the zeolite framework rather than amorphous stacked carbon.

Results and Discussion

As shown in Figure 1B, all three ZTC samples exhibited a type I isotherm [START_REF] Armandi | Study of hydrogen physisorption on nanoporous carbon materials of different origin[END_REF] which is composed of micropores, proving that the parent zeolite's microporous structure has been successfully replicated (Figure 1B). The adsorption-p h f pl 1μ -FAU-ZTC contains a "low-pressure hysteresis". This phenomenon has been reported for various ZTC materials [START_REF] Choi | Large-scale synthesis of high-quality zeolite-templated carbons without depositing external carbon layers[END_REF][START_REF] Lu | Zeolite-templated nanoporous carbon for high-performance supercapacitors[END_REF][START_REF] Liu | Enhancing inward and outward mass transport during chemical vapor deposition of pyrolytic carbon for better synthesis of ZTC[END_REF] , which has been attributed to carbon particles swelling during the nitrogen adsorption measurements, or to the strong interactions of nitrogen molecules with the carbon surface [START_REF] Su | Synthesis and characterization of microporous carbons templated by ammonium-form zeolite Y[END_REF][START_REF] Ma | Synthesis methods for preparing microporous carbons with a structural regularity of zeolite Y[END_REF] . The isotherms of samples 70nm-FAU-ZTC and 10nm-FAU-ZTC differ slightly from that of sample 1μ -FAU-ZTC. An upward turn with a weak hysteresis at high relative pressure was recorded. This hysteresis loop indicates the mesoporosity generated from interparticle space due to the packing of well-defined individual particles with similar sizes [START_REF] Parmentier | Effect of the zeolite crystal size on the structure and properties of carbon replicas made by a nanocasting process[END_REF] . The external surface area of the ZTC increases as the crystal size decreases, but the total surface area and the microporous surface area (S mic ) decrease substantially (Table 1). Apparently, the decrease of S mic of the parent FAU zeolite template leads to the decrease of S mic in the carbon replica. Sample a BET specific surface area was calculated from the adsorption branch in P/P0 = 0.02-0.15.

S BET (m 2 •g -1 ) a S mi (m 2 •g -1 ) b S ex (m 2 •g -1 ) V t (cm 3 •g -1 ) c V mi (cm
b Micropore surface was calculated by the t-plot method.

c Total pore volume was determined at P/P0 = 0.97.

d Micropore volume was calculated from the t-plot method.

Figure 1D shows the RAMAN spectra of the carbon samples, which exhibit D and G peaks at 1380cm -1 and 1600cm -1 corresponding to the A 1g ring respiration mode and the E 2g vibration between adjacent sp 2 hybridized carbon atoms, respectively [START_REF] Zhao | Carbon nanotube formation using zeolite template and applications[END_REF] . The integrity of the carbon structure network, or the level of graphitization, is indicated by the intensity ratio of D and G peaks. The I D /I G ratio of 1µm-FAU-ZTC, 70nm-FAU-ZTC, and 10nm-FAU-ZTC is 0.64, 0.75, and 0.76, respectively. The 1μm-FAU-ZTC has lower I D /I G ratio as the template contains larger zeolite particles that are more robust at high temperature and can provide a stable environment for the growth of the carbon skeleton [START_REF] Nishihara | A possible buckybowl-like structure of zeolite templated carbon[END_REF][START_REF] Cui | Ultrafine Pd nanoparticles supported on zeolite-templated mesocellular graphene network via framework aluminum mediation: An advanced oxygen reduction electrocatalyst[END_REF] . In contrast, samples 70nm-FAU-ZTC and 10nm-FAU-ZTC have higher I D /I G ratio due to the possible collapse of the small zeolite particles at high temperature and presence of defects in the carbon lattice. SEM images of the ZTC samples reveal the progressive decrease in the particle size of the zeolite pl h ZTC (F gu 2). B h 1μ -FA 1μ -FAU-ZTC samples exhibit the classical ortho-octahedral shape, which suggests that the carbon replicates the zeolite morphology cl ly. Th g w h pl 1μ -FAU is possible causing the interconnections in the carbon replica. The 70nm-FAU sample contains individual crystals with a size of 70-150 nm, and less sharp edges along with a change in crystal structure from ortho-octahedral to ellipsoidal. This phenomenon is also observed for sample 70 nm-FAU-ZTC due to the exact replication of the crystal structure of the template. The crystals diameter of sample 10nm-FAU is in the range of 10-40 nm, and these crystals transformed into nearly spherical blocks. Due to the high surface energy and the capillary effect of water during the drying, the small size carbon replicas are easily agglomerated [START_REF] Li | Regulating Pore Structure of Hierarchical Porous Waste Cork-Derived Hard Carbon Anode for Enhanced Na Storage Performance[END_REF] . Therefore no individual particles have been seen for samples 70nm-FAU-ZTC and 10nm-FAU-ZTC by TEM. All three samples exhibited homogeneous microporous structure thus demonstrating the successful replication of the FAU micropores [START_REF] Itoi | Facile synthesis of high surface area zeolite-templated carbons using divinylbenzene and propylene as carbon sources[END_REF][START_REF] Wang | Removing ionic and organic mercury from light hydrocarbon liquids by ordered zeolite-templated carbon doped with sulphur[END_REF] . The sodiation and de-sodiation behavior of the ZTCs electrode during the first five charges and discharges were tested by cyclic voltammograms (CV) (Figure 3a, b, andc). Three peaks located 0.01 V, 0.86 V, 1.22 V c h f g c f pl 1μ -FAU-ZTC. The peak at 0.01 V represents the porefilling of Na + and the adsorption of Na + on defects and edges of ZTC [START_REF] Morikawa | Mechanism of Sodium Storage in Hard Carbon: An X-Ray Scattering Analysis[END_REF][START_REF] Li | Commercial activated carbon as a novel precursor of the amorphous carbon for high-performance sodium-ion batteries anode[END_REF] . The two peaks located at 0.86 V and 1.22 V that disappeared in the subsequent cycles are due to the irreversible reaction between Na + and surface functional groups [START_REF] Xia | Hard Carbon Nanosheets with Uniform Ultramicropores and Accessible Functional Groups Showing High Realistic Capacity and Superior Rate Performance for Sodium-Ion Storage[END_REF] , decomposition of electroyte [START_REF] Mahmood | Ultrafast Sodium/Potassium-Ion Intercalation into Hierarchically Porous Thin Carbon Shells[END_REF] , and formation of SEI layers [START_REF] Yu | Sulfur and nitrogen dual-doped carbon nanofiber with enlarged interlayer distance as a superior anode material for sodium-ion capacitors[END_REF] . The second CV scan did not show any peaks except the one at 0.01 V, indicating that a stable SEI layer was formed and the sodium ion storage behavior of porefilling and the adsorption on the edges and defects sites is predominant [START_REF] Zhao | Tubular MoO2 organized by 2D assemblies for fast and durable alkali-ion storage[END_REF][START_REF] Bobyleva | Unveiling pseudocapacitive behavior of hard carbon anode materials for sodium-ion batteries[END_REF] . The first CV curves for the three samples are compared (Figure S3), the observed peaks were found to become larger with decreasing the particle size of the zeolite template crystals, resulting in enhanced sodium storage. The peaks of surface reaction and SEI layer formation are also enhanced due to the higher external surface area of the template [START_REF] Zhang | Size Effects in Sodium Ion Batteries[END_REF] . For the same reason, the de-sodiation current at 3.0-1.0 V strengthens, indicating that the sodium ions adsorbed on the external surface of ZTCs are being released. In contrast to the other two samples, the subsequent CV curves for sample 10nm-FAU-ZTC revealed a pair of redox peaks at 1.53 V and 1.59 V. Peaks at similar positions have been reported earlier; the formation of pseudometals in the micropores of carbon materials was suggested [START_REF] Morita | Combination of solid state NMR and DFT calculation to elucidate the state of sodium in hard carbon electrodes[END_REF][START_REF] Yang | Boosting sodium storage performance of Mo2C via nitrogen-doped carbon sphere encapsulation and rGO wrapping[END_REF] . This emphasizes the adequacy of sodium ions to enter the pores of sample 10nm-FAU-ZTC. The first five galvanostatic charge-discharge (GCD) profiles of pl 1μ -FAU-ZTC, 70nm-FAU-ZTC, and 10nm-FAU-ZTC at 100 mA•g -1 are shown in Figures 3d, e, andf, respectively. All three samples show a voltage profile of a sloping line without distinct plateaus. The curves all start with a descending slope with a continuous curvature in the high-voltage region, which is the consequence of the quick adsorption of Na + to the defect or edge sites of ZTC [START_REF] Chen | Hard carbon for sodium storage: mechanism and optimization strategies toward commercialization[END_REF][START_REF] Zhong | Mechanistic Insights into Fast Charging and Discharging of the Sodium Metal Battery Anode: A Comparison with Lithium[END_REF] . Unlike the micropore filling that has been reported in hard carbon materials, the filling of the supercages of ZTC does not show a long and distinct plateau, but rather a steady decreasing slope at 1.5-0.1V, which is similar to carbon nanotubes and activated carbon with small interspaces used to store sodium ions [START_REF] Xia | Hard Carbon Nanosheets with Uniform Ultramicropores and Accessible Functional Groups Showing High Realistic Capacity and Superior Rate Performance for Sodium-Ion Storage[END_REF][START_REF] Ren | High-performance sodium-ion storage: multi-channel carbon nanofiber freestanding anode contrived via ingenious solvent-induced phase separation[END_REF][START_REF] Morita | Combination of solid state NMR and DFT calculation to elucidate the state of sodium in hard carbon electrodes[END_REF][START_REF] Kamiyama | MgO-Template Synthesis of Extremely High Capacity Hard Carbon for Na-Ion Battery[END_REF][START_REF] Wang | Encapsulating Silica/Antimony into Porous Electrospun Carbon Nanofibers with Robust Structure Stability for High-Efficiency Lithium Storage[END_REF] . The first discharge capacity and first charge capacity of the three ZTC anodes are presented in Table 2. It can be seen that the smaller particle sizes of samples 70nm-FAU-ZTC and 10nm-FAU-ZTC resulted in a better performance compare to sample 1µm-FAU-ZTC. Among all samples, the 10nm-FAU-ZTC with the smallest particle size exhibited the highest first discharge capacity of 2116 mAh•g -1 due to the highest internal space accessibility. The first discharge capacity contributed by surface reaction and SEI layer formation was calculated using the shoulder region around 1.5V. For samples 1µm-FAU-ZTC, 70nm-FAU-ZTC, and 10nm-FAU-ZTC, the contribution of surface reaction and SEI layer formation has been calculated to be 350, 550, and 850 mAh•g -1 , respectively (calculation details presented in Figure S4). After excluding the capacity contributed by the external adsorbed Na + and SEI formation, the first discharge capacity f pl 1 μ -FAU-ZTC, 70nm-FAU-ZTC, and 10nm-FAU-ZTC is equal to 13.7, 22.6, and 31.0 Na + stored in each supercage, respectively. It can be seen that the short diffusion path of the small particles greatly facilitates the Na + to enter the supercage of ZTC. However, the first charge capacity exhibited a significant loss, indicating that only 1.5, 4.9, and 4.6 Na + were able to be re-x c f h up c g f 1μ -FAU-ZTC, 70nm-FAU-ZTC, and 10nm-FAU-ZTC, respectively. The amount of Na + that can be re-extracted from the supercages of 10 nm-FAU-ZTC is lower than that of 70 nm-FAU-ZTC, which is the result of more SEI layer caused by the large external surface area. Thus 70nm-FAU-ZTC has shown the highest first Coulomb efficiency (23.8%). The shoulder representing the formation of SEI disappears in the subsequent GCD curves for all three samples, implying that a stable SEI layer has been formed during the first discharge. In addition to this, the Na + adsorbed on the external surface contributed mostly to the capacity in the subsequent cycles, for the discharge curve is composed mostly by straight lines with a consistent slope [START_REF] Chao | Array of nanosheets render ultrafast and high-capacity Na-ion storage by tunable pseudocapacitance[END_REF] . However, the contribution of the supercage to the capacity still tends to increase with the increase of the external surface area of the ZTC, which indicates that only the Na + stored in the supercages close to the external surface can be reversibly adsorbed and desorbed. The performance of the three samples after the first 100 cycles was also evaluated and the results are summarized in Figure 4. After several cycles, the ZTC still displays an appreciable reversible capacity and a very high coulombic efficiency. The performance of the samples remained stable even after hundreds of cycles (Figure S5), demonstrating the stability of the ZTC structure. For comparison, the performance of other carbon materials reported earlier are also summarized in Table 2. Commercial activated carbon has a high surface area and complex pore structure, but as a typical anode carbon material, its performance is much lower than that of ZTC. Graphene oxide nanosheets and carbon nanotubes have a polycarbon ring structure similar to that of ZTC, but they do not show significant capacity. However, the first coulumbic efficiency of ZTC samples is still lower than the carbon materials presented in Table 2, which is a main difficult for the ZTC to be considered for practical application. In order to fully explore the sodium storage capacity of ZTC materials and verify the sodium storage mechanism, the samples were discharged with a very low current density (9mA•g -1 ) resulting in more mild surface reaction, SEI layer formation, and more complete porefilling. Thus the sodium storage potential of the internal pores of ZTC can be fully exploited. The first and second discharge capacities of the samples are present in Figure 5. Because of the milder surface reaction and SEI formation, there is no obvious plateau around 1.5 V in the first discharge curve to calculate the resulting capacity contribution. However, it can still be seen that at low current densities, 10nm-FAU-ZTC exhibits the highest first discharge capacity of 3555mAh•g -1 , even though it includes the contribution of more SEI layer formed due to higher external surface area. After the first charge and discharge, the contribution of the SEI layer was excluded, and the second discharge capacity of sample 10nm-FAU-ZTC (642mAh•g -1 ) showed that still 8.07 sodium ions could enter and be stored in each supercage during the second cycle. The sodium ions stored pl 1μ -FAU-ZTC and 10 nm-FAU-ZTC are 4.18 and 6.45, respectively. It can be seen that the utilization efficiency of the supercages increases with decreasing the particle size at low concentration difference polarization due to the very low current density. This confirms that the amount of the SEI layer is an important factor affecting the ZTC capacity at 100mA•g -1 45, 49-52 .

Based on the presented results, the filling process of ZTC by sodium ions is illustrated schematically in Figure 6. First, the sodium ions adsorbed on the external surface of the ZTC [START_REF] Kamiyama | MgO-Template Synthesis of Extremely High Capacity Hard Carbon for Na-Ion Battery[END_REF][START_REF] Zhang | Carbon Anode Materials: A Detailed Comparison between Na-ion and K-ion Batteries[END_REF][START_REF] Licht | Carbon Nanotubes Produced from Ambient Carbon Dioxide for Environmentally Sustainable Lithium-Ion and Sodium-Ion Battery Anodes[END_REF] . Second, the sodium ions start to form an SEI layer and enter the ZTC's supercages and adsorb on the inner surface. Finally, the Na + continues to migrate to the inner part of the ZTC until saturated [START_REF] Morikawa | Mechanism of Sodium Storage in Hard Carbon: An X-Ray Scattering Analysis[END_REF][START_REF] Sun | Understanding of Sodium Storage Mechanism in Hard Carbons: Ongoing Development under Debate[END_REF] . During the sodium removal process, only sodium ions adsorbed on the external layer of the ZTC and in the supercages near the external surface can desorb, while the sodium ions adsorbed in the internal voids of ZTC are retained and cannot leave due to the high adsorption energy and long diffusion path. When the sodiation process occurs on ZTC crystals with large particle size, the pores are unable to be fully filled due to the extended diffusion path. In contrast, the short diffusion path of the small ZTC particles can take full advantage of the supercage, and more sodium ions can be filled into the supercage and stored in a desolvated form, but at the same time, the high external surface area is also responsible for the formation of more SEI layers and this hinders the sodium ions from being extracted from the supercage at high current densities. y, ZTC w h c y l f 1 μ , 70 , 10 w y h z . I w shown that the size of the ZTC particles has a significant effect on their performance as an anode for sodium batteries. On the one hand, increasing the crystal size of the ZTC was beneficial to reduce the contact area with the electrolyte, which suppressed the side reactions and the SEI layer generation occurring on the external surface, resulting in a better first Coulomb efficiency. On the other hand, the large size of the crystals increased the diffusion path of sodium ions, and the pores of ZTC cannot be fully utilized, causing the decrease of both the first discharge capacity and reversible capacity. Therefore, sample 70nm-FAU-ZTC with intermediate particle size was found to exhibit the highest first coulombic efficiency (23.8%) and reversible capacity (180mAh•g -1 at a current density of 100mA•g -1 ), while sample 10nm-FAU-ZTC showed the highest first discharge capacity (2116mAh•g -1 at a current density of 100mA•g -1 ). The potential of ZTC materials to store sodium ions was fully explored at very low current densities (9mA•g -1 ): sample 10nm-FAU-ZTC exhibited the highest first discharge capacity (3555mAh•g -1 ) and second discharge capacity (642mAh•g -1 ). The results suggested that the ZTC materials can still be further optimized to improve their performance as a SIB anode material.

Figure 1

 1 Figure 1 XRD patterns (A), N 2 adsorption-desorption isotherms (B), pore size distribution (C), and RAMA p c (D) f pl 1μ -FAU-ZTC (a), 70nm-FAU-ZTC (b), and 10nm-FAU-ZTC (c).

Figure 2

 2 Figure 2 SEM p c u f pl 1μ -FAU (a), 70nm-FAU (b), 10nm-FAU (c); SEM and TEM p c u f pl 1μ -FAU-ZTC (d, g), 70nm-FAU-ZTC (e, h), and 10nm-FAU-ZTC (f, i).

Figure 3

 3 Figure 3 Cyclic voltammograms (CV) at a scan rate of 0.1 mV•s -1 and galvanostatic charge-discharge profiles (GCD) tested at a current density of 0.1 A•g -1 for samples 1μ -FAU-ZTC (a, d), 70nm-FAU-ZTC (b, e), and 10nm-FAU-ZTC (c, f) at five cycles

Figure 4 1 F

 41 Figure 4 The cycl p f c f pl 1μ -FAU-ZTC (a), 70nm-FAU-ZTC (b), and 10nm-FAU-ZTC (c) at 100mA•g -1

Figure 6

 6 Figure 6 Schematic diagram of ZTC with different particle sizes storing sodium ions on the external surface, internal surface of supercages and porefilling

  7 Al 2 O 3 : 10 SiO 2 : 160 H 2 O (70 nm-FAU) or 9 Na 2 O: 0.7 Al 2 O 3 : 10 SiO 2 : 160 H 2 O (10 nm-FAU) was aged at ambient conditions for 24 h, and then crystallized at 120 °C for 70 min or 50 °C for 45 h,

	p c ly. Th 1 μ FA z l w p p	f	u p	w h	l	f 8 2 O:
	1 Al 2 O 3 : 10 SiO 2 : 400 H 2

Table 1 N

 1 2 adsorption analysis data of ZTC samples.

  [START_REF] Slater | Sodium-Ion Batteries[END_REF]

Table 2

 2 The initial discharge capacity, initial charge capacity, and stable reversible capacity of ZTC samples at 0.1A•g -1 .

	Sample	First discharge capacity (mAh•g -1 )	First charge capacity (mAh•g -1 )	First Coulumbic efficiency (%)	Stable (mAh•g -1 ) a reversible capacity	Reference
	1μ -FAU-ZTC	1718	288	16.8	120	This work
	70nm-FAU-ZTC	1767	420	23.8	180	This work
	10nm-FAU-ZTC	2116	360	17.0	170	This work
	Activated carbon	~820	~150	~18	~100	39
	Graphene oxide nanosheet	~240	~210	~87	~150	40
	Carbon nanotubes	~80	~70	~87	~60	43

ACKNOWLEDGMENT

The authors gratefully acknowledged funding from the National Key Research and Development Program of China of Ministry of Science and Technology (2022YFE116000), the National Natural Science Foundation of China [Grant No. 22175200, No. 21975285], and Innovation Fund of Science and Technology Project of Fujian Province (2022C0021). Thanks to the Yanpeng Li from Advanced Chemical Engineering and Energy Materials Researcher Center in China University of Petroleum for the STEM characterization.


Maximilian R Becker

Roland R Netz
email: rnetz@physik.fu-berlin.de

Philip Loche

Jan Douwe

Bonthuis

Dominique Mouhanna

Hélène Berthoumieux
email: helene.berthoumieux@espci.fr

A H Hansen

L-M Hodgson

D T Liu

A Limmer

P Michaelides

Pedevilla

Dielectric properties of aqueous electrolytes at the nanoscale

Despite the ubiquity of nanoconfined aqueous electrolytes, a theoretical framework that accounts for the nonlinear coupling of water and ion polarization is still missing. We introduce a nonlocal and nonlinear field theory for the nanoscale polarization of ions and water and derive the electrolyte dielectric properties as a function of salt concentration to first order in a loop expansion. Classical molecular dynamics simulations are favorably compared with the calculated dielectric response functions. The theory correctly predicts the dielectric permittivity decrement with rising salt concentration and furthermore shows that salt induces a Debye screening in the longitudinal susceptibility but leaves the short-range water organization remarkably unchanged.

Introduction -The study of nanoconfined electrolytes is exciting both for their ubiquity and for the theoretical challenge they bring [START_REF] Maggs | Simulating nanoscale dielectric response[END_REF][START_REF] Fuentes-Azcatl | Non-Polarizable Force Field of Water Based on the Dielectric Constant: TIP4P[END_REF][3]. The nanometer scale is the typical confinement size of technological and biological devices, the screening length of low concentrated ionic solutions, as well as the range at which water starts to behave as a discrete molecular medium [4][5][6]. The structure of interfacial solutions thus results from a subtle interplay between short-range charge overscreening generated by the solvent [7], ion-ion correlations, water-and ionsurface interactions [8,9]. The balance between these effects encode many aspects of the double layer structure [10,11].

The continuous linear description of water used in Poisson-Boltzmann (PB) model cannot capture this complexity. The well-documented decrement of the bulk permittivity of aqueous ionic solutions induced by an increase of the concentration of the salt [START_REF] Hasted | Dielectric properties of aqueous ionic solutions. part i and ii[END_REF] has been retrieved using nonlinear expansions of a dipolar PB model [START_REF] Levy | Dielectric constant of ionic solutions: A field-theory approach[END_REF]. But a parameterized nonlocal field theory (FT) for structured solvent to fully describe electrolytes at the nanoscale is missing [START_REF] Paillusson | Slits, plates, and poissonboltzmann theory in a local formulation of nonlocal electrostatics[END_REF][START_REF] Ben-Yaakov | Ion-specific hydration effects: Extending the poisson-boltzmann theory[END_REF][START_REF] Blossey | Continuum theories of structured dielectrics[END_REF][START_REF] Blossey | A comprehensive continuum theory of structured liquids[END_REF]. For instance, the simple question "Do electrolytes increase or decrease the range of the longitudinal correlations of water?" is still actively debated [START_REF] Omta | Negligible effect of ions on the hydrogen-bond structure in liquid water[END_REF][START_REF] Chen | Electrolytes induce long-range orientational order and free energy changes in the H-bond network of bulk water[END_REF][START_REF] Zhang | Dissolving salt is not equivalent to applying pressure on water[END_REF]. In contrast, the influence of ionic strength on transverse polarization correlations is largely unexplored although it could be of major impor-FIG. 1. Sketch of the system under study. Water is described as a nonlocal nonlinear continuous dielectric medium and ions as point charges. We evaluate the ionic-strength effect on the water longitudinal P and transverse P ⊥ polarization correlations.

tance for the interactions between solvated objects [START_REF] Schoeger | Universal casimir interaction between350 two dielectric spheres in salted water[END_REF]. The figure 1 shows an illustration of these questions.

Nonlocal model for water -At the molecular scale, water is a correlated dielectric medium, as illustrated by the wave mode dependence of the correlation tensor of the polarization field P [START_REF] Bopp | Static Nonlocal Dielectric Function of Liquid Water[END_REF][START_REF] Bopp | Frequency and wave-vector dependent dielectic function of water: Collective modes and relaxation spectra[END_REF]. In Fourier space, the longitudinal susceptibility of water, χ w (q) = β P (q) • P (-q) / 0 , -with β the inverse temperature, 0 the vacuum permittivity, P the longitudinal part of the polarization, q the wave vector -exhibits an overresponse peak around q=3 Å-1 . This structured response characterizes a charge overscreening at short range, induced by the layering of water molecules organized by the water hydrogen bond network [START_REF] Bopp | Static Nonlocal Dielectric Function of Liquid Water[END_REF]. The transverse susceptibility χ w ⊥ (q) = β P ⊥ (q) • P ⊥ (-q) / 0 , with P ⊥ the transverse part of the polarization -exhibits a Lorentzian decay at low q [START_REF] Bopp | Frequency and wave-vector dependent dielectic function of water: Collective modes and relaxation spectra[END_REF].

Continuum nonlocal electrostatics provides a useful framework to describe correlated fluids [START_REF] Kornyshev | On the non-local electrostatic theory of hydration force[END_REF][START_REF] Basilevsky | Nonlocal continuum solvation model with exponential susceptibility kernels[END_REF][START_REF] Hildebrandt | Novel formulation of nonlocal electrostatics[END_REF]. The electrostatic energy of the fluid U el can be written as a Gaussian functional of P,

U el [P] = 1 2 drdr ∇ • P(r)∇ • P(r ) 4π 0 |r -r | + U G conf [P]. ( 1
)
The first term corresponds to the bare Coulomb interactions between the partial charges -∇ • P(r) of the fluid. The phenomenological configuration energy U G conf [P], where the subscript G stands for Gaussian, can be made explicit using a Landau-Ginzburg expansion [START_REF] Maggs | Simulating nanoscale dielectric response[END_REF] such as,

U G conf [P] = 1 2 0 dr KP(r) 2 + κ l (∇ • P(r)) 2 + κ t (∇ × P(r) 2 + α(∇(∇ • P(r))) 2 (2)
expanded up to the first spatial derivative for the transverse terms (term in κ t ) and up to the second one for the longitudinal ones (terms in κ l and α). Eq. ( 2) has been shown to capture the main features of dielectric properties of water at the nanoscale [6,[START_REF] Berthoumieux | Fluctuation-induced forces governed by the dielectric properties of water-a contribution to the hydrophobic interaction[END_REF]. The polarization susceptibility χ w , obtained by inversion of Eq. ( 1), can be decomposed in Fourier space in a longitudinal χ w and a transverse χ w ⊥ response using isotropy and homogeneity of the system, such that χ ij (q) = χ w (q)q i q j /q 2 + χ w ⊥ (q)(δ ijq i q j /q 2 ), with (i, j) = (x, y, z). Their expressions follow from Eq. (2) as

χ w (q) = 1 1 + K + κ l q 2 + αq 4 , χ w ⊥ (q) = 1 K + κ t q 2 .
(3) χ w is associated with a decay λ l and an oscillating λ o correlation length and χ w ⊥ with a transverse correlation length λ t for which explicit expressions are given in S1.1. The bulk permittivity w of water is defined as [START_REF] Bopp | Static Nonlocal Dielectric Function of Liquid Water[END_REF]. Figure 2 shows χ w and χ w ⊥ (blue line) for which Eq. ( 2) has been parameterized to reproduce MD simulated TIP4p/ water [START_REF] Fuentes-Azcatl | Non-Polarizable Force Field of Water Based on the Dielectric Constant: TIP4P[END_REF]. Details are provided in S1.1 and the values of the Landau-Ginzburg parameters (K, κ l , κ t , α ) are given in the caption of Fig. 2.

w = (1 -χ (0)) -1 = 1 + 1/K
A term scaling like P 4 , unfavorable to large values of the polarization field, can be added to correct the Gaus-sian model. The configuration energy [START_REF] Maggs | Simulating nanoscale dielectric response[END_REF] thus reads:

U conf [P] = 1 2 0 dr γ(P(r) 2 + P 2 0 ) 2 + κ l (∇ • P(r)) 2 + κ t (∇ × P(r)) 2 + α(∇(∇ • P(r))) 2 (4)
and takes into account the saturation of water polarization in the presence of strong fields [START_REF] Alper | Field strength dependence of dielectric saturation in liquid water[END_REF][START_REF] Fedorov | Unravelling the solvent response to neutral and charged solutes[END_REF][START_REF] Paillusson | Dielectric response in the vicinity of an ion: A nonlocal and nonlinear model of the dielectric properties of water[END_REF], in addition to nanometric correlations. Such a medium presents a threshold saturation, P 0 , between a linear response (|P| P 0 ) and a saturation response (|P| P 0 ) regime. We impose P 0 = K/2γ so that the quadratic local terms in Eq. ( 2) -KP 2 -and Eq. ( 4) -2γP 2 0 P 2 become equal. Thus, in the low-field regime, the two models are associated with the same nonlocal properties [START_REF] Paillusson | Dielectric response in the vicinity of an ion: A nonlocal and nonlinear model of the dielectric properties of water[END_REF] and γ remains the only degree of freedom tuning the saturation threshold.

In this letter, we compute the Gibbs free energy and the nonlocal dielectric susceptibility of an aqueous electrolyte for which water is modeled by Eq. ( 2). We perform MD simulations and compare the Gaussian model to a simulated system. We include saturation effects using Eq. ( 4) including the one-loop expansion. We determine the effect of the ionic strength on the bulk permittivity, the longitudinal and transverse water correlations. Finally, we identify the essential building blocks to construct a field theory modeling electrolytes at the nanoscale.

Gaussian model for electrolytes -We consider an electrolyte with N + punctual cations of charge e and N - punctual anions of charge -e solvated in water modeled with Eq. ( 2). The ionic charge density reads ρ(r) = Σ N + i=1 eδ(r-r + i )-Σ N - j=1 eδ(r-r - j ). In the canonical ensemble, the partition function of the system can be written as

Z G = 1 N + ! 1 N -! N + i=1 dr + i   N - j=1 dr - j   × D[P]e -βU G conf [P] e -β 2
drdr ρtot(r)v(r-r )ρtot(r ) (5) with ρ tot (r) = ρ(r) -∇ • P(r) and v(r) = 1/4π 0 |r|. Z G includes the configurational degrees of freedom of the solvent and the Coulomb interactions between free and partial charges. Performing a Hubbard-Stratonivich transformation to get rid of the longrange potential v and switching to the grand canonical ensemble for a more tractable expression, we obtain the grand-canonical ensemble partition function, [START_REF] Levy | Dielectric constant of ionic solutions: A field-theory approach[END_REF]. The action F G u [P, Ψ] is a functional of P and of the electrostatic potential Ψ. Assuming a 1:1 electrolyte and an ionic density n, we find: 8) -as a function of q for water and solutions of increasing salt concentration. The arrow indicates the increase of the peak maximum with an increase of the concentration. (b) Transverse susceptibility χ G ⊥ = χ w ⊥ given by Eq (3) as a function of q, which is identical for water and electrolytic solutions. The parameter values of the FT model given in Eq. ( 2) are K= 1/76, κ l = -0.218 Å-2 , α=0.012 Å-4 and κt=0.013 Å-2 .

Ξ G = D[P] D[Ψ]e -βF G u [P,Ψ]
F G u [P, Ψ] = U G conf [P] -dr 0 2 (∇Ψ) 2 -Ψ∇ • P - 2n β cosh(βΨe) . (6)
See S1.2 for a detailed derivation. The free energy F G of the system follows from F G = -k B T lnΞ G . The associated susceptibility is a 4×4 matrix defined as the inverse of the second functional derivative of F G as

0 χ G χ G P,ψ χ G ψ,P χ G ψ,ψ / 0 (r 1 -r 2 ) =   δ 2 F G u (P ,Ψ) δPi(r1)δPj (r2) δ 2 F G u (P ,Ψ) δPi(r1)δΨ(r1) δ 2 F G u (P ,Ψ) δΨ(r1)δPj (r2) δ 2 F G u (P ,Ψ) δΨ(r1)δΨ(r2)   -1 (7)
where (P , ψ) are the mean fields minimizing the action F G u . They both vanish (See S1.3). As the medium is homogeneous and isotropic, the 2-point susceptibility is a function of the distance r 1r 2 . We focus here on the polarization susceptibility, which is written in Fourier space as χ G ij (q) = χ G (q)q i q j /q 2 + χ G ⊥ (q)(δ ijq i q j /q 2 ). The longitudinal term reads

χ G (q) = w λ 2 D + q 2 w λ 2 D + q 2 (K + κ l q 2 + αq 4 ) + q 2 , (8)
with the Debye length λ D = 0 w /2βne 2 . See S1.4 for details. Fig. (2) a) shows χ G (q) for increasing salt concentration c = n/N a , where N a is the Avogadro number. The Gaussian model predicts an enhancement of the water ordering with an increase of c, as indicated by the magnitude increase of the peak at q=3 Å-1 . This can be understood as follows. In the nonlocal Gaussian framework, an ion, located in r = 0 generates the electrostatic potential, 1/4π 0 w r, and an extra potential decaying on a range λ l and oscillating around zero with a period λ o [START_REF] Vatin | Electrostatic interactions in water: a nonlocal electrostatic approach[END_REF]. This oscillating landscape leads to the organization of the charges and to longer-range correlations that increase with the salt concentration until a nonphysical crystallization of the system occurs, corresponding to a divergence of χ G (q). See S1.4 for details. In contrast, for a very diluted solution λ D → ∞ (see green line in Fig. 2), the longitudinal susceptibility can be approximated as χ G (q) ≈ ( w /λ 2 D + q 2 )/(K( w /λ 2 D + q 2 ) + q 2 ) + χ w (q), the sum of the Debye response function for an homogeneous electrolyte and the pure water spectrum. The Debye length λ D and the water-water correlation lengths λ l and λ o do not couple. The Debye contribution dominates for the low q limit and generates a Lorentzian decay as seen in Fig. [START_REF] Fuentes-Azcatl | Non-Polarizable Force Field of Water Based on the Dielectric Constant: TIP4P[END_REF]. The transverse susceptibility χ G ⊥ (q) is unaffected by the presence of salt and obeys χ G ⊥ (q)=χ w ⊥ (q). Indeed, the coupling between the salt and the solvent occur via the Coulomb interactions and involves only the longitudinal part of the polarization as seen in Eq. ( 5). Finally, we note that the Gaussian model predicts that the dielectric bulk properties of electrolyte solutions χ G (0) = χ G ⊥ (0) = 1/K are independent of the salt concentration.

Comparison with MD simulations -To check the validity of the Gaussian model, we compare its predictions with the dielectric properties of simulated solutions of NaCl in TIP4p/ water for concentrations c up to 1.5 mol.L -1 [START_REF] Fuentes-Azcatl | Non-Polarizable Force Field of Water Based on the Dielectric Constant: TIP4P[END_REF][START_REF] Loche | Transferable Ion Force Fields in Water from a Simultaneous Optimization of Ion Solvation and Ion-Ion Interaction[END_REF]. Simulation details are given in S2. We compute the bulk permittivity and plot (c) in Fig. 3 (a). We observe a linear decay, which is not described by the Gaussian model that predicts a constant permittivity. The susceptibilities χ MD ⊥ (q) and χ MD (q) are shown in Fig 3 (b) and (c) for c=0.15 mol.L -1 -λ D =7.8 Å-, c=0.75 mol.L -1λ D =3.5 Å-, and c=1.5 mol.L -1 -λ D =2.5 Å. The blue markers show the response for pure water. We observe a decay of the bulk value χ MD ⊥ (0) for an increasing concentration, which is not captured by the Gaussian model. In the longitudinal case, we see a very low coupling between the Debye contribution (Lorentzian decay for q ≤2 Å-1 ) and the spectrum of water (q >2 Å-1 ), which remains almost unperturbed. A concentration increase induces a decrease and a flattening of the pseudo resonant peak. The MD results indicate thus an effect of the salt opposite to the Gaussian model predictions that appears at much higher concentration.

A Nonlinear model for the solvent -To get a better agreement between FT and MD simulations, we consider the nonlinear configuration energy U conf in Eq. ( 4). We derive the inverse susceptibility of the system to first order in a loop expansion [START_REF] Levy | Dielectric constant of ionic solutions: A field-theory approach[END_REF]. The grand partition function obeys Ξ = D[P] D[Ψ]e -βFu[P,Ψ] with F u [P, Ψ] obtained by replacing U G conf by U conf in Eq. ( 6). The mean fields (P , ψ) still vanish, see S1.3. The action F u is expanded up to the second order in (P, Ψ) around the mean field solution,

F u [P, Ψ] ≈ F u [P , ψ] + 1/2 drdr (δP (r), δψ(r)) • F (2) u (P , ψ) • (δP (r ), δψ(r )), with F (2) u
the second functional derivative of F u with respect to (P, Ψ) and δP = P -P , δψ = Ψψ. See S1.7 for details. The free energy is then written as F ≈ F u [P , ψ] + Trln(βF

(2) u [P , ψ])/β. The inverse susceptibil- ity follows as χ -1 (r 1 -r 2 ) = χ -1,G (r 1 -r 2 ) + χ -1,1 (r 1 - r 2
), with the one loop correction term χ -1,1 (r 1r 2 ) defined as:

χ -1,1 i,j (r 1 -r 2 ) = 1 2β δ 2 Tr ln βF (2) u δP i (r 1 )δP j (r 2 ) (P , ψ). (9)
χ -1,G is the inverse of the Gaussian susceptibility χ G given in Eq. ( 8). Performing the field derivatives and calculating the trace in Eq. ( 9), one obtains

χ -1,1 i,j (r 1 -r 2 ) = δKδ(r 1 -r 2 )δ ij with δK = 20γ 0 3β χ G (r c ) + 2χ G ⊥ (r c ) . (10)
The first order correction of the susceptibility is purely local and proportional to γ, which tunes the saturation regime of the model, and depends on the Debye length via χ G . The cutoff distance r c has been introduced to get rid of the divergence of χ G in r = 0. See S1.7 for a detailed derivation and expressions of χ G i (r c ), i = , ⊥. It can be written in Fourier space as χ -1,1 ij (q) = δKq i q j /q 2 + δK(δ ijq i q j /q 2 ). The corrected permittivity, = 1+1/(K+δK) depends now on λ D and thus on the salt concentration. δK is expanded linearly in c as δK(c) = δK w + δK c c + τ (c 2 ), where δK w is the one-loop expanded correction to K for pure water and δK c c is the correction induced by the salt at linear order. Their expressions are derived in S1.8. Setting δK w to zero as it is included in the fitted value of K, and performing a linear expansion in c of (c) = 1 + 1/(K + δK c c), we find the permittivity,

(c) = w - δK c K 2 c, with (11)
δK c = 5γN a e 2 ( w -1) 2 24π 2 w (4π 2 λ 2 l + λ 2 o )(4π 2 λ 2 l + 5λ 2 o ) λ l λ 4 o
The permittivity decrement is here expressed as a function of the bulk permittivity and the intrinsic longitudinal correlation lengths of the fluid.

Susceptibility kernels for electrolytes -We now compare the polarization susceptibilities obtained from MD simulations and the nonlinear FT models. Figure 3 (a) shows (c), Eq. ( 11), for a value of γ adjusted to reproduce the MD data and given in the caption of Fig 3.

The one-loop corrected transverse response,

χ FT ⊥ (q) = 1 K + δK c c + κ t q 2 , ( 12
)
is plotted in Fig. 3 (b) for c=0, 0.15, 0.75, 1.5 mol.L -1 and shows very good agreement with simulations. The one-loop expanded longitudinal susceptibility is obtained by replacing K by K + δK c c in Eq. (8). Similarly to the Gaussian model, it foresees an enhancement of the longitudinal correlations (See S1.9) and thus fails to reproduce MD data. 12), (c) MD simulated longitudinal susceptibility χ MD , (d) conjectured expression of the longitudinal susceptibility χ FT given in Eq. [START_REF] Levy | Dielectric constant of ionic solutions: A field-theory approach[END_REF]. Susceptibilities are plotted as function of q for pure water and electrolytes of increasing concentration.

Guided by the absence of coupling between the Debye wavelength and the water-water correlation spectrum illustrated by the simulations, we propose an ad hoc longitudinal susceptibility as follows:

χ FT (q) = w λ 2 D + q 2 (K + δK c c)( w λ 2 D + q 2 ) + q 2 + 1 1 + K + δK c c + κ l q 2 + αq 4 . ( 13
)
The first term corresponds to the susceptibility of an homogeneous electrolyte associated with the corrected permittivity (c) given in Eq. (11). The second term corresponds to the nonlocal susceptibility of pure water associated with this corrected permittivity. We plot χ FT (q) for c=0, 0.15, 0.75, 1.5 mol.L -1 in Fig. 3 (d). It reproduces well MD data as it presents a Lorentzian decay at low q and a flattening of the pseudo-resonant peak at q=3 Å-1 for an increase of the salt concentration. Discussion -In this work, we have derived analytic expressions for the dielectric response functions of electrolytes at the nanoscale. To do so, we have compared susceptibilities calculated from a FT including nonlocal and nonlinear behavior of water and susceptibilities derived with classical MD simulations. We have thus identified the key effects of the salt on the medium organization. For the longitudinal modes, we highlight two length scales that do not couple to each other: at small wave-modes, q ≤2 Å-1 , the medium can be described as homogeneous with a permittivity (c) decaying with the salt concentration. This corresponds to long-range interactions for which the Debye screening occurs. At larger q, the water longitudinal susceptibility is similar to the one of pure water but associated with a corrected permittivity (c). These two decoupled q-domains could indicate two "types" of water molecules that are spatially separated: the one solvating ions polarize in a saturated manner in response to the ionic field, creating an "electrically dead" solvation shell [START_REF] Levy | Dielectric constant of ionic solutions: A field-theory approach[END_REF]. Outside of this hydration shell, the molecular water organization remarkably is unaffected by the ions, as recently predicted by machine learning based simulations [START_REF] Zhang | Dissolving salt is not equivalent to applying pressure on water[END_REF]. We note that local nonlinear models of water also foresee two response regimes (a linear response at low electrostatic excitation and a saturated response at high excitation) with an abrupt switch between the two regimes [START_REF] Berthoumieux | Dipolar Poisson modeels in a dual view[END_REF]. Moreover, our work reveals the absence of coupling between salt screening and transverse polarization modes of water. This surprising result could have important consequences on the interactions between objects immersed in electrolytes [START_REF] Pires | Probing the screening of the casimir interaction with optical tweezers[END_REF] that were assumed to be screened in any circumstances. This study gives a clear picture of the nature and the range effect of the salt on water organization at the nanoscale for unconfined solutions. This paves the way to develop a field theory describing the properties of nanoconfined electrolytes. be written following a Landau-Ginzburg approach. We consider it here up to the second spatial derivative in the polarization P to reproduce the dielectric response of water and set

HB thanks

U G conf [P] = 1 2 0 drdr P(r) • K(r -r ) • P(r ) with P(r)K(r -r )P(r ) = KP(r) 2 + κ l (∇ • P(r)) 2 + κ t (∇ × P(r)) 2 + α(∇(∇ • P(r))) 2 δ(r -r ). (S2)
Note that as we consider bulk water which is homogeneous and isotropic, the two-point kernel K depends simply on the distance rr . We introduce the two-point dielectric susceptibility χ w (rr ) defined as:

U el [P] = 1 2 0 drdr P(r) (χ w (r -r )) -1 P(r ) (S3)
where the superscript w stands for water considered here as a pure solvent. The matrix χ w can be decomposed in Fourier space into a longitudinal χ w and a transverse χ w ⊥ function, such that χ w ij (q) = χ w (q) q i q j q 2 + χ w ⊥ (q)(δ ij -

q i q j q 2 ) χ w (q) = 1 1 + K + κ l q 2 + αq 4 , χ w ⊥ (q) = 1 K + κ t q 2 .
(S4)

The figure S1 shows the MD simulated longitudinal a) and transverse b) response function for the 3 charge water model TIP4p/ [START_REF] Fuentes-Azcatl | Non-Polarizable Force Field of Water Based on the Dielectric Constant: TIP4P[END_REF] (blue markers). MD data show an overesponse of the system χ w (q) >1 around q=3 Å-1 , which corresponds to a negative permittivity, as (q) obeys (q) = 1/(1χ w (q)) [3]. This overscreening phenomenon can be attributed to the H-bond network structuring water at short range [4]. Note that the secondary peak for q >4 Å-1 corresponds to intramolecular correlations and to length scales that are not addressed in this study. The transverse susceptibility is also extracted from MD simulations and shows a spectrum of simpler shape, as it monotonously decays.

The minimal model given in Eq. (S2) can reproduce these main features. For the longitudinal susceptibility, we consider the case (κ l < 0, α > 0). The value of the parameters K, κ l , and α, are adjusted to fit the bulk value, χ w (0), the position and the value of the maximum of χ w (q) of the the MD simulated response. The parameter κ t is chosen to fit the decay of MD data. Fig. S1 shows the FT response functions (black line). The values of the parameters are given in the caption.

The longitudinal susceptibility χ w (q) is associated with two correlation lengths: a longitudinal decay, λ d and an oscillation length λ o , defined as the imaginary and real part of the inverse of the poles of the function. The transverse susceptibility χ w ⊥ (q) is associated with a decay length λ t which is the inverse of its pole. Their expressions obey:

λ d = 2 √ α 2 α(1 + K) + κ l , λ t = κ t K , λ o = 4π √ α 2 α(1 + K) -κ l . ( S5
)
The polarization-polarization correlation function in the real space can be derived by Fourier transforming the susceptibility given in Eq.(S4) and has been studied in a previous work [5]. Using the estimated values of the parameters we get a longitudinal decay length λ d =4.7 Å, an oscillating length λ o =2.1 Å, and a transverse decay length, λ t =1.05 Å.

S1.2. Derivation of Ξ, the partition function in the grand canonical ensemble

In this subsection, we present the main steps of the derivation of the partition function Ξ G in the grand canonical ensemble. We start from the partition function of the canonical ensemble, Z G ,

Z G = 1 N + ! 1 N -! N + i=1 dr + i   N - j=1 dr - j   D[P]e -βU G conf [P] e -β 2 drdr ρtot(r)v(r-r )ρtot(r ) (S6)
with v(rr ) = 1/4π 0 |rr |, the Coulomb potential and ρ tot (r) = ρ(r) -∇ • P(r). We introduce an auxiliary field Φ and perform a Hubbard-Stratonovich transformation using the relation v(rr ) -1 = -0 ∇ 2 δ(rr ) [6]. Dropping the prefactor, we get

Z G = 1 N + ! 1 N -! N + i=1 dr + i   N - j=1 dr - j   D[P]e -βU G conf [P] D[Φ]e -β 2 dr 0(∇Φ) 2 -iβ drΦ(r)(ρ(r)-∇•P(r)) . (S7)
Introducing the expression of the charge density ρ(r

) = Σ N + i=1 eδ(r -r + i ) -Σ N - j=1 eδ(r -r - j ), we get, Z G = D[Φ] 1 N + ! dre -iβeΦ(r) N+ 1 N -! dre iβeΦ(r) N- e -β 2 dr 0(∇Φ(r)) 2 × D[P]e -βU G conf [P] e -iβ drΦ(r)(-∇•P(r)) (S8)
The partition function is brought into a more manageable form by going to the grand canonical ensemble, where we get

Ξ G = ∞ N +=0 e -βµ+N+ N + ! dre -iβeΦ(r) N + × ∞ N -=0 e -βµ-N- N -! dre iβeΦ(r) N - × D[P]e -βU G conf [P] D[Φ]e -β 2 dr 0 (∇Φ(r)) 2 -iβ drΦ(r)(-∇•P(r)) , ( S9
)
with µ i , (i = +, -), the chemical potential for cations and anions. We introduce the field Ψ = iΦ that can be identified as the electrostatic potential [6], we perform the sums over N ± and identify the action

F G u [Ψ, P], such that Ξ G = D[P] D[Ψ]e -βF G u [P,Ψ] . (S10)
Setting e -βµ± = n, the ionic density, we obtain

F G u [P, Ψ] = U G conf [P] -dr 0 2 (∇Ψ) 2 -Ψ∇ • P - 2n β cosh(βΨe) . (S11)
In this section we derive the expressions of the longitudinal χ G (q) and transverse χ G ⊥ (q) polarization susceptibility of an electrolyte. We first define the total susceptibility of the system as

χ tot (r 1 -r 2 ) = 0 χ G χ G P,ψ χ G ψ,P χ G ψ,ψ / 0 (r 1 -r 2 ) =   δ 2 F G u (P ,Ψ) δPi(r1)δPj (r2) δ 2 F G u (P ,Ψ) δPi(r1)δΨ(r1) δ 2 F G u (P ,Ψ) δΨ(r1)δPj (r2) δ 2 F G u (P ,Ψ) δΨ(r1)δΨ(r2)   -1 (S18)
with (P , ψ) the mean fields minimizing the action F G u . Using the expressions given in Eq. (S4), we find for the inverse susceptibility of the system, in Fourier space,

χ -1 tot (q) =     iq x χ -1G 0 (q) iq y iq z -iq x -iq y -iq z -(2ne 2 β + 0 q 2 )     , (S19)
with the matrix χ -1G (q) given by χ -1,G ij (q) = (K + κ l q 2 + αq 4 ) q i q j q 2 + (K + κ t q 2 ) δ ij -q i q j q 2 . (S20) We inverse the matrix χ -1 tot (q) by using a block inversion and obtain

0 1 2 3 4 5 q [ Å-1 ]
χ tot (q) = 0 χ G (q) iQ -iQ χ ψ,ψ (q) 0 (S21)
with Q = (q x , q y , q z ) and

χ G ij (q) = χ G (q) q i q j q 2 + χ G ⊥ (q)(δ ij - q i q j q 2 ) (S22) χ G (q) = w λ 2 D + q 2 w λ 2 D + q 2 (K + κ l q 2 + αq 4 ) + q 2 , χ G ⊥ (q) = 1 K + κ t q 2 (S23) χ ψ,ψ (q) = - K + κ l q 2 + αq 4 w λ 2 D + q 2 (K + κ l q 2 + αq 4 ) + q 2 (S24)
where we have introduced the Debye length λ D = 0 w /2βne 2 . We plot the susceptibility χ G (q) for increasing concentrations (c = n/N a , N a the Avogadro number) in Fig. S2. As we have (κ l <0, α >0) the denominator of χ G (q) will vanish for a small enough Debye length inducing a divergence of the susceptibility. See Fig. S2 , grey dashed line. For the chosen set of parameters, the divergence occurs for c=22 mmol.l -1 . At this concentration, the longitudinal correlations are purely oscillating illustrating an unphysical crystalization of the medium.

S1.5. Susceptibility for Gaussian model in real space

We express the Gaussian susceptibility χ G given in Eq. (S23) in real space. We perform the Fourier transform

χ G ij (r)=F T (χ G ij (q)) defined as χ G ij (r) = 1 2π 3 2π 0 dφ π 0 dθ sin(θ) ∞ 0 dqq 2 e iqr cos(θ) χ G ij (q) (S25)
χ G is here expressed in the intrinsic basis in which the vector r, joining the two correlated points, is aligned with the e z direction of wavemode q basis, as illustrated in Fig. S3. We develop the longitudinal projector q i q j /q 2 (i, j = x, y, z) in the spherical basis, Using Eqs. (S22,S23), χ G ij (r) is splitted into two contributions: χ G ij (r) = F T (χ G (q)q i q j /q 2 ) + F T (χ G ⊥ (q)(δ ijq i q j /q 2 )). We perform the integrals in Eq. (S25) using (S26). We get:

q i q j q 2 =   sin(θ) 2 cos(φ) 2 sin(θ) 2 cos(φ) sin(φ) sin(θ) cos(θ) cos(φ) sin(θ) 2 cos(φ) sin(φ) sin(θ) 2 sin(φ) 2 sin(θ) cos(θ) sin(φ) sin(θ) cos(θ) cos(φ) sin(θ) cos(θ) sin(φ) cos(θ) 2   , (S26)
F T (χ G (q)q i q j /q 2 ) =   1 2 (I 1, (r) -I 2, (r)) 0 0 0 1 2 (I 1, (r) -I 2, (r)) 0 0 0 I 2, (r)   , (S27)
and

F T χ G ⊥ (q)(δ ij -q i q j /q 2 ) =   1 2 (I 1,⊥ (r) + I 2,⊥ (r)) 0 0 1 2 (I 1,⊥ (r) + I 2,⊥ (r)) 0 0 0 I 1,⊥ (r) -I 2,⊥ (r)   . (S28)
We have introduced four elementary functions (I ,1 , I ⊥,1 I ,2 , I ⊥,2 ), defined as follow:

I i,1 (r) = 1 (2π) 3 ∞ 0 dqq 2 π 0 dθ sin(θ) 2π 0 dφ χ G i (q)e iqr cos(θ) , I i,2 (r) = 1 (2π) 3 ∞ 0 dq π 0 dθ sin(θ) 2π 0 dφ χ G i (q)e iqr cos(θ) , (S29)
with i = , ⊥. The Gaussian susceptibility finally reads

χ G (r) =   χ G ⊥ (r) 0 0 0 χ G ⊥ (r) 0 0 0 χ G (r)   (S30) with χ G (r) = I 2, (r) + (I 1,⊥ (r) -I 2,⊥ (r)) (S31) χ G ⊥ (r) = I 1, (r) -I 2, (r) 2 + I 1,⊥ (r) -I 2,⊥ (r) 2 (S32)
in the basis {e i }, i = x, y, z) defined such that r is aligned with e z . See sketch in Fig. S3. The suscpetibility associated with any distance vector r = (x, y, z)) is obtained by performing the following change of basis:

χ G cart = R -1 • χ G (r) • R, (S33)
R being the change-of-basis matrix from Cartesian to spherical coordinates. and

∂ 2 F (2) u (P , ψ) ∂P x (r 1 )∂P y (r 2 ) =    0 4γ/ 0 0 0 4γ/ 0 0 0 0 0 0 0 0 0 0 0 0    δ(r -r 1 )δ(r -r 2 )δ(r -r) (S44)
The other matrices are easily deduced by symmetry. We now have to calculate the trace of matrices in (S43) by integrating over the continuous indices dr dr and summing over the discrete indices. The polarization correlations depend only on the distance u = |rr |. We thus write:

drdr F (2) u -1 • ∂ 2 F (2) u ∂P i (r 1 )∂P j (r 2 ) = dr duu 2 dφ dθ sin(θ) F (2) u -1 ∂ 2 F (2) u ∂P i (r 1 )∂P j (r 2 ) (S45)
Using δ(rr ) = δ(u)/2πu 2 and F

(2) u -1

(P , ψ) = χ G P ,ψ (rr ) given in Eq. (S40), we perform the integral over θ and φ and find for the polarization correlation:

2π 0 dφ π 0 dθ sin(θ)R -1 • χ G ij (u) • R = 4π 3 χ (u) + 2χ ⊥ (u) δ ij . ( S46
)
We now perform the matrix product of matrix given in Eqs. (S43, S46), calculate its trace and find:

χ -1,1 x,x (r 1 -r 2 ) = π β dr ∞ 0 du4γ 0 10 3 χ G (u) + 20 3 χ G ⊥ (u) δ(u) 2π δ(r -r 1 )δ(r -r 2 ) = 20γ 0 3β χ G (0) + 2χ G ⊥ (0) δ(r 1 -r 2 ). ( S47
)
Note that we find the same value for χ -1,1 yy and χ -1,1 zz and that the cross terms are vanishing. We calculate the susceptibility at r = 0, χ G (0), using the elementary functions defined in Eq. (S29) and taken in 0 as follows,

I 1, (0) = 1 2π 2 ∞ 0 dqq 2 w /λ 2 D + q 2 ( w /λ 2
D + q 2 )(K + κ l q 2 + αq 4 ) + q 2 , (S48)

I 1,⊥ (r c ) = 1 2π 2 2π/rc 0 dq q 2 K + κ t q 2 = 1 πKλ 2 t r c (S49) I 2, (0) = I 1, (0) 3 , I 2,⊥ (0) = I 1,⊥ (r c ) 3 (S50)
where we have introduced a cutoff length r c to remove the divergence for I 1/2,⊥ in r = 0. The Gaussian susceptibility tensor reads,

χ G (r c ) = I 2, (0) + I 1,⊥ (r c ) -I 2,⊥ (r c ), χ G ⊥ (r c ) = 1 2 I 1, (0) -I 2, (0) + I 1,⊥ (r c ) -I 2,⊥ (r c ) . (S51)
Using Eq. (S47), we obtain in Fourier space, χ -1,1 (q) = δK q i q j q 2 + δK δ ij -

q i q j q 2 , with δK = 20γ 0 3β χ G (r c ) + 2χ G ⊥ (r c ) . (S52)
Finally, we get the expression for the inverse polarization susceptibility at the first order, χ -1 ij (q) = (K + δK + κ l q 2 + αq 4 )

q i q j q 2 + (K + δK + κ t q 2 ) δ ij - q i q j q 2 . ( S53
)
where we have used the expression of χ G (q) given in Eq. (S20).

S1.8. Linear dependence of χ -1,1 (q) in salt concentration c

The correction δK is now splitted into two contributions,

δK = δK w + δK c c + τ (c 2 ) (S54)
a pure water one δK w and a second one, δK c c depending on the salt concentration and expanded linearly in c.

To get explicit expression of δK w and δK c , we expand linearly in c the functions I i,x (0), i = 1, 2, x = , ⊥ given in Eq. (S48-S50), using w /λ 2 D =c × 2N a e 2 β/ 0 . We get

I i, (0) = I w i, (0) + c × I 1 i, (0) + τ (c 2 ) (S55) I i,⊥ (r c ) = I w i, (r c ), i = 1, 2 (S56)
in which the functions are split into a pure water contribution and a linear correction in c. Note that I i,⊥ functions do not depend of the salt concentration, the associated salt correction is thus vanishing. The functions for pure water obey

I w 1 (0) = 1 2π 2 ∞ 0 dqq 2 1 1 + K + κ l q 2 + αq 4 , I w 2 (0) = I w 1, (rc) 3 , (S57)
I w 1⊥ (0) = 1 2π 2 2π rc 0 dqq 2 1 K + κ t q 2 , I w 2⊥ (r c ) = I w 1,⊥ (r c ) 3 , (S58)
and the linear correction for the parallel function,

I 1 1, (0) = 1 2π 2 0 ∞ 0 dq 2N a e 2 β (1 + K + κ l q 2 + αq 4 ) 2 , I 1 2, (0) = I 1 1, (0) 3 . ( S59
)
The expression of χ w (r c ) and χ w ⊥ (r c ) are obtained by replacing I i, (0) by I w i, (0) and I i,⊥ (0) by I w i,⊥ (0) in Eq. (S47). Finally, we get:

δK w = 20γ 0 3β χ w (r c ) + χ w ⊥ (r c ) . (S60)
We set δK w to zero as it is included in the fitted parameter K. Using the expressions of δK in Eq. ( S52) and (χ G ⊥ , χ G ) in Eq. (S51), we obtain

δK c = 20γ 0 3β I 1, (0). (S61)
We perform the integral in Eq. (S59) and obtain

I 1 1, (0) = βN a e 2 ( w -1) 2 64π 0 w (4π 2 λ 2 d + λ 2 o )(4π 2 λ 2 d + 5λ 2 o ) λ d λ 4 o (S62)
after having expressed K, κ l , α as functions of w , λ o , λ d by inverting Eq. (S5). We deduce that,

δK c = γ 5N a e 2 ( w -1) 2 24π w (4π 2 λ 2 d + λ 2 o )(4π 2 λ 2 d + 5λ 2 o ) λ d λ 4 o (S63)
Expanding linearly the relation (c) = 1 + 1/(K + δK c c ) , we obtain for the permittivity of the electrolytes

(c) = w - δK c K 2 c. (S64)
The inverse susceptibility can thus be written as:

χ -1 (q) = χ -1,G (q) + χ -1,1 (q) = (K + δK c c + κ l q 2 + αq 4 ) q i q j q 2 + (K + δK c c + κ t q 2 ) δ ij - q i q j q 2 . (S65)
that is inverted to get the one-loop corrected susceptibility χ (q) = 1 1 + K + δK c c + κ l q 2 + αq 4 (S66) which we plot for increasing salt concentration in Fig. S4. One sees comparing the susceptibilities obtained for the Gaussian model plotted in Fig. S2 and the one-loop expanded one plotted in Fig ( S5) that the enhancement of the pseudo-resonant peak at q=3 Å-1 is attenuated but not canceled by the one-loop correction.

S2. SUPPLEMENTARY INFORMATION FOR MOLECULAR DYNAMIC SIMULATION

a. Molecular dynamics simulations We simulate a cubic water box of side size L=6.5 nm composed of N w water molecules, N w going from 9033 to 8527 for increasing salt concentration. See a snapshot of the simulated system in Fig. S5. The 0.15 mol.l -1 solution contains 25 ion pairs, the 0.75 mol.l -1 solution contains 124 ion pairs and the 1.5 mol.l -1 solution, 248 ion pairs. Simulations are performed using GROMACS 2021 molecular dynamics simulation package ( [7]), the integration time steps is set to ∆t=2 fs. Simulation boxes are periodically replicated in all directions and long range electrostatics are handled using the smooth particle mesh Ewald (SPME) technique. Lennard-Jones interactions are cut off at a distance r cut =0.9 nm. A potential shift is used at the cut-off distance. All systems are coupled to a heat bath at 300 K using v-rescale thermostat with a time constant of 0.5 ps. We use MDAnanlysis to treat the trajectories. After creating the simulation box, we perform a first energy minimization. We equilibrate the system in the NVT ensemble for 200 ps, and afterwards in the NPT ensemble for another 200 ps using a Berendsen barostat at 1 bar. Production runs are performed in the NVT ensemble for 20 ns.

We performe simulations with TIP4p/ ([2]), a 4 interaction site, three point-charges and one Lennard Jones reference site model. The Lenard-Jones (LJ) center is placed on the oxygen. Charges are placed on the hydrogen atoms and on an additional interaction site, M, carrying the negative charge. The ions (Na + and Cl -) were treated can be written as a sum over the molecular polarization p j (q) of the molecule j which reads as p j (q) = 1 √ V Σ α ez α δr αj iq • δr αj 1e -iq•δrαj (S76)

with δr αj the distance between the charge α and the center of mass of the molecule and V the volume of the simulation box. We then take the transverse part of the polarization P ⊥ (q) = q × P(q)/q and define the transverse susceptibility as χ ⊥ (q) = P ⊥ (q) • P ⊥ (-q) 2k B T 0 .

(S77)

Note that we replace 1e -iq•δrαj /iq • δr αj by 1 for q • δr αj < 10 -5 to prevent numerical errors.

FIG. 2 .

 2 FIG.2. Dielectric susceptibility for electrolytes. (a) Gaussian longitudinal susceptibility χ G -Eq. (8) -as a function of q for water and solutions of increasing salt concentration. The arrow indicates the increase of the peak maximum with an increase of the concentration. (b) Transverse susceptibility χ G ⊥ = χ w ⊥ given by Eq (3) as a function of q, which is identical for water and electrolytic solutions. The parameter values of the FT model given in Eq. (2) are K= 1/76, κ l = -0.218 Å-2 , α=0.012 Å-4 and κt=0.013 Å-2 .

- 1 FIG. 3 .

 13 FIG.3. Comparison of FT derived and MD simulated response functions of electrolytes. (a) Bulk permittivity as a function of c. The field theory expression is given in Eq. (11) and plotted for γ=0.08 Å4 .e -2 , (b) Transverse permittivity -FT expression given in Eq[START_REF] Hasted | Dielectric properties of aqueous ionic solutions. part i and ii[END_REF], (c) MD simulated longitudinal susceptibility χ MD , (d) conjectured expression of the longitudinal susceptibility χ FT given in Eq.[START_REF] Levy | Dielectric constant of ionic solutions: A field-theory approach[END_REF]. Susceptibilities are plotted as function of q for pure water and electrolytes of increasing concentration.

  FIG. S1: MD simulated and FT derived susceptibilities for bulk water. (a) Longitudinal and (b) transverse susceptibilities in Fourier space from FT model Eq. (S4) and from MD simulations (details given in S2). The parameter values of the FT model are K= 1/76, κ l = -0.218 Å-2 , α=0.012 Å-4 and κt=0.013 Å-2 .

- 1 FIG

 1 FIG.S2: Longitudinal susceptibility of electrolytes. The Gaussian susceptibility χ G (q) given in Eq. (S23), for increasing concentration. Parameters are given in caption of Fig.S1.

FIG. S3 :

 S3 FIG.S3: Illustration of the basis considered to calculate the susceptibility tensor χ G .

- 1 FIG

 1 FIG.S4: One-loop expanded longitudinal susceptibility of electrolytes. The susceptibility χ (q) given in Eq. (S66), for increasing concentration with δKc=0.028 mol -1 .L. Other parameters are given in caption of Fig.S1.

We consider a field theory (FT) model describing the dielectric properties of water at the nanoscale. The electrostatic energy of the system is written as a functional of the polarization P as follows [START_REF] Maggs | Simulating nanoscale dielectric response[END_REF]:

The first term corresponds to the bare Coulomb interactions between the partial charges -∇ • P(r) of the fluid and the second term U G conf -where the index G stands for Gaussian -is a phenomenological configuration energy. It can S1.3. Mean field (ψ, P) for the linear and nonlinear models

In this subsection, we derive (ψ, P ), the fields minimizing the action F G u given in Eq. (S11), and the action F u obtained by replacing U G conf in Eq. (S11) by the non Gaussian configuration energy,

The mean fields obey the following equations:

The functional derivative with respect to Ψ gives:

The functional derivative with respect to P i leads to

or the Gaussian action, F G u , and to

for the non Gaussian one, F u :

In the two cases, we obtain

The mean fields are vanishing both for Gaussian and non Gaussian configuration energy.

S1.4. Gaussian susceptibility for an electrolyte S1.6. One-loop free energy expansion

The partition function can be expanded to the second order around the mean field point (P , ψ), as

where we have droped the prefactor. F

u (rr ) is the second functional derivative of the action and can be defined as follows:

u (rr ) = δ 2 Fu(P ,ψ) δPi(r)δPj (r ) δ 2 Fu(P ,ψ) δPi(r)δΨ(r ) δ 2 Fu(P ,ψ) δΨ(r)δPi(r )

We set 2γP 2 0 = K so that the Gaussian asymptotic behavior of the P 4 -model, reached in the low polarization limit, is similar to the Gaussian model. The Gibbs free energy at the first order with a loop expansion can be written as:

S1.7. One-loop correction for the inverse susceptibility

In this subsection, we expand the inverse polarization susceptibility to the first order, χ -1 ≈ χ -1,G + χ -1,1 . Using its expression as a function of the Gibbs free energy of the system,

and the first order expansion of F, (Eq. S39), we get:

The second matrix product in the right hand term is vanishing as δF

(2)

u /δP x (r 2 )(P , ψ) = 0. We thus get:

We use the expression of F

(2) u given in Eq. (S37) and obtain

according to the force field developed in the reference [8]. b. Statistical treatment For the longitudinal and transverse susceptibility, the error bars are derived following the reblocking method [9]. For the bulk permittivity, we cut the trajectory in 5 statistically independent blocks, compute the bulk permittivity of each block, estimate the sample variance σ 2 and define the error bar as σ 2 /5. c. Permittivity Bulk permittivity is calculated from the total system dipole moments [10] M according to:

as implemented in the GROMACS dipoles module. M is the is the volume integral of the polarization as M =

V drP(r) and V is the volume of the box. Note that the ion polarization is not taken into account. d. Susceptibilities To compute the q-dependent susceptibilities, we use the fluctuation-dissipation theorem, relating the response functions to the polarization fluctuations as follows:

One can express the longitudinal susceptibility as a function of the charge structure factor S(q) and gets

The charge structure factor in the Fourier space can be decomposed into a intramolecular and intermolecular part, S(q) = S int (q) + S inter (q) (S70)

with S inter (q) the intermolecular contribution

z is the valency, e being the elementary charge, n w the molecular number density. h IJ is the Fourier transform of g IJ (r) -1, g IJ (r) being the radial distribution function associated with the atoms couple IJ. The intramolecular contribution can be written as

sin(qd HH ) qd HH -4 sin(qd HM )

where d IJ is the intramolecular distance between atom I and J. At low q the precision of this expression of the structure factor become pretty low as the function h IJ (r) is obtained on a finite range imposed by the box size. To solve this problem, we proceed as follows. For q < 2.5 Å-1 , we take into account the periodicity of the system, calculate the charge structure factor for discretized values of the wave length q, q = 2π/L n 2

x + n 2 y + n 2 z . We compute directly the charge structure factor from the charge distribution ρ(q) in the Fourier space, ρ(q) = N w i=1 eze iq•r -2e -iq•rM,i + e -iq•rH1,i + e -iq•rH2,i (S73)

where H 1,i and H 2,i stand for the two hydrogens of the molecule j. the charge structure factor S(q) = ρ(q)ρ(-q) /V , S(q) = 2q 2 H V i,j,j≤i 4 cos(q • d OiOj ) -2 cos(q • d OiH1j ) -2 cos(q • d H1iOj )

-2 cos(q • d OiH2j ) -2 cos(q • d H2iOj ) + cos(q • d H1iH1j ) + cos(q • d H2iH2j )

+ cos(q • d H1iH2j ) + cos(q

where q is a vector and d AiAj stands for r Air Aj . The transverse susceptibility is performed following ( [11]). The polarization of the medium in the Fourier space P(q) = Σ j p j (q)e -iq•rj (S75)


Y Hoarau

M Braza

G Harran

E Longatte

F Baj

T Marcel

Numerical simulation of turbulent flow physic in a tube array

The prediction of fluid-elastic instabilities developing in a tube bundle is of major importance in the design of modern heat exchangers for nuclear reactor cooling to prevent accidents associated with fluid-elastic instabilities leading to flutter, material fatigue, shocks between beams and damage of the solid walls. The fluid-elastic instabilities in tube arrays appear in the laminar regime and persist in the high Reynolds number turbulent flow. They are governed by non-linear interaction between fluid-elastic instability and turbulence. It is of major importance for the design to accurately predict the unsteady loads under turbulent flow and the amplification of the fluid-elastic instabilities. This work focuses on the numerical simulation of turbulent flow in tube bundle and prediction of fluid-elastic instabilities.

Introduction

There is a considerable effort in the state-of-the-art to improve methodologies able to capture the fluid-structure interaction under turbulent flow, and to assess critical thresholds. These aspects are challenging for the research applications involving the nuclear power industry. The present work is part of the BARESAFE project funded by the French National Research Agency (ANR) with LAMSID, Y. Hoarau ICUBE, UDS -CNRS -ENGEES -INSA, France, e-mail: hoarau@unistra.fr M. Braza and G. Harran and T. Marcel Institut de Mécanique des Fluides de Toulouse, UMR CNRS 5502, France E. Longatte Laboratoire de Mécanique des Structures Industrielles Durables, UMR EDF CNRS 2832, France F. Baj Commissariat à l'Energie Atomique CEA / Saclay, France 1 EDF, CEA, IMFT, ICUBE, IRIT and AREVA. The BARESAFE project (Numerical modeling of multi-scale fluid solid coupled systems : advanced hybrid methods for reduction of uncertainty on stability limits and optimization of nuclear safety barrier reliability) is a research project involved in the development of modeling and advanced numerical methods for simulation of large size systems involving multi-physics in the field of mechanics. It addresses the hard issue of stability analysis of dynamical systems submitted to external turbulent flows and aims to establish accurate stability maps applicable for heat exchanger design. The purpose is to provide stability limit dimensionless modeling suitable in a variety of configurations with a maximal accuracy in spite of the large scale of the systems to be considered. The challenge lies in predicting local effects possibly impacting global systems. Therefore the combination of several strategies convenient simultaneously for multi-physics, multi-scale and large size system computation is required. Based on empirical concepts, the heuristic models currently used in the framework of standard stability analysis suffer from a lack of predictive capabilities [Blevins(1974), Blevins(1979), [START_REF] Tanaka | Fluid elastic vibration of tube array in cross flow[END_REF], [START_REF] Tanaka | Flow induced vibration of tube arrays with various pitch-to-diameter ratios[END_REF], Chen(1983), [START_REF] Paidoussis | The mechanisms underlying flow-induced instabilities of cylinder arrays in crossflow[END_REF]], [START_REF] Lever | A theoretical model for fluid-elastic instability in heat exchanger tube bundles[END_REF]]. On the other side, numerical approaches based on fullycoupled fluid solid dynamics system computation remain expensive due to the multiphysics patterns of physics and the large number of degrees of freedom to be involved. In this context, since experimentation is not reachable and numerical simulation is unavoidable but prohibitive, the project proposes an hybrid strategy in order to take advantages from both numerical local and empiric global solutions.

The present study focuses on the improvement of the numerical simulation of the fluid-structure interaction in the cylinders arrays, especially in the parametric range favoring the amplification and negative damping of flutter instability. The main objective is to analyze the physics of the appearance of Movement Induced Vibrations in connection with the instabilities coming from the fluid part in an array bundle. Concerning turbulence modeling methods, among URANS approaches, the Organised Eddy Simulation [START_REF] Braza | Turbulence properties in the cylinder wake at high reynolds number[END_REF], [START_REF] Bourguet | Anisotropic organised eddy simulation for the prediction of non-equilibrium turbulent flows around bodies[END_REF]], OES has been used to capture the coherent flow patterns and non-equilibrium effects associated with the fluid-elastic instabilities. Furthermore, the OES version of the Delayed Detached Eddy Simulation and the Scale Adaptive Simulation has been used for the 3-D simulations.

First we focused on a two-dimensional study in static configuration, in order to assess the predictive ability of 2-D mechanisms. Secondly, the spontaneous motion of a central cylinder, in SDOF (single degree of freedom) vertical displacement is studied, below and beyond critical values for the appearance of flutter instability in the array of cylinders. Then the three-dimensional flow structure is studied by OES-DDES [START_REF] Skopek | Hybrid rans-les modeling of a strongly detached turbulent flow around a tandem cylinders configuration[END_REF]] and SAS [START_REF] Menter | A scale-adaptive simulation model using two-equation models[END_REF], Menter and Bender(2003)].

Only the static approach is presented in this paper.

Problem set-up

The configuration fit an experimental set-up (DIVA) designed by the CEA where 20 cylinders are arranged in 5 columns with a pitch of 1.5D, D being the diameter of the cylinders. The Reynolds number is 20.000 which corresponds to an intertube Reynolds number of 60.000. The grid (Fig. 1) is composed with 392 600 cells distributed on 80 structured blocs. On the entry and on the outlet boundary, the conditions are :

• inlet : far field boundary condition using characteristic variables,

• outlet : subsonic outflow with imposed static pressure based on Riemann invariants. • top-bottom : no-slip walls.

Numerical modeling

The simulations are performed with our parallelized Navier-Stokes Multi-Block (NSMB) solver [START_REF] Vos | Recent advances in aerodynamics inside the NSMB (Navier-Stokes Multi-Block) consortium[END_REF]]. NSMB solves the compressible Navier-Stokes equations using a finite volume formulation on Multi-Block structured grids. Various spatial discretization schemes are available like Jameson's central difference, Roe or AUSM+. The time integration is based on the full matrix implicit LU-SGS (Lower-Upper Symmetric Gauss-Seidel) method [START_REF] Yoon | A Multigrid LU-SSOR Scheme for Approximate Newton Iteration Applied to the Euler Equations[END_REF]] and on the dual-time stepping. In the present work, the artificial compressibility method [Chorin(1968)] is employed. NSMB is parallelized using the Message Passing Interface.

Turbulence modeling

The URANS turbulence models used in this study are the [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF]], the kω models (baseline and SST) of [Menter(1994)], the OES kε and the OES kω baseline models, described in [START_REF] Braza | Turbulence properties in the cylinder wake at high reynolds number[END_REF], [START_REF] Bourguet | Anisotropic organised eddy simulation for the prediction of non-equilibrium turbulent flows around bodies[END_REF]]. The Delayed Detached Eddy Simulation, DDES is employed by using for the RANS part the turbulence length scale from the OES method. The association of the OES modeling in the RANS part of the DES/DDES is reported in a number of studies of the present research group, that have been the object of applications in a variety of test cases in the European research programs DESIDER (Detached Eddy Simulation for Industrial Aerodynamics) and ATAAC (Advanced Turbulence Modeling for Aeronautics Applications Challenges) in aeronautics. The introduction of a turbulence length scale from OES in the Delayed DES, DDES [START_REF] Spalart | A new version of detached-eddy simulation, resistant to ambiguous grid densities[END_REF]] is reported in [Haase and Revell(2009)], [START_REF] Bourguet | Anisotropic organised eddy simulation for the prediction of non-equilibrium turbulent flows around bodies[END_REF]] and in [START_REF] Skopek | Hybrid rans-les modeling of a strongly detached turbulent flow around a tandem cylinders configuration[END_REF]]. The use of this variety of turbulence models allows assessment of their predictive capabilities in respect of the strongly detached unsteady flow character in the array of cylinders.

The OES modeling approach is used to capture the two aspects of fluid structure instability : organized and chaotic. In this approach, the turbulent spectrum is decomposed in a first part regrouping all the coherent structures (resolved part) and in a second part regrouping all the chaotic processes independently on their size (spectrum to be modeled), as presented schematically in Fig. 2. It is recalled that in LES the distinction is done according to the size of the structures which limits its application to moderate Reynolds numbers. The fact that part of the spectrum to be modeled in OES extend from the low to the high frequencies allows the use of statistical turbulence modeling with appropriate modification due the non-equilibrium. In the time-domain, the equations are the phase-averaged NavierStokes equations. Due to the nonlinear interaction between the coherent part and the non-coherent one, there is a slope modification of the inertial part in the spectrum, in the vicinity of the peak. This yields to a reconsideration of the eddy-diffusion coefficient for the class of the two-equation modeling, as well as improved damping functions to attenuate turbulence towards the wall.

Results

The unsteady loads obtained with this four turbulence models are examined and their average values, their minimum and maximum values (Table 1) their frequency distribution (Table 2 We have a good match between the turbulent simulations and the experimental data. The kε-OES and Spalart-Allmaras approaches give the best prediction of the RMS values of the lift coefficient. However, for the extrema values, the kω-BL-OES prediction brings the best results.

Extrema values of the lift coefficient are also important in a study on the fluidstructure interaction. Indeed, it is possible that not only the signal strength (RMS values) but also the extent of the values will trigger the instability. This is why the prediction of unsteady loads have to be accurate for these two aspects. In this context, OES approaches seem to bring better results than conventional URANS approaches.

This can be explained by the fact that this complex flow configuration of the tube bundle at this Reynolds number caused a mixed boundary layer. OES approach allow to have a better capture of the mixed boundary layer than the classical URANS approaches.

Concerning the frequencies, results are presented on Power Spectral Density (PSD) on Fig. 3. We can find also a good match between models and experimental data. This is the first time to our knowledge that unsteady numerical simulations have so good agreement with this test bench of the configuration of a tube bundle under cross flow at a Reynolds number of 20.000.

OES approaches seem to predict a more accurate power level than the URANS classical ones which slightly underestimate this level.

Three predominant frequencies appear on the spectra. After a study of the cartography of separation point, we have linked the values of these frequencies to known physical phenomena : the Von Kármán vortex shedding St V K , the confinement phenomenon St C and the shear layer phenomenon St SL .

We can define a Strouhal number with the free stream velocity, but also a Strouhal number with inter-tube velocity, which we denote St •,it St •,it = f du it = (P * -1P * ) St with P * = Pd

(1)

with P * the reduced step. Values of predominant frequencies are summarized in Tab. 2. The inter-tube value observed for the Von Kármán Strouhal number is nearly that the conventional value known for a Reynolds number of 60.000 for one cylinder (St cyl. ≈ 0.2). Furthermore the four models predict almost the same values for the three predominant frequencies.

As for the study of the cartography of the separation point, we focused only on the central cylinder and we looked at the position of the point of detachment on the wall of the cylinder for a series of snapshots. We found a periodical appearance of grounds which correspond to the predominant frequencies of the spectral study. The frequency which is linked to the Von Kármán vortex shedding phenomenon, whose value is rather low (relative to the other two), is located to the rear of the cylinder.

One that matches the shear layer phenomenon is located just above this area with a higher value. Finally, the third frequency is linked to a phenomenon related to the impact of coherent structures on the zone |α| ≈ 60 • . This zone corresponds to the area where the structures issued from the previous rows of cylinders are crashing. That's why this predominant frequency is linked to a confinement phenomenon, caused by the geometry of the tube bundle.

Conclusions

The present study shows the prediction of unsteady loads in a tube bundle at high Reynolds number in the context of fluid-structure interaction. URANS and OES turbulence modeling comparison on the static configuration is performed. The oscillatory characteristics is well predict, particularly concerning the power spectrum density. OES modeling improves the coherent structure prediction in the cylinder wake, which is under-predicted with the classical URANS approach.

There is a good match with the experimental data, even if the simulation have been computed in two dimensions. The 3-D simulations are in progress in order to find the purely 3-D phenomenon.

There are three phenomena which involve three predominant frequencies in the tube bundle under cross flow at Reynolds number of 20.000: the Von Kármán vortex shedding, the confinement and the shear layer. The demonstration of the existence of this three predominant frequencies is a first step to the numerical study of fluidelastic instability.

The classical URANS and OES approaches are good tools to compute simulations at this Reynolds number. DNS is too expensive in term of computational cost. So there is an obvious interest to use hybrid approaches for 3-D simulation.

Fig. 1

 1 Fig. 1 Grid overview

Fig. 2

 2 Fig.2Dual spectrum splitting : the distinction between the structures to be resolved and those to be modeled is based upon their organized or random character. Part (2) of the non-equilibrium energy spectrum has to be modeled by reconsidering advanced statistical turbulence modeling efficient in high-Re wall flows, due to the inertial-range modulation from equilibrium turbulence, schematically shown on the right.

Fig. 3

 3 Fig. 3 Power spectrum density of the lift coefficient for the models SA, kω-SST, kω-Baseline-OES and kε-OES closures with hanning windows

Table 1

 1 ) are compared. Lift coefficient of unsteady loads

	Spalart-Allmaras 0.44 -0.89 1.04
	k -ω-SST	0.19 -0.68 0.58
	k -ω-BL-OES 0.65 -1.84 1.69
	k -ε-OES	0.46 -1.08 1.11
	DIVA Experiment 0.45 -1.86 2.33

Models

C L,RMS C L,min. C L,max.

Table 2

 2 Values of predominant frequencies

	Models	St V K St V K,it St C St SL
	Spalart-Allmaras 0.34 0.11 0.99 2.25
	k -ω-SST	0.36 0.12 0.99
	k -ω-OES	0.39 0.13 1.01 2.25
	k -ε-OES	0.39 0.13 1.01 2.25

Acknowledgements This study was made using high performance calcul possibility of nationals centers as CINES, TGCC/CCRT and the computing center of the Strasbourg University. This work was conducted in the context of the BARESAFE project (http://www-imfs.ustrasbg.fr/baresafe/index.html), ANR 11-MONU-004.