Structure of electrolyte solutions at non-uniformly charged surfaces on a variety of length scales

The structures of dilute electrolyte solutions close to non-uniformly charged planar substrates are systematically studied within the entire spectrum of microscopic to macroscopic length scales by means of a unified classical density functional theory (DFT) approach. This is in contrast to previous investigations, which are applicable either to short or to long length scales. It turns out that interactions with microscopic ranges, e.g., due to the hard cores of the fluid molecules and ions, have negligible influence on the formation of non-uniform lateral structures of the electrolyte solutions. This partly justifies the Debye-H\"uckel approximation schemes applied in previous studies of that system. In general, a coupling between the lateral and the normal fluid structures leads to the phenomenology that, upon increasing the distance from the substrate, less details of the lateral non-uniformities contribute to the fluid structure, such that ultimately only large-scale surface features remain relevant. It can be expected that this picture also applies to other fluids characterized by several length scales.


I. INTRODUCTION
Historically, the theoretical study of solid-fluid interfaces has naturally started with the investigation of idealized surfaces with laterally uniform properties [1][2][3][4][5][6] instead of realistic models of surfaces with geometrical, chemical, or electrical non-uniformities. This approach was justified, on the one hand, by the initial lack of knowledge about the microscopic structure of real surfaces, and, on the other hand, by the computational advantages gained from exploiting lateral symmetries. However, in particular in the context of electrochemistry and colloidal science, efforts have been made to include surface non-uniformities into the theoretical description. A pioneering contribution is due to Richmond [7,8], who studied the effective interaction of two parallel planar dielectric bodies with non-uniform surface charge distributions mediated by a dilute electrolyte solution in between, assuming that the linearized Poisson-Boltzmann (Debye-Hückel) approximation [9] (see also Refs. [10][11][12]) is applicable. In recent years the issue of electrolyte solutions close to non-uniformly charged substrates within the Debye-Hückel approximation [13][14][15][16], (non-linearized) Poisson-Boltzmann theory [17,18], as well as statistical field theory [19][20][21][22][23] has been addressed intensively (see also the review in Ref. [24]). These studies are focused on large length scales, either by ignoring the microscopic fluid structure of the electrolyte solution or by modelling its long-ranged structure within a square gradient approximation (see Ref. [15]). Moreover, * markus.bier@fhws.de microscopic approaches, e.g., Monte Carlo (MC) simulations [25,26] or classical density functional theory (DFT) [27,28], have been used. But, due to technical reasons, theses studies were limited to rather small systems and special types of surface charge non-uniformities. Thus an approach is missing which exhibits the accuracy of a DFT combined with the efficiency of a Debye-Hückel approximation, in order to span the whole range from microscopic to macroscopic length scales.
The present study suggests a step in this direction. This novel method consists of a quadratic expansion of the density functional not about the bulk profiles (as within the Debye-Hückel approximation), but about the profiles of a planar-symmetric (i.e., quasi onedimensional) system. The surprising observation is, that microscopic hard-core contributions turn out to be quantitatively irrelevant for the formation of the lateral structure. In this respect, disregarding the size of the fluid molecules and ions by using the Debye-Hückel approximation for laterally non-uniform modes, as done in many previous studies, is justified. However, the present investigation suggests, that, in contrast to those previous studies, the Debye-Hückel approximation should not be used for the planar-symmetric contributions, which require more sophisticated descriptions including, e.g., finite size effects.
Our contribution is structured as follows: Section II describes the considered model of an electrolyte solution in contact with a non-uniformly charged substrate and the formalism to infer the structural quantities. Results concerning the number density profiles in the normal and in the lateral directions as well as concerning the interfacial tension as function of the length scales of the lateral non-uniformities are presented and discussed in Sec. III. z x y FIG. 1. A planar, non-uniformly charged substrate with negatively charged (bright, yellow) and charge-neutral (dark, brown) regions is in contact with a dilute univalent electrolyte solution (solvent blue, cations red, anions green). The xy-plane of a three-dimensional Cartesian coordinate system coincides with the substrate surface, whereas the z-direction points in normal direction towards the bulk of the electrolyte solution.
Conclusions about the general structural features of electrolyte solutions in contact with non-uniformly charged substrates are summarized in Sec. IV.

II. MODEL AND FORMALISM
A. Non-uniformly charged substrate We consider a flat substrate with dielectric constant ε s the surface of which coincides with the x-y-plane of a three-dimensional Cartesian coordinate system; the zdirection is pointing towards the fluid at z > 0 (see Fig. 1). The substrate is non-uniformly charged with the surface charge density σ(u) at the lateral position u = (x, y). In the present study periodic surface charge densities of the form are analyzed. The Fourier coefficients σ k ∈ C, which fulfill the constraints σ * k, = σ −k,− for σ(u) ∈ R, and the lateral length scale L > 0 are free parameters. It will turn out that the periodicity of the lateral surface charge distribution is of no physical relevance, but it is technically convenient.

B. Charged hard spheres
The charged substrate is in contact with a dilute univalent electrolyte solution comprising three species of charged hard spheres: the solvent (species i = 0), cations (species i = +), and anions (species i = −). Each species i is characterized by its hard-core radius R i and the valency Z i with Z 0 = 0, Z + = 1, Z − = −1. For simplicity, all radii are chosen to be equal, i.e., R 0 = R + = R − =: R. The bulk number densities of the electrolyte solution are given by 0 and + = − =: I, which is called the ionic strength. This leads to the packing , which is expressed in terms of the thermal energy β −1 = k B T , the elementary charge e, the vacuum electric permittivity ε 0 , and the fluid dielectric constant ε f , one obtains the Debye length κ −1 with κ 2 = 8π B I.

C. Density functional method
Close to the substrate the number density profile i (r) of the fluid species i varies as function of the position r = (x, y, z) = (u, z), whereas i (u, z → ∞) → i . The set of all three number density profiles is abbreviated by := (ρ 0 , ρ + , ρ − ). The equilibrium number density profiles minimize the grand potential density functional Ω[ ] [29][30][31], which, in the present investigation, is approximated by Here and in the following the common convention is in place that a d-dimensional integration runs over R d unless the integration domain is specified. Equation (2) is to be understood as an asymptotic relation in the thermodynamic limit, i.e., first all calculations are performed in a finite domain which is extended to R 3 subsequently. The thermodynamic limit is guaranteed to exist, i.e., βΩ[ ] scales as the volume of the system, because the number density profiles ρ i (r) are bounded due to the imposed lateral periodicity of the system (see Eq. (1)) and due to the bulk limits i (z → ±∞). In Eq. (2), ζ i = Λ −3 i exp(βµ i ), with the thermal wavelength Λ i and the chemical potential µ i , denotes the (bulk) fugacity of species i ∈ {0, +, −}. The hard-wall potential implies that the fluid particles cannot penetrate into the substrate. The hard-core interaction among the fluid particles is described in terms of the White-Bear (mark I) excess free energy [32], which is given by an excess free energy density Φ(n(r)) expressed in terms of ten weighted densities that are indexed by α and that follow from the number density profiles i via the weight functions ω α,i . The electrostatic potential ψ(r) fulfills Gauß's law where Q(r) := e i Z i i (r) and with the boundary conditions for all u ∈ R 2 in lateral direction. In order to guarantee the existence of the thermodynamic limit we consider a globally charge-neutral system; actually Lebowitz and Lieb have shown that slightly weaker but rather artificial conditions would also suffice [33]. Globally charge-neutral systems exhibit the gauge symmetry ψ → ψ + const, which is used to fix the value of the electrostatic potential at z = ∞ by means of the Dirichlet boundary condition (see Eq. (9)). This implies the Neumann boundary condition ∂ψ ∂z (u,z=∞) = 0. For a globally charge-neutral system the electric displacement −ε 0 ε ∂ψ ∂z has to be the same at z = −∞ and at z = ∞, which leads to Eq. (7). Finally, Eq. (8), which is obtained by integrating Eq. (5) over an infinitessimally small box around the point (u, z = 0), describes the discontinuity of the electric displacement at the charged surface z = 0. The equilibrium number density profiles i (r) vanish for z ≤ R i due to the hard wall (see Eq. (3)), whereas for z > R i they fulfill the Euler-Lagrange equations The set of equations (5) and (7)-(10) is technically too demanding to be solvable numerically for an arbitrary lateral length scale L. In order to proceed, Eqs. (5) and (7)-(10) are first solved for the laterally uniform charge distribution σ (1) := σ 00 , which renders the quasi onedimensional number density profiles, the weighted densities, and the electrostatic potential denoted as (1) i (z), n (1) α (z), and ψ (1) (z), respectively. The quadratic expansion of the density functional in Eq. (2) about The equilibrium profiles ∆ρ i (r) fulfill the Euler-Lagrange equations for functions f (u) of the lateral coordinates u ∈ R 2 , from Eq. (12) one obtains Moreover, the lateral Fourier transformation of Gauß's law (see Eq. (5)) leads, due to Eq. (6), to the Helmholtz . Finally, the boundary conditions Eqs. (7)-(9) take the form Here ∆ σ(q) is the lateral Fourier transform of the nonuniform contribution ∆σ(u) := σ(u)−σ (1) to the surface charge density. The Helmholtz equation at z < 0 (see Eq. (17)) and the Neumann boundary condition at z = −∞ (see Eq. (18)) lead to solutions of the form ∆ ψ(q, z) = ∆ ψ(q, 0) exp(|q|z) for z < 0. Then, from Eq. (19) one obtains the Robin boundary condition which, together with the Dirichlet boundary condition at z = ∞ (see Eq. (20)), determines the solution of the Helmholtz equation in Eq. (16). Note that in the first term of Eq. (21) the upper limit of ∂∆ ψ ∂z occurs, because this quantity is discontinuous at the surface z = 0 due to Eq. (19), whereas in the second term of Eq. (21) ∆ ψ can be evaluated at the surface, because the electrosatic potential is continuous everywhere.
In the set of equations (14)-(21) the individual Fourier modes, indicated by q, are decoupled, and the remaining z-coordinate normal to the substrate leads to a quasi one-dimensional problem, which can be efficiently solved numerically.
Moreover, any function with the the Fourier transform which can be non-zero only for lateral wave numbers q = q k := 2π L (k, l) with k, ∈ Z. Therefore, the determination of the (approximate) equilibrium number density profiles i (r) merely requires to calculate the Fourier transforms ∆ i (q, z) as solutions of Eqs. (14)- (21) for q = q k with k, ∈ Z.

D. Interfacial tension
Besides the profiles i (r) = (1) Eqs. (14)-(21) the following discussion also addresses the interfacial tension γ as a common surface quantity.
Here it is defined w.r.t. the geometrical substrate surface at z = 0. If γ (1) is the interfacial tension of a uniformly charged substrate with surface charge density σ (1) , one obtains the deviation ∆γ := γ − γ (1) due to non-uniformities within the quadratic approximation (see Eq. (11)) as where ∆Ω L means integration over (11), i.e., over one lateral periodic image. This expression can be obtained by multiplying Eq. (12) with ∆ i (r), summing over i, integrating w.r.t. r, and inserting the resulting equation into Eq. (11).

E. Parameters
The main focus of the present study is the dependence of the profiles i (r), Q(r), and ψ(r) as well as of the interfacial tension γ on the characteristic length scale L of the lateral charge non-uniformities. The remaining numerous model parameters are fixed to certain realistic values.
As a non-trivial surface structure we choose a twodimensional square lattice with periodicity L > 0 such that the surface charge density takes the constant value σ max for one half of the surface and 0 for the other half. This leads to an average surface charge density σ (1) = σ max 2 and in Eq. (1) to the Fourier coefficients where the sinc function is defined as sinc(t) = sin(t) t for t = 0 and sinc(t) = 1 for t = 0. However, in order to limit the computational demand only Fourier modes with |k|, | | ≤ 5 are used here. The resulting surface The non-uniform surface charge density σ(x, y) comprising 25 Fourier modes (see Eq. (26)) considered in the present investigation is a continuous approximation of a substrate with half of its area being charge-neutral and the other half being made up of charged square patches of side length L/ √ 2 (indicated by the double arrows in the lower left corner) with surface charge density σmax. The restriction to a finite number of Fourier modes gives rise to slight artifacts such as smooth instead of step-like variations as well as undulations (see the apparent substructure in the dark square areas) instead of plateaus. The mean surface charge density is The charged patches are arranged on a (two-dimensional) square lattice with periodicity L, which sets the lateral length scale of this structure. Figure is a continuous approximation of the actually considered step-like structure (see Fig. 2). In addition to the thermal energy β −1 as the energy unit and the elementary charge e as the charge unit, the Debye length κ −1 is chosen as the length unit. Setting the fluid particle radii to be equal, i.e., R 0 = R + = R − = R, the model comprises the following six dimensionless parameters: In the following, the dependence of structural quantities on L * over two decades is discussed, and two values of the parameter σ * ∈ {−1. Given an aqueous electrolyte solution in contact with a uniformly charged surface the saturation surface charge density σ sat = eκ π B denotes the crossover between a weakly charged surface with |σ (1) | < σ sat , for which the linearized Poisson-Boltzmann (i.e., Debye-Hückel) equation is applicable, and a strongly charged surface with |σ (1) | > σ sat , for which the full non-linear Poisson-Boltzmann equation is required [34]. For the aqueous electrolyte solution specified above, the saturation surface charge density is given by σ sat ≈ 2.2 µC cm −2 , which corresponds to a crossover value σ * sat := The two values σ * = −1.1 and −3.3, which will be considered in the following, have been chosen to represent the cases of weakly and strongly charged surfaces, respectively.

III. RESULTS AND DISCUSSION
A. Normal profiles  As the surface charge is negative, the cation number densities + (u, z) close to the substrate surface are larger than in the bulk, whereas the anion number density profiles − (u, z) close to the substrate are smaller than in the bulk. As expected, these trends are particularly pronounced for highly charged surfaces, i.e., large values of |σ * |, and at lateral positions u corresponding to highly charged regions on the substrate.  i (z) for a uniform charge distribution with surface charge density σ (1) is shown (thin black curves). For z → ∞ all profiles approach the corresponding bulk number densities i . Close to the substrate the typical layering due to the hard cores of the fluid particles is clearly visible. Whereas the solvent number density 0(u, z) (see panels (a) and (d)) varies barely as function of the lateral position u or the surface charge σ * , the cation number density +(u, z) (see panels (b) and (e)) and the anion number densities −(u, z) (see panels (c) and (f)) are sensitive to both u and σ * . Note that the surface charge is negative here, i.e., σ * < 0, so that the cations accumulate at and the anions are depleted from lateral positions close to the charged square patches of the substrate (see Fig. 2).
According to Eq. (14) the lateral structure, expressed in terms of ∆ i (q, z), is determined by the electrostatic potential, represented by ∆ ψ(q, z), as well as by the hard-core interaction, given by the third expression in Eq. (14). Upon ignoring the hard-core contribution one obtains approximate lateral number density variations which resemble those within linear Poisson-Boltzmann (i.e., Debye-Hückel) theory. Inverse Fourier transformation leads to so that It turns out, that the approximation i (r) ≈ DH i (r) is reliable to a high degree, i.e., the hard-core contribution as the last term in Eq. (14) can be safely ignored. Whereas the hard-core interaction plays an important role for the number density profiles (1) i (z) close to laterally uniformly charged substrates, it does not influence the lateral structure formation significantly.
Since the hard-core contribution as the last term in Eq. (14) is quantitatively negligible, one ends up with the approximation which is equally valid. Upon inserting Eq. (32) into the Helmholtz equation for z > 0 (see Eq. (16)) one obtains with the abbreviation For z → ∞ the quantity κ(z) approaches the inverse Debye length, i.e., κ(z) → κ, as in this limit (1) ± (z) → I. Figure 5 shows, that κ(z)/κ attains its bulk value 1 already a few particle radii R away from the substrate. Hence beyond a few particle radii R away from the substrate, i.e., at z R, Eq. (33) reduces to with the solution with the normal decay length λ(q) := 1 According to Eq. (29), in the range z R the modes of the lateral structure ∆ i (q, z) decay on the same normal length scale λ(|q|). Whereas the decay length λ(q) is a bulk quantity, the proportionality prefactor of the  (1) − (z) (see Eqs. (33) and (34)).

B. Lateral profiles
The length scale, on which the lateral modes q decay in the normal direction, is given by λ(|q|) (see Eq. (37)). It attains its maximum value κ −1 , i.e., the Debye length, at q = 0. Accordingly, the normal decay length λ(|q|) is not larger than the Debye length κ −1 . Upon increasing |q| the normal decay length λ(|q|) decreases monotonically.
Since σ (1) = σ 00 and i.e., ∆ σ(q = q 00 = 0) = 0 due to Eq. (23), the smallest wave number |q| = |q k | = 2π L k 2 + 2 contributing to a lateral structure is q min = |q ±1,0 | = |q 0,±1 | = 2π L . Hence the lateral structure induced by a non-uniformly charged substrate decays in normal direction on the length scale (39)) crosses over from a linear regime at short length scales L * 2π to λmax κ −1 (i.e., the Debye length) at large length scales L * 2π. Figure 6 displays the dependence of λ max on the length scale parameter L * = κL. At short length scales L * 2π a linear dependence is found, which crosses over to λ max κ −1 (Debye length) at large length scales L * 2π. From the quantitatively reliable approximation i (r) ≈ DH i (r) (see the previous Subsec. III A and in particular Fig. 4) one can infer that the lateral structure of the electrolyte solution, i.e., i (r), is determined by the lateral structure of the electrostatic potential ψ(r) (see Eq. (31)). Accordingly, Fig. 7 displays two sequences of lateral profiles of the electrostatic potential ψ(u, z), i.e., functions of u with z fixed, at the normal positions z = 0, R, λ max , and 2λ max for L * = 0.6 (left column: (a), (c), (e), (g)) and 60 (right column: (b), (d), (f), (h)). Qualitatively, the difference of the electrostatic potential between lateral positions associated with large and with small surface charge densities diminishes with increasing distance from the substrate. However, although the normal decay length λ max is very different for the two cases (κλ max ≈ 0.1 for L * = 0.6 and κλ max ≈ 1 for L * = 60), the decay of the lateral structure of the two as function of z/λ max is similar.

C. Interfacial tension
The findings discussed so far lead to the picture of a surface layer of thickness λ max in which a non-uniform surface charge density, characterized by a lateral length scale L, can be sensed by the electrolyte solution. This thickness λ max is found to increase as function of L as long as L * = κL 2π, whereas it is approximately constant for L * 2π. Therefore, one can expect that the interfacial tension γ (see Subsec. II D) exhibits the same trend. This is indeed the case, as it is shown in Fig. 8 y/L  (31)). The left column corresponds to the case L * = 0.6 (short lateral length scale) and the right column to the case L * = 60 (large lateral length scale). For each case, the lateral structure decays in normal direction on the length scale of λmax (see Eq. (39)). The contrast between charged and neutral parts of the substrate at z = 0 is clearly visible (see panels (a) and (b)). In the first contact layer of the fluid at z = R (see panels (c) and (d)), the contrast is still present, but slightly blurred. At the distance z = λmax (see panels (e) and (f)) the contrast is diminished substantially and even more so at z = 2λmax (see panels (g) and (h)). The decay of the lateral structure as function of z/λmax is similar, irrespective of the lateral length scale L * . Upon increasing z, for L * = 0.6 the rectangular shape of the charged pattern is washed out in favor of circular patterns. For L * = 60 the rectangular shape of the patterns remains even for z = 2λmax.
for the surface charges σ * = −3.3 (red curve) and −1.1 (blue curve). However, the interfacial tension γ of a nonuniformly charged substrate turns out to be limited to at most a few percent above the interfacial tension γ (1) of a uniformly charged substrate with the same mean surface charge density.

IV. SUMMARY, CONCLUSIONS, AND OUTLOOK
The present investigation is devoted to the structure formation in a dilute electrolyte solution close to a nonuniformly charged planar substrate (see Fig. 1). In dilute electrolyte solutions the Debye screening length κ −1 is substantially larger than the size of the fluid molecules R so that, in principle, the spatial region of according thickness κ −1 close to a charged substrate can be sensitive to the surface charge distribution. However, the lateral length scale L of the charge distribution on the substrate turns out to play a role, too. In the present study periodic charge distributions with periodicity L of arbitrary magnitude are considered (see Fig. 2), and the corresponding laterally non-uniform number density profiles of the fluid particles are calculated via expansion about the profiles of a uniform substrate with the same mean surface charge density (see Fig. 3). It is found that the lateral structure is mainly determined by the electrostatic potential, i.e., not by molecular-ranged forces like the hard-core interaction, so that the laterally nonuniform contributions of the number density profiles can be accurately approximated by a Debye-Hückel-like expression (see Eq. (31)), disregarding hard-core contributions (see Fig. 4). As a consequence, for normal distances not too close to the substrate, i.e., at z-coordinates with κ(z)/κ in Fig. 5 close to unity, the lateral contributions of the electrostatic potential, and hence of the number densities, decay on the scale λ max given in Eq. (39) (see Fig. 6). For lateral length scales L with L * = κL 2π, the normal decay length is varying with L according to λ max ≈ L/(2π), whereas for L * 2π it levels off at the value of the Debye length, λ max ≈ κ −1 . As shorter length scales L * 2π decay more rapidly than larger ones, a washing out of fine details at increasing distance from the surface occurs (see Fig. 7). Ultimately only structures at length scales L * 2π contribute to the lateral structure. In terms of the interfacial tension of the non-uniformly charged substrate an increase with L is observed for L * 2π which saturates for L * 2π (see Fig. 8).
Equation (36) in conjunction with Eq. (31) states, that at distances z λ(q) (see Eq. (37)) details of a surface charge distribution with wave number q = |q| become irrelevant for the lateral structure of an adjacent electrolyte solution. Hence, at larger distances from the substrate, only less fine details of a surface charge distribution can be resolved. Ultimately, at distances z κ −1 details with wave numbers q = |q| κ, i.e. with lateral length scales L 2πκ −1 , are washed out so that only surface structures with lateral length scales L 2πκ −1 matter. The strength of the influence of these large-scale structures decay exponentially with a decay length given by the Debye screening length κ −1 . Therefore, when modeling electrolyte solutions with molecular length scale R, one can safely ignore surface non-uniformities at length scales L 2πR, which, for molecular fluids, can be close to a nanometer. Finally, the present study shows that macroscopic descriptions of electrolyte solutions, i.e., on length scales larger than the Debye length κ −1 , are carried out consistently by considering surface details on lateral length scales larger than 2πκ −1 only.
Two main conclusions can be drawn from the present study: (i) Microscopic hard-core interactions have negligible influence on the lateral structure formation of electrolyte solutions close to non-uniformly charged substrates. (ii) Fine details of lateral non-uniformities have negligible influence beyond a certain (short) distance from the substrate. Accordingly, the approach of disregarding the size of molecules and treating them as point particles (see many previous theoretical studies concerning the interaction between non-uniformly charged colloidal particles) can be justified or readily adjusted. Generally, the present study shows that on macroscopic length scales only macroscopically large features of the surface structure are visible. This allows for local descriptions of fluids in terms of partial differential equations, e.g., the Young-Laplace equation in hydrostatics. The dominant correlation length of a fluid, which for a dilute electrolyte solution of a non-critical solvent is the Debye length, separates length scales into macroscopic and microscopic ones. From a microscopic point of view, there is a smooth crossover of the fluid structure from small to large length scales, whereas microsopic details can be safely ignored from a macroscopic point of view.
Several directions of applications of the gained insight are conceivable: The presented approach, i.e., to consider deviations from laterally uniform reference density profiles and to ignore hard-core interactions, could be exploited in various numerical analyses of fluid structures, including computer simulations. This way studies of large laterally non-uniform systems could become feasible. Furthermore, given a certain length scale, the above insight is useful in order to distinguish relevant from irrelevant surface details. This is of importance not only for theoretical considerations or numerical applications, but also for efficiently solving practical problems, such as guiding flows in nanochannels, patterning sur-face structures of catalytic reactors, or designing electrochemical devices. Finally, a common understanding of the small effect microscopic features have on macroscopic length scales (and vice versa) could be helpful for the sci-entific discourse by avoiding confusion when comparing experimental or theoretical results obtained within methods the spatial resolution of which are associated with incompatible length scales.