Electrolyte solutions at heterogeneously charged substrates

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 \delta-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 1,2 and wetting phenomena 3,4 via coating 5 and surface patterning 6,7 to colloid science 8,9 and microfluidics. 10,11 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, 12 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. [13][14][15] 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. [16][17][18][19][20][21][22][23][24][25][26] 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. [27][28][29] 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 ref. 16, 20, 22 and 26, 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. [30][31][32] 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 Section II, selected cases of heterogeneous walls are discussed in Section 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 16,[18][19][20][21] 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 Section IV.

A. Setup
Here, the influence of a chemically and electrically nonuniform wall on the fluid density is studied. In spatial dimension d = 3 the system consists of an impenetrable planar wall for z o 0 and a fluid for z 4 0, both parts being macroscopically large. In the following, the space occupied by the fluid is denoted by V :¼ fr ¼ ðx; y; zÞ 2 A Â Lg; the positions r ¼ ðx; y; zÞ ¼ ðr k ; zÞ are uniquely decomposed into the lateral components r k ¼ ðx; yÞ 2 A & R 2 and the normal component z 2 L ¼ ½0; L relative to the wall surface at z = 0. The size jAj of the wall and the extent L 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 (z = 0) 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.

B. Density functional theory
We use density functional theory [30][31][32] in order to determine the equilibrium number density profiles . ¼ ðR 1 ; R 2 ; R 3 Þ 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 33 is considered: where b = 1/(k B T) is the inverse thermal energy, l ¼ ðm 1 ; m 2 ; m 3 Þ are the chemical potentials of the three species, b 4 0 is a phenomenological parameter with dimension [b] = (length 5 ), which can be inferred from microscopic models (see Section IIIA), e 0 E 8.854 Â 10 À12 A s V À1 m À1 is the vacuum permittivity, 37 e is the relative dielectric constant of the fluid, Cðr; ½.Þ is the electrostatic potential at r 2 V, and h(r 8 ) = (h 1 (r 8 ), h 2 (r 8 ), h 3 (r 8 )) describes the strengths of the nonelectrostatic wall-fluid interactions at r ¼ ðr k ; 0Þ for the three species. Note that for the sake of simplicity, the coupling of number density gradients of different particle types is neglected in eqn (1) (see Section IIIA).
In the present study the bulk state . b ¼ ðR 1;b ; R 2;b ; R 3;b Þ is considered to be thermodynamically far away from any phase transition so that the local contribution bo(.) of the density functional in eqn (1) can be safely expanded around . b up to quadratic order in d.
where the local part of the interactions between different types of particles is captured by the real-valued, symmetric, and positively definite 3 Â 3-matrix M (see, cf., eqn (30)). Furthermore, oð. b ; lÞ ¼ Àp specifies the grand potential density, evaluated for the equilibrium bulk densities . b , which equals minus the bulk pressure p; in the following its value is of no importance. For a given equation of state pð. b ; TÞ the bulk densities Þ are free parameters of the model. Finally, the electrostatic potential Cðr; ½.Þ, which enters into eqn (1) on a mean-field level via the electric field energy density, fulfills the Poisson equation for r 2 V with the boundary conditions for r k 2 A, where s(r 8 ) is the surface charge density at the point r ¼ ðr k ; 0Þ on the wall surface (z = 0), and Z ¼ ðZ 1 ; Z 2 ; Z 3 Þ denotes the valences of the fluid species. The Euler-Lagrange equations, corresponding to the minimum of the density functional specified in eqn (1)-(4), can be written as and À 1 4pl B r 2 beCðrÞ ¼ Z Á d.ðrÞ (6) for r 2 V with the boundary conditions given by eqn (4) and by for r k 2 A, where l B = be 2 /(4pe 0 e) is the Bjerrum length. The linear nature of the Euler-Lagrange eqn (5) and (6) tells that the quadratic (Gaussian) approximation of the underlying density functional in eqn (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 8,34 ) that for sufficiently weak wall-fluid interactions linear response theory provides quantitatively reliable results.

C. Solution of the Euler-Lagrange equations
Instead of solving the Euler-Lagrange equations in eqn (5) and (6) as differential equations for the profiles d. and C as functions of r ¼ ðr k ; zÞ, it is convenient first to perform Fourier transformations with respect to the lateral coordinates r 8 . The resulting transformed profiles andĈ as functions of q k ¼ ðq x ; q y ÞAR 2 and z A R can be combined in the four-component quantity v q k ; z ¼ d. q k ; z ; beĈ q k ; z so that eqn (5) and (6) can be written as (9) where k :  2 . This does not allow for the formation of a scalar product of two vectors of the type v ¼ ðd.; beĈÞ; however, in the following scalar products will not occur.
T is a diagonal matrix with these entries), one obtains Tv 00 q k ; z ¼ T À1 NðkÞT À1 |fflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflffl ffl} ¼: HðkÞ Tv q k ; z : The 4 Â 4-matrix HðkÞ is independent of z and it is symmetric but not real-valued, because the bottom entry of T is imaginary.
HðkÞ is not a normal matrix, i.e., it does not commute with its adjoint matrix HðkÞ y , and hence it does not possess an orthogonal basis composed of eigenvectors. However, the actual structures of the matrix M and of the vector Z used below guarantee the existence of a nonorthogonal basis fK 1 ðkÞ; . . . ; K 4 ðkÞg of eigenvectors of the matrix HðkÞ with respective positive (real-valued) eigenvalues l 1 ðkÞ; . . . ; l 4 ðkÞAð0; NÞ (see Appendix A). Expanding the vector Tv q k ; z in this basis fK 1 ðkÞ; . . . ; K 4 ðkÞg, leads to eqn (10) in the form with the solution where the second boundary conditions in eqn (4) and (7) have been used. Therefore, the solution of eqn (10) can be expressed as Finally, the first boundary conditions in eqn (4) and (7) can be expressed as From eqn (14) and (15) the coefficients g 1 ðq k Þ; . . . ; g 4 ðq k Þ can be determined. Note that according to eqn (14) and (15) the coefficients g 1 ðq k Þ; . . . ; g 4 ðq k Þ and hence the profiles. andĈ depend linearly on the nonelectrostatic wall-fluid interactions ĥ(q 8 ) and the surface charge densityŝ(q 8 ). 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 eqn (1)-(4).

A. Choice of parameters
The present study discusses the influence of the wall-fluid interactions, represented by the nonelectrostatic wall-fluid interactions ĥ(q 8 ) and the surface charge densityŝ(q 8 ), onto the number density profiles . in the adjacent fluid. Applying density functional theory as described in Section II requires knowledge of the bulk number densities . b , the parameter b, and the coupling matrix M all of which are bulk quantities or characterize them. In the bulk local charge neutrality holds, i.e., ZÁ. b = 0. Hence the equilibrium bulk state is determined by the temperature T, the number density R 1,b of the solvent, and the bulk ionic In order to obtain expressions for the parameter b and for the coupling matrix M 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 r is modeled by a square-well pair potential U(r) as displayed in Fig. 1(a). At small distances r o s a hard core repulsion prevents the overlap of two particles. At intermediate distances rAðs; s pot Þ two particles attract each other with an interaction energy Àe pot , and at distances r 4 s pot the nonelectrostatic interaction vanishes. According to the scheme due to Barker and Henderson, 35 the interaction potential U can be decomposed as U = U hc + U pot into the hard core repulsion U hc and the attractive well U pot (see Fig. 1 for the pure solvent (species 1) in the bulk can be approximated by the expression The contribution O mic 1 (here within local density approximation (LDA)) is due to the reference system governed solely by the hard core interaction U hc : Here, contributions of external potentials are neglected, because they do not contribute to the bulk parameters b and M.
The second term on the rhs of eqn (17) is (within random phase approximation (RPA) 30 ) the excess free energy functional due to the square-well attractive interaction U pot : In eqn (18) L 1 is the thermal de Broglie wave length, m 1 denotes the chemical potential of species 1, and f ex,hc (R 1 ) is the excess free energy per volume of the reference system governed by the hard core interaction U hc . Following Cahn and Hilliard, 33 eqn (19) can be approximated by a gradient expansion: with the m-th moment of the pair potential U pot in units of k B T, which, for the present form of U pot , leads to From eqn (20) one obtains the gradient expansion of bO mic 1 [R 1 ] in eqn (17): with the local contribution The comparison of eqn (23) with eqn (1) renders an expression for the parameter b in terms of parameters of the interaction Fig. 1 The nonelectrostatic interaction between fluid particles is modeled by a square-well pair potential U(r) displayed in panel (a), where r denotes the distance between the centers of two spherical particles. For r o s a hard core repulsion prevents the overlap of two particles. For r A (s, s pot ) two particles attract each other with a constant interaction energy Àe pot o 0. At distances r 4 s pot there is no nonelectrostatic interaction. Panel (b) sketches the decomposition U = U hc + U pot of the nonelectrostatic interaction potential U according to the scheme due to Barker and Henderson into the hard core repulsion U hc and the attractive well U pot , which is used in Section IIIA in order to obtain the parameters entering the Cahn-Hilliard square-gradient density functional in eqn (1).
potential U (see Fig. 1 By expanding bo loc 1 ðR 1 ; m 1 Þup to quadratic order in the density deviation dR 1 = R 1 À R 1,b from the equilibrium bulk density R 1,b , which is a solution of the Euler-Lagrange equation one obtains The comparison with eqn (2) leads to the matrix element of the matrix M, where the first term on the rhs stems from the ideal gas contribution of the solvent particles. The argument R 1,b of the second term, which is due to the hard core interaction U hc , is a measure of the packing fraction Z = pR 1,b s 3 /6. The analogue of eqn (23) for the nonelectrostatic contribution of all three particle species is given by the first line of eqn (1) with the local contribution (compare eqn (24)) where R tot = R 1 + R 2 + R 3 denotes the total number density. Note that eqn (29) assumes, that all interactions among the species are the same (see eqn (21)). This implies that the last term in eqn (29) takes the form K 0 2 ðR tot Þ 2 . By expanding bo loc ð.; lÞ up to quadratic order in the density deviations d. = . À . b from the equilibrium bulk densities . b one finally finds (see the steps leading to eqn (28)) where R tot b = R 1,b + R 2,b + R 3,b = R 1,b + 2I denotes the total number density in the bulk. In the present study the hard core excess free energy per volume f ex,hc (R b ) is chosen as the one corresponding to the Carnahan-Starling equation of state: 35 here with the packing fraction Z = pR tot b s 3 /6. Accordingly, from eqn (24) one obtains the following equation of state of the pure solvent: Its derivative with respect to the number density R 1,b , using the relation qp/qR As an example we consider water at room temperature T = 300 K and ambient pressure p = 1 bar (which corresponds to the number density R b,1 = 55.5 M E 33.3 nm À3 and the isothermal compressibility k T = 4.5 Â 10 À10 Pa À1 37 ) with relative dielectric constant e = 80, i.e., with Bjerrum length l B = 0.7 nm, and with a univalent salt of ionic strength I = 1 mM E 6 Â 10 À4 nm À3 . The strength of hydrogen bonds, which generate the dominant attractive interaction contribution, is of the order of e pot E 20 kJ mol À1 E 8k B T. 36,37 Using these data, one obtains from eqn (32) and (33) the bulk packing fraction Z E 0.44 as well as s = 2.9 Å and s pot = 3.4 Å. In the following the Debye length is used as length scale, which, for the present choice of parameters, is 1/k E 10 nm.
In the case of a pure solvent (dR 2 = dR 3 = C = 0), in the bulk the density two-point correlation function Gðr 1 ; r 2 Þ ¼ Gðr 1 À r 2 Þ fulfills an equation similar to eqn (5): Note that the similarity between eqn (5) and (35) is due to the asymptotic proportionality between density deviations and twopoint correlation functions (Yvon equation). 35 From eqn (35), one can readily infer the relation for the solvent bulk correlation length, which characterizes the exponential decay of % G(r). For the present choice of parameters, one has x E 1.3 Å so that kx E 0.013.

B. X-ray scattering
In the following subsections the Fourier transforms d. ¼ Þof the profiles of the density deviations as functions of the lateral wave vector q 8 and of the normal distance z from the wall are discussed in detail. However, from the experimental point of view, it is challenging to directly obtain the z-dependence of the density profiles. One of such direct methods is total internal reflection microscopy (TIRM) 38 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. 39,40 As X-rays are predominantly scattered by the electrons of the fluid molecules one has to consider the electron number density for r k 2 A with the number N j of electrons per molecule of particle species j A {1, 2, 3}, the bulk electron density N j dR j of the electron number density from its bulk value. The X-ray scattering signal for scattering vector q ¼ À the electron number density profile R e in both the lateral and the normal direction. 35 For common specular X-ray reflectivity measurements, i.e., for q 8 = 0, the normalized intensity reflected as function of the normal wave number q z is given by 39,40 where R F (q z ) denotes the Fresnel reflectivity of an ideal, step-like planar interface, 41 and where the notation d b b has been used. Moreover, off-specular diffuse X-ray scattering (q 8 a 0) at grazing incidence (GIXD, Im q z a 0) yields scattering intensities which are proportional to d b b R e q k ; q z 2 . 39 Hence, as the double Fourier transforms d b b R j in eqn (39) of the density deviation profiles dR j are of direct experimental relevance, they will be discussed in the following in parallel to the single Fourier transforms dR j . Note that due to d b b R i 2 C, in Fig. 3 and 4 its absolute value is shown.

C. 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., (1, 0, 0, 0), (0, 1, 0, 0), (0, 0, 1, 0), and (0, 0, 0, 1), which span the four-dimensional space of boundary conditions v 0 À q k ; 0 Á in eqn (15). Therefore, as a first step to study the influence of wall inhomogeneities onto the fluid, these four distinct boundary condition vectors v 0 À q k ; 0 Á are studied. The first one of these vectors is given by which requires (see eqn (15) and which in real space corresponds to the boundary condition This boundary condition corresponds to an attractive, d-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 d b R 1 q k ; z for the solvent and dR 2 q k ; z ¼ dR 3 q k ; z for the ions, as shown in Fig. 2(a) and (b), respectively. Since the boundary condition corresponds to a constant in Fourier space, the density distribution v À q k ; z Á depends on k = |q 8 | only. Actually, the solution v À q k ; z Á for this system is proportional to the first column of the Green's function, which is a 4 Â 4-matrix, of the differential operator corresponding to eqn (5) and (6). Fig. 2 illustrates that for fixed |q 8 | the density deviations from the bulk value increase for smaller normal distances from the wall and, for fixed z, also upon decreasing the absolute value of the lateral wave vector |q 8 |. 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 |q 8 | 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. Fig. 3(a), (c) and (e) show the density profiles for the solvent and Fig. 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. Fig. 3(a) and (b) show the density deviations as function of the normal distance z from the wall for three values of |q 8 |, i.e., Fig. 3(a) and (b) correspond to horizontal cuts through Fig. 2(a) and (b), respectively. For fixed |q 8 |, as in Fig. 2(a) and (b), Fig. 3(a) and (b) clearly show an exponential decay of the density deviation for increasing distances from the wall. In contrast, Fig. 3 Fig. 2(a) and (b), i.e., density profiles as functions of the absolute value of the lateral wave number |q 8 | for three normal distances z from the wall. The dependence of these profiles on the absolute value |q 8 | of the lateral wave vector q 8 implies a laterally isotropic decay of the density deviations in real space. The third pair of graphs, Fig. 3(e) and (f), shows the Fourier transforms of the density profiles of Fig. 3(a) and (b), being additionally Fourier-transformed with respect to the normal direction z, which leads to the Fourier

(c) and (d) show vertical cuts through
in terms of the lateral wave vector q 8 and the normal wave number q z , respectively. All curves in Fig. 3(c)-(f) exhibit a Lorentzian shape as functions of |q 8 | and q z , respectively. These Lorentzian curves in Fourier space correspond to exponential decays in real space in lateral or normal direction. The curves in Fig. 3(c) and (d) show widths of half height which decrease with increasing normal distance z, 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 Fig. 3(e) and (f) exhibit widths of half height which increase with the lateral wave number |q 8 |, 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 Fig. 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 v 0 ¼ À h After having discussed the effects of the boundary condition ĥ a 0 via Fig. 2 and 3, the following second type of boundary condition is analyzed: leading to v 0 q k ; 0 ¼ À bes ð0Þ e 0 e ð0; 0; 0; 1Þ: As before, the physical realization of this boundary condition is a d-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 d-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: This implies that the total ion density deviations vanish dR 2 þ dR 3 ¼ 0. Accordingly, also the density deviation for the solvent vanishes, dR 1 ¼ 0. Fig. 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.  (34)). The plane z = 0 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 v 0 ¼ À h The three panels in Fig. 4 are obtained similarly as the ones in Fig. 3. Fig. 4(a) shows the density profiles d b R 2 q k ; z as functions of the normal distance z from the wall for three values of the lateral wave number |q 8 |. Fig. 4(b) shows the same density profiles d b R 2 q k ; z but as functions of |q 8 | for three distances z from the wall. Fig. 4(c) displays the double Fourier transform d b b R 2 q k ; q z . Compared with the profiles in Fig. 3 for the previously discussed boundary conditions, all present profiles differ significantly from them. Fig. 4(a) reveals a much larger decay length in z-direction, i.e., normal to the wall. Also in Fourier space the decay in lateral direction occurs much more rapidly, i.e., on much longer length scales in real space than in the case of the nonelectrostatic wall-fluid interaction (cf. Fig. 3). This is also indicated by the much narrower peak in the double Fourier transform (see Fig. 4(c)). Furthermore,  (15) and (40)), corresponding to a d-like nonelectrostatic attraction of the solvent particles at the origin of the wall (see Fig. 2 and eqn (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 Section IIIA. Fig. 4(a) shows a variation of the decay length in normal direction as function of |q 8 |. In Fig. 4(b) one observes that the lateral wave numbers |q 8 | at which the profiles d b decay to half of the maximum values decrease upon increasing the distances z from the wall, from which one infers that the lateral decay length in real space increases upon increasing z.
The decay with respect to |q 8 | is much faster than in the previous case (compare Fig. 3(d)), indicating that in real space there is a slower decay in the lateral direction. Moreover, in Fig. 4(b) and (c) the functional form differs from the one shown in Fig. 2 and 3. These differences naturally occur due to the different form of the boundary condition. Since in the case of the boundary condition studied above (see Fig. 2 and 3) the relevant interaction is nonelectrostatic, the length scale dominating the decay is given by the corresponding short-ranged bulk correlation length x (see eqn (36)). In contrast, for the system shown in Fig. 4, due to the electrostatic nature of the corresponding interaction, the dominating length scale is the Debye length 1/k (see eqn (34)), which is much larger than the correlation length x due to the nonelectrostatic interaction, giving rise to the much slower decay in Fig. 4(a) (on the scale of 1/k) and the much faster decay in Fig. 4(b) and (c) (on the scale of k).

D. Circular patch of interaction
Having discussed in Section IIIC actually point-like interactions between the wall and the fluid, as the next step we now study the influence of interaction patterns on the density deviations close to a wall upon broadening the spatial extent of the interaction area. To this end we analyze the influence of a two-dimensional circular interaction patch of radius R centered at the origin (see Fig. 5(a)). Due to the radial symmetry, the spatial structures in Fourier space depend only on the absolute value |q 8 | of the lateral wave vector q 8 . Fig. 6 discusses four distinct configurations.
Alluding to the insights gained in the previous section, Fig. 6(a) and (c) correspond to a homogeneous circular patch of radius R, which interacts with the solvent only, similar to Fig. 2 and 3. This amounts to the boundary condition (see eqn (15)) where the two-dimensional Fourier transform of the Heaviside function Y(R À |r|) is given by In contrast, Fig. 6(b) and (d) refer to a charged circular patch of radius R, similar to Fig. 4: (44)), which corresponds to a d-like surface charge at the origin in real space (see eqn (43)). The profiles for the solvent are not shown, because the deviations linked to the two types of ions cancel out, dR 2 þ dR 3 ¼ 0, 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 1/k (see eqn (34)) instead of on the scale of the much shorter bulk correlation length x (see Fig. 3(b) and eqn (36)) due to the nonelectrostatic interaction. Accordingly, the profiles in (b) and (c) decay on the scale of k more rapidly than their counterparts in Fig. 3(d) and (f). For the remaining relevant parameters see Section IIIA.  (53)). Our results are based on the choice D = 5D peaks (see, cf., Fig. 7). Fig. 6 Density profiles of the solvent (panels (a) and (c)) and of the ions (panels (b) and (d)) as functions of the absolute value |q 8 | of the lateral wave vector for three normal distances z from the wall. The boundary condition corresponds to a circular interaction patch centered at the origin with radius R 4 0. In the panels (a) and (b) the radius of the patch is R = 0.5/k whereas in the panels (c) and (d) the radius is R = 2/k, where 1/k is the Debye length (see eqn (34)). All considered patch radii are much larger than the bulk correlation length, R c x (see eqn (36)). In addition, there are different types of interaction. In panels (a) and (c) the interaction between the wall and the solvent particles is nonelectrostatic (h 1 (r 8 ) = h (0) 1 Y(R À |r 8 |), h 2 (r 8 ) = h 3 (r 8 ) = s(r 8 ) = 0, see eqn (46) as well as Fig. 2 and 3), whereas in panels (b) and (d) the patch contains a constant surface charge and therefore interacts with the ions only (h(r 8 ) = 0, s(r 8 ) = s (0) Y(R À |r 8 |), see eqn (49) and Fig. 4). Besides the profiles, all panels show also the lateral Fourier transform of the boundary condition (eqn (47) and (50)) displayed as a black solid line. In the case of the interaction of the wall with the solvent ((a) and (c)), the decay of the profiles as function of |q 8 | is proportional to the lateral Fourier transform of the boundary condition, which implies, that the density deviations in real space closely follow the shape of the patch. However, in the case of a charged patch at the surface the density distribution of the ions reflects the competition between the length scale R of the radius of the patch and the Debye length 1/k. In the case of small patches (R o 1/k, panel (b)), the Debye length dominates and therefore dictates the decay as function of |q 8 | without noticeable influence of the patch. In contrast, in the case of large patches (R 4 1/k, panel (d)), in which the radius of the patch is the dominating length scale, the shape of the profiles follows the Fourier transform of the charge distribution at the wall, i.e., the patch of radius R. For the remaining relevant parameters see Section IIIA.
leading to the boundary condition v 0 q k ; 0 ¼ À2pR 2 be s ð0Þ e 0 e J 1 q k R q k R ð0; 0; 0; 1Þ: (50) Fig. 6(a) and (b) correspond to the patch size R = 0.5/k whereas Fig. 6(c) and (d) correspond to R = 2/k. In all four panels the black line is given by AJ 1 (|q 8 |R)/(|q 8 |R) with A chosen such that the first maximum of the data for zk = 0 is reproduced. In Fig. 6(a) and (c) only the solvent density profiles are shown, because, due to linearity, Section IIIC indicates, that the ion profiles are proportional to the one of the solvent. Fig. 6(a) and (c) clearly show, that the density deviations are proportional to the Fourier transform v 0 À q k ; 0 Á of the boundary condition (eqn (47), solid black line). This implies that for increasing lateral distances from the center of the patch the decay of the profiles in real space is dominated by the length scale set by the radius R of the patch. This trend holds for both patch sizes. However, as expected, the amplitudes of the density deviations increase for the larger patch size (note the different scales). In contrast, in Fig. 6(b) and (d), where the density profiles of the positive ions are shown and where the profiles for the solvent are omitted for the same reasons as explained in Section IIIC, the profiles do not follow the Fourier transform of the boundary conditions (solid black line). This is particularly pronounced in Fig. 6(b), i.e., for the smaller patch size. In this case the decay as function of |q 8 | is faster than the Fourier transform of the boundary condition, which implies that the profiles decay on a length scale larger than that of the radius R of the patch and also the shape of the decay differs from that of the expression J 1 (|q 8 |R)/(|q 8 |R). This behavior can be understood in terms of two distinct dominating length scales. In Fig. 6(b) and (d), where the effect of electrostatic interactions are shown, the dominating length scale is the Debye length 1/k in contrast to the much smaller correlation length x (x E 1.3 Â 10 À2 k À1 ) induced by the nonelectrostatic interactions characterizing Fig. 6(a) and (c). Since in Fig. 6(b) the radius R of the patch is only half the Debye length 1/k, the dominating length scale is the Debye length 1/k, so that the density deviations decay in real space on a length scale which is larger than the patch radius R. Also, since the profile in Fourier space is not proportional to the Fourier transform of the boundary condition, one can conclude that the shape of the patch has no significant influence on the decay behavior. The competition of the length scales x, 1/k, and R is also borne out in Fig. 6(d), where the patch size R is twice as large as the Debye length 1/k. This case is much more similar to the ones in Fig. 6(a) and (c), because the dominating length scale is set by the radius R, and consequently the profiles follow rather closely the shape (solid black line) dictated by the interaction patch. However, the influence of the smaller Debye length scale is still visible, which is the reason for the deviations from the Fourier transform of the boundary condition (solid black line). In conclusion, as already seen in Section IIIC, the largest length scale sets the decay behavior of the density deviations. In the present case of nonvanishing sizes of the interaction areas, the largest length scale dictates not only the range but also the shape of the density deviations.

E. Periodic distribution of interaction sites
After having discussed the density profiles in the presence of spatially localized, single interaction sites in Sections IIIC and D, here we study the influence of interaction sites forming a regular hexagonal lattice: see Fig. 5(b); the distance between nearest neighbor sites is denoted as D.
Distinct from the previous examples in Sections IIIC and D, the interaction strength around the individual interaction sites r peaks is taken to form a Gaussian distribution, providing either a nonelectrostatic or an electrostatic interaction with equal amplitudes for all interaction sites: and respectively, where D 2 peaks is the variance of the Gaussian interaction. Lateral Fourier transformation leads to the corresponding boundary condition v 0 (see eqn (15) and v 4 0 q k ¼ À be s ð0Þ e 0 e 2pD 2 peaks exp À q k 2 D 2 where jC G j ¼ ð16p 2 =ð3D 2 ÞÞsinð60 Þ is the size of an elementary cell of the corresponding two-dimensional reciprocal lattice G.
Using this boundary condition, we have studied four different systems, as shown in Fig. 7.
The four panels are arranged as in Fig. 6, with the boundary conditions corresponding to an interaction between the wall and solvent only in Fig. 7(a) and (c), and in Fig. 7(b) and (d) boundary conditions corresponding to an interaction between the wall and ions only, i.e., due to a hexagonal lattice of interaction sites with Gaussian intrinsic charge distribution: For the same reason as stated in the context of Fig. 6, in the former case ( Fig. 7(a) and (c)) only the deviations of the solvent density and in the latter case ( Fig. 7(b) and (d)) only the deviations of the ion densities are shown. Fig. 7(a) and (b) correspond to the lattice constant D = 0.5/k whereas Fig. 7(c) and (d) correspond to D = 2/k. The variance D 2 peaks of the peaks is taken as D peaks = D/5, so that D peaks = 0.1/k in Fig. 7(a) and (b) and D peaks = 0.4/k in Fig. 7(c) and (d). Fig. 7(a) and (c) tell that, although different values for q 8 change the amplitude of the profiles in all cases, the solvent profile decays exponentially, upon increasing the normal distance z from the wall, on the scale of the bulk correlation length x. This holds for both values of the lattice constant D. However, the amplitude of the density deviations is slightly increased for the larger lattice constant D, which is in line with the also increased variance D 2 peaks of the interaction sites. In contrast to these findings, Fig. 7(b) and (d) reveal a different picture. In these panels, one still finds an exponential decay of the profiles upon increasing the normal distance z. However, the profiles decay on a much larger length scale than the ones in Fig. 7(a) and (c). Moreover, not only the Fig. 7 Density profiles of the solvent (panels (a) and (c)) and of the ions (panels (b) and (d)) for three lateral wave vectors q 8 = (q x ,q y ) as functions of the normal distance z from the wall in units of the Debye length 1/k. The boundary condition corresponds to a hexagonal lattice of interaction sites with a Gaussian charge distribution characterized by a standard deviation D peaks = D/5. The lattice constant is denoted as D (see Fig. 5). In panels (a) and (b) the lattice constant and the variance are D = 0.5/k and D peaks = 0.1/k, respectively, whereas in panels (c) and (d) the lattice constant and the variance are D = 2/k and D peaks = 0.4/k, respectively, with the Debye length 1/k (see eqn (34)). Panels (a) and (c) correspond to systems with a nonelectrostatic interaction between the wall and the solvent particles (see eqn (56)), whereas panels (b) and (d) correspond to systems with electrostatic interaction sites between wall and ions (see eqn (57)). The insets in (b) and (d) show a magnified version of the respective profiles in the main plot. In all cases, the profiles decay exponentially upon increasing the normal distance z from the wall. However, the decay length differs significantly between the two aforementioned types of interactions. In the case of the nonelectrostatic interaction, the decay length is set by the bulk correlation length x (see eqn (36)) of the fluid, whereas in the case of the electrostatic interaction it is set by the much larger Debye length 1/k. This difference in the decay lengths, both in lateral and in normal direction, which leads to a much faster lateral decay in the case of the electrostatic interactions, is also responsible for the decreasing amplitude of the ion profiles for increased wave vectors (panel (b) and (d)). Another significant difference between the two interaction types is the variation of the decay length as function of the lateral wave vectors. In panels (a) and (c) all profiles decay exponentially on the same decay length x, whereas in panels (b) and (d) the decay length depends significantly on the wave vectors. This effect follows from the dependence of the eigenvalues on |q 8 | as discussed in eqn (A3), corresponding to a lateral decay proportional to exp À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k 2 þ jq k j 2 q z . In principle this occurs for both types of interactions.
However only in the case of the electrostatic interactions it is relevant, which again is due to the difference between the dominating length scales. For the remaining relevant parameters see Section IIIA.
amplitude but also the decay length changes significantly for different values of q 8 . This was already encountered in Fig. 4, where the decay length depends on the value of |q 8 |. This variation of the decay lengths can be inferred from eqn (14), which shows that the eigenvalues and thus the decay length depends on k = |q 8 |. The variation of the decay length can be expressed in terms of the Debye length 1/k, which determines the length scale in case of |q 8 | = 0. From eqn (A3) one finds, that the decay as function of z is proportional to exp À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi k 2 þ jq k j 2 q z . The large differences in the amplitudes of the various profiles in Fig. 7(b) and (d), as well as the pronounced increase of the decay length in comparison to Fig. 7(a) and (c) can be understood in terms of the differences between the dominating length scale. Analogous to the previous sections, for the systems shown in Fig. 7(a) and (c), the dominating length scale in lateral direction is the length scale set by the boundary conditions and the bulk correlation length x, which characterizes the decay of the solvent density in normal direction. However, for the systems shown in Fig. 7(b) and (d), the relevant inherent length scale of the fluid is the Debye length 1/k, which is significantly larger than the bulk correlation length x and thus causes the increase in the length scale of the decay, both in lateral and in normal direction.

IV. Conclusions and summary
In the present analysis the influence of a chemically or electrostatically structured surface on an adjacent fluid has been studied and described in terms of the density profiles of the fluid components. The fluid, which comprises a single solvent species and a single univalent salt far away from bulk and wetting phase transitions, has been investigated within classical density functional theory. [30][31][32] Within this model three examples of heterogeneous walls have been studied. First, single isolated interaction sites are discussed, which interact either nonelectrostatically (between the wall and solvent particles) or electrostatically (between the wall and ions) (see Sections IIIC and D). In the case of a d-like nonelectrostatic interaction, the solvent density increases around the interaction site and decays exponentially on the length scale of the bulk correlation length x. The deviations of the ion number densities from their bulk values are proportional to that of the number density of the solvent (see Fig. 2 and 3). For a d-like electrostatic interaction, within the present model, the solvent does not respond at all, because the deviations induced by the two ion types even out due to symmetries, whereas the density deviations of the ionic particles again decay exponentially. However, the length scale of the latter decay is significantly increased as compared to the former case, because the dominating scale in this case is the Debye length 1/k c x (see Fig. 4). The introduction of another length scale by studying interaction sites of non-vanishing extent (see Section IIID) shows, that the resulting density profiles strongly depend on the dominant length scale (see Fig. 6). If a bulk length scale (bulk correlation length x or Debye length 1/k) dominates, the profiles resemble the ones for d-like interactions. However, if a length scale set by a boundary condition at the wall dominates or is similar to the dominating length scale in the bulk, the decay of the density deviations increasingly reflects the boundary conditions. Finally, the examination of multiple interaction sites, arranged as a regular hexagonal lattice (see Fig. 5(b)), shows, that the size of the interaction sites and the distance between them influence the amplitude and thus the importance of density deviations for large values of the lateral wave number |q 8 | (see Fig. 7).
In summary, the present study provides a flexible framework to determine the influence of various surface inhomogeneities on the density profiles of a fluid in contact with that substrate. The resulting profiles are found to be sensitive to the type of interaction as well as to the size and the distribution of the interaction sites.
This framework is considered as a starting point for extensions into various directions, aiming for the analysis of more sophisticated and realistic models. First, the model used here to describe the fluid is a very simple one, chosen to lay a foundation for further research and to introduce the approach as such. Concerning future work, more realistic descriptions of the fluid and more elaborate density functional descriptions could be used. For instance, the present restriction to low ionic strengths and equal particle sizes can be removed along the lines of ref. 42. Second, for the systems studied here, the fluid is thermodynamically far from any bulk or wetting phase transitions. This is solely done for the sake of simplicity. In future studies of more realistic systems, taking into account the occurrence of phase transitions and their influence on the systems is expected to be rewarding. Third, this study is restricted to linear response theory. Whereas this allows for a broad overview of structure formation in terms of superpositions of only a few elementary patterns, the occurrence of nonlinear structure formation phenomena requires approaches beyond linear response theory. Finally, studying the influence of disordered surface structures within the present framework appears to be very promising.

Conflicts of interest
There are no conflicts to declare.

Appendix A: eigenvectors and eigenvalues of HðkÞ
According to the structure of the matrix M (eqn (2)), with entries given by eqn (30), and of the vector Z ¼ ð0; 1; À1Þ, from eqn (9) and (10) one infers that the matrix HðkÞ has the form HðkÞ ¼