Electronic screening using a virtual Thomas–Fermi fluid for predicting wetting and phase transitions of ionic liquids at metal surfaces

Of relevance to energy storage, electrochemistry and catalysis, ionic and dipolar liquids display unexpected behaviours—especially in confinement. Beyond adsorption, over-screening and crowding effects, experiments have highlighted novel phenomena, such as unconventional screening and the impact of the electronic nature—metallic versus insulating—of the confining surface. Such behaviours, which challenge existing frameworks, highlight the need for tools to fully embrace the properties of confined liquids. Here we introduce a novel approach that involves electronic screening while capturing molecular aspects of interfacial fluids. Although available strategies consider perfect metal or insulator surfaces, we build on the Thomas–Fermi formalism to develop an effective approach that deals with any imperfect metal between these asymptotes. Our approach describes electrostatic interactions within the metal through a ‘virtual’ Thomas–Fermi fluid of charged particles, whose Debye length sets the screening length λ. We show that this method captures the electrostatic interaction decay and electrochemical behaviour on varying λ. By applying this strategy to an ionic liquid, we unveil a wetting transition on switching from insulating to metallic conditions. Ionic and dipolar liquids display unexpected behaviours, especially in confinement, that are relevant to energy storage, electrochemistry and catalysis. An approach that involves electronic screening while capturing molecular aspects of interfacial fluids is now proposed.

level D E F according to λ = ϵ/e 2 D E F (where ε is the dielectric constant and e the elementary charge); E F is directly related to the free electron density n 0 as E F = ℏ 2 (3π 2 n 0 ) 2/3 /2m e , where m e is the electron mass and ℏ = h/2π is Planck's constant (Supplementary Section IIA). Despite this available framework, the development of classical molecular simulation methods to understand the microscopic behaviour of classical fluids in contact with imperfect metals is only nascent. Although insulators are treated using solid atoms with constant charge, metals must be described using an effective screening approach. The concept of image charges can be used for perfectly metallic and planar surfaces 20 , but refined strategies must be implemented for non-planar surfaces, such as variational 21,22 or Gauss's law [23][24][25] approaches to model the induced charge distribution in the metal. A recent proposal 8 builds on our TF framework 7 to propose a computational approach based on variational localized surface charges that account for electrostatic interactions close to imperfect metals.
Here we develop an effective yet robust atom-scale simulation approach that allows us to consider 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 TF fluid made up of light charged particles. Owing to its very fast response, this effective TF fluid mimics metal induction within the confining surfaces on sampling the confined system configurations using Monte Carlo or molecular dynamics (MD) simulations. After a straightforward implementation in existing simulation packages, this strategy provides a means to impose a TF screening length that is directly linked to the equivalent virtual fluid Debye length. By adopting a molecular level description, our approach captures all the specificities inherent to confined and vicinal fluids at the nanoscale (for example, density layering, slip-page and non-viscous effects). This virtual TF model correctly captures electrostatic screening within the confined system on varying the confining host from perfect metal to insulator conditions. The presented molecular approach is also shown to accurately mimic the expected capacitive behaviour, and therefore opens perspectives for the atom-scale simulation of electrochemical devices that involve metals of various screening lengths and/or geometries. To further assess this effective strategy, by considering the freezing of a confined ionic liquid, we showed that the experimental behaviour observed in Comtet et al. 18 can be rationalized as the salt is predicted to melt at a temperature larger than the bulk melting point and that the shift in the melting point should increase with a decreasing screening length. Last, but not least, we used this novel method to consider an ionic liquid between two parallel planes and unravelled a continuous wetting transition as the surfaces were tuned from insulating (non-wetting) to metallic (wetting).
A few remarks are in order. The impact of image forces induced in non-insulating solid surfaces when setting them in contact with charged and/or dipolar molecules was investigated in depth by Vorotyntsev and Kornyshev 27,28 . Using a formal treatment of the Coulomb interactions with appropriate electrostatic boundary conditions, these authors and co-workers considered complex problems that ranged from electric double layers 29,30 to ion and molecule adsorption 31,32 . This robust framework was employed for simple point charges, but also for complex systems, such as ionic liquids (whose chemical structure leads to ion-specific effects beyond electrostatic interactions). In contrast, the TF model-which is at the heart of our MD approach-is based on a simplified description of conductivity in solids. This simple formalism allows us to cover situations from perfect and imperfect metals to semiconductors by varying the screening length λ. The TF model is known to be a rather crude approximation when used to describe the interaction between ions and solid surfaces 33 . In particular, owing to the use of a semiquantum description level of charges in solids, it fails to accurately describe quantum effects, such as Friedel oscillations, the image plane position and motion with respect to metal atoms, and the electron spillover outside the metal (a fortiori, the impact of charge adsorption at the solid surface is even more complicated to account for). As a result, in some cases, such as for transition metals, metallic surfaces will not be appropriately described using the TF approximation.
Yet, despite these drawbacks, as shown here, the TF model provides a simple framework to implement the complex electrostatic interaction-which includes image forces-that arises in the vicinity of metallic surfaces. Such a multicontribution interaction, which results from the theoretical derivation described below (and in more detail in Supplementary Section II), includes time-demanding numerical estimates that cannot be calculated on the fly in MD. As a result, despite its limitations and weaknesses, the TF approach is an approximate scheme to account for electrostatic screening near solid surfaces with properties that range from perfectly metallic to insulating. In fact, in the present molecular approach, rather than formally implementing the equations that correspond to the TF model, we rely on a simple effective procedure in which electrostatic screening induced by a fluid of charges located in the solid material is mapped onto the TF length. As shown in this article, we believe that this matching is a valid approach as the mapped screening length is consistent with the observed capacitance behaviour of the system. Therefore, despite its simplicity, we believe that our effective molecular strategy provides a simple tool to investigate complex surface electrostatic phenomena that occur at solid-liquid interfaces. Even if this approach provides only a first-order picture (in the sense that some of the above-mentioned aspects could be missing), the results reported below suggest that it captures striking phenomena, such as the impact of screening length on capillary freezing. Among its strengths, this method allows us to consider nearly any screening length, but also any electrode shape; indeed, by using a liquid phase of screening charges within the solid material, solid surfaces with a complex morphology can be considered (this contrasts with other approaches that require defining an underlying atomic structure that bears localized polarizable electrons). Moreover, as it can be directly implemented in any standard MD package, our approach can be used for fluids that range from simple charged and dipolar molecules to more complex systems, such as ionic liquids.
Interactions at Thomas-Fermi metal interfaces. Figure 1a depicts point charges in an insulating medium α of relative dielectric constant ε α close to a metal β of TF length λ. As shown in Supplementary Section IID, the electrostatic energy of two charges i and j at distances z i and z j from the dielectric-metal interface and separated by 1/2 with R ij the in-plane distance reads: 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. 1b, U CC is the Coulomb interaction energy between the charges i and j, and 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 ), which corresponds to the interaction with its image in the metal, and a two-body contribution U CI • λ (z i , z j , R ij ), which corresponds to the interaction with the induced charges due to all the other charges j. Analytical expressions exist for U CI • λ and U CI • λ (refs. 7,[9][10][11] ), but U II λ must be estimated numerically from the energy density, that 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 the 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 an expensive integration on the fly as analytical treatment for imperfect metals is only available in closed forms in asymptotic limits 7,14 . Here we model the resulting complex electrostatic interactions between the ions of the liquid thanks to a 'virtual TF fluid' located within the confining solids (Fig. 2a). Our approach builds on the direct analogy between Fig. 1 | Electrostatic interactions in the vicinity of metal surfaces. a, Ionic liquid (IL) in an insulating medium α close to an imperfect metal β having a TF length λ. b, One-and two-body interactions for two point charges i,j at distances z i and z j from the surface and separated by in-plane distance R ij . The induced charge distribution ρ I k r for k = i,j (denoted by half-ellipsoids) within the metal is of opposite sign and decays over λ. The coloured arrows show the different energy contributions U CC , U CI λ and U II λ (equation (1)).
the TF screening of electrons and the Debye-Hückel equation for electrolyte solutions. In the linear TF formalism, the induced electronic charge density in the metal is q TF ρ I (r) = −ε 0 ε β k 2 TF Ψ β (r), where ε 0 is the vacuum permittivity, ε β the relative dielectric constant and k TF = [e 2 D E F /ε 0 ε β ] 1/2 the TF wavevector (Supplementary Section IIA). Combined with the 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 screening length λ between 0 (perfect metal) and ∞ (insulator) can be efficiently mimicked depending on q TF , ρ TF and T TF . This virtual system allows us to simulate the complex electrostatic interactions within the ionic liquid in the vicinity of an imperfect metal.
As mentioned in the previous paragraph, by tuning the different parameters inherent to the TF fluid-namely, the fluid particle charge q TF , temperature T TF and density ρ TF -confining solids with an electrostatic screening length that ranges from metallic to insulating can be mimicked. Yet, equation (2) shows that mapping the fluid of mobile charges onto the TF model only requires us to set a single parameter ρ TF q 2 TF /T TF (fixing ε β = 1). In practice, although this implies that different combinations for these parameters can lead to the same effective screening, there are a few constraints that should be verified. First, as discussed later in this paragraph, to account for the ultrafast dynamic relaxation of charges in the solid compared with that in the confined salt, the TF temperature T TF was chosen to be sufficiently large. Although this is important to capture the order of magnitude difference between these two timescales, it implies that two thermostats have to be used to properly regulate the temperature of the two subsystems. In this regard, we emphasize that the results reported here were obtained using either a Nose-Hoover thermostat or a Langevin thermostat (however, a Langevin thermostat is recommended as it was found that it ensures the two subsystems display uncoupled homogeneous and constant temperatures). Moreover, in contrast with the temperature of the confined charges, T TF should not be seen as a physical temperature, but rather as a parameter that governs the TF fluid dynamics and, hence, an effective screening. Second, as the ionic force in the Debye length is directly the product of the squared charge and density, that is, ρ TF q 2 TF , we can tune the effective screening in the metal by tuning either one or both of these parameters. In practice, after conducting several tests, we found it more efficient to keep the number, and hence the density, of the charges in the TF fluid constant. Indeed, on the one hand, to play with q TF at a constant ρ TF offers more flexibility in tuning the screening length. On the other hand, treating very imperfect metals with constant q TF requires considering a very small number of TF ions, which leads to a poor sampling and/or statistics performance. Before going into more technical details, we emphasize that the exact choice of parameters made to mimic different screening lengths is expected to impact the dynamics and/or kinetics of the observed phenomena (by setting a given temperature T TF , we do impose a relaxation time). However, as in classical thermodynamics, we do not expect this effective relaxation within the virtual TF fluid to affect the equilibrium properties of the charged liquid confined between the metallic surfaces provided it is sufficiently fast.
In more detail, in our MD approach, to ensure that the particles in the effective TF fluid relax fast, their mass and/or temperature were chosen as much smaller or larger than their counterpart in the confined system; typically, m TF ≈ 0.01m and T TF ≈ 10T (which requires typical integration steps of 0.1 and 1 fs, respectively). In practice, as shown in Fig. 2a, 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 TF fluid of density ρ TF . Once ρ TF and T TF are set, λ D is varied by tuning q TF according to equation (2)from q TF = 0 (λ D → ∞) for an insulator to q TF = 1 (λ D = 0.03 nm) for a nearly perfect metal. All the simulations reported in this article were carried out using MD, but they could easily be implemented into Monte Carlo algorithms to perform calculations in other ensembles (Grand Canonical, isothermal and isobaric ensembles, for instance). All the simulation details are provided in Methods. To consider a set of explicit charges with no solvent to model the TF fluid, we used consistently a relative dielectric permittivity equal to one. However, this raises the question of the true dielectric constant of the solid material that we intend to mimic. Describing s-p metals using a TF model is usually done by accounting for the dielectric constant of the underlying ion structure in the metal. Such a dielectric background, which arises from the polarizability of the core electron shells around each metal atom, differs from one metal to another (with values of the order of one to a few times the vacuum dielectric permittivity). In semiconductors, such a dielectric background, in which the outer electrons move, arises from interband transitions with values that can be close to ten. As another example, when semi-metals are modelled, a large TF screening length is used together with a large dielectric constant to account for the low charge-carrier concentration and absence of bandgaps, respectively. In this context, we mention that the dielectric constant in semi-metals can be evaluated by considering the electron density of states, as proposed by Gerischer and co-workers 34,35 . These specific, yet representative, examples illustrate that the nature of the solid material manifests itself in the TF screening length but also in the dielectric constant. Although this duality (screening and dielectric properties) is, of course, of prime importance when modelling the complex physics of metallic and semiconducting surfaces, we simplify the problem in the present approach as we only aim to mimic in an effective but quantitative fashion the electrostatic screening induced by the solid surface on the confined or vicinal charges. To this end, as already introduced above, our proposed molecular approach adopts a different view as it mimics a liquid of charges to produce an effective screening that corresponds to an underlying TF model-this is what is referred to in our article as a 'virtual' TF fluid. In other words, as is shown below, by calculating the induced electrostatic screening length, we can map this pseudo TF medium onto practical situations by establishing a consistent relation between its Debye length λ D and the induced electrostatic screening length λ. Therefore, we are not attempting, per se, to implement the electrostatic screening equations as derived using the TF formalism. In this regard, as is established later, we believe that our mapping procedure is sound and robust as the inferred screening length was found to be consistent with that corresponding to the observed capacitance behaviour of our system.

Two-dimensional crystal at metallic interfaces.
To validate our effective approach, we considered a two-dimensional (2D) square crystal of lattice constant a = 1.475 nm made up of charges of ±1e and located at a distance d from a metal (Fig. 2). Owing to the periodic boundary conditions, a second pore-metal interface is present at a distance d w = 20 nm. Yet, as shown in Supplementary Fig.  9, this second interface does not affect the electrostatic energy as d w is large enough. In the TF framework, the charge density ρ I at a position r in the metal induced by a charge q located in (0,0,d) reads (Supplementary Section IID): where R = (x 2 + y 2 ) 1/2 is the lateral distance to the charge q, J 0 is the Bessel function of the first kind and κ 2 = K 2 + k 2 TF . Figure 2c,d shows the induced charge density ρ I (d,r) as obtained by summing equation (3) for the 2D crystal when d = 0.22 nm and λ D = k −1 TF = 0.25 nm (as discussed in Supplementary Section III), equation (3) must be summed over all the crystal periodic images, but it was found that the sum converges quickly). For comparison, Fig. 2e,f shows ρ I (d,r) as obtained using our effective approach from the local charge density in the metal, that is ρ I = e(ρ + TF − ρ − TF ). In contrast to ρ I (d,r) in the TF model, owing to their finite size, the fluid charges in the simulation cannot approach arbitrarily close to the metal-pore surface. For consistency, the analytical and simulation data were compared by defining z = 0 in the simulation as the position at which the TF fluid density became non-zero. Figure 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 the opposite sign. Moreover, as expected from the TF framework, the induced charge distribution in the effective simulation decays over the typical length λ D .
Our effective approach was assessed quantitatively by probing the energy of the 2D ionic crystal as a function of its distance d to the metal surface for different screening lengths λ. The simulated electrostatic energy U λ (d) consists of all ion-pair contributions in equation (1), as discussed in Supplementary Section V. Figure 3 compares the total energy U λ as a function of the distance d with the numerically evaluated prediction from equation (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 behaviour of the confining medium assuming a screening length λ = c 0 + c 1 λ D + c 2 λ 2 D (with λ D the ion gas Debye length, in our system c 0 = 0.22 nm, c 1 = 0.91 and c 2 = 0.28 nm −1 ). Such values do not simply correspond to fitting parameters that allow us to match the simulated and theoretical energies; as explained in the next paragraph, they were derived so that the capacitance of the virtual TF fluid matches the theoretically expected value C = ε 0 /λ. The fact that the rescaled screening length λ also allows us to recover 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 TF ions, which prevents reaching a screening of λ ≤ σ. This is supported by the fact that c 0 ≈ σ corresponds to the value below which the repulsive interaction potential in the TF fluid becomes larger than k B T. c 1 arises from the non-ideal behaviour of the effective TF fluid, which leads to over-screening compared with that of an ideal gas with the same charge density ρ TF (c 1 = 1 corresponds to the ideal behaviour); c 2 ≠ 0 indicates non-linear effects in electrostatic screening which go beyond the linear approximation used in the TF framework.
To illustrate the ability of our approach to capture the impact of various screening lengths-from insulators to metallic surfaces-on electrostatic interactions, we show in Supplementary Fig. 1a the electrostatic energy that arises from image-charge interactions, U CI (λ), for a molten salt confined between two solid surfaces as a function of the electrostatic screening length λ. To probe the impact of such interactions with induced charges, a fixed liquid configuration was considered as it implies that the direct Coulomb interaction U CC (λ) is constant. As expected, on decreasing λ, the overall electrostatic energy decreased as the interactions with the induced charges in the metallic surfaces become more negative. Moreover, by extrapolating U CI (λ) to perfect metallic conditions (that is, λ → 0), we recovered the expected charge image contribution that corresponded to half of the Coulomb energy. Such data show that our molecular simulation strategy does mimic-albeit in an effective fashion-the electrostatic screening induced by metallic surfaces. In this respect, we emphasize that this approach can be extended to almost any surface geometry and topology. First, this versatile model does not require the input of an underlying atomic structure for the surface (by describing the charges in the metal as a fluid, one does not need to consider an atomic lattice to which charges are linked). Second, any geometry from a simple flat or cylindrical surface to disordered or rough surfaces can be considered, as it simply requires us to encapsulate the screening charges within a mathematically defined region. Correctly accounting for image forces near solid surfaces is crucial to capture the rich and complex behaviour of confined charges. In particular, as theoretically predicted by Kondrat and Kornyshev 36 , it has been observed using molecular simulation that such image forces can lead to superionic states-in which like-charge pairs form-in a metallic nanoconfinement 37 . In this respect, it was found that electron spillover leads to effective pore sizes that are narrower than the nominal pore size. Although the TF model was found to fail to predict this effective pore size reduction 38,39 , we believe thatat least as a first-order approach-the results reported in this Article show that our molecular strategy can be used to capture the impact of electrostatic screening on ions between metallic surfaces. In particular, simulations of ionic liquids at perfect metal surfaces show that when comparing perfect metal conditions (where in the implementation charges are screened on a very short scale) to electrodes at constant surface charge, image charge interactions at the perfect metal electrode do not substantially affect the adsorbed liquid 40 . This highlights the need to consider effective molecular simulation approaches, such as that reported here, to consider imperfect metals and, hence, larger screening lengths, which corresponds with many experimental situations.

Capacitance.
To further establish the validity of our novel molecular approach, an important requirement is to verify that the virtual TF fluid yields the correct capacitance behaviour. With this aim, as shown in Fig. 4a, a simple MD set-up was designed by assembling a composite system made up of two TF fluids that sandwich 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 behaviour). The molten salt (NaCl) was modelled using charged particles, ±1e, that interact via a Born-Mayer-Huggins (BMH) potential 41 (Supplementary Table 1). To prevent mixing of the TF fluid and charged system, a reflective wall of thickness ξ = 0.2 nm was positioned between the two subsystems. The whole composite was placed between two electrodes with an overall charge of +Q and −Q (all the details are in Methods). With such a geometry, the capacitance C = Q/ΔΨ is readily obtained from the potential difference ΔΨ. As shown in Fig. 4a, this molecular simulation set-up mimics in a simplified yet realistic way an experimental electrochemical cell, so we can perform a MD simulation to readily estimate the positive and negative charge-density profiles within the confined salt (ρ + (z) and ρ − (z)). Using the Poisson equation, that is, , ΔΨ(z) is determined by integrating twice the charge density profile ρ(z). Figure 4a shows ΔΨ(z) as a function of the position z within the confined liquid for different screening lengths λ. In practice, two simple situations were considered; the porosity between the two metallic surfaces is either occupied by a vacuum or by a molten salt. As expected, ΔΨ increases with increasing z as the negative electrode is located at z = 0 (positive charge adsorption) and the positive electrode is located at z = d w (negative charge adsorption). Moreover, by considering the datasets for vacuum-filled and liquid-filled pores in Fig. 4a, we observe that the slope of ΔΨ(z) in the pore region is larger for the former than for the latter. Given that C = Q/ΔΨ, this result suggests that, as expected, the sandwiched salt layer has a larger capacitance than the sandwiched vacuum layer. To provide a more quantitative picture of the system capacitance as a function of the screening length λ, we performed a more detailed analysis in which the capacitance of the different elements-confined material and TF fluid-was extracted.
The system considered here simply consists of double-layer capacitors in series so that its capacitance per unit area should verify the combination rule: where the first and second terms on the right correspond to the capacitance of the vacuum slab of width d w and that of the TF fluid (the factor 2 simply accounts for the presence of two TF/vacuum interfaces). As shown in Fig. 4b, 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. 4b, our effective approach captures quantitatively the expected capacitance behaviour of the confining medium on rescaling λ D → λ = c 0 + c 1 λ D + c 2 λ 2 D (see the discussion above). As already mentioned, the vacuum layer in the capacitor was then replaced by a slab of molten salt-see the molecular configuration shown in Fig. 4a. As expected, on inserting such a molten salt, the effective capacitance C drastically increased (that is, the inverse capacitance shown in Fig. 4b decreased). More importantly, as shown in Fig. 4c, the induced capacitance change Δ1/C observed in our simulation data followed the expected behaviour with a λ-independent value: where ε s is the permittivity of the molten salt. Furthermore, as ε s ≫ ε 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 to that of the molten salt).

Capillary freezing and melting in confinement.
In the above, our effective molecular approach was shown to capture the electrostatic energy predicted using the TF formalism for electrostatic screening in metallic materials as well as the capacitive behaviour of a molten salt sandwiched between metallic surfaces with different screening lengths. Yet, in addition to these two validation steps, it should be verified that our simple strategy allows the available experimental data for realistic materials to be reproduced. To do so, we performed calculations using our effective treatment to study the liquid-crystal phase transition in various metallic confinements, as experimentally reported by Comtet et al. 18 By considering the crystallization of an ionic liquid confined between an atomic force microscope tip and a metallic surface, these authors showed that the melting temperature is shifted above the bulk melting point and that the shift in the melting point increased with decreasing screening length. To help rationalize these results, we performed the following molecular simulation study using our effective electrostatic screening strategy in confined charged systems. First, by considering insulating surfaces, we used the direct coexistence method in which a crystalline salt coexists with its molten salt in a slit pore of size d w (various pore widths between 1.5 and 7.1 nm were considered). To mimic a physical system in which the confined phases also interact through simple dispersive and/or repulsive interactions with the surface, a 9-3 Lennard-Jones interaction potential was added between each salt atom and the solid surface. The use of such structureless surfaces to describe the confining solids was made to avoid inducing a peculiar crystalline structure by employing a given molecular periodic lattice. Moreover, to ensure that such dispersive interactions do not impact too much on the melting point in confinement, the potential wall depth chosen was of the order of k B T and, hence, at a value much lower than that of the electrostatic interactions. Using this set-up, MD simulations in the canonical ensemble were performed for different T to determine the melting temperature T m as follows. The crystalline salt melts into the liquid phase for T > T m , whereas the molten salt crystallizes into the crystal phase for T < T m . The insert in Fig. 5a shows the shift in the melting point with respect to the bulk melting point, ΔT m = T m − T m,0 as a function of pore size d w . In agreement with the experimental data 18 , these data show that the salt confined between the insulating surfaces has a melting temperature above the bulk melting point. To assess the impact of the electrostatic screening length λ on capillary freezing, we used the following expression in which the liquid-surface and crystalsurface interfacial tensions between the metallic surface and the crystal (x = c) or the liquid (x = l) is given by its value for the insulating surface corrected for the charge image interactions U CI (λ): γ xw (λ) = γ xw (∞) + ρ x ℓU CI (λ), where ℓ is a scaling length that converts a volume energy U CI into a surface energy. As shown Supplementary  Fig. 1b, such a simple relationship reasonably captures the impact of the electrostatic screening on the liquid-surface interfacial tension, which was assessed using independent simulation through the Irving-Kirkwood formalism: γ(λ) = L z /2〈P N − P T 〉, where the terms in brackets are the average normal and tangential pressures, respectively, L z is the box length in the z direction and the factor 2 accounts for the two interfaces in the slit geometry. Despite the fact that the simple expression γ(λ) ≈ ρU CI (λ) neglects the impact of screening on the entropy of the liquid, it provides an accurate description of the surface tension change induced by electrostatic screening in the metallic surfaces (we note that this simple equation holds even better for the crystalline phase as its entropy is negligible). This allows writing that Δγ(λ) = Δγ(∞) + (ρ l − ρ c )ℓU CI (λ). As shown in Fig. 5a, given that U CI (λ) < 0 becomes more negative on decreasing λ and ρ c > ρ l , this simple scaling predicts that the shift in the melting point ΔT m /T m,0 increases as the surfaces turn from insulating to metallic.
To confront our results with the experimental data on capillary freezing 18 , Fig. 5b shows the impact of the screening length λ on the capillary pore size d s w (λ) = [2Δγ(λ)T m,0 ]/[ρ c ΔHmΔTm] below which salt crystallization is observed. In more detail, to compare quantitatively our data with those obtained experimentally for a room temperature ionic liquid, we plot the shift induced by surface metallicity in this capillary pore size with respect to that for an insulating surface Δd s Fig. 5b by Δd s w (0) = d s w (0) − d s w (∞) allows us to define a quantity in the y axis that varies from 0 for a perfect insulator (λ → ∞) to 1 for a perfect metal (λ = 0). Moreover, in addition to providing a mean to compare with experimental data for any other system, such a normalized quantity provides data that are independent of the specifically chosen value ΔT m /T m,0 . As can be seen in Fig. 5b, our theoretical predictions do capture the experimentally observed behaviour, which indicates that the capillary length d s w increases on decreasing the screening length λ. In this plot, the screening length λ is normalized by a length a, which is a molecular characteristic of the ionic systems under scrutiny (the x axis is plotted as a/λ). As shown in Fig. 5b, a perfect quantitative agreement between our simulated data for a simple molten salt and the experimental data for the ionic liquid is observed. This provides the value a = 1 nm for the simulated molten salt, whereas a slightly smaller value of a = 0.335 nm was used in Comtet et al. 18 in the analysis of the experimental results for the [BMIM][BF4] (1-butyl-3-methylimidazolium tetrafluoroborate) room-temperature ionic liquid (although the simplified modelling in Comtet et al. 18 was used). The parameter a can be seen as a characteristic length that describes the impact of electrostatic screening on freezing. A slightly smaller length a for the room temperature ionic liquids-with a more complex molecular structure-suggests a smaller impact of electrostatic screening on capillary freezing for these complex ions with respect to a simple salt. This is expected given that substantial entropy and molecular packing aspects largely affect the crystallization of room-temperature ionic liquids (in particular, these contributions govern their low melting point). These important results provide a quantitative microscopic picture for this recent experimental finding in which the capillary freezing of an ionic liquid was found to be promoted by metal surfaces. Beyond this important result, this study further suggests that the simple effective approach presented here captures the rich and complex behaviour of charges confined between metallic surfaces.  Fig. 4 | Capacitive behaviour of virtual TF fluids. a, Typical molecular configuration of the MD set-up employed to determine the capacitive behaviour of the virtual TF fluid. A molecular system-here, a molten salt represented by the blue and green spheres-is confined between two virtual TF fluids (yellow and pink spheres). This composite system is sandwiched by two electrodes that have a constant surface charge of ±Q. As shown in the top panel, with this set-up the electrostatic potential profile ΔΨ(z) can be determined by integrating twice the resulting charge density profile ρ(z) = e[ρ + (z) − ρ − (z)]. Such numerical assessments were performed for different screening lengths λ with the confined material either a vacuum layer (dashed lines) or a molten salt (solid lines). b, Reciprocal capacitance 1/C of the empty TF capacitor for different λ D versus the analytical prediction (black line) for two double-layer capacitors in series. Inset: the simulation data collapse onto the same master curve when plotted using the effective screening length λ. c, Reciprocal capacitance 1/C as a function of the screening length λ for a capacitor made up of a vacuum (blue symbols) or a molten salt (red symbols) confined between the TF fluids. As expected, for ε s = ε α ε 0 ≫ ε 0 , the difference between both systems (green symbols) is close to d w /ε 0 .

Wetting transition.
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 modelled using charged particles ±1e). Figure 6a shows the number density profiles ρ n (z) for the salt and TF fluid for different λ. A crossover is observed on decreasing λ; although the salt is depleted at the insulating interface, a marked ion density peak appears under metallic conditions (in contrast, the density profile for the TF fluid is nearly unaffected by λ). This behaviour suggests that the system undergoes a wetting transition on changing the dielectric-metallic nature of the confining medium (perfect wetting or non-wetting for the metal or 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 (introduced above). By considering the salt in its liquid (l) and gas (g) states in contact with the metal (m), we estimated for various λ the gas-metal γ gm (λ) and liquid-metal γ lm (λ) surface tensions. Note that, in MD simulations, the various interfaces (gas-metal, liquid-metal and gas-liquid) are investigated separately 42 ; accordingly, the gas (liquid) phase is metastable when the liquid (gas) phase is stable, i.e. it wets the surface. To investigate the impact of surface metallicity on wetting, we then evaluated the spreading coefficient S from the gas-liquid, surface-liquid and surface-gas interfacial tensions defined as S = γ gm − γ lm − γ lg (refs. 43,44 ). Figure 6b shows the dependence of the spreading coefficient S on the screening length λ. This plot reveals the wetting behaviour of the salt solution on the metallic surfaces under scrutiny, depending on the sign and amplitude of S. As shown in Fig. 6b, the sign of the spreading coefficient S changes from S < 0 to S > 0 for λ ≈ 0.28 nm, that is, when the nature of the surfaces switches from insulating to metallic as the screening length λ decreases. This is the signature of a continuous wetting transition of the liquid salt from partial wetting (S < 0) for large λ (more insulating surfaces) to complete wetting (S > 0) for small λ (more metallic surfaces). In more detail, for λ < 0.28 nm (more-metallic surfaces), S > 0, that is, γ gm > γ lm + γ lg ; this reflects that a wetting film with two interfaces (solid-liquid and liquid-gas) is of a lower surface free energy compared with that of a solid gas interface. As a result, in these conditions, the system is perfectly wetting with a liquid film spreading over the metal surface. However, for λ > 0.28 nm (more-insulating surfaces), S < 0, that is. γ gm ≤ γ lm + γ lg , so that the liquid phase wets the surface incompletely. On a macroscopic surface, this would lead to the formation of a liquid droplet at the solid surface with a contact angle θ related to S according to S = γ lg (cosθ -1) (refs. 43,44 ).
The data in Fig. 6b for partial wetting suggest that cosθ tends to 1 in a linear fashion on decreasing the screening length λ. This scaling suggests that the wetting transition induced by tuning the solid surface from an imperfect to perfect metal is a first-order transition 45 . Moreover, increasing the screening length λ beyond 1 nm (a more and more insulating surface) is expected to lead to complete drying (cosθ = -1). For simple liquids such a drying transition is expected to be a second-order transition, in contrast to the wetting transition discussed above 46 . As shown in the inset of Fig. 6b, the change in Δγ between the insulator and metal is found to scale with the liquidgas density contrast: where ρ l ≫ ρ g was assumed in the second equality. As expected from the TF model, the inset in Fig. 6b shows that α(λ) ≈ U CI λ as the charge interaction with the induced density distributions (including the charge image) dominates the surface energy excess.
Despite the key role of electrostatic interactions-which includes screening induced by metallic surfaces-on the behaviour of charges near the surfaces, we emphasize that surface wetting is also strongly affected by the so-called ion-specific effects. As in bulk electrolytes, these effects, which arise from the ion molecular structure, give rise to a complex physicochemical behaviour interaction of charges and dipolar molecules near the surface. Noteworthy, the classical Frumkin-Damaskin theory describes the relative strength of electrostatic interactions in the vicinity of a charged electrode with respect to interactions responsible for the adsorption of small polar molecules. This model leads to the Frumkin adsorption isotherm, which describes how an electrode polarization increase induces desorption of polar molecules concomitantly with the adsorption of water and ions 47 . In this context, owing to its versatility, our molecular strategy of electrostatic screening between metallic surfaces is suited to account for such molecular and physicochemical effects, as it relies on a general MD approach that can be employed with any available force field. In fact, this is one of the assets of this effective approach that can be used for ionic systems (regardless of the ion-structure complexity), but also for dipolar liquids, which are expected to be affected by electrostatic screening when confined between metallic surfaces.

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

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41563-021-01121-0.

Molecular dynamics simulations.
All the simulations were carried out using the LAMMPS simulation package 48 (stable release 7 August 2019). Electrostatic interactions were calculated using the PPPM method with an accuracy of at least 10 −4 and a real-space cutoff of r c = 12.5 Å. Periodic boundary conditions were used in all dimensions with the non-electrostatic interactions being cut and shifted to zero at r c . For the simulations of the TF fluid and the salt in contact with an insulating-vacuum interface, the interactions between periodic images were not screened, so we employed the slab correction method proposed by Yeh and Berkowitz 49 with a vacuum layer three times the simulation cell height. The non-electrostatic part of the salt-salt interactions are described using the BMH potential, which accurately reproduces the properties of NaCl (either as a crystal or as a molten salt) 41 : The corresponding force-field parameters σ, A, B, C and D are given in Supplementary Table 1. Reflective walls of width ξ = 0.2 nm were used at each metal-dielectric interface so that the TF fluid-charged system did not migrate to the pore space-confining media. The latter implies that, if an atom moves through the wall by a distance δ in a time step, its position is set back to −δ away from the wall and the sign of the corresponding velocity component is flipped.
In all the simulations presented in the main text, the confining media filled with the TF fluid were chosen to have a length d TF = 10 nm. Increasing d TF increased the agreement between the theory and simulations in Fig. 3 due to the decay in the disjoining energy between the two TF surfaces, but at the price of an enhanced numerical cost ( Supplementary Fig. 10). For the TF-TF interaction, a purely repulsive power law of the form U(r) = E/r n was added to avoid numerical infinities when particles overlap. We used n = 8 and E = 10 3 kcal mol -1 Å -8 , but we checked that the detailed form of the interaction potential did not qualitatively influence our simulation results, as shown in Supplementary Fig. 9. The positive and negative TF particles differed only in their partial charge ±q TF and only interact through electrostatic interactions with the salt. The density of the TF fluid was fixed at ρ TF = 57.5 nm −3 at a temperature T TF = 12,000 K and mass m TF = 1 amu to ensure a fast relaxation. The mass of the Na and Cl atoms was set to 22.9898 and 35.446 amu, respectively. Time integration was performed using a Verlet scheme with a timestep of 0.1 fs to allow for a fast relaxation of the TF liquid. The molten salt is simulated at a 2,000 K and the temperature coupling was performed using separate Nose-Hoover thermostats with a characteristic time of 100 timesteps or separate Langevin thermostats for the salt and TF fluids.
Capacitance determination. The capacitive behaviour of our virtual TF 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 with an overall charge of +Q and −Q. As discussed in the main text, two systems were considered to verify the consistency of the obtained results: the virtual TF alone and a composite system made up of a dielectric layer confined by the virtual TF fluid (for the latter, two dielectric materials were considered-either a vacuum layer or a molten salt). The electrodes used for the capacitance measurements consisted of point charges q w = 0.01 arranged on a 1 Å 2D grid (see Supplementary  Fig. 2 for a molecular simulation snapshot), which resulted in a total charge of Q = ± 0.166 C m -2 . It was checked that this value was low enough to ensure that the capacitance response of the system was in the linear response regime so that the capacitance C was readily obtained from the electrostatic potential drop ΔΨ between the two electrodes. The TF fluid was separated from the point charges by 1 Å via a reflecting wall denoted by the grey shaded areas in Fig. 4a (bottom panel). The potential drop was obtained from the Poisson equation by integrating twice the charge density profile, ΔΨ (z) = − ∫ z −∞ dz ′ ∫ z ′ −∞ dz ′′ e(ρ + − ρ − )/ε0, as shown in Supplementary Fig. 2b.

Data availability
All the relevant simulation input scripts are available in this repository: Schlaich, Alexander, 2021, 'Simulation input scripts for "Electronic screening using a virtual Thomas-Fermi fluid for predicting wetting and phase transitions of ionic liquids at metal surfaces"' , https://doi.org/10.18419/darus-2115, DaRUS.