The Role of Energy Scales for the Structure of Ionic Liquids at Electrified Interfaces: A Theory-Based Approach.

Ionic liquids offer unique bulk and interfacial characteristics as battery electrolytes. Our continuum approach naturally describes the electrolyte on a macroscale. An integral formulation for the molecular repulsion, which can be quantitatively determined by both experimental and theoretical methods, models the electrolyte on the nanoscale. In this article, we perform a systematic series expansion of this integral formulation, derive a description of chemical potentials in terms of higher-order concentration gradients, and rationalize the appearance of fourth-order derivative operators in modified Poisson equations, as recently proposed in this context. In this way, we formulate a rigorous multiscale methodology from atomistic quantum chemistry calculations to phenomenological continuum models. We apply our generalized framework to ionic liquids near electrified interfaces and perform analytical asymptotic analysis. Three energy scales describing electrostatic forces between ions, molecular repulsion, and thermal motion determine the shape and width of the long-ranging charged double layer. We classify the charge screening mechanisms dependent on the system parameters as dielectricity, ion size, interaction strength, and temperature. We find that the charge density of electrochemical double layers in ionic liquids either decays exponentially, for negligible molecular repulsion, or oscillates continuously. Charge ordering across several ion diameters occurs if the repulsion between molecules is comparable with thermal energy and Coulomb interactions. Eventually, phase separation of the bulk electrolyte into ionic layers emerges once the molecular repulsion becomes dominant. Our framework predicts the exact phase boundaries among these three phases as a function of temperature, dielectricity, and ion size.

Introduction

Strong electrostatic correlations in crowded environments play an important role in biology, chemistry, and physics.[1−3] For example, in molecular biology, they account for DNA packing,[4] which is crucial for the compaction of genetic material in viruses,[5] impact the cytoskeleton organization,[6] and influence transport in ion channels.[7] Furthermore, such correlations explain the thermodynamic stability of plasmas[8,9] and charged colloidal suspensions.[10,11]

Surprisingly, the complexity of these phenomena can be understood to a large degree by models derived initially for electrolyte solutions.[12] Starting from the fundamental Debye–Hückel theory for dilute solutions,[13] increasingly accurate models for concentrated electrolytes were developed,[14] taking more complex Coulomb correlations into account.

Because ionic liquids (ILs) consist only of positive and negative ions without a neutral solvent, they constitute the extreme limit for the examination of electrostatic correlations in electrolytic solutions. Indeed, ILs possess characteristic properties in the bulk regime,[15,16] but also near electrified interfaces.[17] This makes them highly attractive from both fundamental and applied perspectives.[18−24] The study of interfacial electrochemistry is of wide-ranging interest. For example, the behavior of ILs near electrified interfaces has paramount importance for their performance as battery electrolytes.[25,26]

Theoretical studies of ILs near electrified interfaces discuss the structure of charged electrochemical double layers (EDL) on atomistic/molecular scales. These include classical density functional theory (cDFT) simulations and molecular dynamics (MD) simulations. cDFT gives detailed insights into the arrangement of molecules in the EDL.[27−29] MD resolves the molecular motion and can elucidate the EDL structure.[30−32]

However, cDFT/MD simulations are limited by their computational costs. Simulations at length scales above the nanometer scale are hardly accessible to the atomistic/molecular approach. Thus, continuum theories and mean-field theories (MFT) provide a complementary methodology for the simulation of larger systems, where the microscopic details can be neglected or are used as averaged parameters (e.g., constant dielectric parameters).

Usually, MFTs for electrolytes are based on lattice gas models of ions, as first proposed by Bikermann.[33] Recently, these MFTs have attracted great interest for the study for ILs. As proposed by Santangelo for aqueous systems,[34] the extension of MFTs by higher-order electrostatic correlations is useful for the description of long-range structures emerging in electrolytes. Bazant, Storey, and Kornyshev (BSK) applied this approach to ILs near electrified interfaces.[35] By using a phenomenological model, which is based on a generalized Ginzburg–Landau functional, BSK describes charge oscillations known as overscreening and charge saturation known as crowding. Yochelis et al. rationalized this approach and extended it to bulk properties.[36−38,40] However, MFT models are usually restricted to equilibrium effects of binary ILs with structureless bulk, although rare MFT models, complemented by continuum methods[39,41] and extended to the ternary case,[42] exist.

This highlights the advantage of continuum frameworks which describe dynamic transport processes. In addition, continuum models based on rigorous physical assumptions[43] identify coupled phenomena arising from the interplay of mechanics, thermodynamics, and electromagnetic theory.[44] Furthermore, this approach allows us to develop a unified, thermodynamically consistent framework that provides the common theoretical basis for the description of different electrochemical systems.[45−49] Continuum models are not restricted to binary or ternary systems because they can be formulated for arbitrary many species, charged and uncharged. Thus, they apply to more realistic electrolytes.

Recently, we proposed such a novel continuum transport theory.[45] In this theory, we account for steric effects via the mean volume, which is due to finite molar volumes of the ion species. For this purpose, we impose a volume constraint on the electrolyte. This mechanism stabilizes the bulk structure against Coulomb collapse[50] and leads to charge saturation near electrified interfaces. Thus, our theory resolves the deficiencies of the classical Poisson–Boltzmann (PB) theory, which predicts unrealistically high interface concentrations.[14] Furthermore, the existence of finite molar volumes of the electrolyte species leads to a pressure dependence of the chemical potentials (in accordance with thermodynamic arguments).[51]

However, this bulk framework cannot describe the emergence of long-range structures in ILs near electrified interfaces. Therefore, in a joint experimental/theoretical work, we extended our framework with nonlocal interactions and validated it with results obtained from atomic force microscopy.[52] Thus, we extended the mean volume effect of the bulk theory with molecular volume exclusion due to hardcore repulsion. Our holistic framework allows us to couple dynamic transport processes occurring in the bulk electrolyte with interfacial electrochemical processes. Thus, we provide a continuum model which bridges the length scales from nanometers, e.g., EDL, to millimeters, e.g., battery cells. Moreover, our framework allows us to connect the continuum description with correlation functions generated by MD.

However, the dependence of EDL structures on molecular repulsion, molecular size, temperature, and dielectricity is still unknown. In this article, we derive such an understanding with asymptotic analysis. To this aim, we present our thermodynamically consistent transport theory with an integral formulation of nonlocal interactions in section . These correlations represent atomistic volume exclusion and lead to modified constitutive equations (eqs S-10 to S-13). Moreover, the interactions impose contributions on the stress tensor and thus modify the mechanical coupling to the transport equation (eq S-14). In section , we approximate the interaction functional with a gradient expansion, which facilitates the analytical asymptotic analysis of the EDL structure. In section , we apply our extended framework to study the EDL structure for neat ILs. When we nondimensionalize our dynamic description in section , three competing energy scales describing electrostatic forces among ions, molecular repulsion, and thermal motion appear in the theory. Because our focus lies on the formation of equilibrium structures, we discuss the stationary state in section . In section , we discuss limiting cases of our stationary theory.

We perform numerical simulations and analytical asymptotic analysis to study the interplay and the effect of the competing energy scales on EDL structures. First, in section , we discuss the EDL structure for the mean volume constraint. Second, in section , we incorporate molecular repulsion into our analysis and classify the EDL structure dependence on the relation between competing energy scales.

Theory

Generalized Transport Theory

Recently, we proposed a free-energy functional Fb = ∫ dV ρφH for the dynamical description of ionic liquids in the bulk phase.[45] In this bulk model, the Helmholtz free energy density φH(Υ) = φH(T, c1,..., cN, D, B, κ) is a function of the variables temperature T, concentration cα, dielectric displacement D, magnetic field B, and strain-rate tensor κ. This variable set Υ constitutes material-specific properties of multicomponent, viscous, and polarizable media in the liquid state.

In contrast, models describing nonlocal interactions rely on functionals Fint[Υ], such that the free energy takes the form

This functional approach constitutes a more general description for electrolyte materials and allows the incorporation of nonlocal correlations between field quantities. Such correlations typically arise from microscopic effects occurring on the nanometer scale, e.g., in the vicinity of electrified interfaces. Despite the conceptual difference between the functional approach (eq ) and the canonical bulk approach, the derivation of the resulting transport theory is rather similar to the rationale outlined in great detail in ref (45) for the free energy Fb = ∫ dVρφH. We present a detailed derivation for the functional approach in the Supporting Information (section S-1).

The extension of the free energy according to eq leads to modified constitutive equations for entropy density s, electric field strength , magnetic field , chemical potentials μα, and stress tensor σ in the form of functional derivatives (eqs (S-10) to (S-14) in section S-1.1). We evaluate this framework for the bulk energy of a linear dielectric medium discussed in ref (45) (eq S-17). The resulting forces are supplemented by contributions stemming from the nonlocal correlations (eq S-21).

For the remaining part of this work, we neglect thermal driving forces by setting the temperature equal to constant values and assume the electrostatic limit, B = 0 and . This determines the electric field by the electrostatic potential E = –Φ.

Gradient Expansion of Molecular Interactions

We proposed a model for hardcore-interactions based on a convolution functional for the interaction free energy in ref (52)leading to transport contributions in the form of (section S-1.4)The symmetric potentials determine the correlation length and the magnitude of the interaction. The number of additional parameters describing this interaction depends upon the model for . In a previous publication, we used a Lennard-Jones-type force field for .[52] Such potentials are often used in the literature.[53−57] Furthermore, becausethe potentials determine the direct pair correlation functions used in liquid-state theory.[50]

Experimental results suggest that such interactions typically decay after some ionic diameters.[52] Thus, we focus on potentials ranging over the size of one molecule. Their extent is large compared to the exponential decay of the electric field, i.e., the Debye length,[58,59] yet small compared to the battery cell.

In the SI (section S-1.3), we show that such convolution functionals Fint can be approximated in a power series of concentration gradients when are short-ranged,whereHere, Γ are symmetric perturbation coefficients of dimension [Γ2] = J m3+2 mol–2. We state the complete free-energy functional for IL electrolytes in the SI (eq S-17).

The excess chemical potentials are determined via functional derivatives (eq S-13). In section S-1.4, we show that this ansatz leads to force contributionsThe corresponding electrochemical potentials for ionic liquid electrolytes are

We specify our electrolyte model and assume a one-dimensional Gaussian interaction potential for symmetric ions,Here, denotes the characteristic interaction energy, and a is the extension of ion pairs. This material parameter determines the correlation length of the interaction. We assume that a emerges naturally from the common molar volume ν = ν+ + ν– via ν = NAa3, which is justified below (section and eq ). Thus, only is introduced as a novel independent material parameter. In contrast, potentials of the Lennard-Jones type need at least one more parameter for the well depth. In the case of a binary system, the lowest-order expansion coefficients (eq ) of the interspecies correlations for the Gaussian model in eq are (section S-3.4.1)

Binary IL

In this section, we apply our formalism to binary ILs at electrified interfaces. Thus, we use the extended electrochemical forces (eq ) in our multicomponent framework derived in ref (45).

Scheme of the As-Modeled Setup

The binary IL-electrolyte PYR[1,4]TFSI is subject to the negatively charged interface at the left, which causes the formation of an electrochemical double layer (EDL). Charge ordering diminishes with increasing distance from the interface (towards the right side), and the electrolyte is electroneutral in the bulk.

As discussed in section , binary electrolytes are described with the variables electric potential Φ, charge density ϱ, and center-of-mass convection velocity v. Furthermore, the electrical conductivity κ is the only independent transport parameter in this case.[45] The dynamical transport equations areHere, M± are the molar masses of the ionic species, which sum to MIL, and is the electric current relative to the center-of-mass motion, .

Solutions to eqs to 14 determine the ionic concentrations via ϱ = Fz+(c+ – c–) (charge conservation) and via the Euler equation for the volume (eq S-20).

We restrict our setup to one spatial dimension and assume that the inert electrified interface is located at x = 0. The electroneutral boundary condition ϱ(x → ∞) = 0 implies that the bulk concentration cb = c±(x → ∞) is completely determined by the total partial molar volume ν = ν+ + ν– via cbν = 1. Because binary ILs are electrically neutral, z– = −z+, and we choose z+ > 0.

We neglect viscous forces in our discussion of the EDL (∇τ = 0). Therefore, the Gibbs–Duhem relation (eq S-22) becomes c+∇μ+el + c–∇μ–el = 0, and the expression for the electric flux simplifies to , where we use the chemical potential of the anion species to determine the IL electrolyte (∇μIL = ∇μ–el),Here, we introduced the relative magnitude of the molar volumes γ± = ν±/ν. Thus, the forces given by eq depend upon the model for Fint. Furthermore, by using the Gauss model (eq ), we can either close the forces via the “complete” integral equation (eq ) or we can use the gradient expansion (eq ).

We want to apply the half-cell potential Δϕ. Because the electric potential Φ is continuous across the electrode–electrolyte interface, Φ(0) in the electrolyte is subject to the boundary conditionWithout a loss of generality, we set the electrolyte potential in the bulk to zero, limΦ = 0. Hence, Δϕ = Φ(0) is the potential applied to the electrode.

We perform one-dimensional numerical simulations of this system of eqs and 13 in the completely dissociated state, subject to an inert electrified interface (details in section S-4.1). This electrolyte is part of the IL family composed of TFSI anions and PYR cations. Because of their excellent electrochemical properties, these ILs are widely studied and used for applications in lithium ion batteries.[60] We state the electrolyte parameters in the SI (section S-4.2).

Energy Scales and Dimensions

In this section, we clarify the notation and state the nondimensional form of principal quantities appearing in our theory. For a complete discussion, we refer to the SI (section S-2), where we share the motivation for our choices.

We introduce dimensionless variables for electric potential Φ̃ = ΦFz+/RT, charge density ϱ̃ = ϱνc̃b/Fz+, and concentration c̃α = cανc̃b. As a consequence, the Euler equation for the volume becomes c̃b = γ+c̃+ + γ–c̃–, with the dimensionless molar volumes γ± = ν±/ν and the dimensionless bulk concentration c̃b. The Poisson equation suggests defining the generalized Debye lengthHere, we used F = eNA for the Faraday constant, where e is the elementary charge and NA is Avogadro’s constant, and the model ν = NAa3 for the partial molar volumes introduced above. This Debye length differs from the canonical definition by the asymmetry factors γ±.[61] However, it reproduces the textbook definition for symmetric ions (γ± = 0.5). LD becomes minimal for γ± = 0.5 because the mixing entropy of a binary electrolyte is extreme for equal ion size. Thus, asymmetry increases the Debye screening length.

With this length scale, we nondimensionalize our grid, viz, x̃ = x/LD and ∇̃ = LD·∇ and obtain the dimenionless Poisson equation:

In the SI (section S-2), we nondimensionalize the transport equations (eqs and 14) for binary symmetric ILs. Because we neglect convective effects in our EDL discussion, the complete set of equations consists of the Poisson equation and one transport equation for the charges. By substituting eq into eq , we find for the integral descriptionwhere χ = (γ–M+/MIL – γ+M–/MIL)/c̃b measures the “asymmetry” of the ion species and ∂ = (ε0εR/κ)∂.

The interaction potential is nondimensionalized (see eq S-40) by two energy scales for thermal energy Eth and electrostatic energy Eel,such that In the case of symmetric ions γ± = 0.5, these energy scales take the textbook form for thermal energy and Coulomb energy of charges at distance a. Apparently, both energy scales are coupled by the generalized Debye length LD,

The integral form (eq ) for the transport equation allows us to relate our continuum framework to MD simulations as discussed in section .

In this article, we restrict the gradient expansion of the interaction to the trivial and first nontrivial modes [ and , see eq S-41] and obtain

Stationary State

Because our focus is the formation of equilibrium structures, we discuss the system of equations in the stationary limit. This allows us to integrate the differential equations using electroneutral boundary conditions, which results in a simplified description susceptible to analytical techniques.

Stationarity (∂ϱ = 0) implies that all fluxes are constant. Here, we have no flux conditions , which implies that both species are in equilibrium, μ+el = μ–el, i.e.,

Thus, the stationary state for the binary electrolyte is described by the Poisson equation and eq . Here, we evaluate the equilibrium condition using the gradient description (eq ) in the nondimensionalized form (eq S-41 and eq ) and integrate the result using electroneutral boundary conditions in the bulk,

Apparently, in contrast to the dynamical case where electrolyte momentum is important, the molar masses appearing as parameters in fluxes v and become irrelevant in the stationary limit. Instead, the relative magnitude of the molar volumes γ± enters the system of equations. This highlights the principal role of molar volumes as parameters in the stationary state and is a consequence of the Euler equation for the volume (eq S-20).

For completeness, we state the integral transport equation (eq ) in the stationary limit (section S-3.1)

Small and Large Potentials

Equations and 25 (or eq ) constitute the complete set of equations necessary to describe a binary IL electrolyte in a stationary state. In sections and 4.3, we solve these equations using numerical methods. Our goal is to supplement these numerical methods by an analytical examination of the gradient description. However, the analytical solution of the gradient description is hindered by the higher-order gradients appearing in eq and by the different prefactors of the logarithmic terms (in general, γ+ ≠ γ–). Therefore, we distinguish different limiting cases in our analysis in sections and 4.1. In the SI, we describe the special case of symmetric ion species (sections S-3.2 and S-3.3.3).

In section , we show that the limiting case of small charge densities, ϱ̃ ≪ 1, is useful. In this case, we can expand the logarithmic terms in eqs and 26 around the electroneutral state,such that eqs and 25 becomewhere ε̂R is defined as the dielectric operatorIn the absence of molecular repulsion, , the dielectric operator reduces to the canonical, scalar-valued dielectric parameter ε̂R → 1.

Furthermore, quantities similar to ε̂R also arise in the liquid-state theory of classical statistical mechanics. This expansion corresponds to a small wave vector expansion of the dielectric function expressed as a correlation function of the molecular dipole densities, e.g., ref (50).

In the following sections, we show that the gradient expansion, eqs and 29, allows significant insight into the competing effects of interactions , Eel, and Eth and predicts the EDL structure as a function of the temperature, dielectricity, ion size, and interaction strength.

Mean Steric Effect: Charge Saturation

In this section, we neglect nonlocal interactions, , and discuss the EDL structure of the electrolyte due to bulk effects alone based on Fb (eq ). In this way, we reveal the competition between Coulombic ordering and entropic disordering, i.e., diffusion.

Toward this aim, we consider the system of equations constituted by the Poisson equation (eqs ) and eq subject to ,

First, in section , we solve this system of equations numerically. We supplement this investigation with an analytical analysis and focus on the two limiting regimes of large and small electric potentials. In section , we discuss the case of Φ̃ ≪ 1, and in section , we discuss the case of Φ̃ ≫ 1. For the special case of symmetric ion species (γ± = 0.5), we derive analytical solutions for the electric field Ẽ(Φ̃) and for the charge density ϱ̃(Φ̃) as functions of the electrolyte electric potential in the SI (section S-3.3.3).

Simulations

Figure shows numerical results for the system of eqs and 25. Figure a,b illustrates screening profiles of the electric potential, the charge density, and the ion concentrations for varying electrode potentials Δϕ.

Figure 1

Simulation results of the EDL structure perpendicular to the electrode–electrolyte interface for a binary IL (eqs and 23). If not mentioned otherwise, T = 300 K, εR = 15, and Δϕ = −0.1 V. (a) Profiles of the electric potential and charge density (inset) for different electrode potentials Δϕ. (b) Concentration profiles of the anions and cations for different electrode potentials Δϕ. (c) Concentration profiles for different volume ratios γ+. The inset shows the corresponding electrical potential. (d) Concentration profiles for varying dielectric constants (dashed lines) and temperatures (solid lines).

Apparently, the application of a negative electrode potential (Δϕ < 0) polarizes the electrolyte. The electric potential (Figure a) is continuous across the interface and decays smoothly toward the electroneutral bulk region. The inset of Figure a shows that, for low electrode potentials, the charge density decays exponentially. Similar behavior can be observed in Figure b for the concentrations. The concentration of positive counterions increases toward the interface, whereas negative ions get depleted. Apparently, the electrolyte screens the electrode potential by the accumulation of counterions. However, above Δϕ ≈ −0.05 V, the counterion concentration saturates near the interface. A further increase of Δϕ broadens the EDL.

This behavior can be explained by the mean volume effect. The application of a negative potential Δϕ implies that positive ions accumulate near the interface and negative ions are depleted. However, the Euler equation for the volume, eq S-20, implies the saturation concentration c+sat = 1/ν+. Once the accumulated species reaches this saturation, the screening mechanism transitions from increasing the concentration at the interface to broadening the width of the EDL. The simulated EDL approaches a thickness of 0.6 nm at Δϕ ≈ −0.05 V, which is significantly wider than predicted by the canonical Debye–Hückel theory with the Debye length LD ≈ 0.7 Å (eq ). This phenomenon is typically denoted as “crowding”.[35]

Because the saturation concentration depends upon the molar volume, cαsat = 1/να, the partial molar volumes directly affect the screening behavior. Figure c shows numerical results for the ionic concentrations at different volume ratios γ± = ν±/ν (in which v is kept fixed). The EDL is thinner for smaller γ+ because this allows for a tighter packing of cations.

The effects of temperature T and dielectricity εR on the EDL structure are illustrated in Figure d. The screening is more effective for smaller values of εR, and the EDL width increases with increasing magnitude of εR. This is in qualitative agreement with the screening behavior for dilute solutions as predicted by the Debye–Hückel theory. Likewise, the EDL becomes more diffuse with increasing temperature because of the disordering effect of thermal motion. The observed effects of T and εR highlight the competing interplay between the electrostatic effect of charge ordering and the disordering effect of entropy.

To summarize, the simulations show two distinct regimes of EDL structure. First, for large electrode potentials, where Δϕ̃ ≫ 1, the charge is saturated near the interface. Second, near the electroneutral bulk region, where Δϕ̃ ≪ 1, the charge density decays exponentially. These two distinct EDL structures, charge saturation and exponential decrease, correspond to two disjointed electrolyte regimes, |Φ̃| ≫ 1 and |Φ̃| ≪ 1, respectively.

Asymptotic Analysis

The simulation results in section motivate our procedure for analyzing the EDL structure. First, we study the EDL far away from the electrode close to the electroneutral bulk (large x̃). For this purpose, we expand the stationary equations around the electroneutral bulk for small charge densities ϱ̃ ≪ 1. According to eq , this regime corresponds to small dimensionless potentials |Φ̃| ≪ 1. Note that this coincides with the high-temperature regime because Φ̃ = Φ(Fz+/RT) (section ).

Second, we analyze the behavior close to the electrode (small x̃) at large electrode potentials |Δϕ̃| ≫ 1, where the electrolyte potential satisfies |Φ̃| ≫ 1. This corresponds to low temperatures (section ).

Asymptotic Analysis: Small Potentials |Φ̃| ≪ 1

As outlined above, we begin our analytical examination of the EDL in the limit of small dimensionless potentials |Φ̃| ≪ 1. Our idea is to approach the EDL from the electroneutral bulk region along the direction of decreasing x̃. Thus, we use the expansion of ionic concentrations around the bulk electrolyte in eq , c̃± = c̃b ± ϱ̃γ∓ ≈ c̃b.

For this aim, we insert ϱ̃ from eq into eq for ε̂R = 1, yieldingWith the boundary conditions discussed above, Φ̃(0) = Δϕ̃ = Fz+Δϕ/RT and Φ̃(x̃ → ∞) = 0, we obtain the solutionThe corresponding dimensionalized electrolyte potential,decays exponentially with damping parameter 1/LD. Thus, the decay length in the limit |Φ̃| ≪ 1 is the Debye length LD (eq ).

In Figure , we compare the analytical predictions for this limit (dashed green lines) with the numerical results (blue lines) for different electrode potentials. Apparently, the analytical and numerical results for the electric potential and the charge density are in excellent agreement for small electrode potentials Δϕ̃ ≈ – 0.4 when the condition |Φ̃| ≪ 1 is fulfilled everywhere.

Figure 2

Comparison of the asymptotic analysis of the EDL structure with the numerical results obtained from eqs and 23 (T = 300 K and εR = 15). We consider two different values of the interface potential corresponding to the regimes of small and large potentials at Δϕ = −0.01 V (Δϕ̃ = −0.4) and Δϕ = −0.35 V (Δϕ̃ = −13.5). (a) Profiles of the electrical potentials as predicted by the analytical approximation (eqs and 36). The inset highlights the region close to the interface for the case Δϕ = −0.01 V. (b) Profiles of the charge density as predicted by the asymptotic analysis in eqs and 38. The inset compares numerical and analytical values for the total charge in the EDL (section S-3.3.2 in the SI).

In the SI, we derive the expressions for the total EDL surface charge density and the differential capacitance (section S-3.3.1).

Asymptotic Analysis: Large Potentials |Φ̃| ≫ 1

Next, we discuss the EDL in the limit of large potentials |ϕ̃| ≫ 1. This regime can be found for large electrode potentials ΔΦ̃ ≫ 1 close to the electrode/electrolyte interface. Because Φ̃ = Fz+Δϕ/RT, this analysis is exact in the limit of zero temperature, T = 0.

In this case, the logarithmic terms in eq must compensate for the diverging potential term Φ̃. Because of the mean volume constraint (eq S-20), one of the logarithmic terms is diverging if one species is depleted and the other species saturates, ϱ(x = 0) = Fzα*/να*. Here, we denote the saturating species by the index α*. Because electrical potentials are continuous across interfaces, the saturation species α* is uniquely determined by the sign of the electrode potential, sign(zα*) = −sign(Δϕ̃).

Therefore, ϱ̃sat = −sign(Δϕ̃)c̃sat solves eq , where c̃sat = c̃b/γα*. Upon integration of the Poisson equation (eq ), we findwhere the width of the EDL depends on the electrode potential . Thus, the dimensionalized EDL length isApparently, the decay length increases with ion size because of the mean-volume effect. Also, it scales with the asymmetry ; i.e., it is small for small screening species. Comparison with the Debye screening length shows that LEDL > LD in the limit of small temperatures T or large potentials |Δϕ̃| ≫ 1.

In the limit of vanishing temperature T = 0, the charge profile is box-shaped and is determined by the screening length LEDL,where θ is the Heaviside function. In the SI, we calculate the analytical expressions for the total charge in the EDL and for the differential capacitance in this limit (section S-3.3.2).

In Figure , we compare the analytical predictions for this limit (dashed yellow lines) with the numerical results (solid red lines) for Δϕ̃ = −13.5. The box-shaped charge profile is in good qualitative agreement with the numerical results because it almost predicts the correct width of the EDL. However, the transition from saturation to the bulk state is more diffuse in the numerical profile. This is due to the entropic, thermal influence, which “washes out” the box. Nevertheless, the inset shows that the charge in the EDL, as predicted by the analytical approximation, is quantitatively in good agreement with the numerical results. We note that this profile and its temperature dependence are reminiscent of the Fermi distribution.

Nonlocal Interactions: Charge Oscillations

In this section, we discuss the influence of nonlocal interactions () on the structure of the EDL. As in section , we discuss the two limiting cases of small and large potentials separately.

Static Asymptotic Analysis

Asymptotic Analysis: Large Potentials |Φ̃| ≫ 1

Let us first discuss the regime of diverging electrolyte potentials |Φ̃| → ∞. In this limit, the interaction contribution cannot compensate for the diverging electrolyte potential in eq . The logarithmic terms are diverging if one species is depleted, and we recover the same results as described in section for the case of vanishing interaction contributions .

Asymptotic Analysis: Small Potentials |Φ̃| ≪ 1

In this section, we consider the full theory with molecular repulsion in the regime of small potentials |Φ̃| ≪ 1. As outlined above, we expand the interaction free energy in gradients of the charge density and restrict our analysis to the first two perturbation modes n = 0 and n = 1. With the assumption of small charge densities ϱ̃ ≪ 1, we derived the linear equations above (eqs and 29), which we rephrase in matrix form aswhereWe solve eq via the eigenvalue decomposition with the eigenvalues α̃1,2 of the matrix in eq . These are determined by the relative magnitudes of the three competing energy scales, Eth, Eel, and ,Each eigenvalue α̃ gives rise to a dimensionless wave vector,These determine the general solution of eq together with the eigenvectors ãα̃ = (ã1α̃, ã2α̃) = (−α̃, 1),The expansion coefficients A are determined by boundary conditions and physical arguments. Apparently, the corresponding wave vectors are functions , which determine the structure of the EDL,Thus, the EDL structure depends upon the relative magnitudes of the energies Eth, Eel, and via eq . In particular, the classification of k̃± depends upon the sign of the rootappearing in eq . Thus, the critical values , defined by the condition determine the thresholds for the transition between the phases of the EDL structure,Thus, with eq we can draw the phase diagram for the EDL structures. Because , there are three phases. In the SI, we discuss each case in great detail (section S-3.4.2). Next, we give a short description of each phase.

Phase 1:. In this regime, α̃1,2 ≥ 0, which implies a real-valued wave vector. Thus, this phase corresponds to exponentially damped profiles . A harmonic analysis of the root reveals that (section S-3.4.3)Thus, solutions with the damping parameter k̃1 vanish quickly and are rendered as unphysical, whereas the limit of vanishing interactions for k̃2 reproduces the bulk expansion for Φ̃ ≪ 1 from section (eq ).

Phase 2:. In this regime, , and thus the root (eq ) becomes complex. Therefore, the wave vector has nonvanishing real and imaginary parts, . This corresponds to charge profiles of exponentially damped oscillations,

Phase 3:. In this regime, both eigenvalues are real but negative, α̃1,2 ≤ 0. Therefore, , which corresponds to undamped oscillatory profiles. The limiting case for indefinitely strong interactions yieldsThus, the result for reproduces the experimental findings obtained by AFM measurements that the wavelengths λ ≈ a/2π of the observed oscillations scale with the size of molecules a.[52,62] Apparently, the incompressibility of ions in our model prevents a further decrease of the wavelength.

Thus, the critical values constitute exactly the boundaries between the different EDL phases.

In Figure , we illustrate the phase space of EDL structures as functions of temperature and dielectricity. (In section S-6.1, we also show the phase space as a function of ion size, ion asymmetry, and valency.) Apparently, three distinct phases of EDL structures are present. The exponentially damped EDL phase corresponds to the regions below (red line), whereas the damped-oscillatory EDL phase corresponds to the regions between the blue and red lines. Finally, the undamped oscillatory EDL phase corresponds to the regions above (blue line). Figure a illustrates that temperature T and hardcore interactions are in competition and that the critical interaction strengths increase with increasing temperature T, i.e., thermal energy Eth.

Figure 3

Phase spaces for EDL structure as functions of T and εR for equally sized ions (γ± = 0.5); see eq . (a) Critical interaction energies as a function of temperature (here, εR = 15 and a = 1.3 nm). (b) Critical interaction energies as a function of dielectricity (here, T = 300 K and a = 1.3 nm). Three phases are present: exponentially damped charge density (shaded red), decaying oscillatory charge density (shaded gray), and quasi-crystallinity (shaded blue).

Figure b reveals the influence of dielectricity, i.e., electrostatic forces, on the EDL phases. Apparently, the damped oscillatory phase becomes narrower for ILs with larger dielectricity εR, i.e., smaller electric energy Eel.

We note that the exponentially damped regime for small interaction strengths (small compared to electrostatic and thermal energy) corresponds to the EDL structure found in section in the absence of hardcore interactions. However, as can be inferred from Figure a, this phase is hardly present for reasonable temperatures.

Figure shows the nondimensionalized wave vector as a function of the relative energy scale (where T = 300 K, εR = 15, and a = 1.3 nm). For small interaction energies the wave vector is real, , which corresponds to exponentially damped profiles. In particular, the static profile at reproduces the case described in section where the exponential profile is determined by the scales Eth and Eel alone, i.e., k̃ = 1. Apparently, increases with up to the threshold , beyond which it starts to decrease. Thus, the EDL has minimal extension at . This suggests that the increasing strength of the repulsive ion correlations compresses the screening layer. Once the hardcore potential exceeds , the system overscreens, i.e., the ion layers begin to oscillate. The damping parameter Re(k̃) vanishes exactly when , i.e., when the system transitions into nanosegregation of the ion species. Interestingly, the frequency of the oscillations Im(k̃) exhibits a local maximum and minimum in the regime of damped oscillations. Furthermore, Im(k̃) attains its maximal value in the limit of prevailing interaction strength .

Figure 4

Real and imaginary parts of the nondimensionalized wave vector k̃ = kLD as a function of the relative magnitude of the energy scales and Eth (eqs and 42). Here, T = 300 K, εR = 15, and a = 1.3 nm.

In the SI (sections S-3.4.4 and S-3.4.5), we investigate the influence of individual perturbation modes Γ120 and Γ122 on the phase diagram. As it turns out, neglecting all but the zeroth-order correction Γ120 results in a binary phase diagram comprising only exponentially damped profiles and undamped oscillatory profiles. In contrast, neglecting the zeroth-order correction and taking only the first nontrivial order Γ122 into account results in a binary phase diagram comprising only exponentially damped profiles and damped oscillatory profiles. This is the case for MFTs based on the BSK framework. Thus, for the “complete” set of the three different phases, both perturbation modes Γ120 and Γ122 are necessary.

Interestingly, for the pathologic case of negative interaction strengths , the phase space reduces to the two screening phases of exponentially damped profiles and undamped oscillatory profiles. This follows straightforwardly from eq (and also from the discussions in sections S-3.4.4 and S-3.4.5).

Dynamic Asymptotic Analysis: Linear Stability Analysis

In this section, we complement the static analysis of section by an analytical analysis of the dynamic transport equation in the gradient description (eq ).

For this purpose, we perform a linear stability analysis and consider the limit of small potentials, |Φ̃| ≪ 1. Thus, the logarithmic terms can approximated as in eq , and eq becomesWe expand the electric potential around a uniform bulk state Φ̃b,Here, the equilibrium state is determined by the electroneutral bulk condition Φ̃b = 0 and ϱ̃b = 0. Thus, the first order perturbation takes the formHere, the wavenumber k̃ determines the spatial distribution of the dimensionless perturbation ϵ̃1 ≪ 1, and the parameter s̃ measures the temporal growth rate of this perturbation.

We restrict our analysis to probing the linear stability and substitute eq and the Poisson equation into eq . Next, we collect terms up to the first order in perturbation mode ϵ̃1, which yields a dispersion relation for the growth rate of the perturbation,

The uniform state is stable under perturbation if and only if s̃ < 0. This defines an instability onset k̃c for the wavenumbersThe corresponding stability criterion s̃(k̃1,2c) < 0 determines the phase boundary at which the bulk of the IL electrolyte becomes unstable. This stability threshold exactly equals the phase boundary between the damped oscillatory phase and the nanosegregated phase (eq ),

Thus, for interaction energies the bulk state of the system becomes unstable and phase separation emerges. The initial cause of the structure formation can be driven by external agents or boundary conditions, e.g., by the application of an electrical potential to an IL/electrode interface.

This stability analysis complements the static analysis and rationalizes the emergence of phase separation into ionic layers occurring at interaction energies above .

Validation with Simulation

Our goal in this section is to compare the results of our asymptotic analysis (eqs to 43 in section ), with numerical simulations of the completely coupled system subject to the two theoretical descriptions (integral description, eqs and 19; gradient description, eqs and 23).

We start our numerical investigations with an overview of the screening profiles for the charge density at different values , as obtained from the integral description in eqs and 19 (Figure ). Next, we compare our different EDL descriptions in detail for two different energies (Figures and 7). Finally, we generalize these exemplary findings via a systematic study over the complete phase space of interaction energies (Figure ). This provides a clear illustration of the complete set of phase transitions which the system undergoes and highlights the consistency among the three descriptions.

Figure 5

Screening profiles of the charge density ϱ̃ as obtained from numerical simulations of the integral description (eqs and 19) for different values . The y axis is scaled from −1 to 1, where |ϱ̃| = 1 corresponds to charge saturation.

Figure 6

Results for the electric potential and charge density as obtained from numerical simulations of the integral description (eqs and 19) and of the gradient description (eqs and 23) and as predicted by the analytical predictions (eqs to 43) at .

Figure 7

Screening profile for the charge density ϱ̃ obtained from numerical simulations with respect to the gradient description (eqs and 23) and according to the analytical description (eqs to 43) at .

Figure 8

Meta analysis of the interfacial profiles for four thousand simulations. The dashed vertical yellow lines show the phase boundaries (eq ). The dashed and solid red/blue lines show the peak variance of the complete set of simulations as defined by eq with respect to the integral description (eqs and 19) and with respect to the gradient description (eqs and 23). The left inset shows the onset of the oscillations at small interaction energies. The right inset shows the variance in a nonlogarithmic setting, which highlights the occurrence of phase transitions.

All simulations were performed for a symmetric cell setup, where the IL electrolyte is located within two oppositely charged, blocking interfaces separated by a distance of Lcell = 60 nm. The electrode on the left side is negatively charged with Δϕ = −100 mV, whereas, on the right side the electrode is positively charged with Δϕ = 100 mV. Because charge saturation begins roughly at Δϕ = 70 mV (Figures and 2), the charge distribution can safely be assumed to be saturated adjacent to the interfaces, i.e., |ϱ̃| = 1. The electrolyte is considered to consist of symmetric ions (γ± = 0.5) of size a = 1.3 nm. Hence, the cell geometry allows a “maximal” number of roughly 90 ions. In addition, we assume room temperature, T = 300 K, and εR = 15. The phase boundaries corresponding to these parameters as predicted by our analytical description are and (eq ).

Figure shows the numerical results of the charge density for the integral description (eqs and 19), where takes values across two orders of magnitude. First, at the profile shows charge saturation near the two electrified electrodes, ϱ̃(x = 0) = 1 and ϱ̃(x = Lcell) = −1. Near both electrodes, the profile decays exponentially toward the electroneutral bulk (ϱ̃ = 0). This corresponds to the profiles which we discussed in great detail in section . Because , this is in accordance with the analytical prediction. The next two profiles show results for interaction energies within the intermediate phase, . Both simulations show damped oscillatory profiles, where the long-ranged oscillatory profiles span many nanometers. Apparently, the oscillations in the profile for extend across almost the entire cell. A slight increase of 2 meV to causes the profile to transition to a crystalline phase with undamped oscillatory shape. Note that the amplitudes between the electrodes are smaller than unity, i.e., the bulk region consists of mixed ion layers with one dominant ion species. An increase to enhances the amplitudes of the oscillations further, i.e., enhances the segregation of ion species. The last plot shows the corresponding profile for a significantly enhanced interaction energy (). Here, the amplitudes of the oscillations have reached saturation (|ϱ̃| = 1), and the electrolyte has transitioned into a crystalline phase consisting of alternating pure ion layers. In Figure S-2 (section S-6.2), we highlight that the ionic layers coincide exactly with the ion size a/2. Thus, with increasing energy , the interfacial structure increases into the bulk electrolyte until the bulk itself gets nanostrucured by the layering of the ion species. This phase transition occurs rapidly within a few meV.

Apparently, the numerical results for the integral description confirm the existence of three different screening phases. However, quantitative deviations between our descriptions is present. As we show in the SI (Figure S-4) the phase transitions from exponential decay to damped oscillations occurs roughly at . In addition, as can be inferred from Figure , the transition from damped oscillations to undamped oscillations appears at . Hence, both phase boundaries are slightly shifted to smaller values compared to analytical predictions and (eq ). Thus, the analytical prediction, which is based on the gradient description, slightly underestimates the influence of when compared with Eth and Eel. This can be attributed to the fact that the gradient description is an approximation based on only the first two perturbation modes, whereas the integral description comprises all modes.

Next, we give a quantitative comparison between the numerical results of the two theoretical descriptions and the analytical predictions. Here, we restrict our discussion to the interaction energy , i.e., the intermediate phase of damped oscillations. Figure shows the profiles for the charge density and electrolyte electric potential as obtained from the numerical simulations and as predicted by the analytical description for the first few nanometers of the left half-cell. Figure a illustrates the charge distribution adjacent to the negatively charged electrode. The dashed blue line shows the screening profile obtained from the gradient description, which exhibits a damped oscillatory shape. This confirms the analytical prediction for this interaction energy. The dashed yellow line shows the resulting analytical profile. Note that the analytical prediction in section does not capture charge saturation but determines only the damping parameter and the oscillation frequency of the screening profile. However, in section , we derived an analytical prediction for the saturation width LEDL, which is valid close to the interface (eq ). Hence, to reconstruct the “complete” profile, we supplement the contribution emerging from the bulk (eq ), which is valid far away from the interface, by constant charge saturation ϱ̃ = 1 spanning over the width LEDL. Apparently, the analytical and numerical results of the gradient description are quantitatively in very good agreement. Finally, the solid red line in Figure a shows the numerical results for the integral description. In accordance with the results shown in Figure , these results reproduce the analytically predicted screening phase, but the oscillations are more pronounced. Hence, the influence of the interaction energy is more dominant in the integral description than in the gradient description. Next, in Figure b we show the profiles for the normalized electrolyte electric potential. The dashed blue line shows the profile due to the gradient description. It is in accordance with the charge profile shown in Figure a. Again, we reconstruct the analytical profile by supplementing the profile (eq ), which is valid close to the interface, by the profile (eq ), which is valid toward the electroneutral bulk. Apparently, the analytical results are quantitatively in very good agreement with the results stemming from the gradient description. The red line shows the profile as obtained from the integral description. As in Figure a for the charge density, the oscillations are slightly enhanced when compared with the gradient description. In addition, the brown solid lines show the analytical envelopes for the screening. Apparently, it captures both numerical results qualitatively very well. Interestingly, the differences in electrolyte potential among the three descriptions depend on the electrode potential near the electrode as shown in Figure S-3 (section S-6.3). However, the qualitative agreement is independent of the boundary conditions.

In Figure , we show results for the charge distribution at enhanced interaction energy , i.e., close to the phase boundary . As can be inferred from Figure , the integral description has already transitioned to the phase of undamped oscillations for this interaction energy. Hence, we show only the screening profile as obtained from the numerical simulation of the gradient description (solid blue line) and compare it with the analytical prediction (yellow line). Apparently, in accordance with the analytical prediction, the numerical profile has a damped oscillatory shape, where the oscillations extend over roughly 40 ion sizes. This highlights the influence of the enhanced interaction energy (see also Figure ). Overall, the analytical profile shown here is in nice agreement with the numerical results.

Finally, we conduct a quantitative comparison between the two EDL discriptions and the analytical description across multiple orders of magnitude of . To address this goal, we examine the simulation results of roughly 4000 EDL simulations across the parameter range from 0.1 meV up to 500 meV. As above, we apply Δϕ = ± 100 meV at the electrodes such that we can safely assume charge saturation near the interfaces.

We evaluate the simulation results by extracting two characteristic properties. First, we count the number of peaks appearing in each screening profile. Because of charge saturation, a minimal number of two peaks always occur. At most, roughly 90 ion layers fit into the cell geometry of length Lcell = 60 nm. We present the number of peaks occurring in the full cell as a function of the interaction strengths in the SI (Figure S-4 in section S-6).

However, beyond the number of peaks, we also want to evaluate the peak amplitudes. For this purpose, we investigate the peak variance of the left half-cell, defined byHere, x is the discrete location of the ith peak ϱ̃ = |ϱ̃(x)| (such that 0 ≤ ϱ̃ ≤ 1, where ϱ̃ = 1 corresponds to a saturated peak, i.e., a pure ion layer) appearing in the profile of the charge density. In the SI (section S-6.4), we show analytically that σ converges to Lcell/√3 if the set of simulation energies comprises energies . For such interaction energies, the bulk electrolyte has transitioned to a crystalline phase composed of nanosegregated ion layers (Figure ).

Figure shows the results for variance σ normalized to its maximum Lcell/√3 on a logarithmic scale. In this figure, the vertical dashed yellow lines indicate the phase boundaries , as predicted by the analytical description (eq ). The left inset shows the simulation results for small values , and the right inset comprises the overall results in a nonlogarithmic representation highlighting the transition. The blue dashed line shows the results for σ according to the gradient description (eqs and 23). At small interaction energies , the variance is zero. This corresponds to an exponentially damped, nonoscillatory screening profile. (Note that the only peak, due to charge saturation, is located at x̃ = 0). The variance starts increasing exactly at (left inset). This corresponds to an increasing number of damped oscillations, where the amplitudes of the peaks also increase with . Finally, at , the variance converges to it is constant limiting value Lcell/√3 (right inset). In this energy regime, the bulk electrolyte consists completely of ion layers. Altogether, these results reproduce the phases exactly as predicted by the analytical description.

The red curve shows the results for σ according to the integral description (eqs and 19). In contrast to the gradient description, the variance starts increasing from zero at roughly , i.e., before the analytically predicted phase boundary (left inset). Hence, the phase transition from exponentially damped screening profiles to damped oscillatory screening profiles is slightly shifted to smaller energies. Next, the variance increases exponentially up to roughly , above which it transitions to the constant limiting value Lcell/√3. Altogether, the phase boundaries of the integral description still exhibit qualitatively good agreement with the analytical predictions, although being slightly shifted to smaller values. Apparently, this behavior is due to the cumulative effect of the integral term in eq , which comprises all interaction modes. In contrast, we consider only the first two modes (n = 0 and n = 1) of the gradient expansion in eq .

Multiscale Methodology

In this section, we highlight the relation of our model to theories on smaller and larger length scales. We discuss in section , on the basis of basic concepts from liquid-state theory, how atomistic simulations can directly parametrize our theory. Next, in section we sketch the phenomenologic BSK continuum approach for the description of ILs near electrified interfaces and illustrate its relation to our work. In addition, we state the relation of our framework to AFM experiments in section S-5 (ref (52)).

From Molecular Dynamics to Nonequilibrium Thermodynamics

Here, we explain how the parameters of our continuum theory can be rigorously calculated with quantum chemistry, i.e., DFT and MD.

Ab initio DFT calculations predict the forces between ions and molecules by calculating their electronic structure. The DFT-generated force fields are the focal quantity for MD simulations,[63] which calculate the classical trajectories of ions and molecules. Results from MD simulations are often interpreted via profiles of the radial distribution function g(r).

Liquid-state theory[50] connects this atomistic description to thermodynamic concepts and scattering experiments.[64] On the one hand, the radial distribution function allow a straightforward comparison with the structure factor S from scattering experiments.[65,66] On the other hand, the density distribution function g(r) can be used to calculate different correlation functions. Subtracting its asymptotic value follows the so-called total correlation function used in integral equation theories (IETs), h(r) = g(r) – 1.[3] In IETs, the pairwise total correlation function h relates to the direct correlation function c(2), used in classical density functional theory (cDFT), via the Ornstein–Zernike relation,[67]In cDFT, the direct pair correlation functions c(2) account for pairwise interactions between two ions of species α and β, i.e., the excess free energy due to pairwise ion interactions.[68] Thus, they can be obtained via the 2-fold functional derivative of Fint,[64] i.e., via our interaction potential (eq ),

To summarize, DFT determines force fields for MD, MD determines g(r) for liquid-state theory, g(r) determines c(2) via the Ornstein–Zernike relation, and c(2) determines Fint and generates our nonequilibrium thermodynamic theory.

The dynamic properties of our theory can also be determined from atomistic simulations. These dynamic properties are encoded in the Onsager coefficients,[45] which can be measured experimentally.[69] The Onsager coefficients can be determined by MD simulations (“Green Kubo relations”).[70−72]

From Nonequilibrium Thermodynamics to Phenomenologic BSK Theory

Now we compare our thermodynamically consistent continuum approach with the phenomenologic theory proposed by Bazant, Storey, and Kornyshev (BSK), a seminal MFT approach for ILs near electrified interfaces.[35] In their continuum model of the EDL, BSK incorporates ion correlations using a modified linear dielectric relation between electrostatic fields D̅ and E, where is their dielectric operator. The second-order gradient term in ϱ̅ accounts for nonlocal ion interactions, being effectively short-ranged with correlation length . This Ansatz yields a modified Poisson equation, . The chemical potential connects the electric potential and charge density. Finally,holds in the limit of small potentials Φ̃.

Our model conceptually differs from BSK theory. Because we incorporate electrostatic correlations in the free energy, nonlocal ion interactions enter the set of equations via the chemical potentials. This implies that the MFT quantities appearing in the BSK description, D̅ and ϱ̅ = ∇D̅, differ from the corresponding quantities ϱ and D appearing in our formalism. In contrast to the “mean field charge density” ϱ̅, the charge density ϱ relates to the bulk quantity D, which does not incorporate ion correlations.

Despite these differences, the resulting model equations are very similar. This can be seen as follows. Equations and 29 for the limit of small potentials can be cast into one equation for the electric potential alone,where the dielectric operator ε̂R is defined in eq . Noting the conceptual similarity between the dielectric operators ε̂R and , the similarity between our model and BSK theory becomes apparent. In this way, we give physical meaning to the correlation length in BSK theory and outline its calculation.

Finally, we emphasize that the higher-order gradient terms, which are phenomenologically incorporated in the BSK approach, emerge naturally within our rigorous continuum model. In particular, they merely constitute the limiting case for small potentials of the more fundamental integral formulation (eq ). Furthermore, in contrast to the phenomenological BSK model, our order expansion also comprises a zero-order correction in the dielectric operator (eq ). This mode is mandatory in realizing the complete phase space of interfacial profiles (sections S-3.4.4 and S-3.4.5 in the SI).

Outlook

In this section, we discuss how our framework can be extended to account for additional microscopic IL effects.

In this work, we have supplemented our bulk description for ILs and highly correlated electrolytes, as recently presented in ref (45), with nonlocal interactions. Furthermore, we applied the resulting framework for the case of short-ranged hardcore interactions. However, the generality of our framework based on the modeling of the free energy offers the possibility to incorporate a wide range of nonlocal effects into our framework.

This includes properties such as ion asymmetry, ion geometry, polarization, and charge delocalization, which have a significant influence on the structure of ILs near electrified interfaces.[62,73−76] These effects partially result from the relative orientation between the ions, which makes a one-dimensional approach challenging. Nevertheless, assuming a highly symmetric setup, the one-dimensional description might still capture some basic consequences of these effects.

Similar to the force fields used in atomistic simulations, the short-ranged repulsive interaction can be supplemented by a longer-ranged attractive tail, taking account of higher-order electrostatic effects of the van der Waals type, or larger ions with complex geometry, i.e., long alkyl chains.[53,54,63] Also, by refining the short-ranged repulsive interaction potential, more detailed models for the ion geometry and ion asymmetry can be incorporated into our model. However, the strong influence of these microscopic properties on the EDL structure may lead to some novel features within our framework. For example, the three energy scales which determine the screening profile might transition to field quantities which exhibit spatial variation. Also, the phase space of screening profiles might become higher-dimensional, which can lead to a more complex set of phase boundaries allowing for “mixed” screening types.

Nontrivial polarization effects could be incorporated into our linear constitutive model for the coupling between the electric field and the dielectric displacement. This would result in a spatially varying dielectric function εR(x) appearing on the electrostatic energy scale Eel (eq ) and a direct coupling of the chemical potentials with the ion polarization. For small dielectric perturbations, we hypothesis that the electrostatic energy scale becomes more diffuse, which has the effect that the phase boundaries between the screening profiles wash out. Only in the case of large dielectric variations do we expect the phase space of screening profiles to be altered significantly.

Our dynamical theory offers the possibility to investigate transport processes occurring in electrochemical devices, e.g., the influence of EDL charging on the electrolyte performance, or the influence of the EDL structure on the electrode-transfer kinetics. However, electrochemical devices have some characteristic properties which must be carefully taken into account when they are modeled. For example, overlapping double layers in nanoporous electrodes could be taken into account.[77]

In our description, we assumed an ideally flat surface, which can be a bad approximation for many electrochemical systems.[78] The influence of interface roughness on the EDL structure can be modeled by modifying the entropic contributions in the free-energy functional.[79] In our analysis, this would alter the thermal energy scale Eth (eq ). Depending upon the surface morphology, this would enhance the disordering effect of the thermal energy on the EDL structure. As a result, the formation of crystalline phases might become suppressed at rougher surfaces, similar to increasing the temperature.

Conclusions

In this work, we complement our thermodynamically consistent continuum framework for IL electrolytes by nonlocal molecular repulsion. Our integral formulation can be determined by ab initio MD simulations. Assuming short-ranged interactions, we expand the interaction free energy in concentration gradients and adjust the dynamic equations for transport. The resulting equations connect to the phenomenologic approach of BSK theory. We validate our approach by simulations and find remarkable agreement between the different variants of our theory.

In this way, we develop a predictive multiscale approach to the theory of ILs at electrified interfaces. Atomistic density functional theory calculations parametrize MD simulations, MD simulations yield an integral formulation for molecular repulsion in our thermodynamically consistent transport theory, and our theory can be expanded to give the phenomenological BSK theory.

The expanded continuum approach allows us to perform analytical asymptotic analysis which creates deeper insights into the parameter dependence of the EDL structure as we demonstrate for the example of binary ILs. First, we have neglected molecular repulsion. We can analytically describe both limits: the dilute Debye limit, where charge density is exponentially decaying, and the concentrated crowding limit, where charge is saturated due to steric effects. Second, we have taken into account molecular repulsion. We discuss the structure of the EDL dependence on energy scales for thermal motion, molecular repulsion, and electric Coulomb forces and find three different phases. For small interactions, we recover the dilute Debye limit. For intermediate interactions, a multilayer structure of ions emerges which is washed out over several atom layers. For very large interactions, the analysis predicts the long-ranged, nondecaying crystalline order of the EDL. In simulations of our full theory, we eventually observe charge ordering of quasi-crystalline multilayers in this case.

In summary, we have proposed a thermodynamically consistent description of ILs at electrified interfaces that closes the gap in their multiscale understanding. This makes possible a predictive theoretical approach for tailoring ILs. We prove that the intermolecular forces determine the EDL structure of binary ILs. Future work should be extended to ternary mixtures of ILs and should incorporate the shapes of molecules into the theory.