Modulation of ionic conduction using polarizable surfaces

Switchable conductance is desirable for designing densely interconnected (neuromorphic) systems for information storage, performing complex logic operations, and executing neural network algorithms  [1, 2, 3]. To emulate neural activity, such as voltage spiking [4, 5] and synaptic plasticity [6, 7], researchers use different materials considering switchable ionic or electronic conduction [8, 9]. For example, the gate resistance tunability has been explored for neuromorphic circuits using monolayer MoS₂ multi-terminal memtransistors [9]. To expand functionality and flexibility of device design, integration of ionic and electronic conduction is an attractive option as it may allow for imitating synaptic potentiation, emulating plasticity [10, 11], and achieving neural interfacing. The coupling between ionic and electronic transport is promising for leveraging other applications such as biosensing, energy storage, and responsive materials.

Mixed ionic-electronic conductors are materials that conduct both ions and electronic charge carriers (electrons and/or holes)[12]. Recent developments combine electronic and ionic conductor materials into alternating layers of nanometric dimensions where the ionic and electronic conduction occur simultaneously [13]. The close proximity of ionic and electronic charge carriers in these devices means that the conduction behavior of one influences the other. These materials’ electronic properties can be modified by the stoichiometry [14, 15], gate biasing [9, 16], and structural changes [13]. The coupled ionic and electronic interactions and transport need special understanding beyond the comprehension of their independent behavior.

Despite numerous technological applications envisioned by using mixed ionic-electronic conduction, the lack of fundamental understanding impedes rational materials design. One essential component of mixed conduction is the effect of induced electronic polarization on ionic conduction. Previous work showed that the electronic properties of surfaces affect the nearby ions, specifically via the induced polarization charges due to dielectric mismatch on the material-electrolyte interface [17, 18, 19, 20]. A recent formalism employs the Thomas-Fermi model [21] to consider polarization effects on the ions’ transport in strong confinement by dielectrics and conductors [19]. However, the formalism is derived only for ions constrained to move in two dimensions. Using two different approaches, here we study ion conduction in strong confinement in a slit-like channel and we recover some aspects of the two-dimensional behavior predicted by the Thomas-Fermi model.

Frameworks and models that integrate the interfacial electronic polarization of materials and electrostatic molecular interactions are essential for exploiting the properties of interfaces in nanometer slit confinement [22, 18, 20]. As illustrated in Fig. 1(a), the electrostatic potential generated by a charge in the center of a narrow slit-like channel is highly dependent on the polarizability of the channel. Whereas the range of the potential is increased by dielectric confinement, the potential in strong conducting confinement is so screened that it can be considered short range. Here, we study ionic conduction in strongly confining slit-like channels and show that the ionic conductivity and the ionic adsorption are significantly impacted by changes in the confining material polarization (from dielectric and conductor), the confining distance, and the location of the polarization charges. Results from density functional theory show that the polarization plane location is a material- dependent property [23]. Hence, by adjusting the polarization plane location we consider materials of different electronic properties.

The paper is organized as follows. First, we develop an efficient method for treating ionic interactions in strong confinement by conductors and define the system parameters. We then find a nonlinear ionic conduction response, which is a prerequisite for neuromorphic behavior. Finally, we demonstrate that the location of the polarization charge impacts the ionic distribution and transport through the channel. The placement of the polarization charge offers a mechanism for modeling material-dependent properties of polarizable surfaces.

We consider a system consisting of a 1:1 electrolyte with $N_{+}$ cations and $N_{-}$ anions of diameter $d$ confined between two polarizable surfaces placed parallel to the $x$-$y$ plane and separated by a distance $L$ in the $z$-direction. The (periodic) box has side lengths $L_{x}$ and $L_{y}$ in the $x$ and $y$ directions. Each region has a uniform dielectric constant, which is $\varepsilon_{\rm w}$ and $\varepsilon_{c}$ for the electrolyte and confining material, respectively. To investigate the ions’ transport, we apply an external field E in the direction parallel to the confining walls. The setup is shown in Fig. 1(b).

The ion-ion electrostatic interaction that accounts for the polarization of the confining conducting surfaces is derived from the method developed in Ref. [24]. It circumvents explicit calculation of induced surface charge, applies to strong confinement ($1<L/d<2$), and consists of a single term for $1<L/d\lesssim 1.5$. The method takes advantage of the short-ranged ion-ion electrostatic interaction in strong confinement by conductors, which allows us to use the minimum image convention rather than Ewald summation methods to compute the electrostatic interactions. This significantly accelerates molecular simulations, so we refer to the method as the fast strong conducting confinement (FSCC) method. To derive it, we consider a single confined charge $q_{i}$ at ${\bf r}_{i}=(0,z_{i})$, on the $z$ axis of the cylindrical coordinate system, see Fig. 1(b).

The polarizable infinite planar surfaces are placed at $z=0$ and $z=L$. The Poisson equation was solved for this setup in the context of confined ionic liquids [24]. The electrostatic potential at an arbitrary position ${\bf r}$, generated by the ion $q_{i}$ is given by

where $K_{0}$ is the modified Bessel function of order $0$ and $\Delta r=\sqrt{(x-x_{i})^{2}+(y-y_{i})^{2}}$. For strong confinement, meaning $L$ of order of $d$, just the first term ($n=1$) of the summation in Eq. (1) is necessary because the minimal separation between two ions in the radial direction is $\Delta r_{m}=\sqrt{L(2d-L)}$, see Fig. 1(b). This condition implies that the argument of the modified Bessel function is at least $n\pi\sqrt{(2d/L-1)}$, leading to a fast convergence of sum. Considering $L=1.475d$, this gives $K_{0}(n\pi\sqrt{0.356})\approx 0.13,0.01,0.001$, for $n=1,2,3$, respectively. For the systems studied in this work, we use just the first term, $n=1$, leading to simple expressions, derived below. For larger $L$, $1.5<L/d<2$, one needs to consider more terms in the summation, see Fig. 1(c).

The force between two ions $q_{i}$ and $q_{j}$ at positions $\vec{r}_{i}$ and $\vec{r}_{j}$, converting to Cartesian coordinates, is:

where $F_{x}^{(i,j)}$, $F_{y}^{(i,j)}$ and $F_{z}^{(i,j)}$ are the $x$, $y$, and $z$ components of electrostatic force and $K_{1}$ is the modified Bessel function of order $1$.

The self-electrostatic interaction describes the interaction between an ion and the surface charge it induces. It would be computed by Eq. (1) as $r\rightarrow 0$, but $K_{0}$ diverges in this limit. In this case, we consider the integral based method [24], which gives

The self-force acting on charge $q_{i}$ in the $z$ direction is $F^{self}_{z}=-\frac{q_{i}}{2}\frac{\partial}{\partial z_{i}}\phi_{self}(z_{i})$, which can be written

where $\psi^{(1)}$ is the polygamma function of first order.

To study ionic transport in confinement by conductors, we incorporate the FSCC method into molecular dynamics simulations. Additionally, we consider the case of confining dielectric walls using the method of periodic Green function (PGF) [25]. While both Green function methods (FSCC and PGF) agree for conducting surfaces, the FSCC method is around two orders of magnitude faster.

We investigate a $1:1$ electrolyte under two confinement widths, $L=1.1d$, confining the ions almost exactly to a plane, and a larger value $L=1.475d$. The dielectric constant within the slit is $\varepsilon_{\rm w}=5\varepsilon_{0}$ to represent an organic solvent or strongly confined water [26]. We set $N_{+}=N_{-}=10$, and each ion has a charge $q_{+}=e$ or $q_{-}=-e$ located at its center, where $e$ is the positive elementary charge. The lateral dimensions of the simulation box are $L_{x}=11.9d$ and $Ly=12d$, which were chosen based on surface discretization that will be relevant in Section IV. The ions’ excluded volume is represented with a truncated Lennard-Jones potential (also known as the Weeks-Chandler-Anderson potential) using an energy scale $\epsilon_{\rm LJ}$, with $\sigma_{\rm LJ-ion}=d$. The walls confine the ions via a truncated Lennard-Jones potential, with $\sigma_{\rm LJ}=0.8353d$ and $\epsilon_{\rm LJ}$. The wall-ion interaction is calculated using the Lorentz-Berthelot mixing rule. The centers of the confining walls are located at $z=-0.5\sigma_{\rm LJ}$ and $z=L+0.5\sigma_{\rm LJ}$. The values for the energy, mass, and distance scales are $\epsilon_{\rm LJ}=k_{B}T$ where $T=298$ K, the ion mass set to the mass of sodium $m=22.98$ amu, and $d=0.425$ nm, a typical size of a hydrated ion. We assume symmetry of ion mass and diameter to focus on the effect of polarization charge on ion transport. Effects due to asymmetry of other ion properties will be explored in future work. We study the ionic currents by applying an electric field $\mathbf{E}=E\hat{x}$ tangentially to the surfaces, see Fig. 1(b), and compare the results for conducting ($\varepsilon_{c}\rightarrow\infty$), dielectric ($\varepsilon_{c}=2\varepsilon_{0}$), and non-polarizable surfaces.

We perform molecular dynamics simulations using the Langevin method with a damping parameter of 100 fs and periodic boundary conditions in $x$ and $y$ directions. We consider $100000$ MD steps for equilibration and $100$ MD steps of space between $10000$ uncorrelated samples created for further analysis.

We analyze the ionic density profiles, the current $I$, the radial distributions, and the conductivity of the confined ions. Fig. 2(a) shows that the ions are more adsorbed to conducting surfaces and more repelled from dielectric surfaces than they are from non-polarizable surfaces. This behavior is expected since the ionic interactions with dielectric or conducting surfaces can be understood in terms of equally or oppositely charged images, respectively [27]. Fig. 2(b) shows the pair correlation functions between oppositely charged ions, $g(r)$, which is related to the potential of mean force between ions $w(r)$ as $g(r)=e^{-w(r)/k_{B}T}$. Therefore, Fig. 2(b) shows that the effective interaction between oppositely charged ions is modified by the properties of the confining walls. The ion pairing is seen as a peak in the anion-cation radial distribution function at $r/d\approx 1$. The peak decreases by a factor of $\approx 2.5$ for conducting surfaces compared to dielectric surfaces (see Fig. 2(b)), signifying that the formation of pairs is less favorable between conducting surfaces than dielectrics. The decrease in pair formation in conducting confinement occurs due to the difference in ion-ion interactions near dielectric and conducting surfaces. An ion $q_{i}$ interacts with an oppositely charged ion $q_{j}$ and with $q_{j}$’s equally charged image charges near a dielectric material, but mostly with the $q_{j}$’s opposite image charges near a conducting material. An external electric field E (applied in the channel’s surface parallel direction) decreases the effective anion-cation attraction, but the effect due to the confining walls’ polarization persists (see Fig. 2(c)).

Fig. 3(a) and (b) show that the currents obtained with confining conductive surfaces are much higher than those obtained with dielectric surfaces. The applied electric field and the current are related as $I=A\sigma E$, where $\sigma$ is the conductivity and $A$ is the slit’s cross-sectional area. A field-dependent conductivity (non-linear $I$-$E$ relationship) is essential in memristors [7], which are the basis of memory systems. The three $I$-$E$ curves and their corresponding conductivities in Fig. 3 exhibit diverse regimes, depending on separation between surfaces and polarization. In the field range studied here, for a given value of $L$, the ionic current and conductivity are higher between conductors than between dielectric or non-polarizable surfaces. At the separation distance of $L=1.1d$ between conductors, the current is linear in almost the entire range of studied electric fields. Under such strong confinement, ions are constrained to move nearly on a plane, so the Ohmic behavior in our simulations is consistent with the two-dimensional model between conducting surfaces [19]. The ionic current between conductors at the larger separation distance ($L=1.475d$) shows a non-linear trend for fields below $E\sim 6$ $k_{B}Te^{-1}d^{-1}$. For non-polarizable and dielectric surfaces, the $I$-$E$ curves present a quadratic form, which reflects the linear curves for the conductivity. The exception occurs for dielectric confinement at short separations, which data are best fitted with a $E^{2.6}$ function, which gives a $E^{1.6}$ dependence for conductivity. Combining results using conductor to non-conductor switchable materials suggests a mechanism for inducing hysteresis loops [7] in memristors.

The current nonlinear behavior can be understood in terms of the Onsager’s ion pairing theory [28, 4]. According to this theory, the ions form pairs that last for a time $\tau_{d}$ and stay free for a time $\tau_{a}$. Ion pairs have zero net charge, meaning they do not contribute to the current. Hence, a larger number of ionic pairs between dielectrics than between conductors decreases the overall current. With the increase in electric field the ion pair’s duration time $\tau_{d}$ decreases (reflected as a decrease of the radial distribution peak in Figs. 2(b,c)) and the current increases. Due to the weaker attraction in confinement by conductive surfaces, $\tau_{d}$ is shorter than between dielectric and non-polarizable surfaces. Therefore, in confinement by conductors the ions tend to move dissociated as an electron gas (Ohmic behavior) whereas in confinement by dielectric and non-polarizable surfaces the ion pairs persist in a broader range of electric fields. At high enough $E$-fields all confining materials give an Ohmic response. This variety of behaviors show that the particle-particle interaction strongly influences the ionic current revealing the modulation of surface properties to be a tool for tuning ionic behavior.

The Green function methods apply to planar confinement. For other geometries, such as conical channels suggested for current rectification [29, 30], other methods are necessary. Explicit polarization methods incorporate the surface polarizability into ion-ion interactions by calculating the induced charge in a way that is not restricted to a specific geometry [31, 32]. Conductors are implemented by imposing a uniform surface potential $\psi({\bf s})=const$. In dielectrics with no free surface charges, the electric field boundary conditions at the interface between two media are given by $\varepsilon_{\rm w}\mathbf{E}_{\rm w,n}=\varepsilon_{\rm c}\mathbf{E}_{\rm c,n}$ and $\mathbf{E}_{\rm w,t}=\mathbf{E}_{\rm c,t}$, where the subscripts $w$ and $c$ refer to the ions’ solvent and confining material, respectively, and $n$ and $t$ indicate the perpendicular and tangential electric field components.

The electrostatic boundary conditions define the polarization plane (image plane) which represents the location of the polarization charges. In early density functional theory studies, Lang and Kohn found that the image plane at metallic surfaces is located at the centroid of the induced charge profile [23]. Later work shows that the induced polarization charge peak resides in a range from tenths of an angstrom inside the surface nuclei to a few angstroms outside, depending on electron density of the material, external electric field, and other factors [33]. More recently, it has been suggested that the image plane location affects the adsorption of charged peptide species on gold [34] and shown that the double layer capacitance of silver and graphite surfaces depends on image plane placement [35].

In classical molecular simulations, however, the induced polarization charge is frequently placed at the center of the beads forming the surface [20, 32]. For conducting surfaces, and extending the concept to dielectrics, we study the impact of this choice by placing the polarization charges (1) inside the channel walls ($h=L+\sigma_{\rm LJ}$, Fig. 4(a)) or (2) on the surface ($h=L$, Fig. 4(b)). Because our methods require explicit surface charge to be represented on a discrete mesh, we locate the polarization charges on a hexagonal graphene structure with an average nearest-neighbor distance of $0.332d$ on the $x$-$y$ plane. All the other parameters are the same as described in Section II. We perform molecular dynamics simulations for only the case of $L=1.475d$ using the polarization methods implemented in LAMMPS [36]. In conductors, a constant surface potential is maintained using the Gaussian charge model [37, 38, 31]. For dielectrics, the electrostatic boundary conditions are considered using a boundary element method [32]. Further details of the wall and conductor models are supplied in the Supplemental Material [39].

Fig. 4(c) and (d) show the ion density profiles for conductive and dielectric surfaces, respectively. The figures include the profiles at the polarization charge locations of $h=L$ and $h=L+\sigma_{\rm LJ}$. For comparison, we include the results from the Green function methods (FSCC and PGF) and for non-polarizable systems. Placing the polarization charges on the walls’ surface ($h=L$) leads to conductive surfaces adsorbing the ions and dielectric surfaces repelling them. For both dielectrics and conductors, the effect of surface polarization is significantly dampened when the polarization charges are located inside the surface. Placing the polarization charges inside the surfaces causes the profiles of polarizable systems to converge towards those of non-polarizable systems, while placing the polarization charges on the surface results in profiles that overlap with the Green function methods. Similarly, the ionic current in polarizable models is close to that of non-polarizable system (see Fig. 5(a) and (b)), when the polarization plane is inside the surfaces. When the polarization plane is placed on the walls’ surfaces, the current is similar to that from the Green function methods. We attribute the slight discrepancy between the two methods at high fields to the discrete mesh employed in the explicit polarization models (see Fig. S2 in the Supplemental Material [39]). The pair correlation functions in Figs. 5(c) and (d) explain the current behavior. The systems with surface charge placement indicate decreased ion clustering (lower peak) for conductive surfaces and enhanced ion clustering (higher peak) for dielectric materials, aligning with the Green function methods. Our findings demonstrate that both the surface polarization and the polarization plane placement modulate the ion conductance.

We study the adsorption, interactions, and current of ions confined in slit-like channels made of polarizable surfaces. The ionic current is non-linear and shows different behavior in channels made of conductors, dielectrics, and non-polarizable surfaces. Strong confinement and polarization effects lead to non-Ohmic ionic currents whereas the two-dimensional transport between conductors tends to be linear. The ionic current modulation and non-linear trend (essential for designing memristors) are caused by the screened or enhanced ionic clustering that is different for each type of confining material. We demonstrate that location of the polarization plane and the confining distance significantly affect the confined fluid properties. Considering the surface of materials with different electronic properties is crucial for studying hybrid ionic/electron coupling in numerous fields such as materials with neuromorphic applications, energy harvesting, and water desalination. In quantum-mechanical calculations, however, considering the coupling between ions transport and the electrons of a surface is challenging due to the differences in time and length scales. Here, we employ the polarization and the polarization plane location to overcome that difficulty in a consistent way with density functional theory.