Electrolyte solutions at heterogeneously charged substrates

Electrolyte solutions at heterogeneously charged substrates

Maximilian Mußotter mussotter@is.mpg.de    Markus Bier bier@is.mpg.de    S. Dietrich Max Planck Institute for Intelligent Systems, Heisenbergstr. 3, 70569 Stuttgart, Germany IV Institute for Theoretical Physics, University of Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany
March 28, 2018
Abstract

The influence of a chemically or electrically heterogeneous distribution of interaction sites at a planar substrate on the number density of an adjacent fluid is studied by means of classical density functional theory (DFT). In the case of electrolyte solutions the effect of this heterogeneity is particularly long ranged, because the corresponding relevant length scale is set by the Debye length which is large compared to molecular sizes. The DFT used here takes the solvent particles explicitly into account and thus captures phenomena, inter alia, due to ion-solvent coupling. The present approach provides closed analytic expressions describing the influence of chemically and electrically nonuniform walls. The analysis of isolated -like interactions, isolated interaction patches, and hexagonal periodic distributions of interaction sites reveals a sensitive dependence of the fluid density profiles on the type of the interaction, as well as on the size and the lateral distribution of the interaction sites.

I Introduction

Detailed knowledge of the structure of electrolyte solutions close to solid substrates is of great importance to numerous research areas and fields of application, ranging from electrochemistry Bagotsky2006 (); Schmickler2010 () and wetting phenomena Dietrich1988 (); Schick1990 () via coating Wen2017 () and surface patterning Vogel2012 (); Nee2015 () to colloid science Russel1989 (); Hunter2001 () and microfluidics Lin2011 (); GalindoRosales2018 (). The vast majority of models describing fluids in contact with substrates consider the latter as uniform with respect to the wall-fluid interaction. This approximation is commonly made partly due to a lack of experimental data on the actual local properties of the substrate under consideration and partly for the sake of simplicity. For fluids comprising only electrically neutral constituents and uncharged walls, assuming uniform substrates is typically an acceptable approximation because, in the absence of wetting transitions, heterogeneous substrate properties influence the fluid only on length scales of the order of the bulk correlation length Andelman1991 (), which, not too close to critical points, is of the order of a few molecular diameters. In contrast, nonuniformities of the surface charge density of charged substrates in contact with dilute electrolyte solutions influence the fluid on the scale of the Debye length, which is much larger than the size of the molecules. Furthermore, the charged sites of substrates, such as mineral surfaces and polyelectrolytes, are lateral distances apart which are typically comparable with the Debye length of the surrounding fluid medium Chen2005 (); Chen2006 (); Chen2009 (). Hence, the assumption, that substrates in contact with electrolyte solutions carry a uniform surface charge density, is in general untenable.

In recent years considerable theoretical interest has emerged in the effective interaction between two heterogeneously charged walls (which typically are the surfaces of colloidal particles) mediated by an electrolyte solution Naji2010 (); Ben-Yaakov2013 (); Naji2014 (); Bakhshandeh2015 (); Ghodrat2015a (); Ghodrat2015b (); Adar2016 (); Adar2017a (); Adar2017b (); Ghosal2017 (); Zhou2017 (). In contrast to uniform substrates, this effective interaction can lead to lateral forces, in addition to the common ones in normal direction. However, all the studies cited above model the solvent of the electrolyte solution as a structureless dielectric continuum. This approach precludes coupling effects due to a competition between the solvation and the electrostatic interaction, which are known to occur in bulk electrolyte solutions Onuki2004 (); Bier2012a (); Bier2012b (). In particular, in the presence of ion-solvent coupling and far away from critical points, correlations of the solvent number densities in a dilute electrolyte solution decay asymptotically on the scale of the Debye length. Consequently, under such conditions, nonuniformities of the nonelectrostatic solvent-wall interaction can influence the structure of an electrolyte solution close to a wall and hence the strength and range of the effective interaction between two parallel plates immersed in an electrolyte medium on a length scale much larger than the molecular size. This mechanism differs from the one studied in Refs. Naji2010 (); Ghodrat2015a (); Adar2016 (); Zhou2017 (), in which the walls are locally charged but overall charge neutral.

In the present analysis a first step is taken towards a description of the structure of electrolyte solutions close to chemically and electrically nonuniform walls in terms of all fluid components. The natural framework for obtaining the fluid structure in terms of number density profiles of solvent and ion species is classical density functional theory Evans1979 (); Evans1990 (); Evans1992 (). Here, the most simple case of an electrolyte solution, composed of a single solvent species and a single univalent salt component, is considered far away from bulk or wetting phase transitions. Moreover, the spatial distribution of nonuniformities of the chemical and electrostatic wall-fluid interactions can be arbitrary but their strengths are assumed to be sufficiently weak such that a linear response of the number density deviations from the bulk values is justified. This setup allows for closed analytic expressions which are used to obtain a first overview of the influence of ion-solvent coupling on the structure of electrolyte solutions in contact with chemically or electrically nonuniform walls. This insight will guide future investigations of more general setups within more sophisticated models.

After introducing the formalism in Sec. II, selected cases of heterogeneous walls are discussed in Sec. III. Due to the linear relationship between the wall nonuniformities and the corresponding number density deviations from the bulk values, the latter are given by linear combinations of elementary response features, which are discussed first. Next, two main cases are studied: wall heterogeneities, which are laterally isotropic around a certain center and wall heterogeneities, which possess the symmetry of a two-dimensional lattice; the study of randomly distributed nonuniformities Naji2010 (); Naji2014 (); Bakhshandeh2015 (); Ghodrat2015a (); Ghodrat2015b () is left to future research. For both cases various length scale regimes are discussed, which are provided by the bulk correlation length of the pure solvent, the Debye length, and a characteristic length scale associated with the wall nonuniformities. Conclusions and a summary are given in Sec. IV.

Ii Theoretical foundations

ii.1 Setup

Here, the influence of a chemically and electrically nonuniform wall on the fluid density is studied. In spatial dimension the system consists of an impenetrable planar wall for and a fluid for , both parts being macroscopically large. In the following, the space occupied by the fluid is denoted by ; the positions are uniquely decomposed into the lateral components and the normal component relative to the wall surface at . The size of the wall and the extent of the system in normal direction are both assumed to be macroscopically large. The fluid is an electrolyte solution composed of an uncharged solvent (index “1”), univalent cations (index “2”), and univalent anions (index “3”). Two types of interactions between the fluid and the wall are considered: (i) electric monopoles at the wall surface () and the fluid ions, giving rise to an electrostatic interaction, (ii) all other contributions, in particular those due to nearest-neighbor-like chemical bonds, referred to as nonelectrostatic interactions.

ii.2 Density functional theory

We use density functional theory Evans1979 (); Evans1990 (); Evans1992 () in order to determine the equilibrium number density profiles of the three fluid species. Since we focus on length scales larger than the sizes of the fluid particles and on weak wall-fluid interactions, the following dimensionless density functional within a Cahn-Hilliard-like square-gradient approximation Cahn1958 () is considered:

(1)

where is the inverse thermal energy, are the chemical potentials of the three species, is a phenomenological parameter with dimension , which can be inferred from microscopic models (see Sec. III.1), is the vacuum permittivity Lide1998 (), is the relative dielectric constant of the fluid, is the electrostatic potential at , and describes the strengths of the nonelectrostatic wall-fluid interactions at for the three species. Note that for the sake of simplicity, the coupling of number density gradients of different particle types is neglected in Eq. (1) (see Sec. III.1). In the present study the bulk state is considered to be thermodynamically far away from any phase transition so that the local contribution of the density functional in Eq. (1) can be safely expanded around up to quadratic order in :

(2)

where the local part of the interactions between different types of particles is captured by the real-valued, symmetric, and positively definite -matrix (see, c.f., Eq. (30)). Furthermore, specifies the grand potential density, evaluated for the equilibrium bulk densities , which equals minus the bulk pressure ; in the following its value is of no importance. For a given equation of state the bulk densities are free parameters of the model. Finally, the electrostatic potential , which enters into Eq. (1) on a mean-field level via the electric field energy density, fulfills the Poisson equation

(3)

for with the boundary conditions

(4)

for , where is the surface charge density at the point on the wall surface (), and denotes the valences of the fluid species.

The Euler-Lagrange equations, corresponding to the minimum of the density functional specified in Eqs. (1)–(4), can be written as

(5)

and

(6)

for with the boundary conditions given by Eq. (4) and by

(7)

for , where is the Bjerrum length.

The linear nature of the Euler-Lagrange equations (5) and (6) tells that the quadratic (Gaussian) approximation of the underlying density functional in Eqs. (1) and (2) corresponds to a linear response approach. It is widely assumed and in some cases it can be even quantified (see, e.g., the quantitative agreement between the full and the linearized Poisson-Boltzmann theory in the case that the surface charges are smaller than the saturation value Russel1989 (); Bocquet2002 ()) that for sufficiently weak wall-fluid interactions linear response theory provides quantitatively reliable results.

ii.3 Solution of the Euler-Lagrange equations

Instead of solving the Euler-Lagrange equations in Eqs. (5) and (6) as differential equations for the profiles and as functions of , it is convenient first to perform Fourier transformations with respect to the lateral coordinates . The resulting transformed profiles

(8)

and as functions of and can be combined in the four-component quantity so that Eqs. (5) and (6) can be written as

(9)

where and is the second derivative of with respect to the coordinate normal to the wall. Note that the components of are quantities of different dimensions: and . This does not allow for the formation of a scalar product of two vectors of the type ; however, in the following scalar products will not occur. Writing with (i.e., is a diagonal matrix with these entries), one obtains

(10)

The -matrix is independent of and it is symmetric but not real-valued, because the bottom entry of is imaginary. is not a normal matrix, i.e., it does not commute with its adjoint matrix , and hence it does not possess an orthogonal basis composed of eigenvectors. However, the actual structures of the matrix and of the vector used below guarantee the existence of a nonorthogonal basis of eigenvectors of the matrix with respective positive (real-valued) eigenvalues (see Appendix A). Expanding the vector in this basis ,

(11)

leads to Eq. (10) in the form

(12)

with the solution

(13)

where the second boundary conditions in Eqs. (4) and (7) have been used. Therefore, the solution of Eq. (10) can be expressed as

(14)

Finally, the first boundary conditions in Eqs. (4) and (7) can be expressed as

(15)

with

(16)

as the Fourier transform of with respect to the lateral coordinates and as the Fourier transform of . Note that, as for , the components of are quantities of different dimensions: and . From Eqs. (14) and (15) the coefficients can be determined. Note that according to Eqs. (14) and (15) the coefficients and hence the profiles and depend linearly on the nonelectrostatic wall-fluid interactions and the surface charge density . Such a linear response of the number density profiles inside the fluid to the wall properties requires weak wall-fluid interactions, which is assumed in the present study and which is consistent with the quadratic form of the density functional in Eqs. (1)–(4).

Iii Results and Discussion

iii.1 Choice of parameters

The present study discusses the influence of the wall-fluid interactions, represented by the nonelectrostatic wall-fluid interactions and the surface charge density , onto the number density profiles in the adjacent fluid. Applying density functional theory as described in Sec. II requires knowledge of the bulk number densities , the parameter , and the coupling matrix all of which are bulk quantities or characterize them.

In the bulk local charge neutrality holds, i.e., . Hence the equilibrium bulk state is determined by the temperature , the number density of the solvent, and the bulk ionic strength .

Figure 1: The nonelectrostatic interaction between fluid particles is modeled by a square-well pair potential displayed in panel (a), where denotes the distance between the centers of two spherical particles. For a hard core repulsion prevents the overlap of two particles. For two particles attract each other with a constant interaction energy . At distances there is no nonelectrostatic interaction. Panel (b) sketches the decomposition of the nonelectrostatic interaction potential according to the scheme due to Barker and Henderson into the hard core repulsion and the attractive well , which is used in Sec. III.1 in order to obtain the parameters entering the Cahn-Hilliard square-gradient density functional in Eq. (1).

In order to obtain expressions for the parameter and for the coupling matrix in terms of experimentally accessible quantities, in a first step the pure, ion-free solvent is considered, the particles of which interact only nonelectrostatically. Here this nonelectrostatic interaction between solvent particles at a distance is modeled by a square-well pair potential as displayed in Fig. 1(a). At small distances a hard core repulsion prevents the overlap of two particles. At intermediate distances two particles attract each other with an interaction energy , and at distances the nonelectrostatic interaction vanishes. According to the scheme due to Barker and Henderson Hansen1986 (), the interaction potential can be decomposed as into the hard core repulsion and the attractive well (see Fig. 1(b)). The microscopic density functional for the pure solvent (species ) in the bulk can be approximated by the expression

(17)

The contribution (here within local density approximation (LDA)) is due to the reference system governed solely by the hard core interaction :

(18)

Here, contributions of external potentials are neglected, because they do not contribute to the bulk parameters and . The second term on the rhs of Eq. (17) is (within random phase approximation (RPA) Evans1979 ()) the excess free energy functional due to the square-well attractive interaction :

(19)

In Eq. (18) is the thermal de Broglie wave length, denotes the chemical potential of species 1, and is the excess free energy per volume of the reference system governed by the hard core interaction .

Following Cahn and Hilliard Cahn1958 (), Eq. (19) can be approximated by a gradient expansion:

(20)

with the -th moment of the pair potential in units of ,

(21)

which, for the present form of , leads to

(22)

From Eq. (20) one obtains the gradient expansion of in Eq. (17):

(23)

with the local contribution

(24)

The comparison of Eq. (23) with Eq. (1) renders an expression for the parameter in terms of parameters of the interaction potential (see Fig. 1):

(25)

By expanding up to quadratic order in the density deviation from the equilibrium bulk density , which is a solution of the Euler-Lagrange equation

(26)

one obtains

(27)

The comparison with Eq. (2) leads to the matrix element

(28)

of the matrix , where the first term on the rhs stems from the ideal gas contribution of the solvent particles. The argument of the second term, which is due to the hard core interaction , is a measure of the packing fraction .

The analogue of Eq. (23) for the nonelectrostatic contribution of all three particle species is given by the first line of Eq. (1) with the local contribution (compare Eq. (24))

(29)

where denotes the total number density. Note that Eq. (29) assumes, that all interactions among the species are the same (see Eq. (21)). This implies that the last term in Eq. (29) takes the form . By expanding up to quadratic order in the density deviations from the equilibrium bulk densities one finally finds (see the steps leading to Eq. (28))

(30)

where denotes the total number density in the bulk. In the present study the hard core excess free energy per volume is chosen as the one corresponding to the Carnahan-Starling equation of state Hansen1986 ():

(31)

here with the packing fraction .

Accordingly, from Eq. (24) one obtains the following equation of state of the pure solvent:

(32)

Its derivative with respect to the number density , using the relation with the isothermal compressibility , yields

(33)

As an example we consider water at room temperature and ambient pressure (which corresponds to the number density and the isothermal compressibility Lide1998 ()) with relative dielectric constant , i.e., with Bjerrum length , and with a univalent salt of ionic strength . The strength of hydrogen bonds, which generate the dominant attractive interaction contribution, is of the order of Naik2000 (); Lide1998 (). Using these data, one obtains from Eqs. (32) and (33) the bulk packing fraction as well as and . In the following the Debye length

(34)

is used as length scale, which, for the present choice of parameters, is .

In the case of a pure solvent (), in the bulk the density two-point correlation function fulfills an equation similar to Eq. (5):

(35)

Note that the similarity between Eqs. (5) and (35) is due to the asymptotic proportionality between density deviations and two-point correlation functions (Yvon equation) Hansen1986 (). From Eq. (35), one can readily infer the relation

(36)

for the solvent bulk correlation length, which characterizes the exponential decay of . For the present choice of parameters, one has so that .

iii.2 X-ray scattering

In the following subsections the Fourier transforms of the profiles of the density deviations as functions of the lateral wave vector and of the normal distance from the wall are discussed in detail. However, from the experimental point of view, it is challenging to directly obtain the -dependence of the density profiles. One of such direct methods is total internal reflection microscopy (TIRM) Walz1997 () in the context of the structure of colloidal suspensions close to (optically transparent) substrates. In contrast, for molecular fluids, as the ones considered here, such direct methods are not available and one has to resort to, e.g., X-ray scattering techniques AlsNielsen2001 (); Dietrich1995 (). As X-rays are predominantly scattered by the electrons of the fluid molecules one has to consider the electron number density

(37)

for with the number of electrons per molecule of particle species , the bulk electron density , and the deviation of the electron number density from its bulk value. The X-ray scattering signal for scattering vector is proportional to with the double Fourier transform of the electron number density profile in both the lateral and the normal direction Hansen1986 ().

Figure 2: Density distribution of the solvent (panel (a)) and of the ions (panel (b)) as function of the distance from the wall and of the absolute value of the lateral Fourier wave vector in units of the inverse Debye length (see Eq. (34)). The plane is given by the positions of the fluid particle centers when the surface-to-surface distance between the hard wall and the hard particles vanishes. The data correspond to the boundary condition (see Eqs. (15) and (40)). The physical situation corresponding to this boundary condition is an attraction (see Eq. (42)) of the solvent particles by the wall located at the origin of the wall. Concerning the remaining relevant parameters see Sec. III.1.

For common specular X-ray reflectivity measurements, i.e., for , the normalized intensity reflected as function of the normal wave number is given by AlsNielsen2001 (); Dietrich1995 ()

(38)

where denotes the Fresnel reflectivity of an ideal, step-like planar interface Jackson1999 (), and where the notation with

(39)

has been used. Moreover, off-specular diffuse X-ray scattering () at grazing incidence (GIXD, ) yields scattering intensities which are proportional to Dietrich1995 (). Hence, as the double Fourier transforms in Eq. (39) of the density deviation profiles are of direct experimental relevance, they will be discussed in the following in parallel to the single Fourier transforms . Note that due to , in Figs. 3 and 4 its absolute value is shown.

iii.3 Basis vectors of boundary conditions

As mentioned above, the linear nature of the relationship between wall nonuniformities and the resulting number density deviations leads to the possibility of describing the latter in terms of linear combinations of elementary response patterns. These elementary response patterns correspond to four basis vectors, e.g., , which span the four-dimensional space of boundary conditions in Eq. (15). Therefore, as a first step to study the influence of wall inhomogeneities onto the fluid, these four distinct boundary condition vectors are studied. The first one of these vectors is given by

(40)

which requires (see Eq. (15))

(41)

and which in real space corresponds to the boundary condition

(42)

This boundary condition corresponds to an attractive, -like interaction of the wall with the solvent located at the origin. Solving the Euler-Lagrange equations for this boundary condition, one finds the density distribution for the solvent and for the ions, as shown in Figs. 2 (a) and 2 (b), respectively. Since the boundary condition corresponds to a constant in Fourier space, the density distribution depends on only. Actually, the solution for this system is proportional to the first column of the Green’s function, which is a -matrix, of the differential operator corresponding to Eqs. (5) and (6).

Figure 2 illustrates that for fixed the density deviations from the bulk value increase for smaller normal distances from the wall and, for fixed , also upon decreasing the absolute value of the lateral wave vector . The behavior with respect to the normal distance from the wall can be anticipated, because the effect of the interaction between the wall and the fluid is expected to decrease with increasing distance from the wall. Moreover, also the behavior with respect to is as expected, because a strong attraction at the origin leads to a radially decreasing density deviation, which in Fourier space corresponds to a maximum at the origin. In order to allow for a quantitative analysis of the behavior of the density deviations, Fig. 3 shows various cuts through the data of Fig. 2 along several lines.

Figure 3: Density profiles of the solvent (left column, panels (a), (c), and (e)) and of the ions (right column, panels (b), (d), and (f)) as functions of the normal distance from the wall (top row, panels (a) and (b)), of the absolute value of the lateral wave vector (middle row, panels (c) and (d)), and of the wave number in normal direction (bottom row, panels (e) and (f), with being dimensionless) in corresponding units of the Debye length and the inverse Debye length, respectively (see Eq. (34)). Note that due to , in panels (e) and (f) the absolute values are shown. In each graph, there are three profiles shown corresponding to three values of the other relevant variable. Therefore the profiles correspond to cuts through Figs. 2 (a) and (b) at various positions and in different directions. In this case the boundary condition is (see Eqs. (15) and (40)), corresponding to a -like nonelectrostatic attraction of the solvent particles at the origin of the wall (see Fig. 2 and Eq. (42)). The graphs show, that the density deviations of the ions are proportional to the ones of the solvent, although different in sign. Since only the solvent particles are attracted by the wall, it is favorable for the system to increase their density close to the wall. However, due to the hard core nature of the particles and the equality of the interparticle attraction for all pairs of particles, the increase of solvent particles leads to an extrusion of ionic particles, leading to decreased ion densities at the wall. However, the density deviations of the ions are much weaker. For the remaining relevant parameters see Sec. III.1.

Figures 3 (a), (c), and (e) show the density profiles for the solvent and Figs. 3 (b), (d), and (f) the ones of the positive ions, which in this case are the same as the profiles for the negative ions. This equivalence is due to the nature of the boundary conditions in this special case, which in real space lead primarily to an increased solvent density close to the origin at the wall. The ions, however, react only indirectly via the solvent, with which both ion types interact in the same way. Since the solvent particles get attracted by the wall, it is favorable to increase their density close to the wall. Due to the hard core nature of the particles, the space occupied by the solvent particles is blocked for the ions. Since the solvent is attracted by the wall and the interparticle attraction is the same for all pairs of particles, this leads to an extrusion of the ions in favor of an increased number of solvent particles. Figures 3 (a) and (b) show the density deviations as function of the normal distance from the wall for three values of , i.e., Figs. 3 (a) and (b) correspond to horizontal cuts through Figs. 2 (a) and (b), respectively. For fixed , as in Figs. 2 (a) and (b), Figs. 3 (a) and (b) clearly show an exponential decay of the density deviation for increasing distances from the wall. In contrast, Figs. 3 (c) and (d) show vertical cuts through Figs. 2 (a) and (b), i.e., density profiles as functions of the absolute value of the lateral wave number for three normal distances from the wall. The dependence of these profiles on the absolute value of the lateral wave vector implies a laterally isotropic decay of the density deviations in real space. The third pair of graphs, Figs. 3 (e) and (f), shows the Fourier transforms of the density profiles of Figs. 3 (a) and (b), being additionally Fourier-transformed with respect to the normal direction , which leads to the Fourier transforms in terms of the lateral wave vector and the normal wave number , respectively. All curves in Figs. 3 (c)–(f) exhibit a Lorentzian shape as functions of and , respectively. These Lorentzian curves in Fourier space correspond to exponential decays in real space in lateral or normal direction. The curves in Figs. 3(c) and (d) show widths of half height which decrease with increasing normal distance , i.e., the lateral decay length in real space increases with increasing distance from the wall. This implies that the density distribution broadens upon moving away from the source of the perturbation. The curves in Figs. 3(e) and (f) exhibit widths of half height which increase with the lateral wave number , i.e., the normal decay length in real space decreases with increasing lateral wave number. Consequently, the range of influence of rapidly varying modes of wall heterogeneities onto the fluid is shorter than that of slowly varying modes. This relationship can also be inferred from Figs. 3(a) and (b)). From the above discussions and from Fig. 3 one can conclude, that the response of all species to a simple attraction of nonelectrostatic type is the same up to a proportionality factor. This is confirmed by studying in addition the boundary conditions and ; these results are not shown here.

After having discussed the effects of the boundary condition via Figs. 2 and 3, the following second type of boundary condition is analyzed:

(43)

leading to

(44)

As before, the physical realization of this boundary condition is a -like interaction, with the only difference residing in the type of the basic interaction. Unlike in the previous case, here the interaction is of electrostatic character. Thus the situation corresponds to a -like negative charge distribution placed at the origin of the wall. Since the two ion types respond oppositely, the ion density profiles differ only in sign:

(45)

This implies that the total ion density deviations vanish . Accordingly, also the density deviation for the solvent vanishes, Figure 4 shows the density profiles of the positive ions, which, up to the sign, are the same as the ones for the negative ions. As stated above, for this boundary condition, there is no need to discuss the behavior of the solvent particles.

Figure 4: Density profiles of the ions as functions of the normal distance from the wall (a), of the absolute value of the lateral wave vector (b), and of the wave number in normal direction (c). Note that due to , in panel (c) the absolute value is shown. Each panel shows the profiles for three values of the other relevant variable. These profiles are cuts of the corresponding data (analogous to Fig. 2) along various directions. Here, the boundary condition is given by (see Eqs. (15) and (44)), which corresponds to a -like surface charge at the origin in real space (see Eq. (43)). The profiles for the solvent are not shown, because the deviations linked to the two types of ions cancel out, , leaving the density of the solvent unchanged as if there were no ions. In comparison with Fig. 3, the profiles in (a) decay much slower on the scale of the Debye length (see Eq. (34)) instead of on the scale of the much shorter bulk correlation length (see Fig. 3(b) and Eq. (36)) due to the nonelectrostatic interaction. Accordingly, the profiles in (b) and (c) decay on the scale of more rapidly than their counterparts in Figs. 3(d) and (f). For the remaining relevant parameters see Sec. III.1.

The three panels in Fig. 4 are obtained similarly as the ones in Fig. 3. Figure 4 (a) shows the density profiles