Strong-coupling theory of counterions with hard cores between symmetrically charged walls

By a combination of Monte Carlo simulations and analytical calculations, we investigate the effective interactions between highly charged planar interfaces, neutralized by mobile counterions (salt-free system). While most previous analysis have focused on point-like counterions, we treat them as charged hard spheres. We thus work out the fate of like-charge attraction when steric effects are at work. The analytical approach partitions counterions in two sub-populations, one for each plate, and integrates out one sub-population to derive an effective Hamiltonian for the remaining one. The effective Hamiltonian features plaquette four-particle interactions, and it is worked out by computing a Gibbs-Bogoliubov inequality for the free energy. At the root of the treatment is the fact that under strong electrostatic coupling, the system of charges forms an ordered arrangement, that can be affected by steric interactions. Fluctuations around the reference positions are accounted for. To dominant order at high coupling, it is found that steric effects do not significantly affect the interplate effective pressure, apart at small distances where hard sphere overlap are unavoidable, and thus rule out configurations.


I. INTRODUCTION
The dominant part of colloids release micro-ions of low valence from the surfaces at deionized conditions [1][2][3]. These mobile so-called counterions can be regarded as identical classical particles interacting via the threedimensional 1/r Coulomb potential. The charged surface with the surrounding counterions form in thermal equilibrium a neutral electric double layer [4][5][6]. The geometry of two parallel similarly and uniformly charged walls at distance d with counterions in between provides the simplest setting for studying effective interactions between like-charged macromolecules. It was shown in early experiments [7][8][9][10][11], more recently in membranes, vesicle, or bilayer systems [12][13][14][15], as well as in numerical simulations [16][17][18][19], that like-charged colloid surfaces can attract each other, under the action of Coulombic forces alone. This requires that the coupling be strong enough. In the case of pointlike counterions, the relevant theory involves only one dimensionless thermodynamic parameter, namely the coupling constant Ξ [20]. However, for many systems, one cannot ignore the finite size of the ions. These includes for example systems with bulky counter-ions like ionic liquids [21,22], highly charged surfaces like calcium silicate hydrates where the size of the (hydrated) ions are comparable to the average distance between neighboring surface charges [23], systems with high salt concentrations [24] or with dielectric discontinuities [25]. It is the purpose of the present paper to go beyond the model of point counterions, by accounting for steric effects when these ions feature a finite size.
The weak-coupling (small Ξ) limit at low salt concentrations is well described by the Poisson-Boltzmann mean-field theory [2]. Within a field-theoretic representation of the Coulomb fluid grand-canonical partition function [26], the Poisson-Boltzmann theory is the leading term in a systematic loop-expansion [27]. To describe the opposite strong-coupling (SC) limit, a virial (fugacity) expansion of the grand-canonical partition function in inverse powers of the coupling constant Ξ was proposed [28][29][30][31]. In the case of a single charged surface and to leading virial SC order, each particle moves independently of other particles in the direction perpendicular to the confining surfaces, which was verified by Monte Carlo (MC) numerical simulations [28]. The first SC correction to the particle density profile [29,30] has the right functional form in space, but the wrong dependence of the prefactor on the small parameter 1/Ξ. As concerns the geometry of two parallel equivalently-charged walls, the analytical results for the pressure are accessible only for very small distances d.
Other theoretical attempts to construct a SC theory were based on the ground-state Wigner crystals created by counterions [32][33][34]. In the absence of dielectric wall images, according to Earnshaw's theorem [35] the counterions stick on the wall surfaces at infinite coupling (e.g. temperature goes to zero). For the one-wall geometry, they form a two-dimensional hexagonal (equilateral triangular) Wigner crystal. In the case of two parallel walls, five distinct (staggered) Wigner bilayers I-V were detected as the distance between the walls increases from zero to infinity [36][37][38][39][40][41][42]. Some elusive features of critical properties were revisited in Ref. [43] by using an analytic approach based on an expansion of the energy of the five structures in generalized Misra functions [44]. A SC theory based on the harmonic approximation for particle deviations from their ground-state Wigner positions was proposed in Ref. [45]. The leading order for the density profile and the pressure turns out to be identical to the virial single-particle theory. For the one-wall geometry, the first correction to the particle density profile has the correct functional form in space and the prefactor, pro-portional to 1/ √ Ξ, is in good agreement with MC data. As concerns the two-wall geometry, the harmonic analysis in Ref. [45] is also restricted to very small distances [46].
Taking the harmonic approximation in full [47], i.e. with no restriction to small distances, leads to an effective one-body potential acting on particles at each of the two walls which interpolates correctly between zero for distances d much smaller than the Wigner lattice spacing and the locally linear potential of separated charged walls at asymptotically large inter-wall distances d → ∞. This is why the profile of the particle density and the pressure exerted on the walls are described well also for intermediate distances between the walls comparable with the lattice spacing of the Wigner crystal. The technique was first applied to asymptotically large values of the coupling constant when the system is in a crystal phase. The aspect ratio of the Wigner bilayer structure (around which the harmonic expansion is made) was taken as a free parameter, determined by minimizing the total free energy. For treating the fluid phase present at smaller and more realistic values of the coupling constant Ξ, the Wigner bilayer structure was substituted by a correlation hole, i.e. a depletion region around each particle due to the strong Coulomb repulsion [32,[48][49][50][51][52][53]. It is noteworthy that the details of the structure, crystalline versus strongly modulated liquid, only affect fine details of the effective force, and are rather immaterial.
The aim of this paper is to extend the strong-coupling two-walls analysis of Ref. [47] beyond point-like counterions, treating this species as charged hard spheres of diameter d hc which are impenetrable to other particles as well as to hard walls (primitive model). Theoretical treatment of pure hard-core systems is usually based on a potential of mean force, the so-called depletion potential [54]. The phase diagram of hard spheres (without any charges) between parallel plates was calculated by using MC simulations in a wide range of particle densities and for plate separations ranging from one to two hard-core diameters in Ref. [55]. Besides the standard fluid phase, pure hard spheres freeze into closed-packed versions of the crystal bilayer structures which take place also in the pure Coulomb problem, namely one triangular layer (phase I), the linear buckling structure (phase II), two square layers (phase III), the rhombic structure (phase IV) and two triangular layers (phase V). The primitive model, including both Coulomb and hard-core interactions, was studied mainly numerically by using MC simulations, within an isolated electric double-layer [56][57][58] as well as two-wall geometry [59][60][61], for weak and intermediate values of the coupling constant, up to Ξ 100. Here, we shall assume that the coupling constant Ξ is large, so that the Coulomb interactions dominate in creating the ground state. For small and intermediate interplate distances, the particles are supposed to form basically the Wigner bilayer structure of type I, II or III, their centers being at distance d hc /2 from either plate 1 or 2. We shall look, both analytically and numerically, for steric hard-sphere effects on this structure.
The paper is organized as follows. A short recapitulation of the SC theory for pointlike particles, as developed in Ref. [47], is presented in Sec. II. Subsec. II A reviews the relevant ground-state Wigner bilayers while Subsec. II B deals with the leading SC description of thermodynamics and the density profiles. Sec. III generalizes the SC theory to account for ionic hard core. Subsec. III A brings a list of steric restrictions on the parameters of the Wigner bilayers. Subsec. III B deals with SC thermodynamics of hard spheres. The comparison of the theory with our Monte-Carlo numerical results is given in Sec. IV and Sec. V is for the Conclusion.

II. POINTLIKE PARTICLES
We start with the definition of the model where positions are denoted by r = (x, y, z). There are two parallel plates Σ 1 at z = 0 and Σ 2 at z = d with infinite surfaces |Σ 1 | = |Σ 2 | = S → ∞ along the (x, y) plane. The plates are charged symmetrically by the uniform surface charge density eσ with e being the elementary charge and σ > 0. For this case the resulting electric field vanishes between the plates.
There are N mobile particles between the plates, each with a charge −qe, coined as "counterions". The valency q takes integer values (e.g. q = 1 for Na + ions, q = 2 for Mg 2+ etc.) while e is the electron charge. At this stage we consider classical (i.e. non-quantum) particles to be pointlike. The electro-neutrality of the system is ensured by the equality N q = 2σS. (2.1) The dimensionless distance between the plates is defined as Technically speaking, it is convenient to have a rescaled measure of distance that is temperature independent. This avoids singularities when studying the ground state, that is met under infinite coupling, see below. The particles are immersed in a solution of dielectric constant , the same as that of the walls, so that there are no image forces at work. In Gaussian units, the Coulomb potential at distance r is given by 1/( r). The system of charged particles and plates is in thermal equilibrium.

A. Ground state
In the ground state, corresponding to infinite coupling (Ξ → ∞), our interacting point charges in a slab domain stick to the domain's boundary [35]. In the case of symmetrically charged plates, N/2 particles stick on plate Σ 1 and the remaining N/2 particles stick on plate FIG. 1. a) Geometry for the ground-state structures I, II and III of counterions on two equivalently charged plates, and definition of lattice vectors (a1, a2). Open and filled symbols correspond to particle positions on the opposite surfaces. The ratio |a2|/|a1| defines ∆. b) Side view, with definition of relevant distances. The dimensionless distance η between the plates is defined as d/ √ a1a2. We have d = D − d hc , where d hc is the ionic diameter and D the true distance between the walls; d turns out to be a more relevant quantity than D. Σ 2 . As η goes from 0 to ∞, numerical simulations [36][37][38][39][40][41][42] indicate five distinct bilayer Wigner structures. For small and intermediate values of η studied in this paper, the staggered rectangular structures I, II and III are relevant. As is shown in Fig. 1, a single layer consists in the rectangular lattice with the aspect ratio ∆, defined by the primitive translation vectors a 1 = a(1, 0) and a 2 = a(0, ∆). The lattice spacing a is determined by the electroneutrality requirement that the total surface charge in a rectangle must compensate the charge of just one particle, a = q/(σ∆). The identical rectangular structures on the two plates are shifted with respect to one another by a half period (a 1 + a 2 )/2.
Structure I with ∆ = √ 3 corresponds to a (equilateral) triangular lattice which appears in the monolayer limit η → 0. The aspect ratio is from the interval 1 < ∆ < √ 3 for soft structure II and ∆ = 1 for structure III which is the staggered square bilayer. The phase transformation I-II takes place just at η = 0 [41,43], the phase transition between structures II and III appears at η ∼ 0.263 and phase III provides the lowest energy up to η ∼ 0.621.
Using techniques introduced in Ref. [43], the energy per particle e 0 = E 0 /N is expressible for all three structures I-III in terms of the generalized Misra functions z ν (x, y) = This makes the use of symbolic calculation softwares very efficient. In practice, the infinite series (A4) over (j, k) indices must be truncated at some M . For the well known case of the hexagonal lattice with η = 0 and ∆ = √ 3, the truncation of the series at M = 1, 2, 3, 4 reproduces the Madelung constant up to 2, 5, 10, 17 decimal digits, respectively [43]. To maintain a high accuracy of our results, we truncate all Misra series at M = 6. The calculation of one ground-state energy value takes less than one second of CPU time on a standard PC.
For a given distance η, the value of the rectangular aspect ratio ∆ is determined by the energy minimization condition This condition sets the dependence of the aspect ratio on the dimensionless distance between the plates in the ground-state ∆ 0 (η), see Ref. [43].

B. Crystal phase at strong coupling
The system being in thermal equilibrium at some (inverse) temperature β = 1/(k B T ), there are two relevant length scales. The distance at which two elementary charges interact with thermal energy k B T is the Bjerrum length (2.6) A charge qe at distance z from a wall with the surface charge density eσ has the potential energy 2πe 2 qσz/ . The distance at which the charge qe has the potential energy equal to thermal energy k B T is known as the Gouy-Chapman length The coordinate z, which is perpendicular to the charged surfaces of the walls, will be often expressed in units of µ, z = z/µ. The dimensionless coupling parameter Ξ, quantifying the strength of electrostatic correlations, is defined as the ratio of the two relevant lengths: The strong-coupling (SC) regime Ξ 1 is in practice most conveniently met by increasing the valence q. In doing so, excluded volume effects become prevalent, and the point-like limit of early studies less relevant. Alternatively, the regime of strong coupling corresponds to either low temperatures (a limit that is of little practical interest in view of applications with water, due to the unavoidable freezing of the solvent), or large surface charge densities. The lattice spacing a of the Wigner structure, which is the characteristic length scale in the longitudinal (x, y) plane, is much larger than µ in the SC regime as a/µ ∝ √ Ξ. In the remainder, we take q = 1, without loss of generality, in order not to clutter formulas.
For a single-layer Wigner crystal, experiments [62] and simulations [63] give the estimate Ξ 3×10 4 for the coupling parameter at melting from the ordered crystal to a fluid phase. The coupling parameter at melting of the Wigner bilayer crystal depends on η [38]. Let Ξ be large enough to localize particles near their Wigner-crystal positions. Within the canonical ensemble, the relevant thermodynamic quantities are the partition function Z N and the corresponding (dimensionless) free energy per particle βf = βF/N which are defined, up to some irrelevant constants due to the interaction of surface charge densities with themselves and charged particles, as follows is the Coulomb interaction energy of the particles and λ stands for the thermal de Broglie wavelength. We recall that the electric potential induced by the symmetrically charged plates is constant between the plates. The mean particle number density at point r is defined as ρ(r) = N i=1 δ(r − r i ) , where · · · means the statistical average over the canonical ensemble. It fulfils the conservation condition ρ = N . We here study the (x, y)-averaged density profile ρ(z), which depends only on the perpendicular z-coordinate, ρ(r) ≡ ρ(z), so that With the rescaled particle number density The strong-coupling approach to the counterion system is based on a harmonic expansion of the energy E with respect to particle coordinates around their ground state Wigner bilayer positions [38], where the ground state corresponds to infinite coupling. Numerical simulations in Ref. [47] indicate that at finite although large coupling, the particles form another reference crystal of type I-III with the aspect-ratio parameter ∆ which depends, besides the inter-plate distance η, also on the coupling constant Ξ, i.e. ∆(Ξ, η). We have performed the full harmonic expansion of particle coordinates around this reference crystal and fixed ∆(Ξ, η) of the reference crystal by minimizing the free energy with respect to this parameter. In this paper, we keep only the leading terms linear in z; it turns out that the harmonic deviations in the crystal (x, y) plane as well as quadratic terms in the z-direction (proportional to 1/ √ Ξ) have only minor effects on the results in the SC regime. The neglect of these terms will enable us to include the hard-core interactions in a relatively simple way. The total energy is thus expressed as where the energy change is given by Here, the prefactor to small deviation terms is given by The leading terms are linear in z i for particles sitting in the ground state on plate Σ 1 and in ( d − z i ) for particles i ∈ Σ 2 . The function κ can be viewed as an effective electric one-body field due to the uniform surface charges on the two plates and the particle ground-state layer on the opposite plate. For η → 0 we have κ → 0, i.e. each particle feels the zero field coming from the uniform surface charges on the plates while the effect of the opposite particle layer with the lattice spacing a d is negligible. For η → ∞ we have κ → 1, i.e. each particle feels the field coming from the surface charge at its own plate while the discrete counterion structure on the opposite plate is smeared out and neutralized by the opposite surface charge density on that plate. The function κ thus reflects a continuous interpolation between the two-plate case for small η-values and the one-plate case for large η-values.
The partition function (2.9), with the particle interaction energy given by Eqs. (2.13) and (2.14), reads as Neglecting irrelevant terms which do not depend on η and ∆, the leading SC representation of the free energy per particle is given by The dependence of the aspect ratio ∆ on the coupling constant Ξ and the plate distance η, ∆(Ξ, η), is fixed by the principle of minimum free energy, i.e., This condition is the analogue of the infinite coupling relation (2.5).
The pressure can be obtained via the thermodynamic route as follows The pressure, rescaled as the particle density in (2.11), is given by To assess the consistency of the result, it is appreciable to have an alternative route for computing the pressure. It is offered by the contact theorem [64], that requires the knowledge of the contact ionic density. The particle density profile is derived in appendix B. The contact theorem for planar walls relates the total contact density of particles on the wall and the pressure: The thermodynamic P th and contact P c pressures in general do not coincide in an approximate theory, although they refer to the same quantity. Their difference reveals the accuracy of the approach. It should be kept in mind that in (2.22), the local field κ is distance dependent.

A. Steric restrictions
After having presented the key aspects of the theory for point ions, we now address hard-core effects: each ion is a hard-sphere of diameter d hc . The hard core is impenetrable to other particles (a model referred to as the primitive model) as well as the wall. We shall assume that the coupling constant Ξ is very large, so that the Coulomb interactions dominate and a simple crystal phase (I, II, or III as in the point case with only one ion per lattice cell) is formed, as long as it does not lead to ionic overlap. Scanning only these simple crystal phases was further motivated by viscual inspection of the structures found by our Monte Carlo simulations. The counterions are supposed to be close to the Coulomb bilayer structure of type I-III, their centers being at distance d hc /2 from either of plates 1 or 2 and we shall look for steric hard-sphere effects on this structure.
If D is the true distance between the walls, it is useful to define the reduced distance d via D = d + d hc , where d is the distance available to the center of mass of hardsphere ions; it is equal to 0 in the extreme case when particles touch by their hard-core surfaces simultaneously both plates, see Fig. 1. As above, we use the notation η = d √ σ. It is useful to express lengths in terms of the lattice spacing a b of the hexagonal Wigner bilayer at which compares Coulomb and steric effects in the system. Note that r = 1 when d hc √ σ = 3 −1/4 0.76. When r < 1, the main expectation goes as follows. If Ξ is sufficiently large, the ions strongly repel each other, so that their "in-plane" (xy) motion is essentially frozen: their only possible motion takes place perpendicularly to the plates, along z. It is consequently immaterial to consider point-ions, or hard-sphere ions, as long as r < 1. We then expect that when expressed in terms of the d variable, the pressure curves should be independent of the ionic diameter. This "no-hindrance regime" will be illustrated in section IV.
When r > 1, steric hindrance impinges on the point-like arrangement, and needs to be properly addressed. Supposing that the counterions form basically the Coulomb bilayer structure of type I-III, there are strong steric hard-sphere restrictions on model parameters which have both intra-layer and inter-layer nature. We start by the intra-layer analysis. The existence of the structures I-III is limited by the condition If r ∈ [1, 3 1/4 ], the formula (3.3) yields a restriction on the parameter ∆. For r > 3 1/4 1.316, the bilayer Wigner structures I-III cannot exist at all.
As concerns the inter-layer hard-core restrictions on structures I-III, there exists a minimal distance d min at which the two layers can approach one another. This distance is determined as the one at which two nearestneighbor hard-core particles from the opposite layers 3. The dashed line shows ∆ * defined in (3.6). It discriminates a region where the minimal distance between plates can vanish in the d variable (corresponding to an inter-plate distance equal to the ionic diameter, and thus a monolayer), from another where steric effects preclude this possibility and lead to a non-vanishing minimal distance d, as defined in Fig. 1. The interpretation of the r = 1 threshold is that it corresponds to the maximum hard core size compatible with a possible compaction of the system down to d = 0; the monolayer is then triangular (hexagonal), with ∆ = √ 3. The analysis is here restricted to monolayer or bilayer upon contact, discarding situations with more than three layers, that would be formed in the hatched region.
touch one another (see Fig. 2): Equivalently, For a fixed r, the right-hand-side (rhs) of this equation is a monotonously decreasing function of ∆ (1 < ∆). If  Fig. 4, is also reported. The dashed line corresponds to the ground state optimal configuration when r = 0, i.e. for point-like ions, as derived in [43].
, there is no inter-layer restriction on structures I-III. For r ≥ 1 it holds that η min ≥ 0 for an arbitrary value of ∆ ∈ [1, √ 3], i.e., there is always a hard-core restriction for distances between layers. For r ∈ [0.930605, 1], there is an interval of the aspect ratios ∆ ∈ [1, ∆ * (r)] with η min ≥ 0 and an interval of ∆ ∈ [∆ * (r), √ 3] with η min = 0, ∆ * (r) being given by  Fig. 5 lies in the dashed curve, that shows how the geometry of the ground state problem without hard-core (r = 0) depends on interplate distance η. This curve was obtained analytically in [43]. It lies, although marginally, in the forbidden region of the r = 1 case. Yet, it lies in the acceptable region with r = 0.99. This allows to state that starting from the optimal ground state configuration of point charges, and gradually increasing the radius of hard sphere ions, steric effects will not alter the point-like configuration for r < 0.99. They start to do so for r slightly above 0.99.
For the limiting case η min (∆, r) = 0, the two plates are allowed to touch one another (d = 0). The equality η min ( √ 3, r) = 0 is satisfied for r = 1 which is the threshold beyond which η min is positive. As soon as η min > 0, the pressure is infinite for all inter-plate distances η < η min , since the hard spheres cannot be packed in such a small space.
The crystal state of the counterion system now depends not only on the coupling constant Ξ but also on the r-parameter whose large value can decrease substantially the coupling constant Ξ at which the crystal-fluid phase transition occurs. In the crystal phase, we can treat the hard-core system basically in the same way as the pointlike one in Sec. II, to obtain the effective (dimensionless) potential −κ(η, ∆) z acting on particles at plate 1 and the symmetrically reflected one with respect to the slab center at z = d/2, −κ(η, ∆)( d − z), acting on particles at plate 2. Because of strong Coulomb repulsions in the (x, y) plane, particles move freely along the lines in the perpendicular z-direction defined basically by the ground-state structures I-III. Due to interlayer steric effects, the particles at plate 1 move in a reduced interval z ∈ [0, d − d min ] while those at plate 2 in the interval z ∈ [ d min , d]. In what follows, we shall use the following combination of variables: where d = √ 2πΞ η, with a similar relation between d min and η min .
Finally, previous works [45,47] have shown that even when the ionic system is not coupled enough to be in its crystal phase but rather exhibits a strongly modulated liquid structure, the large-Ξ calculations are nevertheless relevant as an approximate approach. The main reason is that both structures, liquid and solid, do exhibit the common feature of a correlation hole around each ion [30,53]. We thus here develop a theory that is grounded in the large-Ξ regime, the relevance of which at moderate couplings has to be assessed by a direct comparison to numerical simulations. is for a case where one of the 4 plaquette ions (the "upper ion") has moved away from its ground state position, which diminishes the available space for the tagged ion, again materialized by the slab between the two dashed lines. It is assumed here that since the white ion only moves perpendicularly to the plate, the "bottom" black ion does not contribute to the available space.

B. Thermodynamics
To account for steric effects on the Coulomb free energy in the leading SC order (2.18), let us take one of the particles at plate Σ 2 as the reference ion. It has just four nearest neighbors at the corners of one rectangular plaquette of the Wigner crystal at plate Σ 1 ; we denote these particles by 1,2,3,4 and their perpendicular positions respectively by z 1 , z 2 , z 3 , z 4 . Because the particles are supposed to move along the lines determined by Wigner layers in the perpendicular z-direction, from among four positions only the maximal one max( z 1 , z 2 , z 3 , z 4 ) is relevant. The original interval [ d min , d] accessible to the reference particle is thus reduced to [ d min + max( z 1 , z 2 , z 3 , z 4 ), d], see Fig. 6. The contribution of the reference particle on plate 2 to the partition function can be integrated out as follows where h = κ( d − d min ), is a rescaled measure of available space. Performing the above procedure independently for every of N/2 particles on plate 2, the N -particle partition function reduces to the one of N/2 particles at plate 1, where the product is over all i = 1, 2, . . . , N/2 plaquettes of the Wigner rectangular lattice at plate 1, the coordinates of particles localized at the four corners of plaquette i being denoted as z i1 , z i2 , z i3 and z i4 . We see that the elimination of one half of particles implies plaquette four-particle interactions among the remaining half of particles. Finally, making the substitution s i = κ z i one ends up with Through the original variables {s i }, the plaquettes are coupled, which makes the statistical mechanics problem at hand untractable. We shall treat the partition function Q z approximatively by using the Gibbs-Bogoliubov inequality [65]: − ln Tr e −H ≤ Tr (p 0 ln p 0 ) + Tr (p 0 H) , (3.11) where p 0 is any normalized probability distribution, Tr p 0 = 1. Comparing this formula with the studied case (3.10), we identify the (dimensionless) Hamiltonian and Let us choose where α is a free (real) parameter. The reason for this choice is dictated by the observation of ionic density profiles, see below, that appear essentially exponential. In other words, α plays the role of a multiplying factor to the local effective electric field, dressed by steric effects. In absence of hard-core interactions, one would have precisely α = 1, and the treatment of section II would apply. The fact that α = 1 (with an effective field ακ) will be a direct signature of hard-core interactions. Our choice of trial probability p 0 decouples the plaquette, as mean-field treatments do. Since we obtain that Consequently, the free energy of hard spheres with the Coulomb interaction satisfies the inequality βf (η, ∆; r) ≤ √ Ξ 2 3/2 π Σ(η, ∆) + ln κ(η, ∆) + 1 2 ln α − ln 1 − e −αh + 1 2 The free parameter α is chosen to minimize the upper bound for the free energy, i.e., the rhs of this equation. In all the cases studied, the obtained α is in the interval [1, ∞]. The parameter α thus increases the slope of the decay of the particle density from the wall surface; since the density profile is normalized this automatically means the increase of the particle density at the wall as the consequence of the hard-core repulsion from particles close to the opposite wall. This can be thought of as a generalized depletion effect, where ions are pushed to their nominal plate by hard core, that adds to the Coulomb repulsion already at work for point particles. As before, the aspect ratio of the rectangular lattice ∆ is also the minimizer of the free energy, respecting the hard-core restriction (3.3). The (rescaled) thermodynamic pressure is given by formula (2.21). We turn to the density profile. Its part due to the particles in the vicinity of plate 1, ρ 1 (z), is obtained in analogy with the case of pointlike particles by introducing the generating Boltzmann factor w(r), see Appendix B. This means that the dimensionless Hamiltonian (3.12) has to be substituted in the original expression for the whole partition function as follows Within the Gibbs-Bogoliubov formalism, ln w(r i ) appears in the evaluation of Tr(p 0 H) and therefore this is just the trial distribution (3.14) which determines the particle density where θ denotes the Heaviside step function. Given the factorized form taken in Eq. (3.14), this result does not come as a surprise. Although particles at plate 2 have been integrated out within our approach, their contribution to the density profile is determined by its reflection z → d − z symmetry as follows The total density of particles is given by In the analogous formula for pointlike particles (B5), the Coulombic effects are expressed through the function κ ∈ [0, 1] which is coupled to z and d− z in exponentials. Now there is an additional multiplication parameter α ∈ [1, ∞] which reflects the "squeezing" effect of the ionic hard core. The contact version of the pressure follows from the contact theorem (2.22). We have to distinguish between two cases. If d min ≤ 0, the particles from plate 2 can touch plate 1 and therefore If d min > 0, the particles from plate 2 cannot touch plate 1 and therefore

A. Monte-Carlo Simulations
To put to the test the analytic theory, we run Metropolis Monte Carlo simulations of the system composed of two symmetrically charged surfaces with counterions inbetween, at various coupling parameters, separations and hard core radii. For each simulation, we use 512 spherical counterions, which all have their charge located in the center of their hard core. The planar surfaces are modeled as uniformly charged structureless hard walls. Longranged electrostatic interactions are handled by threedimensional Ewald summation techniques with corrections for quasi-2-dimensionality, by adding vacuum slabs on each side of the charged walls (as described elsewhere [47,66,67]). We verified that our vacuum slabs were large enough in order not to influence our results (i.e., pressures and profiles), typically larger than a couple of thousands of Gouy-Chapman lengths defined in (2.7). Besides standard particle trial displacements, we also utilize floppy-box moves at constant box volume at which counterions are displaced conformally: either by shear or coupled biaxial compression-decompressions (where we compress one axis and decompress the other) both in the plane parallel to the surfaces. Trial move parameters were set such to have an acceptance ratio between roughly 20 and 50% for each case. Pressures and profiles at a fixed separation, a given counterion hard core radius, and a given coupling parameter are estimated by first equilibrating for 10 4 Monte Carlo cycles and then sampling over 10 5 subsequent cycles, where a cycle consists of either 512 trial counterion displacements or a trial floppy box move (a fifth of the total cycles). Pressures are evaluated both at the walls (contact theorem) or over the mid-plane by sampling the concentration (entropic contribution), ion-ion correlation (electrostatic energy), and hard-core repulsion (impulse) over the mid-plane [16]. Both measures give the same results within statistical errors. The mid-plane evaluation is usually less noisy and hence all simulation results are reported using this measure. We apply block averaging of ten blocks to estimate the precision in pressures. Starting configurations for our simulations are counterion bilayers of structure I if r < 1 otherwise structure II with ∆ equal to the upper bound of Eq. (3.3), compatible with the minimum separation. Figure 7 shows the numerical results of the equation of state for six different coupling parameters. The two first and lowest ones, Ξ = 0.175 and Ξ = 1.58, yield similar pressure curves. Note that the point-like limit provides a universal equation of state, independent of Ξ provided it is not too large (less than 2). This point-like limit is here in excellent agreement with Poisson-Boltzmann theory results (not shown). Beyond point ions, steric effects result in very similar pressure curves at the two lowest Ξ studied; these effects are responsible for the relevance of the a parameter (or equivalently a b ), for scal- ing out results. These two equations of states are repulsive, irrespective of the hard core radius and separation, with pressure curves increasingly repulsive when increasing hard core radii at a given separation. These low-Ξ results serve as a reference to our strong-coupling analysis, illuminating the importance of increasing electrostatic coupling. Two peaks appear at these low coupling parameters, one at η 0.4 when d hc √ σ = 0.7 (cyan symbols) and the other at η 1.5 when d hc √ σ = 1.
They are fingerprints of the pure hard core system, in this parameter range barely affected by the electric charges. For instance, the change of behaviour for η 0.4 and d hc √ σ = 0.7 is consistent with the confined hard-sphere phase diagram reported by Schmidt and Löwen [55]. Indeed, computing the dimensionless quantities used in [55], we get h 0.56 and ρ H 0.64, which corresponds to the onset of crystallisation, arriving from the fluid sector. Furthermore, for hard core radius d hc √ σ 0.4 (or equivalently, r 0.5), only minor differences are seen compared to the d hc = 0 situation in this low-coupling regime, where electrostatics can be described in a meanfield manner. We will see below when discussing the no-hindrance regime that this insensitivity is even more pronounced in the strong-coupling regime since it holds strictly for r < 1, and also in a sense to be specified for r > 1.
At Ξ = 17.5, and at short separations, one observes a shallow attraction between the two surfaces if the counterions radius is not too large, d hc < 5.25 (or d hc √ σ < 0.5).
Pressure curves start to be influenced by hard core radius as soon as d hc 4. Increasing counterions size makes pressure curves repulsive at all separations for d hc ≈ 5.25, but with local minima for r < 1 ( d hc < 7.98). The peak at η 0.5 seen previously for low couplings and d hc √ σ = 0.7 appears also for the lower d hc √ σ values (0.5 and 0.6) at Ξ = 158 to gradually disappear again at even higher coupling parameters. Ξ = 17.5 is special in the sense that pressures are close to zero around η = 0.5 for the point charge case and hence is sensitive for perturbations (e.g. introducing excluded volume) around this state. The local minimum seen for d hc √ σ = 1 for the low coupling limit persists up to Ξ = 17.5, but vanishes somewhere in the range Ξ = [17.5, 158]. Even though the pressures are in practice zero for the high coupling cases at η = 1.5 we still do not see any effect of the hard core radius (in contradiction to the Ξ = 17.5 case and states around η = 0.5). The relative influence of hard core vs electrostatic interactions is also decreasing with increasing coupling parameter. By increasing Ξ, one turns the r < 1 cases from repulsive at all separations to attractive, except in a narrow range close to zero separation, of extension given by the Gouy-Chapman length. Somewhere around Ξ = 1750, one also turns cases with r > 1 from purely repulsive to having an attractive pressure minimum. For our highest coupling parameter (Ξ = 175000), we see that the r ≥ 1 follows the d hc = 0 curve up to the closest separation for the corresponding r. This can be viewed as an extension to the sector r > 1 of the no-hindrance effect alluded to in section III: it indicates that under such strong couplings, the dominant effect is electrostatics, equivalent to that of point-like ions, while steric effects only matter through the forbidden overlaps. When no overlaps are involved, the Coulombic interactions are largely dominant. The case r = 1 does, however, have a smaller minimum in absolute value compared to the d hc = 0 case even though the closest approaches are the same. This can theoretical be understood as some of the preferred bilayer structures are forbidden due to hard core overlaps (see Fig. 5), leading to a slightly altered and weaker (in terms of attraction) pressure curve.
We come back to this in the next subsection.

B. Comparison with analytic results
We have argued in section III that for d hc < a b , all pressure curves P (d) should collapse onto their pointion limit, provided the coupling parameter Ξ be large enough. This no-hindrance regime is illustrated in Fig.  8. In this figure, the largest value of r reported (respectively 0.92, 0.92 and 0.97 for panels a), b) and c)) is fairly close to 1. Yet, at the largest Ξ, this has no visible effect on the pressure curve, while the quality of the data collapse is altered when decreasing coupling Ξ. In the corresponding equation of state the increasing branch on the right hand side is actually universal, independent of coupling Ξ, when expressed in the proper variable, here η [68,69], as revealed in Fig. 9. The reason behind this universality is that the behaviour is ruled by the infinite coupling attractor of point ions (with a divergent Ξ). The point-ion ground state pressure is indeed shown by the dashed line in Fig. 9. Besides, as hinted at in subsection IV A, the no-hindrance effect extends to case with r > 1, see Fig. 10. For d hc > a b , i.e. r > 1, it is seen that starting from large distances, the pressure curve follows the point counterion limiting curve, down to the smallest distance allowed by non overlap of hard cores. The minimal distance given by Eq. (3.5), matches very well that where P tot diverges in Fig. 10: for the data with r 1.05, we have η min ( √ 3/r 2 , r) 0.29 (see also Fig. 4), while for r 1.18, we get η min ( √ 3/r 2 , r) 0.54. Steric effects are here dichotomic: there are essentially irrelevant due to the strong Coulombic repulsion, or divergent at small d, for those configurations which are not allowed. Quite remarkably, the marginal situation with r = 1 remains close to the point-ion attractor, down to vanishing distances where the two double layers on the opposite walls exactly register.
Before entering into a more precise comparison between analytical and numerical pressure profiles, we test the ansatz underlying our choice of trial exponentialtype density p 0 in Eq. (3.14), by showing the ionic profiles in Fig. 11. The first observation is that they are neatly exponential in the vicinity of the plates, in line with our premises for the trial variational form for p 0 . The profiles significantly deviate from their mean-field, Poisson-Boltzmann counterpart, shown with the dashed lines; on the scale of the figure, the different dashed lines corresponding to different separations are quasisuperimposed, but would depart at larger distances. We have predicted in section III that steric effects make the ionic profiles more peaked at the plates than the equivalent point-ion system, having the same available free space for center-of-mass displacement. This is corroborated by the MC data. This steric-driven enhanced localization is quantified, in the theory, by the parameter α appearing in (3.14) and (3.20). The predicted behaviour of α is shown in Fig. 12-left. When the plates are far away, α → 1, signalling that steric effect do not affect the local field κ, itself distance dependent, that maintains ions in the vicinity of a given plate. On the other hand, decreasing η, it is seen that α increases quite significantly. Besides, it is the product ακ that defines the local field; κ decreases for decreasing η while α shows the opposite trend, and we show in the inset of Fig. 12left the resulting effects for the product, compared to the MC measures. These MC results are extracted from the slopes evidenced in Fig. 11. While the predicted trend seems correct, it is seen that the theory leads to too sharp of a dependence on the distance, while MC yields smoother curves. A similar comment applies to the aspect ratio parameter, displayed in Fig. 12-right. Our variational treatment captures the correct trend, but exaggerates the sharpness of the crossover. We note nevertheless that the agreement between MC and the prediction improves, expectedly, when increasing Ξ. It can be noted that the distance range where the predicted ∆ underestimates the measured one is precisely where the localisation parameter α in Fig. 12-left displays a nonmonotonous local bump. It is also worthwhile here to inspect directly structural features. Fig. 13 shows the projected instantaneous position of ions, which reveals that the arrangement is of type considered in the theo-  [43,68,69]). Indeed, the dashed line displays the corresponding force per unit surface. Symbols and colors as in Fig. 8 with open markers for Ξ = 158, half-filled for Ξ = 1750, and full for Ξ = 175000. retical analysis, with rectangular unit cells. From left to right, the aspect ratio ∆ decreases from 1.53 to a value close to 1, as also shown in Fig. 12-right. This is confirmed by the computation of the in-plane pair correlation function, as displayed in Fig. 14, which also illustrates the relevance of Ξ to maintain in-plane order at large separation. from the analytical theory and numerical Monte Carlo simulations at coupling parameters Ξ = 158, 1750, 175 000 and for various hard core radii of the counterions. We focus here on the r > 1 cases. For Ξ = 158, the analytic theory seems to underestimate the repulsive pressure due to the hard core (or equivalent overestimate the electrostatic attraction). Similar trend is seen at Ξ = 1750 even though the theory yields pressures closer to the numerical results. For Ξ = 175000, we find a good agreement between numerical results and theory. At d hc √ σ = 1 ( d hc = 1050, i.e. r = 3 1/4 ), the numerical Monte Carlo data does however exhibit a significant level of noise. Besides, we have proposed two routes to compute the pressures, a mechanical and a thermodynamical one. Within an exact treatment, both results should coincide. The fact that they yield relatively close results in Fig. 15 assesses the self-consistency of the approach proposed.
To summarize, the agreement between our analytical calculations and the Monte Carlo results is good at the highest coupling studied, where steric effects either do not alter the point-ion pressure, or forbid too close interplate distances, for which the pressure is infinite. Steric effects are thus here dichotomic, all or nothing. At smaller couplings, a crossover sets in, where hard core have a finite and non negligible contribution to the pressure, that we capture semi-quantitatively, see Fig. 15. Difference between the theory and simulation data are due to the mean-field nature of the variational prediction performed.

V. CONCLUSION
We have derived an analytic strong-coupling theory for two like-charged plates, treating counterions as charged hard spheres of diameter d hc . Coulombic coupling is measured by a parameter Ξ that weights electrostatic effects against thermal energy. Starting from point charges and increasing d hc , a regime appears where it is no longer possible to accommodate a layer of counterions between the plates, but a bilayer forms. At this point (corresponding to a ratio of d hc over lattice spacing r = 1), the distant plates nevertheless can accommodate a monolayer of counterions, under larger coulombic couplings. We did not treat the cases of still larger values of d hc , where steric repulsion would lead to more complex arrangements (multilayers), in particular for r > 3 1/4 . Our approach starts from crystalline configuration of counterions that form at large Ξ. These crystals, which are staggered between the two plates, have been assumed to have a rectangular unit cell, which, given the staggering, includes the triangular lattice (often referred to as "hexagonal") that forms at close contact when a monolayer is admissible [70]. By integrating out all ions in the vicinity of plate 2, we obtain a non-trivial effective Hamiltonian ruling the behaviour of ions in the vicinity of plate 1. For the sake of tractability, an upper bound to the corresponding free energy is computed in the Gibbs-Bogoliubov spirit, considering a family of factorized probability distribution for the ions positions that involve a localisation parameter α. By minimizing this bound with respect to α and the lattice aspect ratio, we obtain explicit density profiles and pressures.
Our predictions have been compared to Monte Carlo simulations. At the largest Ξ, the results show a remarkable insensitivity to hard-core diameter, not only for r < 1 but also above 1, provided one works with the shifted distance between the plates (D → D − d hc ).
Only in the small distance range that is ruled out due to unavoidable ionic overlaps is the point-like pressure inapplicable. Decreasing Ξ, packing effects prove more relevant and have a non trivial signature on the equation of state, that our theory captures in a semi-quantitative way (see Fig. 15). We also found numerically that steric effects can completely suppress like-charge attraction within the primitive model, although less efficiently when Ξ increases. Indeed, the larger the coupling, the more ions do repel, and the less relevant their hard core becomes. We believe that this effect (e.g. suppression of like-like attraction) still persists in the low salt concentration cases and that steric effects in general will become more important as the salt concentration is further increased (i.e. increased repulsion in the pressure curves), as has been seen for the Debye screening length [24].
This work paves the way towards a more satisfactory and realistic description of strong-coupling theory, beyond the point ion limit for which it was initially devised. Among interesting perspectives, we mention the study of larger r values, when the ionic diameter would lead to multi-layers at close packing, when the two plates are at closest separation or the effect of the size and structure of the solvent itself. It would also be relevant to address asymmetrically charged walls, together with systems with salt, when micro-ions of both signs are present, not only counterions. In this respect, a promising approach is to extend the analysis of [74], where ions of opposite charges form Bjerrum pairs, i.e. neutral entities that may be in a first approximation discarded from the analysis. This yields an effective salt-free system, as addressed in the present work. Besides, in severely confined configurations, the modification of the solvent (say water) dielectric constant should also be included in the description [13,71].
We would like to thank I. Palaia and J. Zelko for useful discussions. The support of L.Š received from the project  .14) and (3.20), with distance η for Ξ = 1750. From left to right, the curves correspond to d hc = 84, d hc = 94.5 and d hc = 105. The product ακ, shown in the inset, is the local electric field felt by an ion in the vicinity of a given plate. This quantity embodies the influence of steric effects on this field, which does rule the ionic profiles, see Eqs. (3.19) and (3.20). In the inset, the symbols are for the MC measures, as extracted from data like those presented in Fig. 11  Colors are assigned depending on the closest wall for each particle. White dashed lines indicate the main simulation box and black lines correspond to a Voronoi construction for the red particles (say those on plate 1, with neighbors on plate 1). Size of the particles corresponds to the hard-core. Notice that structures become liquid above η 0.6. EXSES APVV-16-0186 and VEGA Grant 2/0003/18 is aknowledged.