for: A screening of results on the decay length in concentrated electro-lytes †

Screening of electrostatic interactions in room-temperature ionic liquids and concentrated electrolytes has recently attracted much attention as surface force balance experiments have suggested the emergence of unanticipated anomalously large screening lengths at high ion concentrations. Termed underscreening, this effect was ascribed to the bulk properties of concentrated ionic systems. However, underscreening under experimentally relevant conditions is not predicted by classical theories and challenges our understanding of electrostatic correlations. Despite the enormous effort in performing large-scale simulations and new theoretical investigations, the origin of the anomalously long-range screening length remains elusive. This contribution briefly summarises the experimental, analytical and simulation results on ionic screening and the scaling behaviour of screening lengths. We then present an atomistic simulation approach that accounts for the solvent and ion exchange with a reservoir. We find that classical density functional theory (DFT) for concentrated electrolytes under confinement reproduces ion adsorption at charged interfaces surprisingly well. With DFT, we study confined electrolytes using implicit and explicit solvent models and the dependence on the solvent's dielectric properties. Our results demonstrate how the absence vs. presence of solvent particles and their discrete nature affect the short and long-range screening in concentrated ionic systems.

of the same kind as the silicon atoms in the mica sheets. The mica sheets are separated from the walls by an interlayer of potassium cations. The system is simulated with 2D periodic boundary conditions for Lennard-Jones interactions and the Yeh-Berkowitz correction 8 for the Ewald sum with an effective vacuum layer between the periodic images of two times the simulation box volume. The initial structure is separated by a distance of D = 15 nm and filled with water, followed by steepest descent minimisations with a maximum of 10000 steps at a 0.01 nm step size and terminated at a maximum force of 5 kJ mol −1 nm −1 and subsequent equilibration in the constant volume NVT ensemble, see Fig. S1(b). Then, a NPT run with anisotropic Berendsen pressure coupling was performed for 10 ns, yielding lateral dimensions (2.88 × 3.00) nm 2 and width in z direction 2.04 nm. In the production runs the distance D is varied and corresponding amounts of water and salt ion pairs are as described below. Note that in the aqueous phase the K surface counter-ions have been replaced by Na to avoid mixing effects of different salts in this study.

S1.2 Simulations at controlled water chemical potential
The computer simulation of interacting surfaces is still posing a challenging task, as the chemical potential of the confined species needs to be accounted for. We here follow the thermodynamic extrapolation method established for a single species (namely water) 9 but generalise this approach to the case of (binary) mixtures. In detail, we imply charge neutrality in the simulation box by treating separately the counter-ions that compensate for the surfaces' charge and the added salt. The number of water molecules, N w , and salt ion pairs, N ion is adjusted for a given, fixed, surface separation D in order to match the chemical potential in bulk at the desired concentration c b in bulk. The latter is determined as described below in a cubic simulation box containing N w = 2111 water molecules, corresponding to a box length of ∼ 4 nm. Using that the concentration of water at atmospheric conditions and T = 300 K is 55 mol/l, the number of salt molecules that need to be added then follows as N w = c ion · N w /(55 mol/l), corresponding to N ion = 4 at the lowest concentration considered, c b = 0.1 mol/l, up to N ion = 77 at the highest concentration considered here, c b = 2.0 mol/l.
In order to measure the chemical potential of species α with an accuracy as high as 0.01 k B T , we split µ into its contributions stemming from the ideal gas, from the Lennard-Jones interactions, and the Coulomb energy, respectively, according to Note that at interfaces or in confinement, the terms appearing in S1 in general depend on the position z normal to the interface. The first term is the ideal gas contribution k B T log ρ α (z)Λ 3 α , where ρ α (z) is the local number density and is the thermal de Broglie wavelength and m α the mass of particles of species α. In thermodynamic equilibrium, the total chemical potential µ α is independent of the position z, therefore it can be evaluated at an arbitrary position, which we choose for convenience in the centre of the water slab between the surfaces. To obtain the excess chemical potentials it has proven convenient to employ a modified Hamiltonian approach 9,10 evaluating the free energy difference between a non-interacting molecule/salt ion pair and a fully coupled one, where the coupling is characterized by a parameter λ ∈ [0, 1] scaling the vdW/electrostatic interactions. The free energy difference between distinct λ states is then evaluated by the MBAR method, minimising the statistical uncertainty. 11 . In our simulations, we employ 38 λ states, cf. Table S1 and evaluate the corresponding free energy changes using PyMBAR 12 . However, a significant reduction of the computation time can be achieved by realising that a TIP4P/ε water molecule consists of a single LJ interaction site only, therefore it is convenient to evaluate µ LJ w via the Widom Test Particle Insertion method (TPI) 13 , rendering the first 20 λ states obsolete for water.  To estimate the dielectric constant of the dumbbell solvent model presented in detail in the main text, we here follow an approach based on a plate capacitor setup, depicted in the inset of Fig. S2. The capacitor consists of two plates with opposite charge density ±Q and is filled with the solvent model at the bulk concentration ρ s = 55.6 mol/l. A potential difference ∆ψ = 2ψ 0 = 2 k B T ≈ 0.05 V is applied over a slit with D = 4 nm. We then determined the apparent dielectric constant of the capacitor, i.e. assuming that it can be described in terms of a homogeneous dielectric background of the same width D based on the surface charge density: for a plate capacitor filled with a homogeneous dielectric, latter is given by Note that the actual dielectric constant in bulk might be underestimated using this approach since interfacial effects (observed in Fig. S2 as deviations from the linear behaviour) need to be taken into account properly, which is typically done either by assuming a serial setup of capacitors 14 or through effective medium theory 15,16 . As observed in the aformentioned studies, the apparent dielectric constant can -depending on the interfacial properties -converge slowly to the bulk dielectric constant with increasing slit width D.
Solving Eq. (S2) for ε r and employing the electrostatic potential profile obtained in our DFT formulation explained in detail in the main text, yields ε r /ε b ≈ 2.4. Since we want to compare solvent-implicit calculations to the explicit solvent case, where all electrostatic interactions are scaled by ε b = 4.1, the corresponding value in the solvent-implicit calculations ε r ≈ 2.4 · 4.1 = 9.8, which is the value given in the main text.

S2.2 Dielectric constant of the dumbbell model from explicit molecular dynamics simulations
In order to validate the estimation of the solvent dielectric constant from the capacitor setup, we additionally performed MD simulations of a pure bulk system of dumbbell molecules using the GROMACS simulation package 1 . Note that we explicitly checked that simulations using ESPResSo 17 and LAMMPS 18 converge to the same value ε r . In detail, we set up a cubic box of 2 nm sidelength containing 268 dumbbell molecules in all simulation packages, corresponding to a density of 55.6 mol/l. We explicitly checked for the absence of finite size effects by performing simulations using the same density in cubic boxes of sidelengths 1.5 − 4 nm using GROMACS, yielding numerically identical results within the statistical uncertainty. The dumbbell charges are scaled by a factor 1/ √ ε b to account for the dielectric background ε b = 4.1 and the hard-sphere interaction is replaced by a Weeks-Chandler-Andersen potential, where u 0 = 12k B T ensures a strong repulsive interaction between the solvent sites of diameter σ s . The bead-bead distance of the solvent model is fixed to 0.2 nm using the LINCS algorithm 19 with 4th order expansion and 1 iteration step. Temperature is controlled using the Bussi-Donadio-Parrinello thermostat 20 with a characteristic timescale of 0.1 ps and electrostatic interactions are treated using the Particle-Mesh Ewald method with the relative strength of the shifted direct potential at the cut-off distance of r Coulomb = 0.4 nm set to 10 −5 . We chose to set the mass of each bead to 9 atomic mass units and a integration timestep of 0.5 fs was employed with positions recorded every 1 ps for a simulation time of 5 ns. After neglecting an initial equilibration time of 100 ps, the dielectric constant is then evaluated from the total dipole moment fluctuations, calculated as M M M = ∑ i q i r r r i , where q i is the charge and r r r i the position of bead i, according to 21