Wetting transition of ionic liquids at metal surfaces: A computational molecular approach to electronic screening using a virtual Thomas–Fermi ﬂuid

Of particular relevance to energy storage, electrochemistry and catalysis, ionic and dipolar liquids display a wealth of unexpected fundamental behaviors — in particular in conﬁnement. Beyond now well-documented adsorption, overscreening and crowding effects 1–3 , recent experiments have highlighted novel phenomena such as unconventional screening 4 and the impact of the electronic nature — metallic versus insulating — of the conﬁning surface on wetting/phase transitions 5,6 . Such behaviors, which challenge existing theoretical and numerical modeling frameworks, point to the need for new powerful tools to embrace the properties of conﬁned ionic/dipolar liquids. Here, we introduce a novel atom-scale approach which allows for a versatile description of electronic screening while capturing all molecular aspects inherent to molecular ﬂuids in nanoconﬁned/interfacial environments. While state of the art molecular simulation strategies only consider perfect metal or insulator surfaces, we build on the Thomas–Fermi formalism for electronic screening to develop an effective approach that allows dealing with any imperfect metal between these asymptotes. The core of our approach is to describe electrostatic interactions within the metal through the behavior of a ‘virtual’ Thomas–Fermi ﬂuid of charged particles, whose Debye length sets the Thomas-Fermi screening length λ in the metal. This easy-to-implement molecular method

while describing very accurately the capacitance behavior -and hence the electrochemical properties -of the metallic confining medium. By applying this strategy to a nanoconfined ionic liquid, we demonstrate an unprecedented wetting transition upon switching the confining medium from insulating to metallic. This novel approach provides a powerful framework to predict the unsual behavior of ionic liquids, in particular inside nanoporous metallic structures, with direct applications for energy storage and electrochemistry.
The fluid/solid interface as encountered in confined liquids is the locus of a broad spectrum of microscopic phenomena such as molecular adsorption, chemical reaction, and interfacial slippage. 7 These molecular mechanisms are key to nanotechnologies where the fluid/solid interaction specificities are harnessed for energy storage, catalysis, lubrication, depollution, etc. From a fundamental viewpoint, the behavior of nanoconfined fluids often challenges existing frameworks even when simple liquids are considered. Ionic systems, either in their liquid or solid state, between charged or neutral surfaces lead to additional ion adsorption, crowding/overscreening, surface transition, and chemical phenomena that are crucial in electrokinetics (e.g. electrowetting) and electrochemistry (e.g. supercapacitors/batteries). 3 Theoretical descriptions of nanoconfined fluids -except rare contributions [8][9][10][11][12][13][14][15] -assume either perfectly metallic or insulating confining surfaces but these asymptotic limits do not fully reflect real materials as they display an intermediate imperfect metal/insulator behavior (only few metals behave perfectly and all insulators are semiconducting to some extent). Yet, the electrostatic boundary condition imposed by the surrounding medium strongly impacts confined dipolar and, even more, charged systems. [16][17][18] For instance, the confinement-induced shift in the freezing of an ionic liquid was found to drastically depend on the surface metallic/insulating nature. 5 Formally, quantum effects leading to electronic screening in the confining metallic walls can be accounted for using the microscopic Thomas-Fermi (TF) model. 5,19 This formalism relies on a local density approximation for the charges in the metal which are treated as a free electron gas (therefore, restricting the electron energy to its kinetic contribution). This simple yet robust model allows considering any real metal -from perfect metal to insulator -through the Thomas-Fermi screening length λ. The latter is defined in terms of the electronic density of state of the metal at the Fermi level D(E F ) according to λ = ε/e 2 D(E F ) (ε is the dielectric constant and e the elementary charge); the Fermi energy is directly related to the free electron density n 0 where m e is the electron mass and = h/2π the Planck constant, see Supplementary Information. Despite this available framework, the development of classical molecular simulation methods to understand the microscopic behavior of classical fluids in contact with imperfect metals is only nascent. While insulators are treated using solid atoms with constant charge, metals must be described using an effective screening approach. The charge image concept can be used for perfectly metallic and planar surfaces 20 but refined strategies must be implemented for non-planar surfaces such as a variational 21,22 or Gauss law [23][24][25] approaches to model the induced charge distribution in the metal. A recent proposal 9 builds on our TF framework 8 to propose a computational approach based on variational localized surface charges that accounts for electrostatic interactions close to imperfect metals.
Here, we develop an effective yet robust atom-scale simulation approach which allows considering the confinement of dipolar or charged fluids between metallic surfaces of any geometry and/or electronic screening length. Following Torrie and Valleau's work for electrolyte interfaces, 26 the electronic screening in the imperfect confining metal is accounted for through the response of a high temperature virtual Thomas-Fermi fluid made up of light charged particles. Due to its very fast response, this effective Thomas-Fermi fluid mimics metal induction within the confining surfaces upon sampling the confined system configurations using Monte Carlo or molecular dynamics simulations. After straightforward implementation in existing simulation packages, this strategy provides a mean to impose an effective electronic Thomas-Fermi screening length that is directly linked to the equivalent virtual fluid Debye length. By adopting a molecular level description, our approach captures all specificities inherent to confined and vicinal fluids at the nanoscale (e.g. density layering, slippage, non-viscous effects). This virtual Thomas-Fermi model correctly captures electrostatic screening within the confined system upon varying the confining host from perfect metal to insulator conditions. This novel molecular approach is also shown to accurately mimic the expected capacitive behavior, therefore opening perspectives for the atom-scale simulation of electrochemical devices involving metals of various screening lengths/geometries. Last but not least, using this novel method, by considering an ionic liquid between two parallel planes, we unravel a continuous wetting transition as the surfaces are tuned from insulating (non-wetting) to metallic (wetting).

Interaction at Thomas-Fermi metal interfaces
where the superscripts C and I refer to the physical charges in the dielectric medium and induced charges within the metal, respectively. As shown in Fig. 1(b), U CC is the Coulomb interaction energy between the charges i and j while U II λ is the interaction energy between the charge densities ρ I i and ρ I j induced in the metal by these two charges. For each ion i, its interaction energy U CI λ with the metal decomposes into a one-body contribution U CI • λ (z i ) -corresponding to the interaction with its image in the metal -and two-body contributions U CI • λ (z i , z j , R ij ) -corresponding to the interaction with the induced charges due to all other charges j. Analytical expressions exist for U CI • λ and U CI • λ 8, 10-12 but U II λ must be estimated numerically from the energy density, i.e. U II λ = drΨ β i (r)ρ I j (r) + Ψ β j (r)ρ I i (r), where Ψ β k and ρ I k are the electrostatic potential and induced charge density in the metal due to the point charge k = i, j. All details are given in the Supplementary Information.

Effective molecular simulation approach
Except for the usual Coulomb energy CC, formal expressions for the CI and II energies cannot be implemented in molecular simulation due to their complexity. In particular, U II λ requires expensive integration on the fly as analytical treatment for imperfect metals is only available in closed forms in asymptotic limits. 8,15 Here we model the resulting complex electrostatic interactions between the ions of the liquid thanks to a 'virtual Thomas-Fermi fluid' located within the confining solids, the Thomas-Fermi wave-vector, and D(E F ) the density of states at the Fermi level (see Supplementary Information). Combined with Poisson equation this leads to the Helmholtz equation for TF screening, ∇ 2 Ψ II = k 2 TF Ψ II , which indeed resembles the Debye-Hückel equation for electrolyte solutions. Accordingly, one can simulate the imperfect metal using a system of virtual (classical) charged particles of charge q TF and mass m TF , with density ρ TF and temperature T TF . The analogous TF screening length λ = k −1 TF can be identified as the equivalent Debye length λ D : Hence, by considering the dynamics of these light ions located in the confining solid, any screen-ing length λ between 0 (perfect metal) and ∞ (insulator) can be efficiently mimicked depending on q TF , ρ TF , and T TF . This virtual system allows simulating properly the complex electrostatic interactions within the ionic liquid in the vicinity of an imperfect metal.
Equation (2) shows that mapping the fluid of mobile charges onto the TF model only requires to set a single parameter ρ TF q 2 TF /T TF (fixing ε β = 1). In our molecular dynamics approach, to ensure that the particles in this effective Thomas-Fermi fluid relax fast, their mass/temperature are chosen much smaller/larger than their counterpart in the confined system; typically m TF ∼ 0.01m and T TF ∼ 10T (requiring typical integration steps of 0.1 fs and 1 fs, respectively). In practice, as shown in Fig. 2(a), the effective simulation strategy consists of sandwiching the charged or dipolar system between two metallic media separated by a distance d w . The confining media of width d TF are filled with the Thomas-Fermi fluid having a density ρ TF . Once ρ TF and T TF are set, λ is varied by tuning q TF according to Eq. (2); from q TF = 0 (λ → ∞) for an insulator to q TF = 1 (λ = 0.03 nm) for a nearly perfect metal. All simulations reported in this article are carried out using molecular dynamics but they could be easily implemented into Monte Carlo algorithms to perform calculations in other ensembles (Grand Canonical and isothermal/isobaric ensembles for instance). All simulation details are provided in the Methods section.

2D crystal at metallic interfaces
To validate our effective approach, we consider a 2D square crystal of lattice constant a = 1.475 nm made up of charges ±1 e and located at a distance d from a metal (Fig. 2). Due to the periodic boundary conditions, a second pore/metal interface is present at a distance d w = 20 nm. Yet, as shown in the Supplementary Information, this second interface does not affect the electrostatic energy as d w is large enough. In the Thomas-Fermi framework, the charge density ρ I at a position r in the metal induced by a charge q located in (0, 0, d) reads (see Supplementary Information): where R = [x 2 + y 2 ] 1/2 is the lateral distance to the charge q, J 0 is Bessel function of the first kind, and κ 2 = K 2 + k 2 TF . Fig. 2(c,d) shows the induced charge density ρ I (d, r) as obtained by  Supplementary Information, Eq. (3) must be summed over all crystal periodic images but it was found that the sum converges quickly). For comparison, Fig. 2(e,f) shows ρ I (d, r) as obtained using our effective approach from the local charge density in the metal, i.e. ρ I = e(ρ + TF − ρ − TF ). In contrast to ρ I (d, r) in the Thomas-Fermi model, due to their finite size, the fluid charges in the simulation cannot approach arbitrarily close to the metal/pore surface. For consistency, the analytical/simulation data were compared by defining z = 0 in the simulation as the position where the Thomas-Fermi fluid density becomes non-zero. Fig. 2 shows that the effective molecular simulation qualitatively captures the predicted density distribution induced in the metal. Each physical charge in the 2D crystal induces in the metal a diffuse charge distribution of opposite sign. Moreover, as expected from the Thomas-Fermi framework, the induced charge distribution in the effective simulation decays over the typical length λ.  Figure 3 compares the total energy U λ as a function of the distance d with the numerically evaluated prediction of the Thomas-Fermi model in Eq. (1). As expected theoretically, the overall energy decays with decreasing λ between boundaries for an insulator (λ → ∞) and a perfect metal (λ → 0). As shown in Fig. 3, our effective approach captures quantitatively the screening behavior of the confining medium assuming a screening length λ = c 0 +c 1 λ+c 2 λ 2 (with λ the ion gas Debye length, c 0 = 0.22 nm, c 1 = 0.91 and c 2 = 0.28 nm −1 in our system). Such values do not simply correspond to fitting parameters that allow matching the simulated and theoretical energies; as explained in the next paragraph, they were derived so that the capacitance of the virtual Thomas-Fermi fluid matches the theoretically expected value C = ε 0 /λ . The fact that the rescaled screening length λ also allows recovering the expected screened interaction energy further supports the physical validity of our effective molecular approach. Moreover, physically, the parameters c 0 , c 1 , and c 2 are not just empirical parameters as they account for the following effects in the screening fluid used in the simulation: c 0 accounts for the finite size σ of the Thomas-Fermi ions which prevents reaching screening λ ≤ σ. This is supported by the fact that c 0 ∼ σ corresponds to the value below which the repulsive interaction potential in the Thomas-Fermi fluid becomes larger than k B T . c 1 arises from the non-ideal behavior of the effective Thomas-Fermi fluid which leads to overcreening compared to an ideal gas having the same charge density ρ TF (c 1 = 1 corresponds to the ideal behavior); c 2 = 0 indicates non-linear effects in electrostatic screening which go beyond the linear approximation used in the Thomas-Fermi framework.

Capacitance
To further establish the validity of our novel molecular approach, an important requirement is to verify that the virtual Thomas-Fermi fluid yields the correct capacitance behavior. With this aim, as shown in Fig. 4(a), a simple molecular dynamics set-up was designed by assembling a composite system made up of two Thomas-Fermi fluids sandwiching a dielectric material of a width d w (either a vacuum layer or a molten salt was considered to verify that the overall capacitance follows the expected physical behavior). The molten salt (NaCl) is modeled using charged particles ±1e that interact via a Born-Mayer-Huggins potential 27 (details can be found in the Supplementary  Information). To prevent mixing of the Thomas-Fermi fluid/charged system, a reflective wall of thickness ξ = 0.2 nm is positioned between the two subsystems. The whole composite is placed between two electrodes having an overall charge +Q and −Q (all details can be found in the methods section). With such a geometry, the capacitance C = Q/∆Ψ is readily obtained from the potential difference ∆Ψ; as illustrated in Fig. 4(a), in our simulations, ∆Ψ is determined using Poisson equation by integrating twice the charge density profile ρ(z) = e[ρ + (z) − ρ − (z)].
The system considered here simply consists of double layer capacitors in series so that its capacitance per unit area should verify the following combination rule: where the first and second terms correspond to the capacitance of the vacuum slab of width d w and that of the Thomas-Fermi fluid (the factor 2 simply accounts for the presence of two Thomas-Fermi/vacuum interfaces). As shown in Fig. 4(b), the simulation data are in reasonable agreement with the prediction from this simple expression with deviations increasing with λ. Interestingly, as shown in the insert in Fig. 4(b), our effective approach captures quantitatively the expected capacitance behavior of the confining medium upon rescaling λ → λ = c 0 + c 1 λ + c 2 λ 2 (see discussion above). As another consistency check, the vacuum layer in the capacitor was replaced by a slab of molten salt -see molecular configuration shown in Fig. 4(a). As expected, upon inserting such a molten salt, the effective capacitance C drastically increases (i.e. the inverse capacitance shown in Fig. 4(b) decreases). More importantly, as shown in Fig. 4(c), the induced capacitance change ∆1/C observed in our simulation data follow the expected behavior with a λ-independent value: where ε is the permittivity of the molten salt. Furthermore, since ε ε 0 , we predict that ∆1/C ∼ d w /ε 0 in very good agreement with the simulation data shown as green circles in Fig. 4(c) (the small deviation is due to the fact that the vacuum permittivity is not completely negligible compared that of the molten salt).
Wetting transition salt (blue symbols) confined between the TF fluids. As expected, for ε = ε α ε 0 ε 0 , the difference between both systems (green symbols) is close to d w /ε 0 .
Having assessed our effective simulation strategy, we now turn to the thermodynamically relevant case of the wetting of an ionic liquid at metal surfaces (as described above, the ionic liquid is taken as a molten salt modeled using charged particles ±1e). Fig. 5(a) shows the number density profiles ρ n (z) for the salt and Thomas-Fermi fluid for different λ (which is modified by tuning q TF ). A crossover is observed upon decreasing λ; while the salt is depleted at the insulating interface, a marked ion density peak appears under metallic conditions (in contrast, the density profile for the Thomas-Fermi fluid is nearly unaffected by λ). This behavior suggests that the system undergoes a wetting transition upon changing the dielectric/metallic nature of the confining medium (perfect wetting/non-wetting for metal/insulator, respectively).
The observed wetting transition was characterized by measuring the surface tension of the liquid salt confined at a constant density within surfaces made of a metallic medium with a screening length λ via the Irving-Kirkwood formula: γ(λ) = L z /2 P N − P T where the terms in bracket are the average normal and tangential pressures, L z is the box length in the z direction and the factor 2 accounts for the two interfaces in the slit geometry. We considered the salt in its liquid (l) and gas (g) states in contact with the metal (m) and estimated for various λ the surface tension difference normalized to the gas-liquid surface tension, ∆γ(λ)/γ gl = [γ gm (λ) − γ lm (λ)]/γ gl . In practice, γ lg was assumed to correspond to the liquid-wall interface at the insulating surface, i.e. γ lg ∼ γ lw (∞) (this approximation does not affect the discussion below as γ lg is used for normalization only). As shown in Fig. 5(b), upon switching the surfaces from insulating to metallic, the confined salt undergoes a continuous transition from non-wetting or partially wetting [∆γ(λ)/γ gl ≤ 1] to perfectly wetting [∆γ(λ)/γ gl > 1]. For imperfect metals with λ 0.1 nm, the spreading parameter S = γ gm − γ lm − γ lg < 0 and the contact angle θ can be inferred from Young equation ∆γ/γ lg = cos θ (with cos θ < 0 and > 0 for non-wetting and partially wetting, respectively). For metals with λ 0.1 nm, S ≥ 0 so that the system becomes perfectly wetting with a liquid film spreading over the metal surface. As shown in the inset of Fig. 5(b), the change in ∆γ between the insulator and metal is found to scale with the liquid/gas density contrast: where ρ l ρ g was assumed in the second equality. As expected from the Thomas-Fermi model, the inset in Fig. 5(b) shows that α(λ) ∼ U CI λ as the charge interaction with the induced density distributions (including the charge image) is dominating the surface energy excess. This important finding provides a microscopic picture for recent experimental results in which capillary freezing and wetting of an ionic liquid was found to be promoted by metal surfaces. 5,6 In particular, these authors showed that the freezing point shift upon varying λ could be rationalized by assuming that the difference between the liquid/wall and crystal/wall surface tensions scales with the density difference between the crystal and liquid. 5

Discussion
We developed a classical molecular simulation strategy that allows considering the confinement within any material ranging from perfect metal to insulator. This approach, which does not require to input any given geometry/molecular structure for the confining material, describes in an effective fashion electrostatic screening within confined/vicinal fluids together with the expected capacitive behavior. After straightforward integration into an existing simulation packages, this method offers a useful framework to investigate the behavior of dipolar and charged fluids in porous materials made up of any material with imperfect dielectric/metal properties. Beyond practical implications, we also unraveled a non-wetting/wetting crossover in nanoconfined liquids as the confining sur-faces vary from insulator to perfect metal. This raises new challenging questions on the complex behavior of charged systems in the vicinity or confined within surfaces with important applications such as electrowetting/switching for energy storage, lubrication, catalysis, etc.

Methods
The corresponding force field parameters are given in Supplementary Table 1. Reflective walls of width ξ = 0.2 nm are used at each metal/dielectric interface to prevent the Thomas-Fermi fluid/charged system to migrate to the pore space/confining media. The latter implies that, if an atom moves through the wall by a distance δ in a timestep, its position is set back to −δ away from the wall and the sign of the corresponding velocity component is flipped.
In all simulations presented in the main text, the confining media filled with the Thomas-Fermi fluid are chosen to have a length d TF = 10 nm. Increasing d TF increases the agreement between theory and simulations in Fig. 3 due to the decay in the disjoining energy between the two TF surfaces but at the price of enhanced numerical cost ( Supplementary Fig. S8). For the TF-TF interaction, a purely repulsive power law of the form U (r) = E/r n is added to avoid numerical infinities when particles overlap. We use n = 8 and E = 10 3 kcal/mol/Å 8 but we checked that the detailed form of the interaction potential does not qualitatively influence our simulation results as shown in Supplementary Fig. S7. Capacitance determination. The capacitive behavior of our virtual Thomas-Fermi fluid was checked as this provides an important benchmark to assess its physical validity. Using a direct measurement approach, the capacitance was estimated using MD simulations in which the system is sandwiched between two electrodes having an overall charge +Q and −Q. As discussed in the main text, two systems were considered to verify the consistency of the obtained results: the virtual Thomas-Fermi alone and a composite system made up of a dielectric layer confined by the virtual Thomas-Fermi fluid (for the latter, two dielectric materials were considered: either a vacuum layer or a molten salt). The electrodes used for the capacitance measurements consist of point charges q w = 0.01 arranged on a 1Å 2D grid (see Supplementary Fig. S12 for a molecular simulation snapshot), resulting in a total charge of Q = ±0.166 C/m 2 . It was checked that this value is low enough to ensure that the capacitance response of the system is in the linear response regime so that the capcacitance C is readily obtained from the electrostatic potential drop ∆Ψ between the two electrodes. The TF fluid is separated from the point charges by 1Å via a reflecting wall denoted by the gray shaded areas in Fig. 4(a). The potential drop is obtained from Poisson equation by integrating twice the charge density profile, ∆Ψ(z) = − z −∞ dz z −∞ dz e(ρ + − ρ − )/ε 0 as shown in Supplementary Fig. S1.