{'fraction_non_alphanumeric': 0.03372126996493962, 'fraction_numerical': 0.012636345929100117, 'mean_word_length': 4.496186270574067, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
One-dimensional counterion gas between charged surfaces: Exact results compared with weak- and strong-coupling analysis

We evaluate exactly the statistical integral for an inhomogeneous one-dimensional counterion-only Coulomb gas between two charged boundaries and from this compute the effective interaction, or disjoining pressure, between the bounding surfaces. Our exact results are compared with the limiting cases of weak and strong coupling which are the same for 1D and 3D systems. For systems with a large number of counterions it is found that the weak coupling (mean-field) approximation for the disjoining pressure works perfectly and that fluctuations around the mean-field in 1D are much smaller than in 3D. In the case of few counterions it works less well and strong coupling approximation performs much better as it takes into account properly the discreteness of the counterion charges.


I. INTRODUCTION
Electrostatic effects are predominant in a variety of soft condensed matter systems such as polymers in solutions, colloidal suspensions and much of the physics of membranes and films [1,2]. The statistical mechanics of these problems cannot in general be solved explicitly; however a number of approaches have been devised, starting with the mean-field Poisson-Boltzmann equation, that allow for an approximate solution to the problem [3]. One can improve the accuracy of mean-field calculations by looking at the fluctuations about the Poisson-Boltzmann solution, and in a number of cases these fluctuation effects are essential to capture the quantitative behavior of the systems under study. The approach of this type can be shown to be valid in the so called weak-coupling regime. More recently the strong-coupling expansion, in its essence similar to the virial expansion, has been developed and successfully applied to a number of interesting problems [3,4]. Despite the success of these approaches many questions remain to be answered about the intermediate regime and how the cross-over between the weak-and strong-coupling limits occurs. Here we show that the one-component 1D inhomogeneous Coulomb fluid can be analytically solved and thus provides an ideal test-bed to study the domains of validity of these commonly used approximation schemes.
For a system of two charged surfaces interacting in an electrolyte solution, the effective coupling strength may depend on the inter-plate separation, and so the calculation of the force between the plates requires an understanding of the both weak and strong coupling regimes and the transition between them. In two important papers by Lenard [5] and Edwards and Lenard [6] a complete solution of the thermodynamics of the two component Coulomb gas in one dimension was given. The second of these papers was based on the fact that the underlying Sine-Gordon theory for this system can be mapped onto quantum mechanics in one dimension and the thermodynamics can be solved in terms of the ground state energy of the corresponding Schrödinger equation. Subsequently the system was studied in the presence of an electric field and it was shown that it always behaves as a dielectric because of a dimerization of the individual charge carriers [7]. When the system is not in an electroneutral state, i.e. when the surface charges corresponding to the applied field are not exactly compensated by the counterions of the system, the physics of the system is subtly different and confinement phenomena between oppositely charged particles appear [8,9]. This confinement is directly related to the fact that the system is dielectric and non-conducting. More recently, the onedimensional Coulomb gas thermodynamics has been studied taking into account additional hard-core and dipolar interactions [10]. The model has also been studied in the finite size setting of a soap film type model where charge regulation mechanism model of the surface, in the presence of salt solution (ions and counterions), is included [11].
In the present paper we study a related model to those outlined above: a one-dimensional model of counterions confined between two plates, carrying fixed but not necessarily equal charges. This situation is quite different to most of those studied above as most of them have concentrated on symmetric electrolyte systems. We place particular emphasis on the computation of the pressure of the system which is the effective interaction between the two plates modified by the presence of counterions. The effect of asymmetric charge distribution on the plates is analyzed in detail. This system, while being somewhat idealized, is a simple model of an array of charged smectic bilayers sandwiched between two fixed parallel charged plates. If the charge on the bilayers is uniform then the use of the one-dimensional Coulomb interaction is justified.
Another reason for studying this simple system is to check the validity of the various approximation schemes alluded to above. The mean-field and strong-coupling approximations in the planar geometry are independent of the dimensionality of the problem. This is due simply to the fact the both Poisson-Boltzmann theory as well as the strong-coupling theory are one-dimensional effective theories since the corresponding fields depend only on the coordinate parallel to the bounding surfaces normals. However, fluctuations about the mean field, or the Poisson-Boltmzann configuration, do depend on the dimensionality of the problem as they correspond to thermal Casimir or zero-frequency van-der-Waals forces. This makes these model calculations particularly appealing since they can clarify the role the fluctuations play in Coulomb systems.
The paper is organized as follows: first we use the Edwards-Lenard path integral formulation for the problem adapted to apply to a finite electroneutral system with surface charges. The exact expression obtained for the force between the two plates is evaluated numerically and compared with the results of Monte-Carlo simulation for a wide range of cases: symmetric/asymmetric surface charges, number of counterions, etc. We find excellent agreement between theory and simulation in all cases as we should expect since the theoretical solution is exact. In the next section we compare the results with the Poisson-Boltzmann mean-field approximation and investigate its domain of validity. We present two methods for calculating the effects of field fluctuations about the mean-field solution and compare with our exact results. The strong-coupling approximation is then applied to the system and a similar comparison with the exact results is made.

II. THE COUNTERION MODEL
We consider a system of charged particles interacting via a Coulomb potential in one dimension. It can also be considered as a system of uniformly charged infinite two dimensional sheets in three dimensions, where the sheets are perpendicular to the direction x, are uniformly charged and can only move in this direction. The Hamiltonian for this system is given by where, e is the unit electric charge, β = 1/k B T , z α is the valence of particle (or sheet) α and x α is its position. We will consider systems in a canonical formulation which are overall charge neutral, so that It is convenient to rewrite the Hamiltonian as where we have used charge neutrality to go from the first to second equation above. We may now write the Hamiltonian as an expectation value over a "Brownian motion" where ψ is Brownian motion with correlation function ψ(x)ψ(y) = βe 2 min(x, y), and started at the value zero at x = 0, so The Boltzmann weight for any configuration can thus be written as In path-integral notation the measure on the Brownian motion ψ as defined by Eq. (5) is given by where L is the overall length of the system. The sum over α has three distinct contributions: that from the counterions with 1 ≤ α ≤ N , for which from now on we assume have valence z α = 1, and those from the surface charge at x 0 = 0, with valence z 0 = −(N + P )/2, and from the surface charge at x N +1 = L, with valence z N +1 = −(N − P )/2. The system is globally electroneutral, consisting of two surface charges and their counterions. For the moment, we will assume that N ± P are even integers, which will be the case if the counterions originate by being released into solution from initially neutral surfaces; i.e., each surface charge is an integer multiple of the unit of charge of the counterions. The generalization of the approach due to relaxing the latter assumption is given in the next section. We then have and the Boltzmann weight for a given configuration is thus The partition function is then obtained as We now introduce an arbitrary fugacity κ having the dimensions of inverse length and write although, clearly, the value of κ should not affect the final physical results since charge neutrality will enforce the constraint N ′ = N . This is ensured by the integral over dλ. If we define a shifted Brownian motion, the starting position of the field φ at x = 0 is simply φ(0) = λ but the integration over the end point remains free and we have The electroneutrality constraint has thus been absorbed into the integration over the starting point of the shifted field φ.
It is convenient, for what comes later, to generalize the result to the case where the counterions have valence z α = q, 1 ≤ α ≤ N . This is easily accomplished by replacing e by qe in the measure for the Brownian motion in Eq. (8). This follows from the forgoing analysis and rescaling the field φ(x) → φ(x)/q. In path integral notation we then have as our final functional integral expression for the partition function.

III. EXACT EVALUATION OF THE PARTITION FUNCTION
We use the Feynman formula, which is the Euclidean version of the equivalence between path integrals and the propagator for the Schrödinger equation. We have that where q is the counterion valence and whereĤ is the complex operator/Hamiltonian given bŷ The partition function can thus be expressed as We should note that the operatorĤ is complex in our problem but in the case of symmetric electrolytes [6] it is real. Because of this we analyse the problem via Fourier analysis rather than on the basis of the eigenfunctions of the operatorĤ.
We define integers M 1 = (N + P )/2 and M 2 = (N − P )/2 and the surface charges are −M 1 qe and −M 2 qe at x = 0 and x = L, respectively, with M 1 + M 2 = N . A quick preliminary check of our approach is to consider the perfect gas limit where e = 0. In this case we find trivially that and thus which is the perfect gas result. In order to proceed further we consider the evaluation of We can write [11] f (z) = n a(n, L) exp(inz) , and the induced evolution equation for the Fourier coefficients is da(n, L) dL = − n 2 q 2 e 2 β 2 a(n, L) + κ a(n − 1, L) .
In our problem we have f (z, M 2 ; 0) = exp (−izM 2 ) and thus the initial condition on the Fourier coefficients a(n, M 2 ; 0) = δ n,−M2 . We define a n = κ n b n and find that with the initial condition b(n, M 2 ; L) = κ n δ n,−M2 and from Eq. (18) we find that we can express the partition function for system of unit charges M 1 and M 2 on the left and right boundaries respectively as but with the initial condition b n (0, M 2 ) = δ n,−M2 ; as we expect, the fugacity κ, which was introduced on dimensional grounds, does not enter in the final physical result. The evolution Eq. (24) can be written as where p M1,M2 is the pressure of the system. Clearly, the last term makes a positive contribution to the pressure as it is the ratio of two partition functions for systems with different surface charges. Equation (26) has a nice physical interpretation. The second term is, up to the factor of β, the average density of counterions at z = 0 because it is the restricted partition function where at least one particle is on the surface z = 0 (thus giving a reduction in the surface charge to M − 1) normalized by the full partition function. We can thus write where is the average value of the counterion density at the surface z = 0 and σ 1 = −M 1 qe the corresponding surface charge. However, this is simply the contact-value theorem for electrostatic systems known to be exact in any dimension [10,12,13,14]. In fact the contact-value theorem can be demonstrated via an extension of the path integral methods used here to higher dimensional systems with (parallel) planar geometries [15]. The average counterion density at the point x can be shown [6,11] to be equal to If we set x = 0 in this formula we obtain in agreement with the previous discussion.

A. The general formalism
The charges on the surfaces at z = 0 and L can be denoted σ 1 and σ 2 respectively. Without loss of generality we may assume that q > 0 and |σ 2 | ≥ |σ 1 |, which allows σ 1 and σ 2 to be of opposite sign, and we can define the charge asymmetry parameter ζ = σ 2 /σ 1 with −1 ≤ ζ ≤ 1. All other cases can be mapped onto this interval with appropriate rescaling of the parameters. Note that the values ζ = 1 and ζ = −1 represent special cases of symmetric system with σ 1 = σ 2 , and antisymmetric system with σ 1 = −σ 2 and no counterions between surfaces (which thus reduces to the trivial case of a planar capacitor).
For the above analysis we see that σ 1 = −M 1 qe and σ 2 = −M 2 qe. The independent parameters are conveniently chosen as σ 1 , ζ, N . With M 1 and M 2 integers, the analysis so far allows only for particular discrete values of ζ. However, the approach can be extended to accommodate all values of ζ and we state the generalization of the algorithm here. Charge neutrality determines that the valence q of the counterion satisfies We are allowing that q need not be an integer.
In addition to the Fourier coefficients b(n, M 2 ; L) we introduce Fourier coefficients c(n, M 1 ; L). Equation (24) gives the evolution of the b(n, M 2 ; L) from z = 0, the surface with assigned charge σ 1 . The coefficients c(n, M 1 ; L) obey a similar evolution equation but now with the surface charge assignment reversed; the charge at z = 0 is now σ 2 . These two sets of coefficients correspond to complementary approaches; the set {b} describes evolution from the surface with charge σ 1 to that with charge σ 2 , and vice-versa for the set {c}.
We define The evolution equations now become The initial conditions are b(n, Then the partition function is given by alternative formulas , and the pressure is given equivalently by The partition function Z(σ 1 , ζ, N ; L) is also given by for any x, 0 ≤ x ≤ L. This acts as a check on the numerics. The counterion number density ρ(x; σ 1 , ζ, N, L) is given by and another check on the numerics is that i.e. the system is electroneutral.

IV. APPROXIMATE EVALUATIONS OF THE PARTITION FUNCTION
A traditional approach to the study of charged (bio)colloidal systems is the mean-field Poisson-Boltzmann (PB) formalism which is applicable for weak surface charges, low counter-ion valency and high temperature [16]. The limitations of this approach are evident when applied to highly-charged systems where counterion-mediated interactions between charged bodies start to deviate substantially from the accepted mean-field wisdom [3,4]. One of the recent fundamental advances in this field has been the systematization of these non-PB effects based on the notions of weak and strong coupling approximations. The latter approach has been pioneered by Rouzina and Bloomfield [17], elaborated later by Shklovskii et al. [18], Levin et al. [19], and brought into final form by Netz et al. [3,4,20,21]. These two approximations allow for an explicit and exact treatment of charged systems at two disjoint limiting conditions whereas the parameter space in-between can be analyzed only approximately [20,21,22,23,24,25] and is mostly accessible solely via computer simulations [3,4,20,21,24,25,26,27,28,29,30,31].
Both the weak-and the strong-coupling approximations are based on a functional integral or field-theoretic representation [32,33] of the grand canonical partition function for a system composed of fixed surface charges with intervening mobile counterions, and depend on the value of a single dimensionless coupling parameter Ξ [20,21]. In three dimensions, Ξ is proportional to the ratio between two relevant length scales, namely, the Bjerrum length and the Gouy-Chapman length. The Bjerrum length, defined as ℓ B = e 2 /(4πεε 0 k B T ), is the distance at which two unit charges interact with thermal energy k B T (in water at room temperature, one has ℓ B ≃ 0.7nm). If the charge valency of the counterions is q then the corresponding length scales as q 2 ℓ B . Similarly, the Gouy-Chapman length, defined as µ = e/(2πqℓ B |σ|), is the distance at which a counterion interacts with a macromolecular surface (of surface charge density σ) with an energy equal to k B T . The 3D electrostatic coupling parameter measures the competition between ion-ion and ion-surface interactions and is given by [20,21] Physically, the weak-coupling (WC) regime Ξ ≪ 1 (appropriate for low valency counterions and/or weakly charged surfaces), is characterized by the fact that the width, µ, of the counterion layer near the surfaces is much larger than the separation between two neighboring counterions in solution, and thus the counterion layer behaves basically as a three-dimensional gas. Each counterion in this case interacts with many others and the collective mean-field approach of the Poisson-Boltzmann (PB) type is completely justified.
On the other hand in the strong-coupling (SC) regime Ξ ≫ 1 (appropriate for high valency counterions and/or highly charged surfaces), the mean distance between counterions, a ⊥ = |σ|/qe, is much larger than the layer width (i.e., a ⊥ /µ ∼ √ Ξ ≫ 1), indicating that the counterions are highly localized laterally and form a strongly correlated quasi-two-dimensional layer next to a charged surface. In this case, the weak-coupling approach breaks down due to strong counterion-surface and counterion-counterion correlations. Since counterions can move almost independently from the others along the direction perpendicular to the surface, the collective many-body effects that enable a mean-field description are absent, necessitating a complementary SC description [20,21].
Formally, the WC limit can be identified with the saddle-point approximation of the field theoretic representation of the grand canonical partition function, and reduces to the mean-field PB theory at the lowest order for Ξ → 0. The quadratic fluctuations around the mean field provide a second-order correction to the mean-field solution for small Ξ < 1 [32,33,34,35,36,37,38,39]. The SC approximation has no PB-like collective mean-field [20,21] since it is formally equivalent to a single particle description obtained from a systematic 1/Ξ expansion in the limit Ξ → ∞, and corresponds to two lowest order terms in the virial expansion of the grand canonical partition function. The consequences and the formalism of these two limits of the Coulomb fluid description have been explored widely and in detail (for reviews, see Refs. [3,4]).
The concept of weak-and strong-coupling electrostatics can be easily generalized to other dimensions. In one dimension (1D), the Coulomb interaction between two unit charges may be written as v 1D (x) = −|x|k B T /ℓ B , where ℓ B = 2k B T /e 2 may be regarded as a 1D Bjerrum length. Likewise, the interaction strength with the boundary charges in an asymmetric system may be characterized by the 1D Gouy-Chapman lengths µ 1 = ℓ B e/q|σ 1 | and µ 2 = ℓ B e/q|σ 2 |. We may proceed by using only µ ≡ µ 1 = ℓ B e/q|σ 1 | in what follows since for any given asymmetry parameter ζ = σ 2 /σ 1 , we have µ 2 = µ/|ζ|. The ratio is then the corresponding electrostatic coupling parameter in 1D. (On a more formal level, this definition for the coupling parameter may be established by looking at the field-theoretical representation (15) and rescaling the coordinates as x i → x i /µ 1 .) Note that in 1D the coupling parameter is related to the asymmetry parameter by virtue of the electroneutrality condition where N is the number of counterions of valency q, as Thus, it may be expected that in the limit N → ∞, contrary to the 3D case, the mean-field theory corresponding to this system becomes exact! On the other hand, the SC description may be expected to follow simply for N = 1. These limiting cases will be discussed further in the forthcoming sections and will be compared with the exact results and Monte-Carlo simulations. One should bear in mind that the 1D system considered here corresponds to a 3D system of mobile charged plates (membranes) confined between two fixed planar charged walls and thus, such a 1D system with even a single "counterion" would have a completely meaningful thermodynamic behavior.
We note also that both the PB theory (Ξ → 0) and the SC theory (Ξ → ∞) for uniformly charged plates are onedimensional theories and should remain valid in 3D as well as in the 1D case that we are studying here. However, even though the mean-field solution depends only on the transverse coordinate, the field fluctuations about the mean-field solution will depend on all coordinates and so it is clear that fluctuations corrections to the mean-field contribution depend on the dimensionality of the system: they will be different for a 1D system than for a 3D system.

A. Weak-coupling limit: Poisson-Boltzmann theory
The expression for the partition function Z N in Eq. (15) is up to multiplicative constants given by the functional integral where the action S is given by In the weak-coupling regime, the leading contribution to the partition function comes from the saddle-point configuration, φ 0 (x), of the action in Eq. (43) [32,33], where the field φ 0 (x) = iψ 0 (x) turns out to be imaginary and is proportional to the mean-field electrostatic potential. The saddle-point configuration can be straightforwardly translated into a solution of the PB equation and corresponds to an exact asymptotic result in the limit Ξ → 0 [20,21,34]. The PB equation for the potential, ψ 0 (z), can be written as [2,16] with boundary conditions stemming from the electroneutrality of the system. Integration of the PB equation gives rise to the first integral of the system of the form where the constant p 0 is nothing but the mean-field PB pressure acting between the bounding surfaces [2] and is the PB number density profile of counterions between the surfaces. As remarked earlier in section III, the choice for κ is arbitrary since it corresponds to a choice of origin for the potential ψ(x). In section III we couched the Fourier solution in terms of Fourier coefficients a n = κ n b n with the choice κ = 1, and so we adopt this choice here, too. The nature of the solution ψ 0 (x) obviously crucially depends on the sign of the pressure p 0 [40,41]. Different forms are obtained for positive and negative pressures, corresponding to repulsion and attraction between the bounding surfaces respectively and were derived previously [42].
Guided by the contact-value theorem Eq. (27), which we emphasize is an exact result, we express the Poisson-Boltzmann, or mean-field, pressure in rescaled units. The repulsive PB pressure is given bỹ and all the other quantities have been defined above. Furthermore,α = αµ and its value is given by the solution of For ζ < 0, and only then, can the PB pressure be attractive. The attractive PB pressure is given bỹ whereα is now given as a solution of The limiting case of zero pressure can be obtained straightforwardly from either limit. The interaction pressure can be obviously computed for any value of the asymmetry parameter, but note that any value of ζ can be mapped onto the interval −1 < ζ < 1 and the pressure only need to be evaluated in that interval of ζ values. Note again that the PB interaction pressure is the same for a 1D as well as for a 3D system.

B. Weak-coupling limit: fluctuations
Fluctuations around the mean field depend on the dimensionality of the system and are different for a 1D than for a 3D system. In order to proceed, one needs to evaluate the appropriate Hessian of the field action in the partition function and study its fluctuation spectrum (see Refs. [32,33,43] for more details). The Hessian of the field about the mean-field value is given by ρ 0 (x) = exp(−ψ 0 (x)) is the mean-field PB counterion density. Hence, The corresponding correction, F 2 , to the free energy of the system comes from the functional integral where so that the corresponding fluctuation part of the free energy can be obtained in the form The functional integral Eq. (54) can be evaluated exactly, in two different ways. The first method is based on the use of the argument principle [36] converting the discrete sum of eigenvalues of the Hessian operator into the logarithm of the secular determinant D λ of the same operator. The second method is based on the Pauli-van Vleck approach to calculating the functional integral of a general harmonic kernel. Both results are the same.
In the first approach the trace-log of the Hessian can be written equivalently in the form that was derived for a 3D case [32,33,36] but can be used in a trivially modified form also for the 1D case under consideration. It gives where D λ is the secular determinant, i.e. the determinant of the coefficients corresponding to the appropriate boundary condition in the solution of the eigenvalue equation, that can be derived from solutions of One can now write A(L) for the quotient (D 1 /D 0 ) −1/2 , since the secular determinant depends explicitly on the value of the inter-surface spacing, L. Using the fact [32,33,43] that the two linearly independent solutions of Eq. (58) for λ = 1 in the repulsive regime are the secular determinant for the symmetric case ζ = 1 comes out as and the corresponding free energy contribution from the quadratic fluctuations is thus given by The fluctuation free energy can be regularized so that all irrelevant constants, i.e. all the terms not depending on the separation between the bounding surfaces, are dropped, amounting to a rescaling F 2 (L) → F 2 (L) − F 2 (L → ∞). This corresponds to a subtraction of the part of the free energy for two separate boundaries at infinite separation from the total free energy.
In the second method we compute the fluctuations about the mean-field path using the Pauli-van Vleck approach [44,45,46]. As before the thickness of the film is L we take the leftmost and rightmost points of the film to be at x = 0 and x = L respectively. Clearly the action of the classical path minimizing s in Eq. (54) is a quadratic function of the initial and final points of the fluctuating field ψ. The generalized Pauli-van Vleck formula tells us that ψ(L)=y where s c is the classical action minimizing s. As the action is quadratic we may write In the case of a symmetric charge distribution we have f = g and thus the fluctuation term is given by (64) Solving the equations of motion we find that .
As we already stated the fluctuation part of the interaction free energy for 3D and 1D are completely different (due to fluctuations in the plane of the film in 3D which are not present in 1D), though the mean-field Poisson-Boltzmann result is exactly the same. In particular, the scaling of the fluctuations contribution to the interaction pressure, say p 2 , with the inter-surface distance L turns out to be very different. In 3D, one finds that [42] whereas in 1D we find for same-sign surfaces (ζ ≥ 0) at sufficiently large separations. Note that the coupling parameter Ξ is defined differently in 1D, Eq. (39), and 3D, Eq. (38). As one goes to the weak-coupling limit in 1D for N → ∞ the fluctuation term obviously makes a vanishing contribution to the total pressure and is thus not particularly important. The logarithmic dependence in 3D is a direct consequence of the contribution from the in-plane modes. In this case, the mean-field pressure scales with the inverse square of the separation as p 0 (L) ∼ 1/L 2 irrespective of the dimensionality of the system. Note that the fluctuations part is always attractive, reflecting the fact that electrostatic correlations mediated by counterions always favor attraction between the charged boundaries. Within the WC analysis the fluctuations are always assumed to be small as compared with the leading order mean-field contribution. Thus, the total pressure, p 0 + p 2 , is dominated by the mean-field contribution (except at the equilibrium point where p 0 = 0 [42]) and in particular for same-sign surface charges, it still remains repulsive.
It is also interesting to note that the magnitude of fluctuation pressure relative to the mean-field pressure decreases with the separation distance L in 3D, while in 1D, it increases with L. For same-sign surfaces, the ratio |p 2 /p 0 | scales as ∼ ln L/L in 3D, while it scales as ∼ L in 1D. This indicates qualitatively different distance-dependent behaviors for the fluctuations and thus, qualitatively different regimes of validity for the loop-expansion approach in 1D and 3D. The latter is determined by assuming that |p 2 /p 0 | ≪ 1. Hence, the weak-coupling validity criterion for same-sign surfaces in 3D reads This means that at a given coupling parameter, the WC analysis becomes increasingly more accurate at larger separations, while as the surfaces get closer a smaller coupling parameter needs to be chosen. In 1D, one has which indicates the opposite trend. In particular, for given ζ and noting from Eq. (41) that Ξ 1D = (1 + ζ)/N , the WC analysis becomes increasingly more accurate as the separation decreases for fixed N , or as N increases for fixed separation.

C. Strong-coupling limit
The strong-coupling approximation coincides with the lowest order of a non-trivial expansion of the partition function in terms of the fugacities of the counterions. This expansion may be expressed as a 1/Ξ series expansion [20,21], whose leading order term (Ξ → ∞) corresponds to the SC theory. We will not delve into the strong-coupling expansion in more detail since it has been exhaustively reviewed in the literature [3,4,20,21].
At leading order, the SC free energy is obtained as where W 0 is electrostatic interaction energy of charged surfaces with S representing surface area of each plate, and W 1 and W 2 are electrostatic interaction energies between a single counterion and individual charged surfaces, i.e.
Since in the strong-coupling regime the free energy is given via simple quadratures, it is much simpler to evaluate it than at the weak-coupling level. Defining the rescaled free energỹ we obtainF whereL = L/µ. Differentiating the free energy with respect to the surface-surface distanceL we get the corresponding pressure acting between the bounding surfaces where, following the discussion in section IV A, we have definedp(L) = 2p(L)/σ 2 1 . Note that the SC pressure can become attractive for both like-charged and oppositely charged surfaces which contrasts with the mean-field theory that does not allow attraction between like-charged surfaces. This is because of the strong electrostatic correlations mediated by counterions between the charged surfaces for Ξ ≫ 1 and has been investigated thoroughly before for equally charged surfaces in [20,21] and for asymmetric surfaces in [42].

V. NUMERICAL SIMULATIONS
We next consider Monte-Carlo simulations of the 1D system of counterions treated in the preceding sections using both exact as well as approximate (limiting) analytical approaches. Monte-Carlo simulations enable us to access the parameter space inaccessible to the limiting WC and SC theories, and thus may be compared directly with the exact solution presented in Section III.
We proceed by simulating a system of N counterions in a finite interval x ∈ [0, L] in the canonical ensemble by applying a standard Metropolis algorithm. The Hamiltonian of the system is defined in Eq. (1) and the electroneutrality condition is imposed via Eq. (40). We run the simulations by up to 10 8 − 10 9 Monte-Carlo steps per particle with 10 7 − 10 8 steps used for relaxation purposes.
The interaction pressure between the two charged boundaries is calculated via the contact-value theorem (see also Eq. (27)) where ρ(0) and ρ(L) represent the counterion density at contact with the first (σ 1 ) and the second (σ 2 ) plate, respectively. The resulting pressure computed from the contact density at either surface is found to be the same within the numerical errorbars. Simulations were conducted in rescaled units for various number of counterions N = 1, . . . , 20 and for the asymmetry parameters ζ = −0.5, 0.5 and 1.0. The coupling parameter in each case follows from Eq. (41). The simulated pressure is shown in rescaled units in Fig. 1 (symbols) for the rescaled pressurep = βpℓ B /(σ 1 /e) 2 as a function of the rescaled distance,L = L/µ, between the charged boundaries. As seen, the interaction pressure decays monotonically with separation distance L for all values of ζ and N . Also in all cases, the simulation data (symbols) and the exact results (thin solid curves) are nicely bracketed by the mean-field PB result (thick solid curve) and the SC result (dashed curve).
For the symmetric case ζ = 1, the simulation results are spot on the SC curve (dashed line) for N = 1 and approach slowly the attractive asymptotic SC pressure at large separationsL → ∞, that isp SC ∞ = −ζ 2 . As the number of counterions increases the pressure becomes less attractive and eventually for sufficiently large N , the simulation data tend to the PB curve (thick solid curve) exhibiting repulsive pressure at all separations. This is because for large N the system effectively splits into two nearly electroneutral halves, where each bounding surface is neutralized by its corresponding layer of counterions. This also indicates that the entropic contribution from counterions which favors repulsion becomes important at large N . For finite N , the simulation data is described neither by the PB theory nor by the SC theory. In this case, an excellent agreement is found between the data and the exact results given by Eq. (34) (shown by thin solid lines).
Similar trends are observed for asymmetric systems as shown for ζ = 0.5 and ζ = −0.5 in the figure. In the case of oppositely charged boundaries (ζ < 0), it turns out that the PB and SC curves and hence the exact results roughly coincide and become less distinguishable. A reasonable explanation for this would be that for oppositely charged surfaces the counterions mostly feel the effect of the strong uniform external field provided by the surface charges, which acts similarly in the strong as well as the weak-coupling limit. Thus the mean-field and the strong-coupling approaches should converge. For small external field, as in the case of similarly charged surfaces, the mean-field theory depends more on the local counterion density whereas the strongly coupled counterions still feel mostly the external field. Thus the difference between the WC and the SC frameworks in the ζ < 0 and ζ > 0 cases. We also emphasize that the above discussion holds in rescaled representation as the pressures are plotted here in rescaled units; in actual units, Fig. 1 corresponds to different ranges of separation, L, for the WC and SC regimes as the Gouy-Chapman length, µ, is typically very different between the two limits, i.e., it would be small at high couplings and large at small couplings as may be realized, e.g., by changing the counterion valency at fixed surface charge densities and Bjerrum length.
Note that for oppositely charged boundaries, both the PB and the SC pressure can be attractive at large separations, and hence the pressure at any finite N can be attractive. It is also notable that for charged surfaces of opposite sign, the PB analysis in general performs much better than for the surfaces of equal sign and that the asymptotic SC value is approached more quickly.
In all cases considered here the theoretical and the simulated values of the interaction pressure converge for very small surface-surface separations. In fact, the SC and PB results coincide in the leading order as the rescaled distance, L, tends to zero, since both are dominated by the osmotic pressure of counterions. One thus finds which is of course nothing but the ideal-gas osmotic pressure of counterion confined between the two plates (i.e., p = N k B T /L in actual units).

VI. CONCLUSIONS
We have analyzed the statistical physics of an overall electroneutral system, of one-dimensional counterions confined between two charged surfaces. This can be seen as a simple model for charged lipid multilayers or stiff mica sheets neutralized by counterions and confined between external charged plates. The thermodynamics can be solved exactly and a widely varying behavior of effective interaction between the confining plates can be discerned.
Apart from performing exact statistical mechanics and MC simulations, we have also analyzed these systems within the mean-field (Poisson-Boltzmann) or weak-coupling approximation and within the strong-coupling approximation. These approximations are in fact independent of the dimensionality of the problem and it is therefore interesting to see how they compare with our exact results in one dimension. For a large number of particles it is found that the mean-field approximation works well. This is physically understandable as the counterion distribution can be reasonably assumed to have a continuous profile as predicted by the mean-field theory and also the inter-particle interaction in one dimension is long range. In the case of few counterions it works less well as the effect of correlations becomes important. In this regime the SC approximation will perform better as it takes into account properly the discreteness of the counterions charges. Indeed, as the strong-coupling expansion is basically a form of the virial expansion it is to be expected that it works better in the case of few counterions.
On the level of the approximate evaluations of the partition function, only the fluctuation contribution to the weakcoupling limit depends on the dimensionality of the problem because of fluctuations in the plane of the film in 3D that are not present in 1D for which the summation over the in-plane modes is absent. In particular, the fluctuation contribution to the interaction pressure scales with the distance L as p 2 ∼ −Ξ 1D /L in 1D and as p 2 ∼ −Ξ 3D ln L/L 3 in 3D and for same-sign surfaces (ζ ≥ 0) at sufficiently large separations. The fluctuational contribution to the pressure is always attractive, reflecting the fact that electrostatic correlations mediated by counterions always favor attraction between the charged boundaries. Since the coupling parameter in 1D, Ξ = Ξ 1D , depends inversely on the number of counterions, the fluctuation contribution to the interaction pressure that scales linearly with Ξ is in this case vanishingly small, a situation clearly confirmed by exact and MC evaluation of the partition function. We can solve for p 2 (L) exactly in 1D and 3D and evaluate the relative importance of the contributions from fluctuations and from mean field theory. For example, in 1D with N = 20, ζ = 1 (Ξ 1D = 0.1), we find p 2 /p 0 0.1 forL < 5; whereas in 3D we have p 2 /p 0 0.1 forL > 0.05 when Ξ 3D = 0.1, andL > 20 when Ξ 3D = 1.0. This confirms the results inferred from the asymptotic behaviour of p 2 (L).
In the case of symmetric surface charges the mean-field predictions always give a positive (repulsive) pressure between the two plates. However the SC and exact results do predict a possible attraction at large inter-plate separations. This effect is basically due to strong correlations arising when the counterions cannot neutralize both surface charges simultaneously and thus the residual surface charges are both attracted to the residual counterion charge in between the two plates. As the number of counterions increases this effect disappears and the system exhibits a more mean-field-like behavior.
As regards future work it would be interesting to see if the one-dimensional theory can be adapted to model a full three dimensional system. There is some hope that a proper reformulation of the statistical mechanics of this system, by isolating explicitly the dependence of the local fields in the transverse as opposed to the longitudinal direction with respect to the bounding surface normals, might open up new prospects for approximations that would work also in the gray area where the WC and the SC fail. This approach may work due to the fact that the mean-field and strong-coupling approximations are dimensionality independent. Whether or not it can be expected to work will depend on whether charge correlations in the direction perpendicular to the film dominate the physics or whether lateral correlations (such as the formation of locally two dimensional Wigner-crystal-like structures in the plane of the film) dominate.