Effect of solute size and polarizability on the equilibrium of like-charged macromolecules: interplay between enhanced repulsion and depletion attraction

The complex internal structure of the solute molecules involved in biomolecular processes necessitates the characterization of the forces between charged macromolecules beyond the point-ion description. In this article, we develop a functional integral theory of like-charged membrane interactions in liquids composed of structured ions modeled as charged Drude oscillators. The electrostatic mean-field (MF) regime of the system is characterized by an integro-differential Poisson-Boltzmann (PB) equation and a generalized contact value identity. The approximate analytical and exact numerical evaluation of these equations indicates that the average solute size and polarizability give rise to two distinct membrane coupling regimes. In the ideal gas regime of short intermembrane distances and dilute salt, the finite solute size induces a twofold enhancement of the PB-level repulsion mediated by point-like ions. At intermediate separation distances, the short-range enhanced repulsion is followed by a regime of depletion attraction between the like-charged membranes. This radical alteration of the intermembrane interactions by the intramolecular salt structure identifies the solute specificity as a key ingredient of the thermodynamic stability in colloidal systems.


I. INTRODUCTION
The internal structure of the microscopic particles invisible to the naked eye is often discernible by the macroscopic behavior mediated by them [1]. The salt specificity of protein stability [2] and the surface tension of electrolyte solutions [3,4] obeying the Hofmeister series are notorious examples of this upward causation. The largescale macromolecular aggregation induced by depletion interactions is another illustration of the macroscopic behavior regulated by intramolecular solute structure. The depletion forces induced by the excluded volume of the solute particles are frequently used in industrial applications requiring phase separation [5,6], such as paint making [7] and water purification by flocculation [8].
The first formulation of depletion interactions has been introduced by Asakura and Oosawa in 1954 [9]. The corresponding Asakura-Oosawa (A-O) theory predicts that if two uncharged macromolecules suspended in a solution of smaller particles approach each other by a distance shorter than the size of the particles, the sterically induced particle exclusion from the intermolecular region leads to an inward osmotic pressure gradient driving the macromolecules towards each other. Subsequent works extended the A-O model to interacting charged macromolecules immersed in a solution of spherical particles [10], formulated the problem within the framework of the density functional theory [11], and investigated the depletion interactions induced by polymeric depletants [12,13]. The depletion attraction effect has been equally studied experimentally via atomic force microscope and light scattering techniques [14][15][16][17]. The his- * email: buyukdagli@fen.bilkent.edu.tr tory of the discovery and the current applications of this entropic mechanism are elaborated in Refs. [5][6][7].
The electromagnetic Casimir interactions between two plates located in a vacuum environment bear some qualitative similarities with the aforementioned depletion interactions. In a parallel plate configuration, the exclusion of the long-wavelength modes from the interplate zone leads to the mutual attraction of the plates even in the absence of any excess charge on their surface [18][19][20][21]. The thermal van der Waals (vdW) interactions occurring between solvated macromolecules correspond precisely to the classical counterpart of the Casimir interactions [22][23][24][25][26][27]. The interplay between the short-ranged vdW attraction and the longer-ranged double layer repulsion [28,29] is the key ingredient of the Derjaguin-Landau-Verywey-Overbeek (DLVO) theory describing the electrostatic equilibrium of macromolecules in salty solutions.
The primitive model at the basis of the DLVO formalism is limited by its implicit solvent framework and the underlying point-ion assumption. Explicit solvent structure and ionic polarizability have been introduced into the primitive model by point-dipole approaches [31][32][33][34]. In order to relax the point-dipole approximation of these formalisms, we have previously developed a nonlocal electrostatic theory of finite-size solvent molecules and polarizable ions in contact with single interfaces [35]. Within this non-local PB (NLPB) framework, we also characterized the hydration-induced alteration of ion polarizabilities in bulk ionic liquids [36]. Recently, we reconsidered this model in the limit of point-ions with vanishing polarizability to investigate the effect of explicit solvent on macromolecular interactions [37]. Additional works focused on the electrostatics of rigid particles to explain the ordering transition in bulk liquids [38], and the bridging attraction between like-charged macromolecules [39][40][41][42][43].

arXiv:2209.11818v1 [cond-mat.soft] 23 Sep 2022
In this article, we characterize the effect of the extended ion structure on like-charge macromolecular interactions. To this aim, we develop an implicit solvent theory of structured electrolytes in contact with charged macromolecules. In our formalism, the primitive ion model is extended by the incorporation of the average solute size and polarizability. The schematic depiction of this electrolyte model and the inner solute structure are displayed in Figs. 1(a)-(b). The novelty of our formalism with respect to the aforementioned works is the full treatment of the nuclear and electronic charge separation without any multipolar expansion. The resulting alteration of the interfacial solute partition leads to two competing effects on macromolecular interactions.
Our manuscript is organized as follows. In Section II, we derive the field theory of structured liquids composed of finite-size solute molecules modeled as charged Drude oscillators [44]. In the present work, we limit ourselves to the electrostatic MF regime of weak macromolecular charges, and consider the specific case where the macromolecules consists of two negatively charged membranes. In the corresponding configuration, the intermembrane interactions are characterized by a one dimensional integro-differential PB equation and a generalized contact value identity. In Section III, via the approximate analytical and exact numerical evaluation of these equations, we show that the average solute size and polarizability give rise to a twofold enhancement of the shortrange PB-level repulsion, followed by a longer-range depletion attraction between the like-charged membranes. The corresponding pressure profile is schematically represented in Fig. 1(c). We also carry out a detailed investigation of the physical mechanisms driving these effects. In Conclusions, we summarize our findings, and discuss the limitations and potential extensions of our formalism.

A. Field-theoretic partition function
The geometry of the inhomogeneous liquid is displayed in Fig. 1(a). The system is composed of two anionic membranes of surface charge density −σ < 0. The membranes of separation distance d are immersed in an electrolyte solution. The electrolyte is composed of an implicit solvent modeled as a dielectric continuum of relative permittivity ε w = 78, and s species of polarizable molecules. The liquid confined between the membrane walls at z = 0 and z = d is in chemical equilibrium with a bulk reservoir of the same species. The entire system is at the ambient temperature T = 300 K. In the liquid, the solute molecules of the species i have fugacity Λ i and bulk concentration ρ ib . The nuclear and electronic charges of the polarizable ions illustrated in Fig. 1(b) are n i > 0 and −e i < 0, respectively. This corresponds to the net ion valency q i = n i − e i .
The electronic cloud fluctuations are taken into ac- count with a Drude oscillator model characterized by the harmonic potential with the inverse thermal energy β = 1/(k B T ) including the Boltzmann constant k B [44]. The electronic cloud radius b fluctuates around its isotropic average a i corresponding to the permanent dipole length. The variance α i of these fluctuations set the polarizability of the molecules [46]. The grand-canonical partition function of the polarizable charges reads where r ij and b ij stand for the center-of-mass position corresponding to the location of the nuclear charge, and the distance vector between the nuclear and electronic charges of the molecule j in the species i, respectively. Moreover, the factor c i normalizing the probability distribution characterizing the electronic cloud fluctuations is defined as where we used the error function erf(x) [47]. One notes that in the limit of a vanishing average size (a i → 0), the coefficient (3) tends to unity (c i → 1).
In Eq. (2), the interaction energy between the charges has been defined as (4) The first term of Eq. (4) corresponds to pairwise Coulombic interactions between the charges of total density where σ s (r) is the fixed charge density of the macromolecules corresponding here to the interacting membranes. The Coulomb potential is defined in terms of its inverse as v −1 c (r, r ) = −(k B T /e 2 )∇ · ε(r)∇δ(r − r ), where e is the electron charge, and ε(r) stands for the dielectric permittivity profile. In bulk liquids characterized by a uniform dielectric permittivity ε(r) = ε 0 ε w , with the vacuum permittivity ε 0 , the Coulomb potential reads v c (r) = B /r, where B = e 2 /(4πε 0 ε w k B T ) ≈ 7Å is the Bjerrum length. Finally, the second term in Eq. (4) takes into account the steric potential w i (r, b) acting on the polarizable solute molecules. This potential will allow to derive the ion densities and to impose the confinement of the molecules to the intermembrane region.
Applying to Eq. (2) a standard Hubbard-Stratonovich transformation, the grand-canonical partition function takes the form of a functional integral over the fluctuating electrostatic potential φ(r), i.e. Z G =´Dφ e −βH [φ] , with the Hamiltonian functional In Eq. (6), the first term corresponds to the electrostatic free energy of the implicit solvent, the second term is the contribution from the fixed macromolecular charge sources, and the third term takes into account the presence of the mobile polarizable charges.

B. Electrostatic equation of state
In the present study, we limit ourselves to the electrostatic MF regime of low-valency ions and weak membrane charges. At the corresponding MF-level, the equation of state satisfied by the average potential profile follows from the saddle-point condition δ(βH)/δφ(r) = 0. Passing from the complex to the real potential via the transformation φ(r) → iφ(r), one obtains an extended PB equation of the form k B T e 2 ∇ · ε(r)∇φ(r) + σ s (r) (7) with the potential energies experienced by the nuclear (n) and electronic (e) charges on the polarizable ions, respectively. From now on, we will refer to Eq. (7) as the polarizability-augmented PB (PPB) equation.
In order to derive the number densities, we specify the steric ion potential as w i (r, b) = w i,n (r) + w i,e (r + b), where the first and second terms on the r.h.s. are the steric potentials acting on the positive nuclear and negative electronic charges, respectively. At the MF-level, the number densities of these charges follow from the thermodynamic identities ρ i,n (r) = −δ ln Z G /δw i,n (r) and ρ i,e (r) = −δ ln Z G /δw i,e (r). This yields Hence, the PPB Eq. (6) takes the form of an integrodifferential Poisson equation, [n i ρ i,n (r) − e i ρ i,e (r)] = 0.
(12) In the bulk reservoir located far from charge sources and system boundaries, the average potential φ(r), and the steric potentials w i,n (r) and w i,e (r) vanish. Hence, taking the bulk limit of Eqs. (10)- (11), and carrying out the remaining gaussian integrals over the electronic cloud radius, one finds that the fugacity of each molecular species equals its bulk concentration, i.e. λ i = ρ ib .
At this point, we take into account the planar distribution of the fixed membrane charges, i.e.
Then, we assume the uniformity of the dielectric permittivity, i.e. ε(r) = ε w /ε 0 . Moreover, in order to account for the confinement of the charges to the intermembrane region, we choose their steric potentials as where θ(z) stands for the Heaviside step function. The resulting translational invariance along the membrane surfaces implies the presence of an exclusively perpendicular average field, i.e. φ(r) = φ(z). Thus, carrying out the integral in Eqs. (10)-(11) over the transverse component b = b xûx + b yûy of the electronic radius b = b + b zûz , and imposing the constraint (14) to the gaussian integral over the variable b z , the PPB Eq. (12) finally becomes with the nuclear and electronic charge densities In Eqs. (16)- (17), the conditional probability for the electronic cloud fluctuations reads (18) In the case of a vanishing polarizability, the rigid rod limit of Refs. [40,41] follows from Eq. (18) as Then, if the fixed dipole moment vanishes, the conditional probability (18) reduces to a gaussian probability distribution, i.e.
Integrating Eq. (15) in the vicinity of the membrane surfaces at z = 0 and z = d, the boundary conditions required for the solution of the PPB Eq. follow in terms of the Gouy-Chapman length µ = 1/(2π B σ) as Finally, in the bulk limit where the steric and electrostatic potentials in Eq. (15) vanish, the electroneutrality condition follows directly as C. Intermembrane pressure from the polarizability-augmented contact value theorem Within the framework of the PB formalism for point charges, the contact value identity relating the surface density of the ions and the pressure P (PB) net is given by where ρ i (z) = ρ ib e −qiφ (PB) (z) is the ion density, and φ (PB) (z) is the potential solving the PB equation [1]. In this part, we derive a generalized contact value theorem extending Eq. (23) to finite-size polarizable ions. First, we note that the MF grand potential of the system corresponds to the Hamiltonian (6) evaluated with the solution of the PPB Eq. (7). Thus, the dimensionless grand potential per surface defined as f = βH/S reads The membrane pressure follows from the variation of the grand potential density (24) with respect to the intermembrane distance d, i.e. βP = −δf /δd. This yields where we defined the integral J i =´d 0 dz∂ d ρ i,n (z). First, plugging into this integral the nuclear charge density in Eq. (16), using the Leibniz integral rule, and introducing the change of variable z → d + b z , one finds that the integral equals the contact density of the electronic charge, i.e. J i = ρ i,e (d). Then, we note that as the average potential satisfies the PPB Eq. (7) obtained from the minimization of the grand potential density (24), the fifth term on the r.h.s. of Eq. (25) vanishes. Taking into account these simplifications in Eq. (25), and subtracting the bulk osmotic pressure βP b = s i=1 ρ ib acting on the outer boundary of the interacting membranes, the net intermembrane pressure P net = P − P b follows as, The comparison of the generalized contact value theorem (26) with its PB counterpart (23) indicates that the former involves separately the contact density of the electronic and nuclear charges located on the polarizable species. In Section III B, we show that this feature has a major impact on the magnitude of the short-range intermembrane repulsion.

D. Calculation of the nematic order parameter
The orientational tendency of the polarizable molecules determining the shape of their electronic cloud can be characterized by the nematic order parameter where θ stands for the angle between the z axis and the dipolar length b. Vanishingly small values of this order parameter indicate random orientations of the dipolar vector b. Then, the parallel and perpendicular orientational tendencies of the ions with respect to the membrane walls correspond to the intervals −1/2 ≤ S i (z) < 0 and 0 < S i (z) ≤ 1, respectively. In Eq. (27), the orientational average is defined with the density function of the nuclear charge in Eq. (16) as In order to reduce the dimension of the integrals D i and N i located in the numerator and the denominator of Eq. (28), we will evaluate the order parameter (27) for vanishing average size. Thus, we set here a i = 0. First, we take into account the geometric identity y the radius component parallel with the membrane surface. Next, we use the corresponding cylindrical coordinate system to express the integrals D i and N i over the fluctuations of the electronic cloud, account for the planar symmetry of the system implying φ(r) = φ(z), and impose the steric constraints in Eq. (14). One obtains Evaluating in Eq. (30) the integral over the variable b , the average in Eq. (28) finally reduces to where we used the exponential integral function Ei(x) [47].

E. Solution schemes
1. Monopolar approximation to the PPB Eq. (15) The evaluation of the intermembrane pressure (26) requires the solution of the integro-differential PPB Eq. (15). Here, we introduce an approximate solution that will be used in Sec. II E 2 for the exact numerical solution of Eq. (15). This consists of the expansion of the potential functions in the exponentials of Eq. (16)- (17) in terms of the dipolar length b z . At the monopolar order O b 0 z , the charge densities become where we defined the monopolar potential φ 0 (z), and the steric exclusion function . Evaluating the latter integral over the electronic cloud fluctuations, the exclusion function takes the closed form Plugging the ion densities in Eq. (32) into Eq. (15), the monopolar approximation to the PPB equation emerges as a generalized PB equation with a non-uniform screening parameter, with the DH screening parameter per species The monopolar potential profile φ 0 (z) satisfying the ODE (34) can be numerically obtained via standard algorithms such as the 4th order Runge-Kutta method. In Fig. 6(c), we show that the monopolar potential φ 0 (z) (squares) is a highly accurate approximation of the potential φ(z) satisfying the PPB Eq. (15) (solid curves). This indicates that the multipolar moments of the polarizable charges bring a negligible contribution to the net MF potential. Hence, the major impact of the molecular polarizability on the MF electrostatics is the steric depletion effect embodied by the exclusion function (33).
2. Recursive solution of the PPB Eq. (15) In this part, we explain the recursive solution of the PPB Eq. (15) using the numerical solution of Eq. (34). First, we recast the PPB Eq. (15) such that its l.h.s. has the form of Eq. (34), where we omitted the BCs on the r.h.s. of Eq. (34), and introduced the corrective source function −n i k i,n (z) + e i k i,e (z)} .
In Eq. (37), the charge partition functions are defined as In this work, Eq. (36) was solved iteratively via a standard Runge-Kutta algorithm. This iterative solution is based on the treatment of the source function W [φ(z); z] as an external function. This has been achieved by including in Eq. (36) the corresponding source function evaluated with the potential profile of the preceding iterative step. To this aim, we modify Eq. (36) as where the index α denotes the number of the iterative step. The recursive cycle was initialized at α = 1 by solving first Eq. (34) to obtain the potential profile φ 0 (z). Then, this monopolar potential was used in the source term of Eq. (38) to obtain from the Runge-Kutta algorithm the updated potential φ 1 (z). At the next iterative step α = 2, Eq. (38) was solved by using the updated potential φ 1 (z) as the new input potential on the r.h.s., and this cycle was continued until numerical convergence is achieved. For the model parameters used in this work, the corresponding scheme enabled the rapid convergence of the solution cycle in less than five iterative steps.

III. RESULTS
In this part, we characterize the effect of solute size and polarizability on macromolecular interactions. In the remainder, we consider a symmetric 1:1 electrolyte and take s = 2. Thus, the valencies of the nuclear and electronic charges are set to n − = 1, e − = −2, n + = 2, and e + = −1, where the indices i = + and i = − refer to the cationic and anionic species, respectively. Together with the electroneutrality condition (22), this yields q ± = ±1 and ρ ±b = ρ b . Finally, with the aim to simplify the analysis, both ionic species are assumed to have the same average size and polarizability, i.e. a i = a and α ± = α.
A. Neutral membrane interactions

Ion-size induced depletion attraction between neutral membranes
In the case of neutral membranes with vanishing surface charge density σ = 0, the average potential satisfying the PPB Eq. (15) vanishes, i.e. φ(z) = 0. As a result, the charge densities in Eqs. (16)-(17) simplify as where the steric exclusion coefficient γ i (z) is given by Eq. (33). Consequently, the pressure (26) becomes where erfc(x) is the complementary error function [47].
The pressure (40) is illustrated in Fig. (2)(a). One sees that the finite extension of the electronic cloud induces a purely attractive interaction force. At vanishing polarizability α → 0 corresponding to rigid molecules (black curve), the range of the linearly decaying pressure is set solely by the average solute size a. Then, in the presence of a finite polarizability α a, the intermembrane pressure acquiring an enhanced magnitude and range decays over the characteristic separation distance d ∼ 2α. Finally, in the specific case of a finite polarizability and a vanishing average size, the pressure (40) simplifies to At short separation distances, the interaction force (40) decays linearly for a < d α as and for α d < a as Then, at large distances a < α d, the net pressure decreases according to a gaussian law, (45) In Fig. 2(a), the asymptotic laws (43) and (45) are displayed by the dashed red curves. These equations indicate that the presence of an average ion size a enhances the characteristic range of the interaction force (42) set by the molecular polarizability α.

Depletion mechanism behind the membrane attraction
In Fig. 2(b), the mechanism driving the attraction between the neutral membranes is illustrated in terms of the terminal charge density in Eq. (39). In the simplest case of purely rigid ions with vanishing polarizability α = 0 (black curve), the terminal charge density reducing to exhibits interfacial solute exclusion solely over the distance a. However, the plot shows that in the presence of a substantial polarizability α d/2, the entropic constraint limiting the variations of the electronic cloud radius to the intermembrane distance causes the partial exclusion of the solute from the entire intermembrane region, i.e. α ↑ ρ i (z) ↓. As a result, the contact density of the terminal charges drops below the bulk ion density, i.e. γ i (d) < 1/2. This point can be analytically shown for vanishing ion size where the density (33) reduces to From Eq. (47), the solute contact density follows as Thus, in the presence of a finite polarizability α, due to the imbalance of the contact ion densities between the reservoir and intermembrane sides, the outer osmotic pressure on the membrane walls exceeds the inner pressure, driving the net force (26) into the attraction regime of Eq. (42). We found that the contact density of the molecules drops significantly in the regime d 2α where the r.h.s. coincides nearly with the variance of the electronic cloud fluctuations. This feature can be equally illustrated via the pore-averaged ion density t i =´d 0 dzγ i (z)/d given by Fig. 2(c) shows that upon the increase of the polarizability into the regime α d/2, the partition function t i (solid curves) decays rapidly from the rigid rod limit (dotted lines) towards vanishingly small values. This confirms that in Fig. 2(a), the range of the pressure is set by the critical intermembrane distance below which the electronic cloud overlaps with the membrane walls.

Distortion of the electronic cloud by confinement
Before investigating the case of charged membranes, we characterize the effect of ion confinement on the shape of the electronic cloud. For vanishing potential φ(z) = 0, the evaluation of the nematic order parameter (31) yields where we used the exponential integral function E n=1 (x) [47]. Fig. 2(d) shows that as one moves from the regime α d/2 of large intermembrane distances or weakly polarizable ions to the opposite regime α d/2 of strongly polarizable molecules or short separation distances, the order parameter (50) drops from vanishingly small values corresponding to random dipole fluctuations down to the limiting value of S i (z) = −1/2 indicating parallel molecular orientation with the membrane surfaces. Thus, as the membrane walls get closer, the partial salt exclusion from the intermembrane region is accompanied with the alteration of the electronic cloud of the confined species from a spherical to an ellipsoidal shape. In the case of neutral membranes, we found that the ionic size and polarizability triggering salt depletion from the intermembrane region induces the purely attractive force (40). In order to consider the case of charged membranes, in Figs. 3(a)-(b), we plotted the pressure (26) versus the distanced = κd rescaled by the inverse DH length κ = √ 8π B ρ ib at the dimensionless parameter s = κµ = 15. These plots show that the extended inner structure of the ions confined to charged membranes gives rise to the coexistence of a large-distance attractive regime and a short-distance repulsive force of enhanced magnitude. Namely, Fig. 3(a) indicates that upon the rise of the dimensionless polarizabilityᾱ = κα at vanishing average sizeā = κa = 0, the large distance branch of the PB pressure (dashed curve) turns from weakly repulsive to attractive. However, at short distances, the repulsive PB pressure rises rapidly with the polarizability nearly by a factor of two, and saturates atᾱ d . Fig. 3(b) shows that the increase of the fixed ion sizeā at vanishing polarizabilityᾱ = 0 causes an identical amplification of the short-range repulsion, and the switching of the large distance pressure from repulsive to attractive. However, one notes that the fixed ion size leads to shorter range attractive interactions than the intrinsic polarizability, and also triggers the emergence of a cusp at the characteristic distance d = a above which the rotational entropy of the confined ions rises significantly.
Finally, with the aim to identify the effect of the salt concentration ρ b and the membrane charge density σ on the intermembrane force, in Fig. 3(c), we display the pressure at various values of the dimensionless parameter s ∝ √ ρ b /σ. One sees that upon the reduction of this parameter, the short-range repulsive pressure regime overrides gradually the attractive longer-range branch. Below the characteristic value s ∼ 10, the attractive pressure branch disappears, and the interaction force be-comes purely repulsive. Hence, salt decrement or equivalently the rise of the membrane charge strength drives the system from the depletion-induced attraction to the polarizability-enhanced pure repulsion regime.

Analytical pressure formula from the Donnan approach
In order to provide physical insight into the ion size effects identified above, we derive here a closed-form expression for the interaction pressure (26). To this aim, we introduce a constant Donnan potential approximation that neglects the spatial variations of the electrostatic potential. Thus, we set φ(z) ≈ φ D . Substituting this approximation into the nuclear charge density (16), the grand potential density (24) reduces to with the steric partition coefficient t i given by Eq. (48). Evaluating the saddle-point condition ∂f /∂φ D = 0 with Eq. (51), one obtains the second order algebraic equation e −φD − e φD − 2r = 0 yielding the Donnan potential as where we defined the ratio of the fixed membrane and mobile charge densities per surface as r = σ/(dρ b t i ). Evaluating now the inner pressure βP = −δf /δd with Eq. (51), and subtracting the bulk osmotic pressure, the net pressure follows as with the density renormalization factor of entropic origin, We note that for neutral membranes (r → 0), Eq. (53) yields the exact MF pressure (40). Finally, taking the point charge limit a → 0 and α → 0 of Eq. (53), the Donnan approximation for the PB pressure follows as In Fig. 3(c), the comparison of the circles and the solid curves shows that the Donnan pressure (53) provides a highly accurate approximation of the exact MF pressure (26). Thus, below, we use Eq. (53) to characterize the asymptotic behavior of the ion size effects on the charged membrane interactions.

Asymptotic behavior of the intermembrane pressure
For large salt concentrations and substantial separation distances σ/(ρ b d) 1, the Taylor expansion of Eq. (53) for d α a leads to Then, the expansion of Eq. (53) for d a α yields Finally, in the same large-distance regime σ/(ρ b d) 1, the Donnan approximation of the PB pressure for point charges follows from the Taylor expansion of Eq. (55) as  (53) and (55), and the circles are from the numerical evaluation of the contact value identities (23) and (26).
The comparison of Eqs. (56)-(57) with Eq. (58) shows that ion size and polarizability bring to the PB pressure an attractive contribution of linear dependence on the salt concentration, and a membrane charge-dependent repulsive contribution scaling inversely with the salt concentration. In Fig. 3(c), this mixed effect of the ion size is precisely responsible for the transition of the pressure from the polarizability-enhanced pure repulsion to the ionic depletion-driven pure attraction via salt increment.
In order to compare the individual effects of the fixed solute size and polarizability, Figs. 4(a)-(b) display this transition for solely polarizable ions (a = 0) and purely rigid charges (α = 0) together with the long-range asymptotic laws (56)-(57) (circles). The comparison of the left and right panels indicates that due to the longer range effect of the intrinsic polarizability on the intermembrane force, the switching from repulsion to attraction at finite polarizability occurs at lower salt densities than in the presence of a fixed ion size of the same value.
The opposite short-distance regime corresponding to a counterion-only liquid is characterized by the inequalities σ/(ρ b d) 1, d a , and d α. The expansion of Eqs. (53) and (55) for short separation distances yields , the corresponding effect is illustrated against the ion polarizability and size, respectively. One sees that as these parameters increase and exceed the intermembrane distance, the pressure quickly rises above the PB value and saturates at a plateau located at P net 2P (PB) net . Then, as the membrane separation distance is reduced, this plateau approaches from above the value of two, i.e. P net /P (PB) net → 2 + as d → 0. According to the previous MC simulations and perturbative theories of point ions, the electrostatic MF and strong-coupling (SC) regimes of the short-distance macromolecular interactions are mutually governed by the ideal gas law corresponding to the first term of Eq. (60) [48][49][50]. Considering this universality, the ion size-induced twofold enhancement of the short-distance pressure in Eq. (59) is a key prediction of our work.

Mechanism driving the short-range enhanced repulsion
Here, we probe the underlying mechanism behind the enhancement of the intermembrane repulsion by the extended solute structure. Above, we showed that the fixed ion size and polarizability affect the pressure in a qualitatively similar fashion. Thus, with the aim to simplify the analysis, we set a = 0. Figs. 6(a)-(c) display the nuclear charge densities of the salt species characterizing their C.M. distribution, and the average potential at various polarizability valuesᾱ. One sees that due to the overlap of the electronic cloud with the membrane walls, together with the intramolecular charge separation amplifying the coupling of the nuclear and membrane charges, the rise of the polarizability reduces the anion density in the entire intermembrane region, and also decreases the cation density in the middle of the slit. As a result of the electroneutrality constraint, the corresponding salt exclusion amplifies the pore potential by a factor of two, enhancing the surface cation density in Fig. 6(b). Thus, the first effect of the finite solute size consists in strengthening the anion rejection and the cation adsorption by the charged membranes, driving the system closer to the counterion-only ideal gas regime.
To illustrate the main effect of the finite polarizability leading to the twofold amplification of the short-distance PB-level repulsion, in Fig. 6(d), we plotted the electronic charge density of the interfacially adsorbed cations. The comparison of Figs. 6(b) and (d) shows that despite the similar sign of the electronic and membrane surface charges, the contact density of the nuclear and electronic charges of a given cation are significantly close.
To explain this feature, in Fig. 6(e), we reported the nematic order parameter (31) of the cation species. Similar to the neutral membrane case displayed in Fig. 2(d), the order parameter is seen to be located in the negative regime −1/2 < S + (z) < 0. Hence, as the compression of the electronic cloud by the membrane walls dominates its electrostatic repulsion by the membrane charges, the electronic charge distribution adopts an ellipsoidal shape indicating solute orientation parallel with the membrane surfaces. In the vicinity of the membrane substrate, this orientation leads to an electronic charge density as localized as the nuclear charge density, causing the nearly identical momentum exchange of the electronic and nuclear charges with the membrane walls. As a result, the net force exerted by an extended solute molecule on the membrane is twice as strong as that induced by a pointlike ion. Consequently, Fig. 6(f) shows that the rise of the polarizability into the regime α d amplifies the total contact density of the PB theory (dashed curve) and the corresponding osmotic pressure by a factor of two.

Derivation of the ideal gas limits from electroneutrality
For an analytical insight into the aforementioned mechanism driving the twofold amplification of the PB-level repulsion, we derive the ideal gas terms in Eqs. (59)-(60) directly from the electroneutrality condition. In the well-known case of point charges, the electrolyte configuration at short intermembrane distances corresponds to a thin counterion-only layer parallel with the membrane surfaces [48,49]. Within this picture, the membrane walls of area S and total surface charge Q s = −2Sσe confine a total amount of Q + = Sdρ + (d)q + e counterion charges. Accounting for the electroneutrality constraint Q s + Q + = 0, the cationic contact density follows as ρ + (d) = 2σ/(dq + ). Substituting the latter into the contact value identity (23), one obtains the PB-pressure βP (PB) net ≈ 2σ/(q + d) − 2ρ b corresponding to the first two terms of Eq. (60).
We follow now the same lines for the more general case of extended ions. First, we note that the total cation charge is Q + = Sde [ρ +,n (d)n + − ρ +,e (d)e + ]. Thus, the electroneutrality condition Q s + Q + = 0 gives Taking now into account the comparable contact density of the nuclear and electronic charges displayed in Figs. 6(b) and (d), Eq. (61) yields ρ +,n (d) ≈ ρ +,e (d) ≈ 2σ/(q + d). Finally, injecting into the extended contact value identity (26) the above equalities implying equal osmotic pressure contributions from the electronic and nuclear charges, the dominant terms in Eq. (59) follows as βP net ≈ 4σ/(q + d)−2ρ b , confirming the twofold enhancement of the short-range PB repulsion calculated above.

IV. CONCLUSIONS
The complex charge structure of the solute molecules ubiquitous in living systems requires the characterization of the macromolecular interactions governing biological processes beyond the point-ion approximation. Motivated by this need, we developed a field-theoretic formalism of intermolecular interactions incorporating the extended charge structure as well as the polarizability of the solute molecules via a Drude oscillator model. The novelty of our formalism with respect to the prior works is the full tratment of the finite solute size and polarizability, and the resulting entropic effects usually bypassed by the point-dipole approaches [31,33,34].
In the electrostatic MF regime of our model, the average potential and the macromolecular interaction force follow from the integro-differential PPB Eq. (15) and the generalized contact value identity (26). Via the numerical evaluation of these equations, and the use of the analytical Donnan pressure formula (55) of significant accuracy, we showed that the extended charge partition on the solute molecules generates a short-range enhanced repulsion, and a longer-range depletion attraction absent in the PB formalism for point-like ions. The physical mechanisms driving these effects are summarized below.
At short separation distances characterized by a counterion-only liquid, due to the entropically driven solute orientation parallel with the anionic membranes, the contact density of the electronic charges located on the adsorbed counterions is nearly as high as the contact density of the cationic nuclear charges. As a result, the momentum transfer from a finite-size ion to the membranes is twice as strong as the momentum exchange induced by a point charge of the same valency. The resulting twofold enhancement of the short-range repulsion by the ion size and polarizability is the main prediction of our work.
In the subsequent regime of intermediate membrane distances where the anionic and cationic salt species are mutually present in the intermembrane region, the restriction of the electronic cloud size by the membrane walls causes a partial solute rejection into the reservoir. This entropic effect reduces the inner osmotic pressure below the bulk pressure, generating a depletion attraction between the like-charged membranes.
We also characterized the effect of salt on the interplay between these opposing forces. We found that at dilute salt concentrations, the membrane interactions are governed by the polarizability-enhanced pure repulsion. Upon salt increment, the emerging depletion attraction overrides gradually the short-range repulsion. Beyond a characteristic salt concentration, the like-charge membrane coupling becomes purely attractive.
We emphasize that our predictions are not limited to polarizable ions; from polyamines to non-uniformly charged proteins and viruses, the finite molecular size and flexibility are common features of a large variety of solute molecules involved in biological processes [51,52]. Therefore, our model and findings may help to comprehend qualitatively the role and behavior of various force mediating molecules under nanoconfinement. Moreover, we note that via the inclusion of radial charge distributions obtained from X-ray spectroscopy [53] and quantum mechanical calculations [54], the present formalism can be generalized to characterize the effect of more specific ion structures on macromolecular interactions. The coupling of our electrostatic theory with hydrodynamic transport equations would also enable to incorporate salt specificity into ion transport formalisms used in nanofiltration and nanofluidic technologies [55,56]. We plan to consider these relevant extensions in future works.
In the present study, we investigated exclusively the electrostatic MF regime known to be valid for low valency ions and weak membrane charges [49]. In the future, we would like to characterize the effect of the finite solute size and polarizability in the opposite strong-coupling regime of macromolecular interactions [48]. We note in passing that in the corresponding regime characterized by the emergence of ionic self-energies, the dielectric cavities associated with the inner ion structure should be also taken into account [4]. An additional limitation of our work is its implicit solvent nature neglecting the electrostatics of water molecules. It should be however noted that according to our recent results on explicit solvent liquids including point-like salt ions, the cancellation of various solvent effects of large magnitude but opposite sign leads to a nearly vanishing contribution from explicit solvent to intermolecular forces [37]. The inclusion of solvent molecules into the present model will be necessary to verify the validity of this cancellation mechanism in the more general case of ions with extended charge structure.