Electric Double Layer and Orientational Ordering of Water Dipoles in Narrow Channels within a Modified Langevin Poisson-Boltzmann Model

Abstract

The electric double layer (EDL) is an important phenomenon that arises in systems where a charged surface comes into contact with an electrolyte solution. In this work we describe the generalization of classic Poisson-Boltzmann (PB) theory for point-like ions by taking into account orientational ordering of water molecules. The modified Langevin Poisson-Boltzmann (LPB) model of EDL is derived by minimizing the corresponding Helmholtz free energy functional, which includes also orientational entropy contribution of water dipoles. The formation of EDL is important in many artificial and biological systems bound by a cylindrical geometry. We therefore numerically solve the modified LPB equation in cylindrical coordinates, determining the spatial dependencies of electric potential, relative permittivity and average orientations of water dipoles within charged tubes of different radii. Results show that for tubes of a large radius, macroscopic (net) volume charge density of coions and counterions is zero at the geometrical axis. This is attributed to effective electrolyte charge screening in the vicinity of the inner charged surface of the tube. For tubes of small radii, the screening region extends into the whole inner space of the tube, leading to non-zero net volume charge density and non-zero orientational ordering of water dipoles near the axis.

1. Introduction

The electric double layer (EDL) is a central phenomenon found at the boundary between a charged surface and an electrolyte solution [1,2,3,4,5,6,7,8]. The counterions are accumulated close to the charged surface and the coions are depleted from this region, resulting in a non-homogeneous distribution of ions. The physical properties of the EDL are crucial in understanding colloidal systems, transport of charged molecules across biological membrane channels or binding of charged proteins to biological surfaces.

Recently, much attention is being devoted to inorganic and organic hollow cylindrical structures in the nanometer range due to their potential benefit in technology, biology and medicine [9]. Potential applications range from microelectronics to microfluidics [10]. Ion channels or pores in biological membranes and blood capillaries are also examples for cylindrical nanotubes.

In some biological systems, the walls of organic nanotubes are charged and in contact with electrolyte solution, where the primary agents of interaction are electrostatic forces, both between charged particles and polar water molecules. Due to the surface charge of the walls, counterions and coions of the electrolyte are, respectively, accumulated and depleted near the walls. At the internal surfaces concave electrical double layers of cylindrical geometry are formed [11].

Furthermore, when bound to a cylindrical geometry, the effect of curvature on EDL properties is significant on small enough scales. Such biological cylindrical channels, where EDL interactions are important, encompass axons or tunneling nanotubes [12]. When artificially made channels, for example, those found in nanoporous materials, are used in the manufacture of electrochemical nanocapacitors, their power and energy densities are dependent on EDL characteristics such as capacitance [13,14,15].

EDL was first modeled by Helmholtz who assumed that the charged surface attracts the surrounding point-like counterions and a single layer is formed to screen the charge [16,17]. Later, these ions have been described by a Boltzmann distribution, forming a diffuse layer extending into the bulk [18,19]. The finite size has been incorporated by Stern with the so-called distance of closest approach [20] and later developed further by numerous authors [3,21,22,23,24,25,26]. In recent decades, EDL has been the subject of numerous analytical and numerical studies from Monte-Carlo methods, DFT theories and lattice models [3,7,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44]. Additionally, interest in nanostructured materials [45,46,47,48] requires that theoretical models of EDL are revisited [49,50,51], also by taking into account the possible quantum effects [52,53].

It has been shown that close to the charged surface, orientational ordering and depletion of water molecules may result in a strong decrease in the local permittivity of the electrolyte solution [54,55,56,57,58,59,60,61]. Considering the orientational ordering of water and finite size of molecules, Outhwaite and collaborators developed a modified Poisson-Boltzmann’s (PB) theory of EDL composed of a mixture of hard spheres with point-like dipoles and finite-sized ions [54,62]. Later, Szalai et al. [63] published a mean spherical approximation-based theory [64] that can reproduce simulation results for the electric field dependence of the dielectric permittivity of a dipolar fluid in a saturation regime. The problem was also considered within a discrete lattice statistics model taking into account the asymmetric size of ions and orientational ordering of water dipoles [44]. Recently, ion-ion and ion-water correlations were also considered in a mean-field approach [65,66].

In the present paper, we first discuss the relative permittivity of water molecules within a cavity field model. We then go on to the derivation of a modified Langevin Poisson-Boltzmann (LPB) equation for point-like ions and water dipoles for planar geometry and then generalize the equations for arbitrary geometry. In derivation of modified LPB equation we construct a Helmholtz free energy functional and minimize it to derive the analytical expressions for ion distributions and spatial dependence of statistical averages orientations of water dipoles. The free energy expression also includes contributions from configurational entropy of ions and rotational entropy of water dipoles. In the second part of the paper the modified LPB equation and the corresponding boundary conditions, generalized for an arbitrary geometry, are utilized to present the numerical solution for a cylindrical geometry with special emphasis given to very narrow cylindrical channels (Figure 1).

2. Relative Permittivity of Water

The dipole moment of an isolated water molecule is around 1.85 D (Debye is $3.336\times {10}^{-30}$ Cm). In a solution, the dipole moment of a single water molecule differs from an isolated one since each molecule is also polarized by the electric field of the neighboring water molecules, creating an effective value of the dipole moment around 2.4 D–2.6 D [67,68]. The body of literature dealing with the dielectric permittivity of water is voluminous and comprehensive, from analytic models detailing the state of bound water molecules and water in charged crevices [69,70] to molecular dynamics simulations with nonlinear response to external electric fields [71,72].

The effect of a polarizing environment can be reproduced in the most simple way by introduction of the cavity field [61,73,74,75]. Cavity field is derived by solving the Poisson’s equation of a model water molecule placed in an outside homogeneous electric field (for a detailed derivation, see Reference [76]). The present section deals with polarization of water dipoles that follows directly from the cavity field.

The water molecules are described within the modified Kirkwood approach [75] as point-like dipoles $\mathbf{p}$ with magnitude $\left|\mathbf{p}\right|=p$ at the centres of finite sized spheres, embedded in a medium with electric permittivity representing the ion-water solution ${\epsilon }_r$ (Figure 2) [7,61]. Within this medium, a spatially homogeneous electric field, $\mathbf{E}$, is present. Due to the built up charge at the interface between the inside and outside of the sphere, the dipole experiences the so called cavity field ${\mathbf{E}}_c$. The relative permittivity of water is given by ${\epsilon }_r=1+P_{\mathrm{tot}}/\left({\epsilon }_{0}E\right)$, where $P_{\mathrm{tot}}$ is the total polarization of water dipoles, E is the magnitude of the spatially homogeneous electric field and ${\epsilon }_0$ is the permittivity of vacuum. The total polarization is the sum of electronic polarization, $P_e$, and orientational polarization due to the permanent water dipoles P, so that $P_{\mathrm{tot}}=P_e+P$. The electronic polarization determines the refractive index of water [51,61] $n^2=1+P_e/\left({\epsilon }_{0}E\right)\approx 1.8$ and ${\epsilon }_r$ can be expressed as

$${\epsilon }_r=n^2+\frac{P}{{\epsilon }_{0}E}.$$

(1)

To find the expression for P we must take into account the constant number density of water $n_w$ and the statistical-average orientation of water molecules in the solution [7]:

$$P=n_{w}p\langle \cos \theta \rangle .$$

(2)

Here, $\theta$ is the angle between $\mathbf{p}$ and the cavity field ${\mathbf{E}}_c$ acting on it (see Figure 3). Statistical averaging is labeled by $\langle \dots ..\rangle$. To estimate $\langle \cos \theta \rangle$, we must first find the expression for ${\mathbf{E}}_c$. This involves solving the Poisson equation for a sphere with electric permittivity $n^2$ embedded in a medium with a relative permittivity ${\epsilon }_r$ described in detail in Reference [76]. Neglecting the short range interactions between dipoles, the local electric field strength at the centre of the sphere at the location of the permanent point-like dipole (Figure 2) can be expressed as [7,76]

$${\mathbf{E}}_c=\frac{3{\epsilon }_r}{n^2+2{\epsilon }_r}\mathbf{E}.$$

(3)

When the surrounding medium has a relative permittivity much larger than the refractive index of water ${\epsilon }_r\gg n^2$, it follows that

$$E_c\approx \frac{3}{2}E\to {\mathbf{E}}_c\approx \frac{3}{2}\mathbf{E}.$$

(4)

So far we have neglected the reaction field, which is the field of the point-like dipole at the center of the cavity itself. This reaction field is directly proportional to the strength of dipole ${\mathbf{E}}_{\mathrm{react}}\propto \mathbf{p}$. In vacuum, in the case of a single isolated water molecule, the external dipole moment is also the experimentally measured dipole moment of a single water molecule ${\mathbf{p}}_e$ given by [7,76]:

$${\mathbf{p}}_e=\frac{3}{n^2+2}\mathbf{p}.$$

(5)

Here, $\mathbf{p}$ is the permanent point-like internal water dipole at the center of the sphere. The energy of an internal point-like dipole $\mathbf{p}$ in a local field ${\mathbf{E}}_c$ is [61]

$$W_e=-\mathbf{p}·{\mathbf{E}}_c.$$

(6)

Substituting from Equation (4) and Equation (5), we can express the dipole energy as [61]

$$\begin{array}{ccc}\hfill W_e& =& -\frac{3}{2}\left(\frac{2+n^2}{3}\right)p_{0}Ecos\theta ,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill W_e& =& \gamma p_{0}Ecos\omega .\hfill \end{array}$$

(8)

Here, $p_0$ is the magnitude of ${\mathbf{p}}_e$ and $\omega$ is supplementary to $\theta$, as shown in Figure 3. The constant $\gamma$ equals [7,61] (see Equations (7) and (8))

$$\gamma =\frac{3}{2}\left(\frac{2+n^2}{3}\right).$$

(9)

With this in mind, the ensemble average in Equation (2) can be calculated as:

$$\langle cos\omega \rangle =\frac{\int cos\omega e^{-(\beta \gamma p_{0}Ecos\omega )}d\mathsf{\Omega }}{\int e^{-(\beta \gamma p_{0}Ecos\omega )}d\mathsf{\Omega }}=-\mathcal{L}\left(\beta \gamma p_{0}E\right).$$

(10)

Here, $\beta$ is the Boltzmann’s factor equal to $\beta =1/kT$, where $kT$ is the thermal energy. The element of solid angle is $d\mathsf{\Omega }=2\pi sin\omega d\omega$, meaning that the integral runs from 0 to $\pi$ with assumed azimuthal symmetry. The Langevin function is defined as $\mathcal{L}\left(u\right)=cothu-1/u$. By taking into account Equations (1), (2), (5) and (10), we can express the relative water permittivity as [7]:

$${\epsilon }_r=n^2+\frac{n_w{p}_0}{{\epsilon }_0}\left(\frac{2+n^2}{3}\right)\frac{\mathcal{L}\left(\beta \gamma p_{0}E\right)}{E}.$$

(11)

In the limit of vanishing electric field strength ($E\to 0$), the above expression for the relative permittivity of water yields to the Onsager limit [7]

$${\epsilon }_r=n^2+\frac{n_w{p}_0^2\beta }{2{\epsilon }_0}{\left(\frac{2+n^2}{3}\right)}^2.$$

(12)

For $p_0=3.1$ D and $n_w/N_A=55$ mol/L, [7,44], where $N_A$ is the Avogadro number, Equation (12) yields the value ${\epsilon }_r=78.5$ at room temperature. Returning to Equation (2), we can write the final result for the orientational polarization of water dipoles P, which will be needed for our Helmholtz free energy minimization in the following section:

$$P=-n_w{p}_0\left(\frac{2+n^2}{3}\right)\mathcal{L}\left(\beta \gamma p_{0}E\right).$$

(13)

3. Derivation of the Modified LPB Equation by Minimization of Helmholtz Free Energy

Our model assumes the electrolyte solution is a mixture of point-like monovalent co- and counterions and permanent water dipoles, representing the water molecules. The expression for the spatial dependence of the solution permittivity ${\epsilon }_r\left(x\right)$, arising as a direct consequence of the spontaneous ordering of water dipoles, can be obtained by using the minimization of the Helmholtz free energy in a one-dimensional setting with the charged planar surface located at $x=0$. In the minimization procedure, the local electric field at the positions of the hydrated point-like ions in the electrolyte solution is denoted by $E\left(x\right)$, while the local cavity field at the positions of the water internal point-like dipoles is denoted by $E_c\left(x\right)$. We can write the Helmholtz free energy of the system F as (see also Reference [58]):

$$\begin{array}{cc}\hfill F& =\underset{F_1}{\underbrace{\frac{{\epsilon }_0{n}^2}{2}\int E_c^2\left(x\right)dV}}+\underset{F_2}{\underbrace{\frac{{\epsilon }_0{n}^2}{2}\int E^2\left(x\right)dV}}+kT[\underset{F_3}{\underbrace{\int (n_{+}\left(x\right)ln\frac{n_{+}\left(x\right)}{n_0}-(n_{+}\left(x\right)-n_0))dV}}+\hfill \\ \hfill & +\underset{F_4}{\underbrace{\int \left(n_{-}\left(x\right)ln\frac{n_{-}\left(x\right)}{n_0}-(n_{-}\left(x\right)-n_0)\right)dV}}+\underset{F_5}{\underbrace{\int ({\lambda }_{+}{n}_{+}\left(x\right)+{\lambda }_{-}{n}_{-}\left(x\right))dV}}+\hfill \\ \hfill & +\underset{F_6}{\underbrace{\int n_w{\langle \mathcal{P}(x,\omega )ln\mathcal{P}(x,\omega )\rangle }_{\omega }dV}}+\underset{F_7}{\underbrace{\int n_w\eta \left(x\right)\left({\langle \mathcal{P}(x,\omega )\rangle }_{\omega }-1\right)dV}}].\hfill \end{array}$$

(14)

The thermal energy is given by $kT$, while n is the refractive index. For greater clarity, we split the particular contributions to the free energy as marked by the underbraces in Equation (14). The mean field created by coions and counterions and the water dipoles polarization contribution are given by terms $F_1$ and $F_2$, respectively. Mixing entropy free energy contributions of point-like counterions and coions are accounted for in terms $F_3$ and $F_4$. The constraint of a constant number of ions in the system is given in $F_5$, where ${\lambda }_{+}$ and ${\lambda }_{-}$ are the global Lagrange’s multipliers for counterions and coions. The free energy that corresponds to orientational entropy of permanent water dipoles is given by $F_6$, while the last term, $F_7$, gives the local constraint for orientation of dipoles. $\mathcal{P}(x,\omega )$ is the probability that a permanent water dipole located at x is oriented at angle $\omega$ with respect to the normal to the charged surface (Figure 3). The brackets ${\langle \dots \rangle }_{\omega }$ denote the average:

$${\langle \mathcal{F}(x,\omega )\rangle }_{\omega }=\frac{1}{4\pi }{\int }_0^{\pi }\mathcal{F}(x,\omega )2\pi sin\omega d\omega .$$

(15)

Here, $\omega$ is the angle between the internal dipole moment vector, $\mathbf{p}$, and ${\mathbf{n}}_{\varphi }=\nabla {\varphi }_c/|\nabla {\varphi }_c|$ (see Figure 3). We perform variation on the Helmholtz free energy, F, in Equation (14), so that $\delta F=0$. Let us deal with the variational approach of every contribution in Equation (14) particularly, beginning with $F_1$ and $F_2$. For clarity of notation, direct spatial dependence will sometimes be omitted, so that for example, $n_{+}\left(x\right)\equiv n_{+}$.

3.1. Variation Procedure

3.1.1. Electric Fields ($F_1$ and $F_2$)

Since there are no time dependent magnetic fields, we can express the electric fields as potentials $E\left(x\right)=-{\varphi }^{\prime }\left(x\right)$, $E_c\left(x\right)=-{\varphi }_c^{\prime }\left(x\right)$, where the prime labels the spatial derivative, and perform a variation on the electrostatic term pertaining to water dipoles.

$$\delta \left(\frac{{\epsilon }_0{n}^2}{2}\int {\varphi }_c^{{}^{\prime }2}dV\right)=\frac{{\epsilon }_0{n}^2}{2}\int 2{\varphi }_c^{\prime }\delta \left({\varphi }_c^{\prime }\right)dV.$$

(16)

We can rearrange this term, if we consider the rules of differentiating a function product

$$\begin{array}{cc}\hfill {\left({\varphi }_c\delta {\varphi }_c^{\prime }\right)}^{\prime }& ={\varphi }_c^{\prime }\delta {\varphi }_c^{\prime }+{\varphi }_c\delta {\varphi }_c^{″},\hfill \\ \hfill {\varphi }_c^{\prime }\delta {\varphi }_c^{\prime }& ={\left({\varphi }_c\delta {\varphi }_c^{\prime }\right)}^{\prime }-{\varphi }_c\delta {\varphi }_c^{″}.\hfill \end{array}$$

(17)

The integral in Equation (16) can be rewritten as,

$$\begin{array}{cc}\hfill {\epsilon }_0{n}^2& \int {\varphi }_c^{\prime }\delta \left({\varphi }_c^{\prime }\right)dV={\epsilon }_0{n}^2(\underset{=0}{\underbrace{{\varphi }_c\delta {\varphi }_c^{\prime }{|}_0^{\infty }}}-\hfill \\ \hfill & -\int {\varphi }_c\delta \left({\varphi }_c^{″}\right)dV),\hfill \end{array}$$

(18)

where the first term on the right-hand side equals 0 at infinity, since we impose the electric potential there to be constant and equal to 0. Taking into account the Poisson’s equation for the water dipoles, namely ${\varphi }_c^{″}\left(x\right)=\nabla ·\mathbf{P}/{\epsilon }_0{n}^2$, where $\mathbf{P}$ represents the net polarization of the permanent water dipoles, we get

$$-{\epsilon }_0{n}^2\int {\varphi }_c\delta \left({\varphi }_c^{″}\right)dV=\int {\varphi }_c\delta {\rho }_{c}dV.$$

(19)

Here, ${\rho }_c$ corresponds to the bound charge density due to the dipoles’ polarizations, which is related to net polarization ${\rho }_c=-\nabla ·\mathbf{P}$. We observe that $\delta (\nabla ·\mathbf{P})=\nabla ·\delta \mathbf{P}$. The integral in (Equation (19)) can now be rewritten:

$$\int {\varphi }_c\delta {\rho }_{c}dV=-\int {\varphi }_c(\nabla ·\delta \mathbf{P})dV.$$

(20)

The product rule for divergence can be used $\nabla ·\left({\varphi }_c\delta \mathbf{P}\right)=(\nabla {\varphi }_c)·\delta \mathbf{P}+{\varphi }_c(\nabla ·\delta \mathbf{P})$, so that the integral of Equation (20) can now be written differently again:

$$\int {\varphi }_c(\nabla ·\delta \mathbf{P})dV=\underset{=0}{\underbrace{\int \nabla ·\left({\varphi }_c\delta \mathbf{P}\right)dV}}-\int (\nabla {\varphi }_c)·\delta \mathbf{P}dV.$$

(21)

Here, the first integral on the right hand side vanishes by virtue of the divergence theorem; since the potential far away from the plates is constant and equal to zero. We therefore arrive at the final result

$$\delta F_1=\int (\nabla {\varphi }_c)·\delta \mathbf{P}dV.$$

(22)

The polarization, $\mathbf{P}$, is related to the average orientation of all water dipoles (Equation (2)):

$$\mathbf{P}\left(x\right)=n_w{\langle \mathcal{P}(x,\omega )\rangle }_{\omega }p{\mathbf{n}}_{\varphi }.$$

(23)

Here, $n_w$ is the number density of water molecules in the solution, $p=\left|\mathbf{p}\right|$ is the internal point-like dipole magnitude, ${\mathbf{n}}_{\varphi }=\nabla {\varphi }_c/|\nabla {\varphi }_c|$ is the unit vector away from the charged plate and ${\langle \mathcal{P}(x,\omega )\rangle }_{\omega }$ is defined by Equation (15) (see Figure 3). Since our case deals with a negatively charged surface ($\sigma <0$), $\mathbf{P}$ points in the direction opposite to the direction of the x-axis and is thus negative (for details see Reference [76]). Since the variation of $\mathbf{P}$ can be written $\delta \mathbf{P}\left(x\right)={\langle n_w\mathbf{p}\delta \mathcal{P}(x,\omega )\rangle }_{\omega }$, we arrive at the variation of $F_1$:

$$\delta F_1=n_w\int {\langle \delta \mathcal{P}(x,\omega )(\nabla {\varphi }_c)·\mathbf{p}\rangle }_{\omega }dV.$$

(24)

Similarly, for $F_2$ by taking into account Equation (17), we get

$$\delta \left(\frac{{\epsilon }_0{n}^2}{2}\int {\varphi }^{{}^{\prime }2}dV\right)=\int \varphi \delta {\rho }_{\mathrm{free}}dV.$$

(25)

The Poisson equation is different for free charges (ions): ${\varphi }^{″}\left(x\right)=-{\rho }_{\mathrm{free}}/{\epsilon }_0{n}^2=e_0(n_{+}\left(x\right)-n_{-}\left(x\right))$. The variation by ${\varphi }^{″}\left(x\right)$ in Equation (25) can be written with macroscopic net volume charge density ${\rho }_{\mathrm{free}}\left(x\right)$, which in turn is the sum of the contributions of the local net ion charges. Performing the variation on ion charge distribution ${\rho }_{\mathrm{free}}\left(x\right)$ gives

$$\delta {\rho }_{\mathrm{free}}=e_0(\delta n_{+}-\delta n_{-}),$$

(26)

finishing the variation of the term $F_2$:

$$\delta F_2=\int e_0\varphi (\delta n_{+}-\delta n_{-})dV.$$

(27)

3.1.2. Ion Mixing ($F_3$, $F_4$ and $F_5$)

It makes sense to perform the variation of the ion mixing terms ($F_3$ and $F_4$), together with their Lagrange multipliers ($F_5$), since the variation $\delta n_{+}$ and $\delta n_{-}$ will be a common term for positive and negative ions, respectively. It is easily shown from Equation (14) that

$$\delta F_3+\delta F_4+\delta F_5=kT\int \delta n_{+}({\lambda }_{+}+ln\frac{n_{+}}{n_0})dV+kT\int \delta n_{-}({\lambda }_{-}+ln\frac{n_{-}}{n_0})dV.$$

(28)

3.1.3. Dipole Mixing ($F_6$ and $F_7$)

Variation of the terms $F_6$ and $F_7$ is straightforward. Since the bulk water number density, $n_w$, is taken to be constant, the variation of $F_6$ is

$$\delta F_6=kT{n}_w\int \left({\langle \delta \mathcal{P}(x,\omega )ln\mathcal{P}(x,\omega )+\delta \mathcal{P}(x,\omega )\rangle }_{\omega }\right)dV.$$

(29)

The last variation of $F_7$ is performed over the probability, $\mathcal{P}(x,\omega )$, and the Lagrange multiplier, $\eta \left(x\right)$. Expanding and applying the product rule, we find that

$$\delta F_7=kT{n}_w\int (\delta \eta \left(x\right){\langle \mathcal{P}(x,\omega )\rangle }_{\omega }+\eta \left(x\right){\langle \delta \mathcal{P}(x,\omega )\rangle }_{\omega }-\delta \eta \left(x\right))dV.$$

(30)

3.2. Euler-Lagrange Equations

Combining the variations of all the integrals given in Equation (14), their sum $\delta F$ gives us the variation of Helmholtz free energy. Factoring all the variation terms with respect to $n_{+}\left(x\right),n_{-}\left(x\right),\mathcal{P}(x,\omega )$ and $\eta \left(x\right)$ gives

$$\begin{array}{cc}\hfill & \delta F=\int dV\delta n_{+}\left(x\right)\left[kT\left(ln\frac{n_{+}\left(x\right)}{n_0}+{\lambda }_{+}\right)+\varphi e_0\right]+\int dV\delta n_{-}\left(x\right)\left[kT\left(ln\frac{n_{-}\left(x\right)}{n_0}+{\lambda }_{-}\right)-\varphi e_0\right]+\hfill \\ \hfill & +\int dV{n}_w{\langle \delta \mathcal{P}(x,\omega )\left(\nabla {\varphi }_c·\mathbf{p}+\frac{ln\mathcal{P}(x,\omega )+\eta \left(x\right)+1}{\beta }\right)\rangle }_{\omega }+kT\int dV{n}_w\delta \eta \left(x\right)\left({\langle \mathcal{P}(x,\omega )\rangle }_{\omega }-1\right).\hfill \end{array}$$

(31)

The volume differentials in a planar geometry are $dV=Sdx$. Since the minimization condition demands $\delta F=0$, the expressions multiplied by $\delta n_{+}\left(x\right)$, $\delta n_{-}\left(x\right)$, $\delta \mathcal{P}(x,\omega )$ and $\delta \eta \left(x\right)$ in the last equation must equal zero, resulting in a system

$$\begin{array}{c}\hfill kT\left(ln\frac{n_{+}\left(x\right)}{n_0}+{\lambda }_{+}\right)+\varphi e_0=0,\end{array}$$

(32)

$$\begin{array}{c}\hfill kT\left(ln\frac{n_{-}\left(x\right)}{n_0}+{\lambda }_{-}\right)-\varphi e_0=0,\end{array}$$

(33)

$$\begin{array}{c}\hfill E_{c}pcos\omega +\frac{ln\mathcal{P}(x,\omega )+\eta \left(x\right)+1}{\beta }=0,\end{array}$$

(34)

$$\begin{array}{c}\hfill {\langle \mathcal{P}(x,\omega )\rangle }_{\omega }-1=0.\end{array}$$

(35)

Here, we write $\beta =1/kT$ and expand the dot product $\nabla {\varphi }_c·\mathbf{p}=E_{c}pcos\omega$ (see Figure 3). Solving Equations (32) and (33), we obtain the ion spatial distributions

$$n_{+}\left(x\right)=n_{0}exp\left(-\beta e_0\varphi -{\lambda }_{+}\right),$$

(36)

$$n_{-}\left(x\right)=n_{0}exp\left(\beta e_0\varphi -{\lambda }_{-}\right).$$

(37)

The boundary conditions state that $\varphi (x\to \infty )=0$ and $n_{+,-}(x\to \infty )=n_0$, which renders ${\lambda }_{+}={\lambda }_{-}=0$. We may now turn our attention to the variation of permanent water dipoles orientation. Solving for $\mathcal{P}(x,\omega )$, Equation (34) gives

$$\mathcal{P}(x,\omega )=\mathsf{\Lambda }\left(x\right)exp\left(-\beta E_{c}pcos\omega \right),$$

(38)

where $\mathsf{\Lambda }\left(x\right)=exp(-\eta (x)-1)$. Substituting the cavity field $E_c$ by E (Equation (4)) and dipole moment magnitude p by $p_0$ (Equation (5)) gives

$$\mathcal{P}(x,\omega )=\mathsf{\Lambda }\left(x\right)exp\left(-\beta \frac{3E}{2}\left(\frac{2+n^2}{3}\right)p_{0}cos\omega \right),$$

(39)

where $p_0$ is the magnitude of ${\mathbf{p}}_e$. The final result is expressed using the constant $\gamma$ defined in Equation (9):

$$\mathcal{P}(x,\omega )=\mathsf{\Lambda }\left(x\right)exp\left(-\beta \gamma E{p}_{0}cos\omega \right).$$

(40)

We can now evaluate the average internal dipole moment by integrating over mean orientations (considering Equation (23)),

$$\begin{array}{cc}\hfill & p\langle cos\omega \rangle =p_0\left(\frac{2+n^2}{3}\right)\langle cos\omega \rangle \hfill \\ \hfill & =\frac{{\int }_0^{\pi }\left(\frac{2+n^2}{3}\right)p_{0}cos\omega exp(-\beta \gamma E{p}_{0}cos\omega )d\mathsf{\Omega }}{{\int }_0^{\pi }exp(-\beta \gamma E{p}_{0}cos\omega )d\mathsf{\Omega }}\hfill & \\ \hfill & =-p_0\left(\frac{2+n^2}{3}\right)\mathcal{L}\left(\beta \gamma E{p}_0\right).\hfill \end{array}$$

(41)

The orientational polarization of water is thus (see Equations (2) and (5)):

$$\begin{array}{ccc}\hfill P\left(x\right)& =& n_{w}p\langle cos\omega \rangle \hfill \\ & =& -n_w{p}_0\left(\frac{2+n^2}{3}\right)\mathcal{L}\left(\beta \gamma E\left(x\right)p_0\right).\hfill \end{array}$$

(42)

If we insert the above result and the ion distribution functions (Equations (36) and (37)) into the average microscopic charge density equation $\rho \left(x\right)={\rho }_{\mathrm{free}}\left(x\right)-dP/dx$ [61,77], where ${\rho }_{\mathrm{free}}$ is the contribution of the net ion charges Equations (26), (36) and (37) and $P\left(x\right)$ is the polarization due to partially oriented water dipoles, we get the expression for $\rho \left(x\right)$ in a one-dimensional case:

$$\rho \left(x\right)=-2e_0{n}_{0}sinh\left(\beta e_0\varphi \left(x\right)\right)+n_w{p}_0\left(\frac{2+n^2}{3}\right)\frac{d}{dx}\left(\mathcal{L}\left(\beta \gamma E\left(x\right)p_0\right)\right).$$

(43)

Inserting the above expression for average microscopic volume charge density $\rho \left(x\right)$ into the Poisson’s equation,

$${\varphi }^{″}\left(x\right)=-\frac{\rho \left(x\right)}{n^2{\epsilon }_0},$$

(44)

we get the modified LPB differential equation for the electric potential $\varphi \left(x\right)$:

$${\varphi }^{″}\left(x\right)=\frac{1}{n^2{\epsilon }_0}\left[2e_0{n}_{0}sinh(\beta e_0\varphi \left(x\right))-n_w{p}_0\left(\frac{2+n^2}{3}\right)\frac{d}{dx}\left(\mathcal{L}\left(\beta \gamma E\left(x\right)p_0\right)\right)\right],$$

(45)

where ${\varphi }^{″}\left(x\right)$ is the second derivative of the electric potential $\varphi \left(x\right)$ with respect to x and $E\left(x\right)=-{\varphi }^{\prime }\left(x\right)$. Equation (45) can be factorized via a product rule if we take into account that the Langevin function is odd and its derivative is ${\mathcal{L}}^{\prime }\left(u\right)=1/u^2-1/{sinh}^{2}u$ in the following form [7]:

$$\frac{d}{dx}\left[{\epsilon }_0{\epsilon }_r\left(x\right){\varphi }^{\prime }\left(x\right)\right]=2e_0{n}_{0}sinh\left(\beta e_0\varphi \left(x\right)\right),$$

(46)

$${\epsilon }_r\left(x\right)=n^2+n_w\frac{p_0}{{\epsilon }_0}\left(\frac{2+n^2}{3}\right)\frac{\mathcal{L}\left(\beta \gamma E\left(x\right)p_0\right)}{E\left(x\right)},$$

(47)

where ${\epsilon }_r\left(x\right)$ is the relative permittivity (Equation (11)). This modified Langevin Poisson-Boltzmann (LPB) differential equation (Equation (46)) is subject to two boundary conditions. The first boundary condition arises from the electro-neutrality of the system, which assumes that the total net charge of the system is zero, hence

$$\int {\rho }_{\mathrm{free}}\left(x\right)dV-\sigma S=0,$$

(48)

where $\sigma$ is the negative membrane surface charge density, S is the total membrane surface area and ${\rho }_{\mathrm{free}}\left(x\right)=-2e_0{n}_{0}sinh\left(\beta e_0\varphi \left(x\right)\right)$ is the macroscopic (net) volume charge density of coions and counterions. Since the macroscopic volume charge density is only dependent on x (Equation (43)) and the differential $dV=Sdx$, Equation (48) may be rewritten

$${\int }_0^{\infty }2e_0{n}_{0}sinh\left(\beta e_0\varphi \left(x\right)\right)dx=-\sigma .$$

(49)

If we integrate Equation (45) once over the whole system, we get

$${\varphi }^{\prime }(x=0)=-\frac{1}{n^2{\epsilon }_0}\left[\sigma +n_w{p}_0\left(\frac{2+n^2}{3}\right)·\mathcal{L}\left(\beta \gamma E{p}_0{|}_{x=0}\right)\right].$$

(50)

The second boundary condition states that the electric potential far away from the charged surface (in the bulk) is constant ${\varphi }^{\prime }(x\to \infty )=0$, rendering $\mathcal{L}\left(\beta \gamma E{p}_0{|}_{x\to \infty }\right)=0$. The modified LPB equation (Equation (46)) was derived in one dimension, but can be rewritten in a more general form to apply to an arbitrary three-dimensional geometry. In three dimensions, the steps are analogous and discussed in detail in a previous work [58], where a three-dimensional Lagrangian was derived for a model of finite-sized ions. With this in mind, the modified LPB equation (Equation (46)) can be rewritten:

$$\nabla ·\left[{\epsilon }_0{n}^2\nabla \varphi \left(\mathbf{r}\right)\right]+n_w{p}_0\left(\frac{2+n^2}{3}\right)\nabla ·\left({\mathbf{n}}_{\varphi }\mathcal{L}\left(\beta \gamma E{p}_0\right)\right)=2e_0{n}_{0}sinh(\beta e_0\varphi \left(\mathbf{r}\right)),$$

(51)

where ${\mathbf{n}}_{\varphi }=\nabla \varphi /|\nabla \varphi |=\nabla \varphi /E$. We may factor the last equation, so that

$$\nabla ·\left[{\epsilon }_0\left(n^2+\frac{n_w{p}_0}{{\epsilon }_0}\left(\frac{2+n^2}{3}\right)\frac{\mathcal{L}\left(\beta \gamma E{p}_0\right)}{E}\right)\nabla \varphi \left(\mathbf{r}\right)\right]=2e_0{n}_{0}sinh(\beta e_0\varphi \left(\mathbf{r}\right)).$$

(52)

This modified LPB equation can be written even more compactly, considering the definition of spatially dependent permittivity ${\epsilon }_r\left(\mathbf{r}\right)$ given by Equation (47) (for details, see Reference [58]):

$$\nabla ·[{\epsilon }_0{\epsilon }_r\left(\mathbf{r}\right)\nabla \varphi \left(\mathbf{r}\right)]=2e_0{n}_{0}sinh\left(\beta e_0\varphi \left(\mathbf{r}\right)\right),$$

(53)

$${\epsilon }_r\left(\mathbf{r}\right)=n^2+n_w\frac{p_0}{{\epsilon }_0}\left(\frac{2+n^2}{3}\right)\frac{\mathcal{L}\left(\beta \gamma E\left(\mathbf{r}\right)p_0\right)}{E\left(\mathbf{r}\right)}.$$

(54)

Here, ${\rho }_{\mathrm{free}}\left(\mathbf{r}\right)$ is the macroscopic (net) volume charge density of coions and counterions. A corresponding three-dimensional variant of the boundary condition (Equation (50)) is

$$\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})=-\frac{1}{n^2{\epsilon }_0}\left[\sigma {\mathbf{n}}_{\varphi }+{\mathbf{n}}_{\varphi }{n}_w{p}_0\left(\frac{2+n^2}{3}\right)·\mathcal{L}\left(\beta \gamma E\left(\mathbf{r}\right)p_0\left(\mathbf{r}\right){|}_{\mathbf{r}=\mathrm{surf}}\right)\right].$$

(55)

Rearranging, it follows that

$$\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})\left[1+\frac{{\mathbf{n}}_{\varphi }}{\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})}\frac{n_w{p}_0}{n^2{\epsilon }_0}\left(\frac{2+n^2}{3}\right)·\mathcal{L}\left(\beta \gamma E\left(\mathbf{r}\right)p_0\left(\mathbf{r}\right){|}_{\mathbf{r}=\mathrm{surf}}\right)\right]=-\frac{\sigma }{n^2{\epsilon }_0}{\mathbf{n}}_{\varphi }.$$

(56)

Evaluating the second expression on the left hand side of the last equation gives

$$\frac{{\mathbf{n}}_{\varphi }}{\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})}=\frac{\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})}{|\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}}\left)\right|}\frac{1}{\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})}=\frac{1}{E(\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})}.$$

(57)

Combining this simplification with Equation (42), Equation (55) becomes

$$\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}}){\epsilon }_r(\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})=-\frac{\sigma {\mathbf{n}}_{\varphi }}{{\epsilon }_0}.$$

(58)

Here we also take into account the expression for ${\epsilon }_r$ (Equation (47)). We see that the term inside the square brackets on the left-hand side of Equation (56) is precisely the definition of the relative permittivity on the surface of charged membrane ${\epsilon }_r(\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})$ (Equation (54)), yielding the general result

$$\nabla \varphi (\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})=-\frac{\sigma {\mathbf{n}}_{\varphi }}{{\epsilon }_0{\epsilon }_r(\mathbf{r}={\mathbf{r}}_{\mathrm{surf}})}.$$

(59)

4. Results

Figure 4 shows the dependency of the calculated macroscopic (net) volume charge density of the electrolyte solution inside the nanotubes (${\rho }_{\mathrm{free}}\left(r\right)$) on the radial distance from the geometrical axis of the tube. It can be seen in the figure that for larger radii of the inner cross-sections of the nanotubes (R), the value ${\rho }_{\mathrm{free}}$ at the geometrical axis of the tube is zero, which means that the number densities of counterions and coions there are equal and the electric potential is constant, that is, zero in our case (see the right panel in Figure 5).

On the contrary, for smaller values of the nanotube radius R, the value of ${\rho }_{\mathrm{free}}$ at geometrical axis of the nanotube is not zero (Figure 4). Accordingly, for small values of the radius of the inner nanotube also the gradient of the electric field (Figure 6) and the electric potential at the nanotube geometrical axis are not zero (left panel in Figure 5). Hence, the bulk condition of the equal number densities of counterions and coions is fulfilled outside the interior of the nanotube.

Figure 7 shows the dependency of the average orientation ${\langle cos\left(\omega \right)\rangle }_{\omega }$ and the relative permittivity ${\epsilon }_r$ on the radial distance from the geometrical axis of tube (r), calculated for four different values of nanotube inner radius R. It can be seen that for small radii, R, the average orientational of water dipoles is relatively strong also in the vicinity of geometrical axis of the tube, while for larger R the average orientation of water dipoles is strong only in the region near the charged inner surface of the tube.

5. Discussion and Conclusions

In this paper, we derived a modified Langevin Poisson-Boltzmann (LPB) model and the modified LPB equation to theoretically describe the electric double layer (EDL) for a monovalent electrolyte solution inside very narrow nanotubes with a negatively charged inner surface. In the modified LPB approach, the electronic polarization of the water is taken into account by assuming a permanent dipole embedded in the center of the sphere with a volume equal to the average volume of a water molecule. The effect of a polarizing environment is reproduced by introduction of the cavity field in the saturation regime [7,61,76]. In past EDL studies, treatments of cavity fields and structural correlations between water dipoles were limited to cases of relatively small electric field strengths, far away from the saturation limit of polarization and orientational ordering of water molecules [73,74,75]. High magnitudes of electric field strength were later added in several works [44,61,63,78]. A commonly oversimplified assumption when theoretically describing the EDL is the assumption of a surface charge density-independent relative permittivity in the inner (Stern) layer. Due to orientational ordering of water dipoles, the relative permittivity of the Stern layer depends on the electric field strength, that is, on the surface charge density ($\sigma$) of the electrode [51,79,80,81,82]. Furthermore, fitting the model curves with a range of free parameters to the experimental points [83] cannot prove that the Stern layer capacitance and permittivity is $\sigma$-independent. The decrease in the relative permittivity close to the charged surface (electrode) is obviously partially the consequence of orientational ordering of water dipoles close to the saturation regime or in the saturation regime as shown theoretically in References [6,27,44,54,58,59,61,62,63,64,80,82].

Within a recently presented phenomenological approach it is claimed that close to the charged surface, almost all water molecules belong to water shells around the ions, while the free water molecules are excluded [83]. The results of simulations clearly refute this fact [84] by showing increased water ordering in the direction towards the charged surface (including the region close to the charged surface) (Figure 7, upper panel) even for high salt concentrations [84], in quantitative agreement with mean-field theoretical predictions [7,82]. For example, for a magnitude of $0.16{\mathrm{As}/\mathrm{m}}^2$ surface charge density, there is practically no difference in the orientational ordering and space distribution of water dipoles close to the charged surface between water with and without NaCl (of concentration 500 mmol/L) [84]. In general, for magnitudes of surface charge density up to around $0.3\mathrm{As}/{\mathrm{m}}^2$, where the mean-field approach can still be justified [7,82], there is only a weak quantitative influence of salt on the profile of orientational ordering of water dipoles in Stern and diffuse layers, but not qualitative [84]. Note that the multi-layering of water predicted in simulations [84] cannot be predicted within our mean-field approach [44,61] as well also not in the oversimplified phenomenological models [83].

Besides the saturation in polarization/water dipole orientation at high magnitudes of the electric field strength, the important thing to consider in the EDL studies is also the saturation in the counterion concentration near the charged surface due to the finite size of ions. These steric effects were first predicted in the Wicke-Eigen’s model (also called the Bikerman’s model) and their modifications [3,5,22,25,27,35,85,86]. For finite sized ions, the dielectric permittivity profile in the vicinity of a charged surface is modulated by the depletion of water dipoles at the charged surface due to accumulated counterions [58,82]. In the modified LPB model [7,59], described in the present paper, these steric effects were not taken into account.

The described decrease in the relative permittivity relative to its bulk value in the present paper is the consequence of strong orientational ordering of the water dipoles in the vicinity of the charged surface (Figure 7). Contrary to our results it is claimed in Reference [87] that the relative permittivity is increased in direction to the charged surface. As pointed out in publications of different authors the predicted increase of relative permittivity near the charged surface in Reference [87] is unphysical [6,59] and defies the common wisdom in electrochemistry [56]. In addition, the experiments report just the converse as predicted in Reference [87], that is, the experiments indicated the decrease of relative permittivity near the charged surface [88,89]. The predicted substantial increase of relative permittivity in the inner part of the double layer near the charged surface in Reference [87] is due to arise in the dipole density near the surface as pointed out in Reference [56]. This unphysical result [6] is the consequence of inconsistency of so-called dipolar PB theory presented in Reference [87] as indicated in Reference [59]. Namely, it was shown [59] that the dipolar PB theory for point-like ions in Reference [87] assumes an orientationally averaged Boltzmann factor for spatial distribution function for water dipoles, which is however not compatible with the assumption of point-like ions. Energy dependent spatial distribution of water dipoles cannot be taken into account simultaneously with the assumption of point-like ions, but only if the finite size of molecules in the electrolyte solution are taken into account [35,61]. This means that the dipolar PB model presented in Reference [87] is not a self-consistent model and consequently predicts unphysical results which are not compatible with experimental results even qualitatively, as noticed in References [6,56,59] and other publications.

The other important difference between our modified LPB model and the theoretical model presented in Reference [87] is that our value for (external) water dipole moment 3.1 D [7,44,51,61] is considerably smaller than the corresponding value 4.86 D used in Reference [87]. The value 3.1 D is closer to the experimental values of the effective dipole moment of water molecules in clusters (2.7 D) and in bulk solution (2.4–2.6 D) (see for example Reference [68]). The value 4.86 D is so large in order to compensate for the cavity field [6,61,74,75,78] that is not taken into account in Reference [87], as noticed also in Reference [6], but is considered in the present modified LPB model. The model value 3.1 D can be additionally decreased by taking into account structural correlations between water dipoles [60,78]. The ion-ion and ion-water correlations were taken into account also in the mean-field models of References [8,65,66].

It has been shown that for finite-sized ions the drop in the number density of water near a charged surface results in an additional decrease of permittivity [7,58]. A further generalization of the modified LPB model with steric effects taken into account within a lattice-statistics model of a modified LPB is found in References [44,51,82]. By taking into account asymmetric finite size of ions the modified LPB equation was generalized to (modified Langevin Eigen-Wicke model) [44,51,82]:

$$\frac{d}{dx}\left[{\epsilon }_0{\epsilon }_r\left(x\right)\frac{d\varphi }{dx}\right]=2e_0{n}_s{n}_0\frac{sinh\left(\beta e_0\varphi \right)}{{\mathcal{D}}_A(\varphi ,E)},$$

(60)

where ${\epsilon }_r\left(x\right)$ is the spatial dependence of relative permittivity:

$${\epsilon }_r\left(x\right)=n^2+n_{0w}{n}_s\frac{p_0}{{\epsilon }_0}\left(\frac{2+n^2}{3}\right)\left(\frac{\mathcal{F}\left(\gamma p_{0}E\beta \right)}{{\mathcal{D}}_A(\varphi ,E)E}\right)$$

(61)

and

$${\mathcal{D}}_A\left(\varphi \right)={\alpha }_{+}{n}_0{e}^{-e_0\varphi \beta }+{\alpha }_{-}{n}_0{e}^{+e_0\varphi \beta }+\frac{n_{0w}}{\gamma p_{0}E\beta }sinh\left(\gamma p_{0}E\beta \right).$$

(62)

Here, the parameters ${\alpha }_{+}$ and ${\alpha }_{-}$ are the number of lattice sites occupied by a single positive and negative hydrated ion, respectively, where a single water molecule is assumed to occupy just one lattice site. The reduced number density of lattice sites $n_s/N_A=55\mathrm{mol}/\mathrm{L}$ is equal to the concentration of pure water [44,51,82]. The symbol $n_{0w}$ stands for the bulk number density of water molecules. The function $\mathcal{F}\left(u\right)$ is defined as $\mathcal{F}\left(u\right)=\mathcal{L}\left(u\right)sinh\left(u\right)/u$, where $\mathcal{L}\left(u\right)$ is the Langevin function.

The results of the present paper are important when considering electric fields within artificial as well as biological channels containing an electrolyte. Much attention has recently been given to understanding tunneling nanotubes (TNTs), small tubular structures that drive cell communication and spreading of pathogens [12]. Not yet fully understood, it is thought that these tubular structures initiate from local membrane bending facilitated by laterally distributed proteins or anisotropic membrane nanodomains. Further research is needed to clarify the role of EDL in the inception of these structures, since cytoplasmatic proteins and other elements are electrically charged. When such motor proteins are complemented by protruding cytoskeletal forces provided by the polymerization of f-actin, TNT formation is crucial in determining cell morphology, sometimes even leading to endovesiculation of the red blood cell membrane [90,91,92]. Recently, within a molecular mean-field approach and taking into account the asymmetric size of ions, polarization of water, and ion-ion and ion-water correlations, the ionic and water flows through biological ion channels was theoretically considered [65,66].

To conclude, in the present paper we started from a mean-field Helmholtz free energy functional, presented a thorough derivation of the modified LPB equation and model by minimization of the system free energy for the case of planar geometry. A special emphasis was devoted to orientational ordering of water dipoles, taken into account in the expression for the free energy by rotational entropy. Our approach provides a distinct analytical description of the interplay between mean-field electrostatic and entropic effects arising from the mixing entropy of ions and rotational entropy of water dipoles in EDL.

The derived modified LPB equation in planar geometry is then generalized for arbitrary geometry and then used to calculate numerically the average orientation of water dipoles, relative permittivity ${\epsilon }_r$, magnitude of electric field strength, electric potential and the macroscopic (net) volume charge density of coions and counterions for a cylindrical geometry (in dependence on radial distance from the center of the tube).

Among other things it is indicated that in the saturation regime close to the charged surface, where the magnitude of electric field is very large (Figure 6), strong orientational water dipole ordering (Figure 7, upper panel) may result in a strong local decrease of permittivity (Figure 7, lower panel). The relative permittivity of the electrolyte solution decreases with increasing magnitude of the electric field strength.

Most interesting, we have shown that in the case of very narrow nanotubes the macroscopic (net) volume charge density of coions and counterions (${\rho }_{\mathrm{free}}$) at geometrical axis of the nanotube is not zero (Figure 4). In addition, in narrow nanotubes the water dipoles are partially oriented also close to the axis of the nanotube (Figure 7, upper panel), as schematically shown in (Figure 8). The potential importance of this phenomena for the transport through the narrow channels with the charged inner surface, specific only for very narrow nanotubes, should be investigated in the future. The channels in biological membranes can be an interesting example of such systems.

6. Materials and Methods

To solve Equation (52), a partial differential equation, we have used Comsol Multiphysics and its electrostatics stationary solver. The mesh consists of 4946 elements, the boundary condition (Equation (59)) was applied on the 2D cross-section of the nanotube and the geometry was solved for 10293 DoFs. The numerical results were solved using a stationary nonlinear solver (Automatic (Newton)), which implements a damped Newton’s approach, with a minimum damping factor of ${10}^{-6}$.