1. Introduction

Charged systems on an aqueous medium are necessary to understand a range of physical, chemical and biological processes. For example, in biological systems: the DNA condensation and packaging inside viral shells^1,2, the self-assembly of DNA into cationic liposomes³, the concentration of charged ions near a membrane channel and the interaction of proteins with membranes^4,5. In soft matter systems, charges systems are critical to the stabilization of colloidal dispersions, emulsions and help us to describe phenomena as wetting⁶.

One of the basic models for studying the properties of charged systems has been the system of two (or one) flat charged surface in an aqueous medium of dielectric constant ϵ. This type of system contains the main statistical information for an understanding of the charged surfaces and their interaction with a solution of charged particles, and allows to generalize its study to other systems (new geometries, new interactions). For this model different frameworks have been proposed with various levels of complexity to study its physical and chemical properties. For example, mean field models, such as Poisson-Boltzmann (PB) equation, that does not take into account the correlations between the different elements of the system have been widely used with good results only for systems with a low level of charges^6,7.

Theories that incorporate correlation effects between different elements of the system have been proposed at the level of integral equations, such as Ornstein-Zernike (OZ) models^8-11. For example, O. de la Cruz using the so-called anisotropic hypernetted-chain (AHNC) approximation for the Ornstein-Zernike found distinct ion-induced force in aqueous solution and introduce the concept of soft-structure to visualize the deformation of the local environment around the ions¹². For high charge systems a strong coupling theory has been proposed, see for example the work of R. Nezt^13,14 and R. Podgornik¹⁵. Where using different limits on the integral representation of the partition function, the model of weak coupling (WC) and strong coupling (SC) can be derived for a Coulombic fluid and its interaction with charged surfaces. In addition, there has been a variety of different simulations to study the system of charged surfaces. For instance, the method of molecular dynamics was used by A. Travesset to determine three regimes for distribution of ions and counterions (plasma, binding and uniform regimes) in a system of discrete charged surfaces¹⁶. Also, some new convergence techniques have been implemented to increase efficiencies in the evaluation of the electrostatic potential in Monte Carlo methods^17,18. In all of these studies, a constant concern has been to determine the ion concentration profiles generated between the charged surfaces and then to determine the system pressure.

Recently others effects have been incorporated into these theories and models. For example, charge image has been included in the calculation of particles profiles^19,20. Also, the effects of discrete charges on surfaces havestudied^16,17,21, and a charged lipid membrane with headgroups²² and models that take into account the effect of the dielectric constant of the medium^23-26. Finally, new membrane geometries (cylindrical, spherical) have been studied too^27,28.

In this paper, we consider a system of two surfaces introducing a new degree of freedom: the thickness of the surface charge distribution dm. Where ions can be freely distributed, see, Fig. 1. The two surfaces are separated a distance d 𝐿 . The total length is 2dm+dL. The model has three regions where charged particles can be distributed, all parts have the same dielectric constant 𝜖, but a new system with different values for each region is in progress. Using the method developed by Kjellander and Marčelja^9,10 we calculate the correlation between each particle and the particles profiles for each region. The concentration profile allows us to find the contact density and then using the value of the bulk concentration, the pressure between the surfaces can be obtained with the contact theorem⁶. In this work, we only calculate the contact density instead of the pressure. For systems with low surface charge the PB models and the integral theories seem to agree⁹, the problem arises when effects of high charge densities or multicomponent systems are studied²⁹. In previous studies, it has been shown that the pressure between the charged surfaces can become attractive. However, our new results suggest that this effect could decrease since the contact concentration profiles are of the same order as those found in the PB solution (see Fig. 3, for example). Our model does not permit direct comparison with other systems that work with different dielectric constant and salt, but we are working on a new system that allows to include more interactions between the elements which form the model.

Figure 1 Diagram showing three uniform surfaces where particles can be distributed. In the figure, we can see the distances (d_m, d_L) and the positions (Z₀, Z_m, Z_h, Z_L) used to calculate system properties. The charge of the particles in the region 1 and 3 are negative and positive in region 2. The distance D set an exclusion zone between regions and the particles cannot move from one region to another. All regions have the same dielectric constant ∈. 

Figure 3 The contact concentration profile ρ_L(Z_LL) as a function of the distance d_L using three charge densities σ₁ = 0:0938, C/m², σ2 = 0:267 C/m² and σ₃ = 0:348 C/m². The black line is our model with dm = 6:25 Å ; the red line is a model with d_m = 0:0 Å , and the blue line is the PB solution, Eq. 6. The insert shows ρ_L(Z_L) now as a function the distance dm for a charge density σ₂. The line with symbol -■- is for d_L = 9:25 Å , the symbol -●- is for d_L = 14:25 Å and symbol -▲- for d_L = 19:25 Å . The solid lines are the solution of PB equation. 

2. Model system

Figure 1 shows the system, consisting of three regions 1, 𝟐, and 𝟑, where particles are free to move and get separated from a core region of length D. In this paper, we considered that the three regions have the same dielectric constant ϵ and assume that there are no interactions of charge images. The surfaces are infinite in the radial direction 𝑟 and the concentration profiles ρ can vary along the 𝑧 coordinate. In the figure, we show some characteristics distances, Z_m and Z_L for contact profiles, Z₀ and Z_h for profiles at the midpoint distances. The system has positive particles in region 𝟐 and negative for regions 𝟏 and 𝟑. Thus the total net charge is then zero. We study the equilibrium properties of concentration profile 𝜌 as a function of the separation distance d_m, d_L and the charge density σ in the region 𝟏 (or 𝟑). In a previous paper¹¹, we study a system with a fixed charge density in the surface dm=0, in this new model the particles have a profile distribution in the region 𝟏 or 𝟑 and they interact electrostatically with particles from other regions. It is worth mentioning that the theory allows determining the correlation function between particles of different regions.

3. Theoretical Framework

We study the equilibrium distribution of particles with the use of the anisotropic HNC theory. The theory was originally proposed by Kjellander and Marčelja^9,10 and has been widely used to study the thermodynamic properties of interfaces and surfaces on charged system in a planar confinement^30,11,8. Here we only review the principal ideas of the theory. For more details see, for examples, Refs. ³¹,³² and ¹¹. In this theoretical scheme, the particle distribution 𝜌(𝐫) is calculated from

where β ψ(r) is the average external electrostatic potential, μ(r) is the excess chemical potential of the particles. The total correlation function h(r1,r2)=g(r1,r2)-1 is determined from the Ornstein-Zernike (OZ) equation

where c(r1,r2) is the direct correlation function. This integral equation is solved with the Hyper-Netted-Chain (HNC) closure approximation

where β=1/kBT, with kB the Boltzmann’s constant and T=298 K the absolute temperature. The set ofEqs.(1-3) is solved iteratively. The correlation function h(r1,r2) and c(r1,r2) are determined with Eqs. (3) and (2) and used as inputs for correcting the new profile given byEq.1, the process is repeated until self-consistency is achieved for two successive correlation function, with a small numerical error of about ±0.001 from each other solution. However, the convergence of the system of equations was poor for large distances dL where the bulk density is defined. The charged particles interact via a pairwise Coulomb potential

where 𝑒 is the elementary electric charge and r3D is the center to center distance of separation of two particles, and it is given in cylindrical coordinates and 𝑎 is the diameter of particles. In the model, the interaction with all images has been removed. A cut of for the long range tails of all correlation function due to the Coulomb interactions were performed as described in the Ref. ³². The total charge of the system must be zero. Therefore the concentration profiles in each region (see Fig. 1) satisfy the condition

where z is the normal direction to the surfaces. The density of surface charge 𝜎 is a parameter that matches with the integral σ=∫ρm(z) dz in the region 1. We compare the results of our model with those found from the Poisson-Boltzmann equation^7,6

when the charged particles of region 1 and 3 are continuous and located only in one surface. Also, we considered a model without structure dm=0) in region 1 and 3 as we previously studied¹¹.

4. Results

The concentration profile ρ(z) in each region (1, 2 and 3) is shown in Fig. 2(A) for three different charge densities (symbol -■- for σ1=0.0938 C/m², symbol -○- =0.267 C/m², and symbol -▲- σ3=0.348 C/m²). The distance of separation between the two surface is equal to dL=14.25 Å and dm=6.25 Å. The concentration profile in region 𝟐 increases as the charge density increases at σ3 the particle profile shows a local maximum at Z_h and two local minima, while for smaller charge densities 𝜎 1 and 𝜎 2 , there is only a local minimum. This structure is produced by pair correlations between different particles of the region 1, 2 and 3. The blue dashed line shows the concentration profile for a model with zero distance in the region 1, dm=0 Å, and in this case, there is no structure in the profile. For comparison, the continuous line shows the concentration profile found with the solution of PB Eq. (6), where distances were adjusted to consider point particles. Overall, for the concentration profile ρ(z), there is good agreement between the solution of PB and the results of our model with a thickness of the surface charge of dm=6.25 Å, while for the system with a thickness dm=0 Å there is a markedly different, see Fig. 2(A) and (B). However, at the contact profile ρ(zL), a significant difference is present between the three models studied. These differences are more pronounced for systems with small charge density, as we can see in Fig. 2(B) that show the concentration profiles for a σ1=0.0938 C/m². For this system, the PB solution and our model for a thickness of dm=6.25 Å are almost equal, and different from the profile with a thickness of dm=0 Å. The concentration profiles are symmetric in each region, but their contact values are not equal, this produces a concentration gradient which can induce instability in the membrane. This effect can be compensated by considering regions with different dielectric constant value.

Figure 2 (A). Concentration profile ρ_L(z) for region 2 and ρ_m(z) for the regions 1 and 3, with a distance d_L = 10:0 Å and d_m = 6:25 Å, using three charge densities σ1 = 0:0938 C/m2 (symbol -■-), σ2 = 0:267 C/m2 (symbol -○-) and σ3 = 0:348 C/m2 (symbol -▲-). The green line is the PB solution equation 6 and the dotted line is the answer for a system with d_m = 0 Å, see the Ref. ¹¹ (B) shows the concentration profile ρL(z) only in the region 2 for a surface charge density σ1, with the same parameters. 

The concentration profile in the contact position ρL(z=ZL) is shown in Fig. 3 as a function of the distance d_L in region 2 and three different values of charge density 𝜎 in region 1 ( or 3). The figure shows that when we have a large distance dL>20 Å of surface separation and small charge densities, there is not a difference in contact profile between PB solution, the model with a thickness of the surface charge of dm=6.25 Å and dm=0 Å. For distances dL<10 Å the profiles ρL(z=ZL) are the same for the PB solution and our model with dm=6.25 Å, however for the model with a thickness of dm=0.0 Å the solution overestimates the concentration at contact. When we increase the charge density σ2=0.267 C/m₂ and σ3=0.348 C/m₂ (for the model with a thickness of dm=0 Å and 𝜎 3 it was not possible to find a convergent solution of the system of equations) a significant difference appears in the contact profile betweenthe three model systems. The concentration profile has local

maxima and minima near the surface for a model with a thickness of dm=0Å and a charge density of σ2=0.267 C/m₂ , while the concentration profile of PB and the model with a thickness of dm=6.25 Å shows a continuous decay concentration profile. The figure also shows the contact profile, for σ3=0.348 C/m₂, but now with a significant difference between our model with a thickness of dm=6.25 Å and the PB solution. These results, show that the internal structure of the membrane in the region 𝟏 is critical in determining the interaction between the particle at the contact layers. By the contact theorem, this effect is essential for calculating the pressure between the two surfaces. In the Ref. ¹¹ was shown that the pressure could be negative to high surface charge densities. However, this effect might change if we consider the internal structure of the membrane. In contrast, the box in Fig. 3 shows the concentration profile at the contact position. In this case as a function of distances 𝑑 𝑚 in region 1 and for three different values of separation 𝑑 𝐿 . The PB solution that is always constant(solid line) was compared with a model of a thickness dm=6.25 Å (symbol line). We can notice a minimum difference in the contact profile as a function of 𝑑 𝑚 and a significant difference with the PB solution. The difference increased as we decrease the distance d_L, this means that the size of the region 1 is not a major factor in calculating the contact concentration profile in region 2.

Finally, in Fig. 4 we have the concentration profile at the contact points ρm(z=Zm) and ρL(z=ZL) as a function of the charge density 𝜎 in region 1 and a distance dL=14.25Å. For the system with dm=13.25 Å the contact concentration profiles ρm(z=Zm) is less than ρL(z=ZL) (symbol −▫−), ρL(ZL) > ρm(Zm). The PB solution (red line) match with ρL(z=ZL) only for small charge density σ<0.125 C/m 2 ), and therefore correlation effects are not important at these charge densities. However, when the charge density increases the PB solution shows differences with the model dm=13.25 Å and hence the effects on the system structure are important. A similar situation is present in a thin membrane, a distance dm=6.25Å , with the important difference that now the profiles satisfy the relation ρL(ZL) < ρm(Zm). This is expected because the particles have a small amount of space to spread, and also the correlation effects are present even at low densities. A remarkable fact is that the concentration profiles 𝜌 𝐿 (𝑧= 𝑍 𝐿 ) are the same for the two distances dm=13.25 Å and dm=6.25 Å). The box in Fig. 4 shows the concentration profile now at the middle positions z=Z0 and z=Zh of the regions 1 and 2 for the same set of parameters. Once again for the distance dm=13.25 Å we have ρL(Zh)>ρm(Z0) and for the distance dm=6.25 Å, ρL(Zh)<ρm(Z0), but in this case the contact concentration profiles ρm(z) and ρL(z) exhibit more structure (local minima and maxima appear). The concentration profiles in these positions ZL and Zh are important, because they are necessary for the calculation of the pressure between the surface^6,7, in a new work we are calculating the pressures for this system and other complex systems.

Figure 4 The contact concentration profile ρ_L(z_L) (blue line) and ρm(zm) (black line) as a function of the charge density σ. The lines with symbols -▲- and -■- are for a distance d_m = 6:25 Å and the lines with symbols -▲- and -■- are for d_m = 13:25 Å. The solid red line is the solution of PB equation 6. The insert shows the concentration profile for distances z = Z₀ in the region 1 and z = Z_h in the region 2, (see Fig. 1), using the same parameters. 

5. Conclusions

We have calculated the concentration profile and the correlation functions via the formalism of AHNC for a membrane system with internal structure (model with a thickness dm). We showed that the concentration profiles calculated by our model are in the same order of magnitude to the solution of the PB equation for small charge densities, and are different from a membrane system without structure¹¹ (model with a thickness dm=0 Å). For the case with high charge densities, important differences appear in the concentration profile between the three models (PB solution, a model with dm=0 Å and dm=6.25Å). We found that the correlation effects between the particles of regions 𝟏, 𝟐 and 𝟑 are important for small distances 𝑑 𝐿 and high charge densities σ>0.2 C/m²). These effects could be significant for the evaluation of the net pressure between membranes and may generate positive pressures, on systems that have been shown present attractions.