On the nature of screening in Voorn-Overbeek type theories

While it was originally formulated in the context of phase separation between charged polyions, different variants of the Voorn-Overbeek (VO) theory [1] still represent the basic conceptual underpinning of diverse phenomena ranging from the theory of phase separation in solutions of weakly charged polyelectrolytes [2, 3], swelling of charged gels [4, 5], polyelectrolyte complex coacervates [6, 7, 8], as well as liquid-liquid phase coexistence phenomena in biology [9, 10].

Formally, the VO theory is anchored in the competition between the configurational entropy of charged polyions and electrostatic correlation attraction between them, the former being evaluated within the Flory-Huggins polymer mixing framework [11], while the latter is based on the Poisson-Boltzmann (PB) electrolyte theory in the Debye-Hückel (DH) approximation [12]. The system is modelled as consisting of three components: water and two types of polyions, denoted by $p^{+}$ and $p^{-}$, of charge $\pm eN_{p^{\pm}}$, and degree of polymerization $N_{p^{\pm}}$. The monomers and water molecules are all assumed to have the same molecular volume, $v=a^{3}$.

The components of the VO free energy density, following the local density approximation where the inhomogeneities described only by the coordinate dependence of the densities, are then assumed to have the following forms:

(i) The Flory-Huggins polymer solution free energy is given by [13, 11]

$\displaystyle\frac{f_{FH}(\phi_{p_{+}},\phi_{p_{-}})~{}a^{3}}{k_{B}T}={\frac{\phi_{p_{+}}}{N^{+}}}\ln\phi_{p_{+}}+{\frac{\phi_{p_{-}}}{N^{-}}}\ln\phi_{p_{-}}+{(1-\phi_{p_{+}}-\phi_{p_{-}})}\ln(1-\phi_{p_{+}}-\phi_{p_{-}})+{\textstyle\frac{1}{2}}\sum_{j,k}\chi_{jk}\phi_{j}\phi_{k},$

(1)

where $f_{FH}$ is the free energy density, $k_{B}T$ is the thermal energy and $\phi_{\pm}$ the volume fraction of $p^{\pm}$ polymers.

We included also the $\chi_{jk}$ interaction term, describing the short range interactions of non-electrostatic nature such as the van der Waals interactions, which is also frequently included in later formulations of the theory [14]. The connection between the volume fractions $\phi_{p^{\pm}}$ and the concentrations $c_{p^{\pm}}$ is given by $\phi_{p^{\pm}}=a^{3}c_{p^{\pm}}N_{p^{\pm}}$.

(ii) The DH dilute electrolyte correlation free energy can be obtained in different ways [15, 16, 17] and is given by

	$\displaystyle\frac{f_{DH}(\phi_{p_{+}},\phi_{p_{-}})~{}a^{3}}{k_{B}T}=-{\frac{(\kappa_{D}a)^{3}}{12\pi}}=$
	$\displaystyle~{}~{}~{}~{}~{}~{}-{\textstyle\frac{2}{3}}\sqrt{\frac{\pi}{a^{3}}}~{}\ell_{B}^{3/2}\Big{(}{N_{p^{+}}}\phi_{p^{+}}+{N_{p^{-}}}\phi_{p^{-}}\Big{)}^{3/2},$		(2)

where again $f_{DH}$ is the free energy density, and the inverse Debye screening length, $\kappa_{D}^{-1}$, is given by $\kappa_{D}^{2}=4\pi\ell_{B}\left({c_{p^{+}}{N_{p^{+}}}^{2}}+{c_{p^{-}}{N_{p^{-}}}^{2}}\right)$, where $\ell_{B}$ is the Bjerrum length, $\ell_{B}=e^{2}/(4\pi\varepsilon k_{B}T)$, $\varepsilon=\epsilon\epsilon_{0}$ with $\epsilon$ the relative dielectric permittivity and $e$ the elementary charge. The assumption in the above formula is that the valency of polyions coincides with the number of monomers.

Often the DH correlation free energy is taken in the form corresponding to a finite charge radius (see below) that of course reduces to the above expression for vanishing ionic size [15, 16, 12]. The DH correlation free energy is generally obtained from the thermal fluctuations around the mean-field (saddle-point) PB theory [18], either by a coupling constant integration of the electrostatic potential of a single ion [17], or by integrating out the quadratic electrostatic potential fluctuations around the PB saddle-point [19]. The two methods are completely equivalent and either way, the inverse Debye screening length can be obtained as [20]

$\displaystyle\kappa_{D}^{2}=\varepsilon^{-1}\frac{\partial\rho(\psi)}{\partial\psi}{\Big{|}}_{\psi=0}=\frac{{\cal C}}{\varepsilon},$

(3)

where $\rho$ is the mobile charge density, $\psi$ the mean-field PB electrostatic potential, and ${\cal C}$ is the capacitance density at thermal equilibrium. The Debye screening length can be interpreted as the thickness of an equivalent parallel plate condenser whose surface charges result from an imposed potential difference. By the general statistical mechanical relationship between fluctuations and response functions, it can be expressed also in terms of the thermal fluctuations of the electrostatic potential by the Einstein formula [21].

The Debye screening length is of course straightforward to calculate in the Debye-Hückel framework, where the entropy of the ions is given by the ideal gas expression, and - as we will see shortly - the Legendre transform of the free energy density and its second derivatives can be calculated explicitly and analytically. It is less clear how to approach the screening problem for a general free energy which is not in the form of the ideal gas entropy, as is usually the case for the VO type theories.

The correlation free energy Eq. I thus quantifies the electrostatic potential fluctuations in the solution, whatever its composition and whatever its other degrees of freedom are. The importance of electrostatic potential fluctuations is the fundamental insight of the VO theory. In some sense it can be seen also as a special case of the van der Waals theory, that recognizes the importance of the electrodynamic fluctuations via the dispersion interaction [22]. Similarly to the van der Waals theory, the VO theory displays features of a mean-field theory [6] even if the attractive Debye-Hückel term corresponds to fluctuations around a zero potential mean-field state [18].

The characterizing features of the theory can be and were criticized on different levels and as a consequence it has been generalized/ammended to include either better approximations, better models or both (see the discussion in the recent excellent review Ref. [7]), with some of the major new directions emerging from the implementations of the polymer field theory [23, 24], with more detailed consideration of connectivity and excluded volume effects, the scaling or ‘blob’ theory approaches [25], based on the competition of thermal and electrostatic degrees of freedom that set the key length scales, as well as explicit inclusion of counterion condensation and/or localization mechanisms that includes additional degrees of freedom such as non-ionic interactions in the description of the polyelectrolyte complexation [26]. Nevertheless, the most widely used point of departure for further developments remains the Voorn-Overbeek theory [14] that endures also as the baseline for comparison [7].

Here, however, our ambition is more technical in nature and can be formulated as follows: since the VO free energy is not of the ideal gas form, as assumed in the DH theory that is consistent with the Debye screening length, what is then the correct form of the screening length that is consistent with the VO or related theories? The methodology employed is based on standard properties of the Legendre transform as reviewed by Zia et al. [27] as well as a generalization of the screening length for a general free energy Ansatz as developed in Ref. [28]. We will write down the general equations for the screening length, based on the local curvature properties of the Legendre transform, apply them specifically to the VO model free energies and finally comment on when the screening actually reduces to the simplified and standardly assumed Debye form in the context of the polyelectrolyte condensation. The approach advocated here is rather similar in essence to the mean-field analysis of the steric [29, 28], correlation [30, 31] and structural [32, 33] effects that result in modified Poisson-Boltzmann equations with different types of local nonlinearities which cover steric effects, higher-order derivative terms that describe the correlation effects or both for the structural effects. While these approaches certainly address the main features of the phenomena and are computationally tame, more sophisticated theories are needed to address the combined effects of packing constraints, correlation interaction and solvent structure [34, 35].

Our aim is thus not to generalize or improve on the VO theory. What we set out to do has a more modest but nevertheless a fundamental goal: we will address the problem of the form of the screening length, avoiding the standard definition of the Debye screening length, and making it consistent with the other terms in the VO free energy Ansatz.

We start with the free energy, or rather free energy density, of an uncharged system composed of ${\cal N}$ components: $f(c_{j})=f(c_{1},c_{2},\dots c_{\cal N})$. In an analysis proceeding from the Legendre transform and the local thermodynamics approximation, Maggs and Podgornik [28] have recently shown that the thermodynamic potential of a charged system, where each component $j$ has a charge $e_{j}N_{j}$, can be written as

$${\cal F}[\mu_{j},\psi]=-\int_{V}d^{3}{\bf r}\Big{(}{\textstyle\frac{1}{2}}{\varepsilon}(\mbox{\boldmath$\nabla$}\psi)^{2}+p(\mu_{j}+e_{j}N_{j}\psi)\Big{)},$$

(4)

where $p(\mu_{j})=p(\mu_{1},\mu_{2},\dots\mu_{\cal N})$ is the thermodynamic pressure, or the equation of state, defined as the Legendre transform of the free energy [36], $f(c_{j})-\sum_{j}\mu_{j}c_{j}=-p(\mu_{j})$, while $\mu_{j}$ is the chemical potential of the $j$-th species. $\psi({\bf r})$ is the Legendre multiplier field, identified as the electrostatic potential, that ensures the local imposition of the Gauss’ law. Clearly the above result reduces to the standard PB form since for a two component uni-univalent electrolyte the equation of state has the ideal gas form $p(\mu\pm e\psi)=\rho_{+}+\rho_{-}=2e^{\mu}\cosh{e\psi}$ [37].

While the whole derivation proceeded entirely on the mean-field level, it can be extended to the case when the Coulomb interactions are included exactly and the mean potential becomes the fluctuating local potential in a functional integral representation of the partition function as derived by Wiegel [38].

In fact the field theory representation of the partition function is formally obtained with a field action at imaginary values of the electrostatic potential, ${\cal F}[\psi]\longrightarrow{\cal F}[i\psi]$ (for details see Ref. [39]). The mean field or the PB solution then corresponds to the saddle-point of the field action with $\psi\longrightarrow\psi_{0}$ and free energy $F_{PB}(\psi_{0})$, while the Gaussian fluctuations around the mean field can be integrated out yielding the one-loop fluctuation corrected free energy density, $F=F_{PB}(\psi_{0})+F_{c}$, where the correlation free energy is given by

$\displaystyle F_{c}=-{\textstyle\frac{1}{2}}k_{B}T~{}{\rm Tr}\ln{\frac{\delta^{2}{\cal F}[i\psi]}{\delta\psi({\bf r})\delta\psi({\bf r}^{\prime})}}{\Big{|}}_{\psi_{0}}.$

(5)

with field Hessian

$\displaystyle\frac{\delta^{2}{\cal F}[i\psi]}{\delta\psi({\bf r})\delta\psi({\bf r}^{\prime})}=\left(-\varepsilon\nabla^{2}+\frac{\partial^{2}p(\mu_{j}+ieN_{j}\psi)}{\partial\psi^{2}}{\Big{|}}_{\psi_{0}}\right).~{}$

(6)

If the bulk system is electroneutral then the second derivative of $p(\mu_{j}+ieN_{j}\psi)$ is evaluated at $\psi_{0}=0$, otherwise it has to be evaluated at the value of the Donnan potential, $\psi_{0}=\psi_{D}$. Evaluating the last term Eq. 5 in the Fourier space, subtracting the single ion contribution and assuming a vanishing size of ions, yields after some rearrangements the first identity in Eq. I with the inverse screening length defined as

$$\kappa^{2}\equiv\frac{\partial^{2}p(\mu_{j}+ieN_{j}\psi)}{\varepsilon~{}\partial\psi^{2}}{\Big{|}}_{\psi_{0}=0}\!\!\!=\frac{e^{2}}{\varepsilon}\sum_{j,k}N_{j}N_{k}\frac{\partial^{2}p(\mu_{j})}{\partial\mu_{j}\partial\mu_{k}},$$

(7)

where we assumed an electroneutral bulk system with the potential $\psi_{0}=0$. Invoking furthermore the Gibbs-Duhem relation, $c_{j}={\partial p(\mu_{j})}/{\partial\mu_{j}}$, we can rewrite the inverse screening length with the density derivatives, which are of course the relevant response functions. Equation 7 is then the proper generalization of the Debye definition of the screening length Eq. 3.

While the point charge form of the DH correlation free energy Eq. I is commonly used in the VO theory, it can be straightforwardly generalized to include the finite size of the ions. In this case the electrostatic correlation free energy change can be written as

$\displaystyle\Delta F_{c}=-{\textstyle\frac{1}{2}}k_{B}T~{}{\rm Tr}\ln{\frac{-\nabla^{2}+\kappa^{2}}{-\nabla^{2}}}$

(8)

with $\kappa^{2}$ given by Eq. 7. This can be furthermore rewritten as a coupling constant integral of the resolvent operator

$\displaystyle\Delta F_{c}=-{\textstyle\frac{1}{2}}k_{B}T~{}{\rm Tr}\int_{0}^{1}d\lambda~{}{\cal G}_{\lambda}({\bf r},{\bf r}^{\prime})$

(9)

with the Green’s function ${\cal G}_{\lambda}({\bf r},{\bf r}^{\prime})$ given as a solution of

$\displaystyle\left(-\nabla^{2}+\lambda\kappa^{2}\right){\cal G}_{\lambda}({\bf r},{\bf r}^{\prime})=\delta({\bf r}-{\bf r}^{\prime})$

(10)

where one needs to subtract the contribution of the $\kappa^{2}=0$ solution, see Eq. 8. This can be obtained straightforwardly as [15, 16, 17]

$\displaystyle{\cal G}_{\lambda}({\bf r},{\bf r})=\frac{\lambda^{1/2}\kappa}{4\pi~{}(1+\lambda^{1/2}\kappa a)}$

(11)

since in the evaluation of $\rm Tr$ only the value at identical arguments is needed. Insering back into Eq. 9 one derives the expression commonly used in the free energy Ansatz of most current VO type theories [7, 8, 5]

$\displaystyle\Delta F_{c}=-\frac{V}{4\pi a^{3}}\left(\ln{(1+\kappa a)}-(\kappa a)+{\textstyle\frac{1}{2}}(\kappa a)^{2}\right),$

(12)

with ${\rm Tr}"1"=V$, the volume of the system. In the limit $a\longrightarrow 0$ this reduces back to Eq. I. Notably in the above derivation the Debye charging process [15, 16, 17] is recognized simply as the ”coupling constant integration” of the resolvent operator.

The form of the inverse Debye screening length as given in Eq. 7 together with the curvature duality of the Legendre transform [27], leads to its straightforward and elegant calculation for any form of the free energy even when the equation of state $p(\mu_{j})$ is not explicitly available [28].

We can test these expressions on the original two component Debye-Hückel theory with $N_{1}=N_{2}=1$ and volume fractions $\phi_{1,2}=a^{3}c_{1,2}$, where the entropy of the uncharged systems is assumed to be just the ideal gas entropy

$\displaystyle\frac{f_{D}(\phi_{1},\phi_{2})~{}a^{3}}{k_{B}T}={\phi_{1}}\left(\log{\phi_{1}}-1\right)+{\phi_{2}}\left(\log{\phi_{2}}-1\right),$

(13)

the Legendre transform of which can be obtained analytically in the simple form of

$$p_{D}(\mu_{1},\mu_{2})=\frac{k_{B}T}{a^{3}}\left(e^{\beta\mu_{1}}+e^{\beta\mu_{2}}\right)=p(\mu_{1})+p(\mu_{2}),$$

(14)

clearly additive in the two components. The inverse Debye length from Eq. 7, or equivalently from the expansion of the PB equation, then follows as

$$\kappa_{D}^{2}=\frac{4\pi\ell_{B}}{a^{3}}\left(N_{1}\phi_{1}+N_{2}\phi_{2}\right)=4\pi\ell_{B}(N_{1}^{2}c_{1}+N_{2}^{2}c_{2}),$$

(15)

where $\ell_{B}$ is again the Bjerrum length. This is of course nothing but the standard Debye screening length.

In what follows we will analyze the screening length of two models, see Fig. 1: a model of polyions (model 1) that reduces to the Flory–Huggins theory for uncharged polymers and a model of polyions (model 2) with monovalent salt that is a combination of the Flory-Huggins theory for polyions and DH theory for simple salt ions [40, 8]. The ensuing screening length will in general differ substantially from a simple Debye screening form as we elucidate next.

We analyze a couple of non-trivial generalizations of the Debye screening for models that are standardly used in the description of complexation in mixtures of different polyelectrolytes [41, 7, 42, 40, 14, 8].

Model 1. We start with a polydisperse polyion mixture, model 1 , composed of species “1” at concentration $c_{1}$, itself composed of $N_{1}$ monomers, each of charge $e$, and species “2” at concentration $c_{2}$, itself composed of $N_{2}$ monomers, each of charge $-e$, in an aqueous solvent of (water) molecules of diameter $a$. This model has been used for the description of complexation in mixtures of oppositely charged polyelectrolytes [41, 7, 42].

Expressed in terms of the volume fractions $\phi_{1},\phi_{2}$, defined as $\phi_{1,2}=a^{3}c_{1,2}N_{1,2}$, the Flory-Huggins lattice level free energy density is [11, 40]

	$\displaystyle\frac{f(\phi_{1},\phi_{2})~{}a^{3}}{k_{B}T}=\frac{\phi_{1}}{N_{1}}\log{\phi_{1}}+\frac{\phi_{2}}{N_{2}}\log{\phi_{2}}$
	$\displaystyle+(1-\phi_{1}-\phi_{2})\log{(1-\phi_{1}-\phi_{2})}.$		(16)

We omit the $\chi$ interaction term as it is irrelevant for subsequent derivations. The Hessian of the pressure $p(\mu_{1},\mu_{2})$ is positive definite, and according to Eq. 7, the inverse Debye screening length is given by

$$\kappa^{2}=\frac{e^{2}}{\varepsilon}\left(N_{1}^{2}\frac{\partial^{2}p}{\partial\mu_{1}^{2}}-2N_{1}N_{2}\frac{\partial^{2}p}{\partial\mu_{1}\mu_{2}}+N_{2}^{2}\frac{\partial^{2}p}{\partial\mu_{2}^{2}}\right).$$

(17)

Contrary to the original two component Debye theory, derived in the previous section, the chemical potential derivatives cannot be evaluated explicitly as the lattice entropy does not posses a simple analytical Legendre transform.

At this point, however, one can take recourse to the general properties of the Legendre transform, specifically to its curvature duality property. One of the fundamental properties of the Legendre transform is that the Hessian of the Legendre transform is the inverse of the Hessian of the function itself [27], so that one can write

$$\sum_{m=1}^{2}\frac{\partial^{2}p(\mu_{1},\mu_{2})}{\partial\mu_{j}\partial\mu_{m}}\frac{\partial^{2}f(c_{1},c_{2})}{\partial c_{m}\partial c_{k}}=\delta_{jk}$$

(18)

where all the matrices are $2\times 2$. This implies that the local curvature of the Legendre transform is inverse to the local curvature of the original function in a manner reminiscent of the uncertainty relation, as observed by Zia et al. [27]. The above relation remains valid only for strictly convex functions so that neither derivative ever vanishes.

The matrix of partial derivatives $\frac{\partial^{2}f(c_{1},c_{2})}{\partial c_{m}\partial c_{k}}$ from the free energy Eq. 16 can be calculated straightforwardly and from Eq. 18 one then obtains the derivatives $\frac{\partial^{2}p(\mu_{1},\mu_{2})}{\partial\mu_{j}\partial\mu_{m}}$ which then yield $\kappa^{2}$ from the combination in Eq. 17, so that

	$\displaystyle\kappa^{2}(N_{1},N_{2},\phi_{1},\phi_{2})$
	$\displaystyle~{}~{}~{}~{}~{}~{}=\frac{4\pi\ell_{B}}{a^{3}}\frac{u\left(N_{1}\phi_{1}+N_{2}\phi_{2}\right)+4N_{1}N_{2}\phi_{1}\phi_{2}}{\left(N_{1}\phi_{1}+N_{2}\phi_{2}+u\right)},$		(19)

where we used the abbreviation $u=(1-\phi_{1}-\phi_{2})$ and the Bjerrum length $\ell_{B}$ was defined before.

In general the inverse square of the screening length is therefore not a linear function of the concentrations as is the case for the Debye screening length. The bulk electroneutrality furthermore restricts the concentrations of components to $N_{1}c_{1}-N_{2}c_{2}=0$ or equivalently

$\displaystyle\phi_{1}-\phi_{2}=0,$

(20)

while both $N_{1},N_{2}$ can still remain arbitrary. Denoting $\phi_{1,2}=\phi_{M}$, Eq. 19 can be recast in the form

$\displaystyle\kappa^{2}(N_{1},N_{2},\phi_{M})=\frac{4\pi\ell_{B}}{a^{3}}\frac{\left(N_{1}+N_{2}\right)\phi_{M}+2\big{(}2N_{1}N_{2}-(N_{1}+N_{2})\big{)}\phi_{M}^{2}}{\big{(}1+({N_{1}}+{N_{2}}-2)\phi_{M}\big{)}}.$

(21)

In the case of $N_{1}=N_{2}=1$, as well as for any symmetric case, $N_{1}=N_{2}$, this obviously reverts back to the standard Debye form . However, for any other case the screening length is a much more complicated function of the volume fractions, or concentrations of the species.

On Fig. 2 we show the dependence of the ratio

$$\tilde{\kappa}^{2}(N_{1},N_{2},\phi_{M})=\frac{\kappa^{2}(N_{1},N_{2},\phi_{M})}{(N_{1}+N_{2})\phi_{M}}$$

(22)

for an electroneutral system. The denominator would be the standard inverse square of the Debye screening length expected for point ions. Clearly the dependence of $\tilde{\kappa}^{2}$ is in general not a linear function of $N_{1},N_{2}$ and can in addition show strong non-monotonic behavior with a minimum at the standard Debye screening value corresponding to $\tilde{\kappa}^{2}=1$. This non-monotonic regime in the vicinity of $N_{1}=N_{2}$ state, at which the screening length exhibits a minimum value, is the more pronounced the larger is the volume fraction of the polyions. Obviously the non-symmetric systems are then under-screened with the screening length exceeding the Debye screening length expectations.

Model 2. We now proceed to a more complicated system, model 2, composed of a uni-univalent salt as well as a polydisperse polyion mixture that we already discussed above. We refer to this model as the Voorn-Overbeek-type model, which has been recently used to model complexation in mixtures of oppositely charged polyelectrolytes with explicit salt ions [40, 14, 8].

Positively “3” and negatively “4” charged polymers have $N_{3}$ and $N_{4}$ monomers, while the salt is composed of univalent positively “1” and negatively “2” charged simple salt ions. The free energy can be taken in a form generalizing the Flory-Huggins lattice entropy as

	$\displaystyle\frac{f(\phi_{1},\phi_{2},\phi_{3},\phi_{4})~{}a^{3}}{k_{B}T}$	$\displaystyle=$	$\displaystyle{\phi_{1}}\log{\phi_{1}}+{\phi_{2}}\log{\phi_{2}}+\frac{\phi_{3}}{N_{3}}\log{\phi_{3}}+\frac{\phi_{4}}{N_{4}}\log{\phi_{4}}$		(23)
		$\displaystyle+$	$\displaystyle(1-\phi_{1}-\phi_{2}-\phi_{3}-\phi_{4})\log{(1-\phi_{1}-\phi_{2}-\phi_{3}-\phi_{4})}.$

The first two terms describe the simple monovalent salt, the next two terms correspond to polyions, while the last term is the solvent entropy. The above free energy is clearly a straightforward generalization of Eq. 16.

We now proceed in the same way as before, except that now the analysis is a bit more involved since we have a four component system: salt and polyions, so that both the free energy, $f(\phi_{1},\phi_{2},\phi_{3},\phi_{4})$, as well as the pressure, $p(\mu_{1},\mu_{2},\mu_{3},\mu_{4})$, are functions of four variables and consequently the Hessian matrices of derivatives will now be $4\times 4$, instead of $2\times 2$, with the algebra correspondingly more difficult, but not unmanageable analytically.

The Debye length in this case is obtained with complete analogy to Eq. 17 as

$$\kappa^{2}=\frac{e^{2}}{\varepsilon}\left({\partial_{\mu_{1}}}\!\!-\!\!{\partial_{\mu_{2}}}\!\!+\!\!N_{3}{\partial_{\mu_{3}}}\!\!-\!\!N_{4}{\partial_{\mu_{4}}}\right)^{2}p(\mu_{1},\mu_{2},\mu_{3},\mu_{4}).$$

(24)

In order to evaluate this, we need to invert the Hessian of the original function and calculate the Hessian of its Legendre transform through the curvature duality relation of the Legendre transform, Eq. 18, which in this case assumes the following form

	$\displaystyle\sum_{m=1}^{4}\frac{\partial^{2}p(\mu_{1},\mu_{2},\mu_{3},\mu_{4})}{\partial\mu_{j}\partial\mu_{m}}\frac{\partial^{2}f(c_{1},c_{2},c_{3},c_{4})}{\partial c_{m}\partial c_{k}}=$
	$\displaystyle~{}~{}~{}~{}~{}~{}=\sum_{m}p_{j,m}f_{m,k}=\delta_{jk},$		(25)

where we have introduced the following notation for the double derivatives of free energy density and pressure

$$f_{m,k}=\partial_{c_{m}}\partial_{c_{k}}f(c_{1},c_{2},c_{3},c_{4})$$

(26)

with $f(\phi_{1},\phi_{2},\phi_{3},\phi_{4})$ given by Eq. 23, and

$$p_{m,k}=\partial_{\mu_{m}}\partial_{\mu_{k}}p(\mu_{1},\mu_{2},\mu_{3},\mu_{4})$$

(27)

Again, Eq. 25 implies that the local curvatures of the Legendre transforms are inverse to each other as observed before [27].

Expressing the derivatives of the equation of state $\partial_{\mu_{j}}\partial_{\mu_{m}}p(\mu_{1},\mu_{2},\mu_{3},\mu_{4})$ with the derivatives of the free energy $\partial_{c_{m}}\partial_{c_{k}}f(c_{1},c_{2},c_{3},c_{4})$ and inserting the result into Eq. 24, we can then obtain the inverse square of the screening length in terms of the derivatives of the equation of state in the form analogous to Eq. 17 as

$\displaystyle\kappa^{2}=\frac{e^{2}}{\varepsilon}\Big{(}p_{1,1}-2p_{1,2}+p_{2,2}+2N_{3}p_{1,3}-2N_{3}p_{2,3}-2N_{4}p_{1,4}+2N_{4}p_{2,4}+N_{3}^{2}p_{3,3}-2N_{3}N_{4}p_{3,4}$

$\displaystyle+$

$\displaystyle N_{4}^{2}p_{4,4}\Big{)}.$

(28)

Evaluating the matrix inversion of $f_{m,k}$ explicitly from Eq. 25 we furthermore obtain

	$\displaystyle\kappa^{2}$	$\displaystyle=$	$\displaystyle\frac{4\pi\ell_{B}}{a^{3}~{}{\cal D}et}\Bigg{(}\phi_{1}\Big{(}1-\phi_{2}+(N_{3}-1)\phi_{3}+(N_{4}-1)\phi_{4}\Big{)}+\phi_{2}\Big{(}1-\phi_{1}+(N_{3}-1)\phi_{3}+(N_{4}-1)\phi_{4}\Big{)}$		(29)
		$\displaystyle+$	$\displaystyle\Big{(}2\phi_{1}\phi_{2}-2N_{3}\phi_{1}\phi_{3}+2N_{3}\phi_{2}\phi_{3}+2N_{4}\phi_{1}\phi_{4}-2N_{4}\phi_{2}\phi_{4}+2\phi_{3}\phi_{4}\Big{)}+N_{3}\phi_{3}\Big{(}1-\phi_{3}+(N_{4}-1)\phi_{4}\Big{)}$
		$\displaystyle+$	$\displaystyle N_{4}\phi_{4}\Big{(}1-\phi_{4}+(N_{3}-1)\phi_{3}\Big{)}\Bigg{)}$

where ${\cal D}et$ is given by

$${\cal D}et=\big{(}1+(N_{3}-1)\phi_{3}+(N_{4}-1)\phi_{4}\big{)}.$$

As in the case of Model 1 in general the inverse square of the screening length is not a linear function of the concentrations.

In order to investigate the case of a bulk electroneutral system we need to evaluate the above expression with electroneutrality condition $c_{1}-c_{2}+N_{3}c_{3}-N_{4}c_{4}=0$ or equivalently

$\displaystyle\phi_{1}-\phi_{2}+\phi_{3}-\phi_{4}=0.$

(30)

The electroneutrality condition, involving now four variables, exerts less of a constraint on the values of the different volume fractions as in the case of Model 1, where we had a system with only two components.

Let us first consider the limit of $\phi_{1}=\phi_{2}=\phi_{s}$ and $\phi_{3}=\phi_{4}=\phi_{M}$, i.e., the salt and the polyions are electroneutral separately, being just a particular case of the general electroneutrality condition.

In this case we obtain for the inverse square of the screening length the expression

	$\displaystyle\kappa^{2}$	$\displaystyle=$	$\displaystyle\frac{4\pi\ell_{B}}{a^{3}~{}{\cal D}et}\Bigg{(}2\phi_{s}\Big{(}1+(N_{3}\!\!+\!\!N_{4}\!\!-\!\!2)\phi_{M}\Big{)}+2N_{3}N_{4}\phi_{M}^{2}$
		$\displaystyle+$	$\displaystyle N_{3}\phi_{M}\Big{(}1+(N_{4}\!\!-\!\!2)\phi_{M}\Big{)}\!\!+\!\!N_{4}\phi_{M}\Big{(}1+(N_{3}\!\!-\!\!2)\phi_{M}\Big{)}\Bigg{)}~{}$
					(31)

with ${\cal D}et=1+(N_{3}+N_{4}-2)\phi_{M}$. The screening length in this case is obviously much more complicated then in the case of the Debye screening length and is partitioned jointly between the simple salt and the polyions.

On Fig. 3 we show the dependence of the screening ratio

$$\tilde{\kappa}^{2}(N_{3},N_{4},\phi_{s},\phi_{M})=\frac{\kappa^{2}(N_{3},N_{4},\phi_{s},\phi_{M})}{(2\phi_{s}+(N_{3}+N_{4})\phi_{M})}$$

(32)

for an electroneutral system. The denominator is again the expected “naive” Debye screening length. The dependence of the screening length is strongly non-monotonic with a minimum at the “naive” Debye value for $N_{1}=N_{2}$ and exhibits pronounced under-screening with $\tilde{\kappa}^{2}\leq 1$. The under-screening depends on the number of monomers the broader is the regime displaying this non-monotonicity.

Figure 4 is the same as Fig. 3 except that we show the dependence of the screening ratio as a function $\phi_{s}$ and $\phi_{M}$. The dependence of the screening ratio on both $\phi_{s}$ and $\phi_{M}$ is monotonic, but is an increasing function in the former and a decreasing function in the latter case, but nevertheless exhibiting under-screening in both cases. This behavior remains valid for different values of the number of charged monomers of the polyions, $N_{3}$ and $N_{4}$. The dependence on $\phi_{M}$ furthermore indicates that the “naive” Debye screening length represents a minimal screening length, and the the system is strongly overscreened.

The Debye-Hückel electrostatic correlation (fluctuation) free energy that enters the Voorn-Overbeek type theories is the only electrostatic contribution to the free energy because the system is considered in the bulk and the homogeneity eliminates any other mean-field electrostastic contribution(s). The fluctuational part of electrostatics, usually calculated in some variant of the one-loop approximation, depends on the screening length which was identified as the standard Debye screening length, based on the assumption that the polyions act as point particles and that the conformational degrees of freedom do not contribute to screening, so that the simple ions as well as the polyions contribute proportionately to the screening. The major assumptions inherent in this approach are the point particle approximation for the polyions and omission of the conformational fluctuations of the polyions, both of which have been discussed and generalized by many later developments (for a recent review see [7]).

There are two separate issues with this approach: one is the expression for the screening length consistent with the other terms in the free energy, which we address in this paper, and the other one is the accuracy of the free energy Ansatz, based on the local density approximation and the lattice gas entropy, to describe the effects of the hard core repulsions, which we do not scrutinize. In this work, our modest goal was not to introduce any new approach or model, but to remain well within the Voorn-Overbeek paradigm, except for the nature of screening and the associated screening length, which we derive self-consistently not by assuming a priori a standard Debye form, but by deriving it explicitly in such a way that it is consistent with all the other terms in the Voorn-Overbeek type free energy Ansatz. This was accomplished by studying the Legendre transform of the Voorn-Overbeek free energy and based on this, calculating the Hessian of electrostatic fluctuations around the mean field as described in detail in Ref. [28]. The fluctuations around the mean field lead to the standard Coulomb fluid correlation free energy, except that the screening length in general differs from the Debye form. The difference can be characterized broadly as under-screening, $\kappa\leq\kappa_{D}$, and thus implies a reduces correlation contribution to the free energy relative to the Voorn-Overbeek theory. This becomes particularly relevant for non-symmetric systems where the polycations and polyanions components are composed of a different number of monomers.

Unlike the case of the ion screening in simple electrolytes, the screening properties of the polyelectrolyte solutions, which are crucial for the onset of complexation phenomenology [43], depend importantly on the polyion chain conformation, so that the ionic screening and the polymer chain conformations self-consistently determine the nature of screening [44, 45]. The naive approach that would simply treat the polymers as finite size particles misses the crucial point of the connection between polymer conformations and screening, whose experimental fingerprint still remains poorly understood [46]. This connection, epitomized by Muthukumar as ”double screening”, couples excluded volume and electrostatic interactions in a polyelectrolyte solution [47] and has been later elaborated by different attempts, most notably the field-theoretic renormalized Gaussian fluctuation theory that accounts self-consistently for the coupling between polyelectrolyte chain conformations and electrostatic interactions in the screening cloud of the polyions [48, 49].

While our contribution is mostly methodological, its usefulness was already demonstrated before in the case of an asymmetric lattice gas [28, 50], and in the case of structural interactions in ionic liquids [32]. We are convinced that it should have important consequences also for various generalizations of the Voorn-Overbeek theory [40, 14, 41, 7, 42, 8] that are often based on complicated expressions, way beyond the ideal gas form, and the corresponding screening lengths cannot be assumed to be of the Debye form. We do not and cannot claim that our methodology can improve on the more fundamental approaches [51, 45, 48], but it does introduce a measure of self-consistency and a possibility for a simple and fast generalization of the electrostatic correlation free energy in the cases where screening is vastly different from the case of the simple salts.

A general conclusion that emerges from this analysis is that the Debye form of the screening length is incompatible with any theory that is not based on an ideal gas entropy term. There are limiting cases, however, such as completely symmetric systems, where even more complicated free energies lead back to the Debye screening length. We give explicit general formulas for the screening length that are valid for any form of the free energy, including the ideal gas free energy, the Flory-Huggins free energy or the Voorn-Overbeek free energy and find that the screening length exhibits a non-monotonic under-screening behavior as a function of the number of monomers of the polyions, and displays a different functional dependence on the volume fractions of the components then the standard Debye screening length. We believe the importance of our methodology is not only to correctly evaluate the proper screening length but also to underline the consistency one needs to strive for in defining it.

The authors would like to acknowledge the support of the 1000-Talents Program of the Chinese Foreign Experts Bureau, as well as the support of the School of physical sciences, University of the Chinese Academy of Sciences, Beijing and the Institute of the physics, Chinese Academy of Sciences, Beijing. RP would like to thank Tony Maggs for his illuminating comments on an earlier version of this manuscript.