Charge and electric field distributions in the interelectrode region of an inhomogeneous solid electrolyte

A solid ionic conductor with cation conductivity in the interelectrode region is studied. Due to their large size, the anions are considered fixed and form a homogeneous neutralizing electric background. The model can be used to describe properties of ceramic conductors. For a statistical mechanical description of such systems, which are characterized by short-range Van der Waals interactions and long-range Coulomb interactions, an approach combining the collective variables method and the method of mean cell potentials is used. This formalism was applied in our previous work [Bokun G., Kravtsiv I., Holovko M., Vikhrenko V., Di Caprio D., Condens. Matter Phys., 2019, 29, 3351] to a homogeneous state and in the present work is extended to an inhomogeneous case induced by an external electric field. As a result, mean cell potentials become functionals of the density field and can be described by a closed system of integral equations. We investigate the solution of this problem in the lattice approximation and study charge and electric field distributions in the interelectrode region as functions of plate electrode charges. The differential electric capacitance is subsequently calculated and discussed.


Introduction
We would like to dedicate this paper to the memory of our good friend and colleague Vyacheslav Vikhrenko who passed away nearly a year ago. Together, we had an opportunity to work and discuss many aspects of the theory of solid electrolytes including the charge and electric field distributions in the interelectrode region. We should note that we also discussed with Vyacheslav the general ideas of the present paper at the initial stage of these studies.
Solid ceramic materials remain to be an area of intensive scientific activities. They are widely used in electrochemical sensors, rechargeable batteries, supercapacitors, memory devices, fuel cells and other technological applications [1][2][3]. These materials are characterized by high defect crystal structure and inhomogeneities inducing the emergence of electrical double layers close to intergrain boundaries [4][5][6][7][8].
In consequence, atomic level description of charge transfer in such systems proves to be a challenging problem. In this paper we restrict out treatment of this problem to the case of inhomogeneities caused by the electrical field of confining electrodes and neglect the adsorption phenomena on the surface of the solid electrolyte connected with specific interactions between ions and electrodes.
A peculiar feature of the materials under consideration is the presence of anions and cations that are drastically different in size. This fact allows us to apply the adiabatic approximation and reduce the problem to a one-component model. In the framework of this description, we consider thermal movement of cations against a homogeneous negatively charged background of fixed anions confined by the crystal lattice [1][2][3][4][5][6][7][8]. To this end, we apply a cell model in the framework of lattice description whereby we consider only the case of cation leaps between lattice sites.
Out of the theoretical approaches fit to study the ionic systems, notably the density functional theory [9][10][11][12][13][14][15], field theoretical methods [16][17][18][19], the collective variables method [20][21][22], and mean field theories, the latter remains the most popular one. The approach employed in this paper is self-consistent and takes into account pair correlations and long-range interactions as well as single particle potentials forming the mean field.
Solid ionic systems with conductivity induced by ions of the same sign [23] are characterized by short-range Van der Waals attraction and long-range Coulomb repulsion. They can also be viewed as systems interacting via a so-called SALR (short-range attractive and long-range repulsive) potential [24][25][26][27][28]. A special feature of our model is that the repulsive part of the SALR potential is of Coulomb nature. Solid electrolytes are characterized by high polarizability and high values of permittivity. As a result, the energy of short-range Van der Waals interaction is comparable to the energy of Coulomb interaction [4][5][6][7][8]. To address this fact, we combine two methods that had previously been used to describe systems of neutral and charged particles, namely the methods used in [23,[29][30][31][32][33] with the collective variables approach capable of taking into account long-range interactions [20,21].
In our previous paper [23], the homogeneous case was considered. Such a case appears in the absence of an electric field when local charge compensation takes place. Due to an external electric field, the cations in the interelectrode region cause an inhomogeneous distribution of mobile charges and the electric field. To describe the solid nature of the objects under study, the partition function is expanded with respect to mean potentials [29,30]. The approach [23] applied to homogeneous systems here is extended to inhomogeneous states emerging under an external field induced by the electrodes. In the present work the results [32,33] are generalized by taking into account the effect of correlations. In contrast to other approaches such as the density functional theory, where a density field functional is employed in combination with correlation functions, in this work we solve integral equations for distribution functions at a given density. This approach allows us to avoid the procedure of interpolating the distribution function for two different densities by two functions for respective densities of homogeneous states. This means that for an inhomogeneous system, the distribution functions are functionals of the density field. Hence, in this paper we develop a procedure for the construction and calculation of such functionals. Furthermore, we calculate thermodynamic potentials as functionals of a given density field. Unlike the approaches aimed at a fast direct calculation of experimentally measurable quantities based on the model functionals, the method proposed herein involves solving a system of integral equations. To this end, we resort to a lattice approximation which reduces the equations under consideration to a system of algebraic equations coupled only to the nearest neighbors of a given site.
The paper is arranged as follows. The generalization of the cell theory for the description of an inhomogeneous system is presented in section 2. In section 3 the reference distribution for the inhomogeneous case is considered. In section 4 the solution of the system of equations for the potentials and the lower order correlation functions in the lattice aproximation are discussed. In section 6 the charge and electric field distributions in the region between two plate electrodes are considered. In section 7 the results obtained are applied for the calculation of the electric capacitance. Finally, we conclude the paper by a few remarks in section 8.

Cell theory for the description of inhomogeneous systems
We partition a system of volume containing particles into cells ( ) of volumes = / . Each cell can be vacant or occupied by one or more particles. The state of some cell is characterized by the occupation numbers = 0 , 1 , 2 , ... subject to the condition where is the mean concentration of particles in the system.

23501-2
In contrast to our previous work [23], here we consider the case of an inhomogeneous system characterized by an inhomogeneous distribution of the density field. To model such a distribution, we take the mean occupation number to be different and extend the number of possible occupation numbers by describing the systems in terms of quasiparticles. We consider the case of an inhomogeneous medium, whereby the concentration is a function of the cell position and is denoted by . The distribution of particles inside a micro cell is characterized by a generalized coordinate that determines the potential of interaction between quasiparticles corresponding to different occupation numbers of the cell. Hence, each occupation number corresponds to a certain sort of quasiparticles. Particles of the first sort with = 0 corresponding to vacancies neither interact with one another nor with the particles of other sorts. Particles of the second sort with = 1 interact with one another and with other particles as real particles. The case of = 2 therefore corresponds to two real particles such that interaction of a particle of the third sort ( = 2) with a particle of the second sort ( = 1) corresponds to a system consisting of two real particles with coordinates 1 and 1 interacting with a particle of the Cartesian coordinate 1 . Such a description means that we consider the cell as a box with three different states: = 0 -no particles in the cell, = 1 -only one particle in the cell and = 2 -two particles in the cell. As a result, particles of the third sort have the energy of self-interaction ( 2 ) equal to the energy of interaction Φ( 1 , 1 ) of two real particles located in cell in the states 1 and 1 , respectively.
The potential energy of the system characterized by parameters that describe the distribution of particles across the system can be written in the form where ( ) = Φ( ) + ext ( ), and ext ( ) is the external field potential.
Now, we assume that the energy consists of the short-range and the long-range parts where Φ ℎ is the Van der Waals interaction which can be described, for instance, by the Lennard-Jones potential; Φ is the long-range Coulomb potential. Thermodynamic properties of the system can be found by expanding the configuration integral on the states of the reference system, which includes the long-range component of the energy, the energy of the external field, and the field of mean potentials corresponding to short-range interaction responsible for forming the crystal state of the system.

Reference distribution of the inhonomogeneous system
The Hamiltonian of the reference system can be written in the form where ( ) denotes summation over the cells surrounding cell .

23501-3
The potentials ( ) present in equation (3.1) have the meaning of the potential of the singleparticle mean force acting by the particles distributed in cell on a quasiparticle of sort fixed in cell in the state . In contrast to the case we had investigated before, whereby the potentials ( ) are functions of the mean density of particles, here ( ) are functionals of the density field described by the values of parameters . These quantities have the meaning of mean densities of quasiparticles of sort in cell . Due to the normalization condition, we have At a given density field, the distribution of quasiparticles across cell is characterized by the singlet distribution function Similar to [23], the pair distribution function can be presented in the form where ( , ) is the screening potential [33]. The correction factors −1 ( ) can be determined from the system of equations Extending the derivations performed in [23] to the present case, we see that the mean field potentials ( ) forming the crystal state can be found from the system of equations The relationships (3.6)-(3.8) make up a closed system of integral equations for the potentials ( ) and the normalization factors ( ). The latter makes it possible to find the singlet and the pair distribution functions of the reference system and subsequently calculate the internal energy and the Helmholtz free energy of the system.
Following the derivations presented in [23], one can write the free energy of the system as (3.10)

The solution of the system of equations for the potentials and the lower order correlation functions in the lattice approximation
In the framework of the lattice approximation, we take into account only the locations of the particles corresponding to the minimum of the potential energy of the system. In the crystal state this corresponds to the case when particles are located in the cell sites at = 0, 1 and between the sites at = 2 ( , * ), where and * correspond to the coordinates of the two particles which can be positioned on the opposite edges of the cell. Under this condition, due to the relation (3.4), the integral equations (3.6) reduce to algebraic equations in which and are the dummy variables and the indices and are fixed Unlike our previous work [23], in which only one equation was considered (due to all , being equal), here in the inhomogeneous case we have a system of equations with a given density distribution across the entire system. Given the fact that the screened potential is symmetric with respect to the densities and and depends only on the densities at sites and , equations (4.1)-(4.2) form a closed system of equations with respect to , and , . As a result, we obtain a closed system of six equations.
Notably, in the case when = 0, 1, this system is simplified and takes on the form Equations (4.4) and (4.5) can be closed by the similar equations with indices and swapped. In the equation (4.5) we account for the fact that due to the absence of interaction between vacancies, the relation (4.3) is not equal to one only when = 1 and = 1. We, therefore, simplify the notation: Due to definitions (4.7) and (4.8), the system (4.4) and (4.5) takes on the form The solution of equations (4.9)-(4.10) can be written as From equation (4.7) we also have that = . (4.14) From equation (4.14) one can see that in the lattice approximation the pair distribution function does not depend on single-particle mean potentials , ( ). Therefore, in this approximation the systems (3.6) and (3.8) become partially decoupled allowing to find the solution of equation (3.6) followed by equation (3.8). In the lattice approximation, the latter can be rewritten as * , Since the value of , is different from one, when = 1 and = 1, we put * = 1 ,1 allowing the system (4.15) to take on the form (4.4) and present the solution in the form (4.11).

Free energy of the inhomogeneous system
In the inhomogeneous case, the free energy can be expressed as the sum of local contributions. Each term in this series is a functional calculated in accordance with the algorithm outlined for a given density field. In our previous model [23], the free energy was a function of one variable, namely the mean density in the system. The equilibrium distribution of density across the system is found from the condition of an extremum of the functional. We note that in this case the number of fluctuating parameters is equal to the number of flat monomolecular layers which is proportional to the cubic root of the number of particles in the system. This contribution is important in the critical region, which we do not consider in this work.
The notations (4.8) and (4.14) allow us to write the short-range part of expression (3.9) in the form ℎ = 1 ln 1 + 0 ln 0 = 1 ln 1 0 + ln 0 For the long-range part of (3.9), considering that Φ ≠ 0 when = 1, we have where ℎ 2 ( ) is the correlation function and is the distance between sites and . When the long-range interaction corresponds to Coulomb interaction, we have where = 1/ is the reference temperature measured in units of the Lennard-Jones well depth, = ( 2 ) (4π 0 ) is the Bjerrum radius measured in units of cell parameter . The distances in equation (5.5) and herein are, therefore, measured in units of .
According to [23], from equations (4.11) and (4.14) we can write From equations (5.5) and (5.6) we see that it is convenient to measure the parameter in units of , i.e., * = / . As a result, due to equation (5.5), the quantity in equation (5.3) takes on the form We should note that in the exponent for in equation (5.6) we have a typical Debye-Hückel (DH) form for the screening potential of point ions [34]. However, there is a considerable difference for the dependence of the inverse Debye radius on the concentration of mobile ions. In the DH theory, this radius is proportional to the square root of concentration while in the case considered here it is proportional to the square root of the product (1 − ) [33]. The results of the two theories coincide only when ion concentrations are low. Another important result in the framework of the theory developed was obtained in [33] from the condition of equilibrium between bulk and surface properties of mobile ions. As a result, it was shown that the distribution of ion concentration can be presented in the form which looks like a Fermi-Dirac distribution corrected for the interparticle interaction contribution and concentration dependent multiplier in front of the exponent. In equation (5.9), is the electric potential at lattice , the subscript indicates the position of the corresponding lattice site, is the elementary electric charge, = 1/ , is the temperature, is the Boltzmann constant, is the surface excess part of Coulomb interaction and is the short range part of the interpaticle interaction contribution. As we will show in the next two chapters, both of the discussed effects are very important for the description of electrophysical properties of a double layer between a solid electrolyte and a hard charged electrode.

Charge and electric field distributions in the region between two plate electrodes
We consider an electrolyte confined between two parallel electrodes and divide the interelectrode region into layers parallel to the plates. As a result, the variables 1 and 1 have the meanings of the average concentrations of particles in layers and , respectively. The total number of layers is , i.e., = 1, 2, ..., . The long-range component of the free energy, given by the formula (5.2) consists of correlated and non-correlated parts The correlated part is expressed in terms of the correlation function while the non-correlated part corresponds to the energy of Coulomb interaction between capacitor plates modelling molecular layers. It is convenient to calculate this energy via the Poisson equation [33]. As a result, we have where * = / , = − , and is the average concentration of cations in the homogeneous electrolyte. Equation (6.2) holds for the case when electrolyte layers are numbered consecutively from = 2 to = − 1. The first and the last layers model the electrode plates, whose charges are used as boundary conditions. Hence, the energy of the external electric field at location is given by the following expression The plates have opposite charges ( 1 + = 0), is the index of the middle layer. The model, therefore, is analogous to a parallel-plate capacitor.
As a result, the electric potential induced by the inner electrolyte charges and the electrodes is described by the expressions According to equations (5.2)-(5.4), the correlated part of the free energy is expressed in terms of the short-range correlation function. As the first approximation for , we take into account only the first neighbors. Due to equation (5.8), we can write for the cubic lattice where , is the distance between the nearest neighbors in layers and .

23501-8
We consider the case of small deviations of charge carrier concentrations from the average value . Expansion of the function (5.4) on leads to the equations =˜( −1 + 4 2 + +1 ). (6.8) The linear terms in expansion (6.8) vanish since in the expression for the total free energy such terms cancel due to the condition of electroneutrality of the system An expression similar to (6.8) also appears in the expansion (5.1). Since in the approximation chosen, the Van der Waals and the correlated Coulomb parts have the same structure, they can be combined by introducing an expansion coefficient˜. Differentiation of the free energy with respect to the variable leads to an expression for the chemical potential˜in a local layer . It is convenient to use non-dimensional chemical potential such that =ã s well as non-dimensional potentials Ψ =Ψ and =˜. This potential consists of a combinatorial part a part corresponding to non-correlated Coulomb interaction between charges (i.e., an internal electric potential) Charge distribution in the interelectrode electrolyte film is presented in figures 1-3 in the right-hand panels, while the left-hand panels correspond to the case of non-interacting ions. Due to the concentration dependence of interparticle correlation term , the equation (5.9) has a more complex concentration distribution compared with the traditional Fermi-Dirac distribution for non-interacting particles when = 0. In the considered case, this term is taken into account by the usual iterative procedure. Although the condition of electroneutrality is always satisfied, for the strong charge 2 = 0.5, the distribution becomes more asymmetric. For the non-interacting particles ( figure 3 on the left), the cations are strongly repelled toward the negatively charged electrode and the non-adjacent electrolyte layers have a nearly constant charge value determined by the average ion concentration = 0.1. When interaction is taken into account, the particles are redistributed in a much more symmetric manner. In the latter case, one can also observe the tendency for the formation of a neutral plateau in the middle of the interelectrode region.   Electric field distributions induced by internal charges in the interelectrode region are shown in figure 4 for different numbers of layers as well as different plate charges 2 . One can see that the potential profile becomes more linear for wider slabs as well as higher values of electrode plate charges.
The chemical potentials for different values of are presented in figure 5. In each case, we show only the curve for one particular value of 2 , noting that a similar quasi-constant behavior takes place for

Electric capacitance
The differential electric capacitance is defined as where is the variation of the charge due to a variation of the potential. We can find the charge in the bulk region by summing up the charges in all the layers The potential Ψ is, therefore, the difference of the potentials of the boundary layers To calculate the quantity (7.1), we find the variations (7.2) and (7.3) caused by the variation of the charge of the electrode layers 2 = ℎ. The quantities and Ψ can be determined as (7.5) Equations (7.2) and (7.3) tell us that in the case of many layers the quantities (7.4) are the differences of large numbers, which requires high precision when determining these differences from the condition of a constant chemical potential (see figure 6). For this reason, the accuracy of calculations sharply diminishes at high potentials and at large numbers of variables. We demonstrate the potential dependence of the differential capacitance for the case when the total error of the calculations is approximately five percent ( figure 6). To get the values of at large and Ψ, we need to address the problem of a decreasing accuracy with increasing values of these parameters. In figure 6, we present the diffetential capacitance as a function of the potential at different mean concentrations of the cations. At small values of the potential, the differential capacitance differs very little from the constant value corresponding to the linear theory.
Smoother changes of the differential capacitance as a function of the potential with increasing mean concentrations is the consequence of equation (6.10) at = 0.5 containing the inflection point, which widens the applicability region for the linear description. This also explains the shift of maxima to the right as the mean concentration increases. A fall of the differential capacitance at high potentials is the result of saturation in the distribution of charges in the regions close to the electrodes. In order to increase the density of cations in one region and decrease in another, one should, therefore, increase the potential as the density of charges in the sample increases. When the potential exceeds some critical value, the charge can no longer go up and the differential capacitance falls to zero as in figure 6. It is also worth noting that the curves obtained can be extended to the negative values of the electric potential. For the  case under consideration, these curves will be mirror symmetric with respect to those presented for the positive values of the electric potential. Such a shape of the differential capacitance as a function of the electric potential is known in the literature as camel-like shape and has recently been widely discussed in connection with the room temperature ionic liquids [35,36].

Conclusions
The paper provides a statistical mechanical description of an inhomogeneous equilibrium state of a solid electrolyte by extending the approach constructed earlier for a homogenous state [23]. As in our previous work, we consider the case of mobile cations in a neutralizing field of fixed anions. The Van der Waals interaction is accounted for by means of mean cell potentials while the long-range Coulomb interaction is described in the framework of the collective variables formalism. In addition, the possibility of taking into account particle fluctuations not only around lattice sites but also between the sites is shown. An expression for the free energy as a functional of the particle density field is derived. Similar to [32], the electric potential distribution and the charge density distribution in the interelectrode region of the electrolyte are computed though with additional, more accurate considerations. Notably, the expression for the chemical potential used in [32], which is characteristic of ideal solutions, is replaced by the calculation of the free energy functional developed in [23] and extended here for inhomogeneous systems. This enabled us to take into account correlations and interparticle interactions more accurately using mean potentials. Furthermore, in contrast to [32], the equilibrium states are calculated based on numerical minimization of the free energy and on the enforcement of the condition that the total chemical potential across the system should be constant. Such an approach allowed us to avoid the problems that arose in [32] from the direct solution of the system of equations, which proved to be poorly defined. The results obtained are in qualitative agreement with those reported in [25][26][27][28].
An important aspect of the present paper is the use of the results of our previous work [33], namely the expression for the screened potential. This potential is of the same form as that in the Debye-Hückel theory but with a different screening radius. In contrast to the Debye-Hückel theory, in which the inverse screening radius is proportional to the square root of the concentration of mobile ions , in the case considered it is proportional to the square root of the product (1 − ) [33]. The results of the two theories coincide only when the concentrations of mobile ions are small. Another principal result of this work was also obtained in [33] and comes down to the equation (5.9) for ionic concentration distribution. This expression is presented in the form of the Fermi-Dirac distribution corrected for the contribution from interparticle interaction and a concentration-dependent factor in front of the respective exponent. Both effects proved to be highly instrumental in the description of electrophysical properties of the electric double layer for the model under consideration. Notably, we showed that the electric differential capacitance as a function of the electric potential calculated employing these effects has a typical camel-like shape. At small values of the electric potential, the capacitance gradually increases, reaches a maximum, and then falls sharply as the potential is further increased, which is consistent with saturation of the electric charge in the near-electrode region in accordance with the distribution (5.9). A more gradual change of the capacitance at higher concentrations of mobile ions is due to the fact that the relation (6.10) at = 0.5 contains an inflection point that widens the region of linear description. This is also the reason for the shift of the curve maxima and the decrease of their values in the direction of higher potentials.