text:
# Analytical and numerical solution of coupled KdV-MKdV system

## 1 Introduction

There are two complementary approaches to integrable systems: analytical and numerical ones to be developed. Even most profound analytical IST method cannot give explicit solution of general Cauchy problem while numeric can, but is rather compulsory in use: calculations could need powerful computers. May be most transparent of analytical methods are based on algebraic structures associated with a problem. To such structures belongs Darboux transformations covariance of Lax representation of nonlinear equations that yields a powerful tool for explicit solutions production. We investigate applications of special kind of such discrete symmetry - to be called elementary ones . Its elementarity simply means that a product of such transformations generate the standard one . Here we study combinations of such transforms that do not coincide with binary ones  and hence are not so known.

The main ideas of numerical integration of such integrable systems go up to the famous properties of the equations as the Lax pair and infinite series of conservation laws existence . From a point of view of general theory of such systems some hopes are concerned with a development of the finite-difference or other approximations of the systems. Namely if one could prove a convergence and stability theorems for such difference systems (existence of solutions is implied), a way to existence and uniqueness of solutions is opened .

The coupled KdV-MKdV system arises in many problems of mathematical physics. Some integrable systems are associated with a polynomial spectral problem and have Virasoro symmetry algebras are considered . A dispersive system describing a vector multiplet interacting with the KdV field is a member of a bi-Hamiltonian integrable hierarchy . Recently a multisymplectic numerical twelve points scheme was produced. This scheme is equivalent to the multisymplectic Preissmann scheme and is applied to solitary waves over long time interval . The coupled KdV-MKdV system is also connected to other physical applications .

The general system we consider in this work have the following form

$$\theta _t^n+\underset{m,k}{}\left(g_{m,k}^{n,1}\theta ^m\theta _x^k+g_{m,k}^{n,2}(\theta ^m)^2\theta _x^k+g_{m,k}^{n,3}\theta _x^m\theta _x^k+g_{m,k}^{n,4}\theta ^m\theta _{xx}^k+g_{m,k}^{n,5}\theta ^m(\theta ^k)_x^2\right)+d_n\theta _{xxx}^n=0,$$
(1)
where n, m, k =1,2,…N are the dependent variables numbers. Nonlinear coefficients are $`g_{m,k}^{n,l},l`$ =1,2,..5 and $`d_n`$ are dispersion coefficients.

In particular, for the system under consideration ($`N=3`$) the variables $`\theta ^1,\theta ^2,\theta ^3`$ are denoted by f, u, v to have the integrable system

$$f_t+\frac{1}{2}f_{xxx}+\frac{3}{2}(uf)_x\frac{3}{4}f_xf^2=0$$

$$u_t\frac{1}{4}u_{xxx}\frac{3}{2}u_xu+3vv_x+\frac{3}{4}u_xf^2\frac{3}{2}(f_xv)_x=0$$

$$v_t+\frac{1}{2}v_{xxx}+\frac{3}{2}v_xu\frac{3}{4}(vf^2)_x+\frac{3}{4}u_{xx}f+\frac{3}{2}u_xf_x=0$$
(2)
The Lax pair is given in . The system exhibits two integrable reductions having explicit solutions, Hirota-Satsuma  and a two components KdV-MKdV system . Krishnan  showed that a generalized KdV-MKdV system have solitary wave solutions and investigate the effects of increasing the nonlinearity of one variable on the existence of solitary waves. Some form of KdV-MKdV system have explicit solutions in terms of Jacobi elliptic functions .

In this work we present explicit solutions for a system of three equations (2) that have not been specified in . We study this two-parameter explicit solutions and show effects of choosing these parameters on the solutions. We demonstrate the use of two arbitrary elementary DTs  and its special choice that holds a hereditary of the reduction to built explicit solutions to the KdV-MKdV system (2).

Also we modify a numerical method  for solution of system (1). It is a difference scheme for Cauchy problems for arbitrary number of equations with constants coefficients. The scheme preserves two conservation laws for the KdV type equations and the order of error of the difference formulas is improved . The convergence is proved and stability is analyzed giving the conditions taken in account in choosing time and space step sizes  .

The present work is organized as follows. Section 2 introduces the matrix spectral equation of the second order with 2 x 2 matrix coefficients and two elementary DTs. We select the second equation of the Lax pair and derive the compatibility conditions. The product of these two transformations yields the standard DT . Section 3 illustrates how the, first, elementary DT is used to produce solution to the KdV equation as well as the general evolution equation generated by the compatibility conditions of the Lax pair. Explicit solutions are introduced for the case of zero initial potentials of the matrix problem. In section 4 we consider a reduction constraints on the potential of the matrix spectral equation. This reduction gives an automorphism that relates two pairs of solution of the spectral equation for two spectral parameters. We use this results in the compound elementary DTs to produce an explicit solution to a coupled KdV-MKdV system that results from the compatibility conditions of Lax pair under this reduction. The effects of these parameters on the solution (reality, singularity) is analyzed. Section 5 introduces a numerical method for solving coupled KdV-MKdV system (1). We produce a difference scheme for a Cauchy problem with initial condition rapidly decreasing at both infinities. The main steps of the scheme convergence and stability analysis is shown while the details are explained in appendix A and B. The scheme is tested by applying it to integrable coupled KdV-MKdV system and the numerical results are compared with explicit formulas obtained in section 4.

## 2 Lax pair spectral equations and the elementary DTs

Consider a matrix spectral equation of the second order with spectral parameter $`\lambda `$ and $`2\times 2`$ matrix coefficients.

$$\mathrm{\Psi }_{xx}+F\mathrm{\Psi }_x+U\mathrm{\Psi }=\lambda \sigma _3\mathrm{\Psi }$$
(3)
where the vector $`\mathrm{\Psi }=(\mathrm{\Psi }_1,\mathrm{\Psi }_2)^T`$and the matrix potentials are $`U=\{u_{ij}\},F=\{f_{ij},f_{ii}=0\},i=1,2`$ while $`\sigma _3=diag(1,1)`$ is the Pauli matrix.

For equation (3) we perform two elementary Darboux transforms . The first one is

$$\stackrel{}{\mathrm{\Psi }}_1=\mathrm{\Psi }_{1x}+ϵ_{11}\mathrm{\Psi }_1+ϵ_{12}\mathrm{\Psi }_2,ϵ_{11}=\left(\phi _{1x}+\frac{1}{2}f_{12}\phi _2\right)/\phi _1,ϵ_{12}=f_{12}/2,$$

$$\stackrel{}{\mathrm{\Psi }}_2=\mathrm{\Psi }_2+ϵ_{21}\mathrm{\Psi }_1,ϵ_{21}=\phi _2/\phi _1,$$

$$\stackrel{}{\phi }_1=\left(_x+ϵ_{11}\right)\phi _3+ϵ_{12}\phi _4,\stackrel{}{\phi }_2=\phi _4+ϵ_{21}\phi _3.$$
(4)
where $`(\phi _1,\phi _2)^T`$ and $`(\phi _3,\phi _4)^T`$ are two solutions of (3) corresponding to different spectral parameters.

Substituting the above expressions for $`\stackrel{}{\mathrm{\Psi }}_1,\stackrel{}{\mathrm{\Psi }}_2`$into (3) and collecting the coefficients of $`\mathrm{\Psi }_1`$, $`\mathrm{\Psi }_2`$and their derivatives we obtain the expressions for new potentials as

$$\stackrel{}{f}_{12}=u_{12}+f_{12}ϵ_{11},$$

$$\stackrel{}{f}_{21}=2ϵ_{21},$$

$$\stackrel{}{u}_{11}=u_{11}2ϵ_{11x}\stackrel{}{f}_{12}ϵ_{21}f_{21}ϵ_{12},$$

$$\stackrel{}{u}_{12}=u_{12x}ϵ_{12xx}+ϵ_{11}u_{12}ϵ_{12}\left(\stackrel{}{u}_{11}+u_{22}\right),$$

$$\stackrel{}{u}_{21}=f_{21}2ϵ_{21x}\stackrel{}{f}_{21}ϵ_{11},$$

$$\stackrel{}{u}_{22}=u_{22}ϵ_{21}u_{12}\stackrel{}{u}_{21}ϵ_{12}\stackrel{}{f}_{21}ϵ_{12x}.$$
(5)
The second elementary DT is performed after the first one and can be obtained by reversing the indices 1$`2`$and 2$`1`$to get, for example, the following potentials

$$\stackrel{}{f}_{21}=\stackrel{}{u}_{21}+\stackrel{}{f}_{21}\stackrel{}{ϵ}_{22},\stackrel{}{ϵ}_{22}=\left(\stackrel{}{\phi }_{2x}+\frac{1}{2}\stackrel{}{f}_{21}\stackrel{}{\phi }_1\right)/\stackrel{}{\phi }_2,$$

$$\stackrel{}{u}_{22}=\stackrel{}{u}_{22}2\stackrel{}{ϵ}_{22x}\stackrel{}{f}_{21}\stackrel{}{ϵ}_{12}\stackrel{}{f}_{12}\stackrel{}{ϵ}_{21},\stackrel{}{ϵ}_{12}=\stackrel{}{\phi }_1/\stackrel{}{\phi }_2,\stackrel{}{ϵ}_{21}=\stackrel{}{f}_{21}/\mathrm{\hspace{0.17em}2}.$$
(6)
The spectral equation (3) is considered as the first equation of the Lax pair, take the second as

$$\mathrm{\Psi }_t=\mathrm{\Psi }_{xxx}+B\mathrm{\Psi }_x+C\mathrm{\Psi },$$
(7)
where $`B=\frac{3}{2}diagU+\frac{3}{2}F_x+\frac{3}{4}F^2`$
and $`C=\frac{3}{2}U_x\frac{3}{4}diagU_x\frac{3}{4}(f_{12}u_{21}+f_{21}u_{12})I+\frac{3}{8}(f_{12,x}f_{21}f_{12}f_{21,x})\sigma _3+\frac{3}{4}(u_{11}u_{22})\sigma _3F.`$
Equation (7) is also covariant under transformations (4), (5). The compatibility conditions have the following form

$$F_tF_{3x}+B_{2x}3U_{2x}+2C_x+FB_x\sigma _3B\sigma _3F_x+UB$$

$$\sigma _3B\sigma _3U+FC\sigma _3C\sigma _3F=0,$$

$$U_tU_{3x}+C_{2x}+UC\sigma _3C\sigma _3U+FC_x\sigma _3B\sigma _3U_x=0.$$
(8)
and the transformations (4), (5) determine a discrete symmetry of (2)

## 3 Solution of two coupled KdV-MKdV equations and KdV equation via the first elementary DT

For a spectral parameter $`\lambda `$ and a seed potential F, U we obtain the solutions $`\phi _1,\phi _2`$ to the pair (3), (7). Then performing the first elementary DT to obtain the new potentials $`\stackrel{}{F},\stackrel{}{U}`$ which are solutions to the system (2). For the case of zero seed potential the solutions, $`\phi _1`$ and $`\phi _2`$ of the system $`(\text{3}),(\text{7})`$ have the form

$$\phi _1=c_1e^{ax+a^3t}+c_2e^{(ax+a^3t)},$$

$$\phi _2=d_1e^{iax+(ia)^3t}+d_2e^{(iax+(ia)^3t)}.$$
(9)
where $`c_1,c_2,d_1,d_2`$are arbitrary constants, $`a=\sqrt{\lambda }`$ and i is the imaginary unit. System (2) reduced (for the only nonzero elements ) to the following

$$f_{21t}+\frac{1}{2}f_{21xxx}+\frac{3}{4}f_{21}u_{11x}=\frac{3}{2}u_{11}u_{21},$$

$$u_{11t}\frac{1}{4}u_{11xxx}\frac{3}{2}u_{11}u_{11x}=0,$$

$$u_{21t}+\frac{1}{2}u_{21xxx}+\frac{3}{4}u_{21}u_{11x}+\frac{3}{2}u_{21x}u_{11}=\frac{3}{4}u_{11}f_{21xx}+\frac{3}{4}u_{11}^2f_{21}.$$
(10)
where $`f_{12}=0,u_{12}=0,u_{22}=0`$ and tildes are omitted for simplicity. This system with explicit solution obtained from (4), (5) as

$$f_{21}=\left(2e^{\left(1i\right)a\left(a^2t+x\right)}\left(d_2e^{2ia^3t}+d_1e^{2iax}\right)\right)/\left(c_2+c_1e^{2a\left(a^2t+x\right)}\right),$$

$$u_{11}=\left(8a^2c_1c_2e^{2a\left(a^2t+x\right)}\right)/\left(c_2+c_1e^{2a\left(a^2t+x\right)}\right)^2,$$

$$u_{21}=\left(2iae^{\left(1i\right)a\left(a^2t+x\right)}\left(d_2e^{2ia^3t}d_1e^{2iax}\right)\right)/\left(c_2+c_1e^{2a\left(a^2t+x\right)}\right).$$
(11)
where $`c_1,c_2,d_1,d_2`$are arbitrary constants and $`a=\sqrt{\lambda }`$.

The second equation in (10) is the KdV equation while the remaining are a two components coupled KdV-MKdV system that was solved by elementary DT.

## 4 Solution of three coupled KdV-MKdV equations via the compound elementary DTs

Existence of different kinds of automorphism causes special constraints . Multiplying (3) by $`\sigma _1=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)`$ to have

$$\sigma _1\mathrm{\Psi }_{xx}+\sigma _1F\mathrm{\Psi }_x+\sigma _1U\mathrm{\Psi }=\lambda \sigma _1\sigma _3\mathrm{\Psi }$$
(12)
but $`\sigma _1\sigma _3=\sigma _3\sigma _1`$ and consider the conditions $`\sigma _1F=F\sigma _1`$ and $`\sigma _1U=U\sigma _1`$ that means

$$f_{12}=f_{21}=f,u_{11}=u_{22}=u,u_{12}=u_{21}=v.$$
(13)
So (12) becomes

$$\left(\sigma _1\mathrm{\Psi }\right)_{xx}+F\left(\sigma _1\mathrm{\Psi }\right)_x+U\left(\sigma _1\mathrm{\Psi }\right)=\lambda \sigma _3\left(\sigma _1\mathrm{\Psi }\right)$$

The above automorphism $`\mathrm{\Psi }\left(\lambda \right)\sigma _1\mathrm{\Psi }\left(\lambda \right)`$ relates two pair of solutions $`(\phi _1,\phi _2)`$ and $`(\phi _3,\phi _4)`$ of (3) corresponding to different values of spectral parameter $`\lambda ,\lambda `$ as

$$\left(\begin{array}{c}\phi _3(\lambda )\\ \phi _4(\lambda )\end{array}\right)=\sigma _1\left(\begin{array}{c}\phi _1(\lambda )\\ \phi _2(\lambda )\end{array}\right)=\left(\begin{array}{c}\phi _2(\lambda )\\ \phi _1(\lambda )\end{array}\right)$$

Using this result in the elementary DTs (4), (5) and (6) to obtain the expressions for the new potentials f, u, v. In the case of zero initial potentials these new potentials have the following forms

$$f=2\frac{\phi _1(\phi _2)_x\phi _2(\phi _1)_x}{(\phi _1)^2(\phi _2)^2},$$

$$u=\left(\frac{(\phi _1^2)_x(\phi _2^2)_x}{(\phi _1)^2(\phi _2)^2}\right)_x+2\left(\frac{\phi _1(\phi _2)_x\phi _2(\phi _1)_x}{(\phi _1)^2(\phi _2)^2}\right)^2,$$

$$v=2\left(\frac{\phi _1(\phi _2)_x\phi _2(\phi _1)_x}{(\phi _1)^2(\phi _2)^2}\right)_x+\frac{\left(\phi _1(\phi _2)_x\phi _2(\phi _1)_x\right)\left((\phi _1^2)_x(\phi _2^2)_x\right)}{\left((\phi _1)^2(\phi _2)^2\right)^2}.$$
(14)
where $`\phi _1,\phi _2`$ are as in (9) with $`c_1,c_2,d_1,d_2`$are arbitrary constants and $`a=\sqrt{\lambda }`$.
The above expressions are solutions of system (2) that reduced under the reduction conditions (13) to system (2).

The choice of the arbitrary constants ($`c_1,c_2,d_1,d_2`$) affects on the behavior of the solution in formula (4). For example choosing equal constants $`c_1=c_2=d_1=d_2=0.5`$ (we choose the value to be $`0.5`$ to simplify the resulting formula but the idea valid for any value) the solutions have the form

$$f=2a(sin\eta _1cosh\eta _2cos\eta _1sinh\eta _2)/(cosh^2\eta _2cos^2\eta _1),$$

$$u=2a^2(sin\eta _1cosh\eta _2+cos\eta _1sinh\eta _2)^2/(cosh^2\eta _2cos^2\eta _1)^2,$$

$$v=2a^2(cos3\eta _1cosh\eta _22sin\eta _1sinh\eta _2(cos2\eta _1+cosh2\eta _2+2)cos\eta _1cosh3\eta _2)$$

$$/(cosh^2\eta _2cos^2\eta _1)^2.$$
(15)
where $`\eta _1=a^3tax,\eta _2=a^3t+ax,a=\sqrt{\lambda }`$ is real.

We see that the above expression (15) is singular at $`\eta _2=0,\eta _1=n\pi ,n=0,1,2,\mathrm{}`$. Hence we have singularity at $`(x=\frac{n\pi }{2a},t=\frac{n\pi }{2a^3})`$.

To obtain continuous solutions we can choose $`c_1=c_2,d_1=d_2=r.c_1`$, r is real constant. We again choose $`c_1=0.5`$ following the previous concept. So (4) have the form

$$f=2ar\left(\mathrm{cosh}\eta _2\mathrm{sin}\eta _1\mathrm{cos}\eta _1\mathrm{sinh}\eta _2\right)/\left(\mathrm{cosh}^2\eta _2r^2\mathrm{cos}^2\eta _1\right),$$

$$u=a^2\left(1r^4r^4\mathrm{cos}2\eta _1+\mathrm{cosh}2\eta _2+r^2\mathrm{sin}2\eta _1\mathrm{sinh}2\eta _2\right)/\left(\mathrm{cosh}^2\eta _2r^2\mathrm{cos}^2\eta _1\right)^2,$$

$$v=2a^2r(((7+6r^2+2r^2\mathrm{cos}2\eta _1)\mathrm{cos}\eta _1\mathrm{cosh}\eta _2\mathrm{cos}\eta _1\mathrm{cosh}3\eta _2)$$

$$2(1+r^2+r^2\mathrm{cos}2\eta _1+\mathrm{cosh}2\eta _2)\mathrm{sin}\eta _1\mathrm{sinh}\eta _2))/(1+r^2+r^2\mathrm{cos}2\eta _1\mathrm{cosh}2\eta _2)^2$$
(16)
where $`a,\eta _1,\eta _2`$ as in (15)

Choosing this parameter ($`r`$) to be $`r`$ $`<1`$ gives real nonsingular solutions. The above formula (4) is built from elliptic and periodic functions so it does not preserve its symmetry but its localized as shown in figures (1.a,b) below.

Fig.(1.a) Non-singular solutions, f, u and v (r=0.5) , a=2, t=0.

Fig.(1.b) Propagation of solutions, f, u and v (r=0.5) , a=2, t=1.

Fig.(1) The solutions in (4) does not preserve its symmetry but its localized.

Choosing the parameter r to be $`r`$ $`>1`$ in formula (4) gives singular solutions as shown figure (2) below.

Fig.(2) Singular solutions, f, u and v (r=2), a=2, t=0.

Moreover the choice of these arbitrary constants ($`c_1,c_2,d_1,d_2`$) as well as the spectral parameter $`\lambda `$ affects the reality of the resulting solution. As example for $`\lambda =2im^2`$, m is real and choosing $`c_1=c_2=d_1=d_2=0.5`$, we get real solution

$$f=m(cos2\zeta _1sinh\zeta _2sinh\zeta _2sin\zeta _1cosh2\zeta _2+sin\zeta _1)/(0.25cosh2\zeta _20.25)(1cos2\zeta _1)$$

where $`\zeta _1=2mx+4m^3t,\zeta _2=2mx4m^3t`$, while choosing $`c_1=c_2=1,d_1=d_2=2`$ give the following complex solution

$$f=m(8(5sinh\zeta _2cos2\zeta _1+5sinh\zeta _2+5sin\zeta _1cosh2\zeta _25sin\zeta _1)8i(6sinh\zeta _2$$

$$+3sin2\zeta _1cosh\zeta _2+3cos\zeta _1sinh2\zeta _2+6sin\zeta _1))/(17cosh2\zeta _2+10+36cos\zeta _1cosh\zeta _2$$

$$8cos2\zeta _1cosh2\zeta _2+17cos2\zeta _1)$$

## 5 The numerical method

### 5.1 The difference scheme

For the coupled KdV-MKdV system (1) we introduce a numerical (finite difference) method of solution . This scheme is valid for arbitrary number of equations with constants coefficients and of the form

$$\frac{\theta _i^{n,j+1}\theta _i^{n,j}}{\tau }+\underset{m,k}{}(g_{m,k}^{n,1}\theta _i^{m,j}\frac{\theta _{i+1}^{k,j}\theta _{i1}^{k,j}}{2h}+g_{m,k}^{n,2}(\theta _i^{m,j})^2\frac{\theta _{i+1}^{k,j}\theta _{i1}^{k,j}}{2h}$$

$$+g_{m,k}^{n,3}\frac{\theta _{i+1}^{m,j}\theta _{i1}^{m,j}}{2h}\frac{\theta _{i+1}^{k,j}\theta _{i1}^{k,j}}{2h}+g_{m,k}^{n,4}\theta _i^{m,j}\frac{\theta _{i+1}^{k,j}2\theta _i^{k,j}+\theta _{i1}^{k,j}}{2h}$$

$$+g_{m,k}^{n,5}\theta _i^{m,j}\theta _i^{k,j}\frac{\theta _{i+1}^{k,j}\theta _{i1}^{k,j}}{2h})+d_n\frac{\theta _{i+2}^{n,j}2\theta _{i+1}^{n,j}+2\theta _{i1}^{n,j}\theta _{i2}^{n,j}}{2h^3}=0$$
(17)
where i and j are the discrete space and time respectively. The time step is denoted by $`\tau `$while h denotes spatial step.

### 5.2 Stability analysis of the scheme

We prove stability with respect to small perturbations of initial conditions . It is the boundness of the discrete solution with respect to small perturbation of the initial data. We give here the main steps while the details are presented in Appendix A. We can write

$$d\theta _i^{n,j+1}=T_{i,r}^{n,j+1}d\theta _r^{n,j}=T_{i,r}^{n,j+1}T_{i,r}^{n,j}d\theta _r^{n,j1}=\underset{𝑟}{\mathrm{\Pi }}\left(T_{i,r}^n\right)^rd\theta _r^{n,o}$$

where $`d\theta _i^{n,j+1}`$ is perturbations of the discrete solution, $`d\theta _r^{n,o}`$ small perturbation of the initial data and $`T_{i,r}^{n,j+1}`$ is a differentiable operator. Stability required the boundedness of $`\underset{𝑟}{\mathrm{\Pi }}\left(T_{i,r}^n\right)^r`$i.e $`T^r`$is bounded.

We found that

$$T^{j+1}^2e^{a(\tau ,h)\tau },$$

$$a(\tau ,h)=2\underset{l,n,m,k}{\mathrm{max}}\left|g_{m,k}^{n,l}\right|\underset{i,m,k}{\mathrm{max}}\left(\right|\theta _{x,i}^{m,j}\left|\right|\theta _{x,i}^{k,j}\left|\right)+\tau [\underset{l,n,m,k}{\mathrm{max}}|g_{m,k}^{n,l}|\underset{i,m,k}{\mathrm{max}}\left(\right|\theta _{x,i}^{m,j}\left|\right|\theta _{x,i}^{k,j}\left|\right)$$

$$+\frac{1}{h}\underset{l,n,m,k}{\mathrm{max}}\left|g_{m,k}^{n,l}\right|\underset{i,m,k}{\mathrm{max}}\left(\right|\theta _i^{m,j}\left|\right|\theta _i^{k,j}\left|\right)+\frac{3}{h^3}\underset{𝑛}{\mathrm{max}}\left|d_n\right|]^2$$
(18)
The scheme is stable if $`a(\tau ,h)const`$. We have here a conditional stability. That is we require that $`\tau 0`$ more faster than $`h0`$. Namely we need

$$\tau \left(\text{constant}\right).h^6$$

### 5.3 Convergence proof for the scheme

We prove that the solution of (17) converges to the solution of (1) if the exact solution is continuously differentiable one . We introduce here the main points for the scheme convergence and give the details in Appendix B.

$`\theta _i^j`$ is the difference solution of (17), u$`{}_{}{}^{j}{}_{i}{}^{}`$ is the exact solution. Hence the error v$`{}_{}{}^{j}{}_{i}{}^{}`$ is given by $`v_i^j=\theta _i^ju_i^j`$. Introducing L$`_2`$norm defined by $`V^j=\left(\underset{𝑖}{}\underset{𝑛}{}\left(v_i^{n,j}\right)^2h\right)^{1/2}`$

The scheme converges when the norm of that error $`V^j0`$ as $`\left(\tau ,\text{ }h\text{ }0\right)`$
We found that $`V^{j+1}P(M)O\left(\tau +h^2\right)`$, where $`P(M)`$is a polynomial in the bounded constant $`M=\tau \frac{e^{a\tau j}1}{e^{a\tau }1}`$ and $`a`$ as in (18). Hence the convergence proved.

### 5.4 Numerical calculations and test

The coupled KdV-MKdV system (2) is solved numerically using scheme (17) with initial condition from (4) at $`t=0`$ and the results are compared with the explicit formulas (4). The percentage errors are shown in the following plots.

Fig.3 percentage errors of the numerical solutions relative to the explicit solutions.

The results of the test confirms the validity of the numerical scheme we propose. It also illustrates the errors of evaluation that could be estimated by the resulting inequalities of the scheme convergence proof.

## 6 Conclusion

Darboux transformations covariance of Lax representation of nonlinear equations is a powerful tool for explicit solutions production. Here we investigate applications of special kind of such discrete symmetry - to be called elementary ones. We use these elementary DT to produce explicit solutions for coupled KdV-MKdV system. The iteration of DT can be formulated in form of determinants representations . A numerical method for general coupled KdV-MKdV system is introduced. It is a difference scheme for Cauchy problems for arbitrary number of equations with constants coefficients. We analyze stability and prove the convergence of the scheme. The scheme keeps two conservation laws chosen in analogy with KdV type equations. Analyzing stability and proving the convergence beside comparing the numerical results with explicit formulas allow us to use the numerical scheme to systems with arbitrary coefficients that is presumably non-integrable. Obviously the coupled KdV systems are successfully treated by our scheme .

Acknowledgment  We thank S. B. Kshevetskii for useful discussions about numerical scheme for the problem under consideration.

## 7 Appendix A <br>Stability analysis of the scheme

We prove stability with respect to small perturbations (because we consider nonlinear equations) of initial conditions. Strictly speaking it is the boundness of the discrete solution in terms of small perturbation of the initial data. Consider the differential

$`T_{i,r}^{n,j+1}=\left\{\theta _i^{n,j+1}/\theta _r^{n,j}\right\},d\theta _r^{n,j}=\left\{\theta _{i2}^{n,j}\theta _{i1}^{n,j}\theta _i^{n,j}\theta _{i+1}^{n,j}\theta _{i+2}^{n,j}\right\}^t`$ and define the norm $`d\theta ^j=\left(\underset{𝑟}{}\underset{𝑛}{}\left(d\theta _r^{n,j}\right)^2h\right)^{1/2}`$

We can write $`d\theta _i^{n,j+1}=T_{i,r}^{n,j+1}d\theta _r^{n,j}=T_{i,r}^{n,j+1}T_{i,r}^{n,j}d\theta _r^{n,j1}=\underset{𝑟}{\mathrm{\Pi }}\left(T_{i,r}^n\right)^rd\theta _r^{n,o}`$

where $`d\theta _i^{n,j+1}`$ is perturbations of the discrete solution, $`d\theta _r^{n,o}`$ small perturbation of the initial data. Stability required the boundedness of $`\underset{𝑟}{\mathrm{\Pi }}\left(T_{i,r}^n\right)^r`$i.e $`T^r`$is bounded. We calculate $`T`$from (17) as follow

$$T_{i,r}^{n,j+1}=\delta _{i,r}\tau \underset{m,k}{}(\frac{g_{m,k}^{n,1}}{2h}[\theta _i^{m,j}(\delta _{i+1,r}\delta _{i1,r})+\delta _{i,r}(\theta _{i+1}^{k,j}\theta _{i1}^{k,j})]$$

$$+\frac{g_{m,k}^{n,2}}{2h}\left[\left(\theta _i^{m,j}\right)^2\left(\delta _{i+1,r}\delta _{i1,r}\right)+2\theta _i^{m,j}\delta _{i,r}\left(\theta _{i+1}^{k,j}\theta _{i1}^{k,j}\right)\right]$$

$$+\frac{g_{m,k}^{n,3}}{2h}\left[\left(\theta _{i+1}^{m,j}\theta _{i1}^{m,j}\right)\left(\delta _{i+1,r}\delta _{i1,r}\right)+\left(\delta _{i+1,r}\delta _{i1,r}\right)\left(\theta _{i+1}^{k,j}\theta _{i1}^{k,j}\right)\right]$$

$$+\frac{g_{m,k}^{n,4}}{2h}\left[\theta _i^{m,j}\left(\delta _{i+1,r}2\delta _{i,r}+\delta _{i1,r}\right)+\delta _{i,r}\left(\theta _{i+1}^{k,j}2\theta _i^{k,j}+\theta _{i1}^{k,j}\right)\right]$$

$$+\frac{g_{m,k}^{n,5}}{2h}[\theta _i^{m,j}\theta _i^{k,j}(\delta _{i+1,r}\delta _{i1,r})+\theta _i^{m,j}(\theta _{i+1}^{k,j}\theta _{i1}^{k,j})\delta _{i,r}+\theta _i^{k,j}(\theta _{i+1}^{k,j}\theta _{i1}^{k,j})\delta _{i,r}])$$

$$\frac{\tau d_n}{2h^3}\left[\delta _{i+2,r}2\delta _{i+1,r}+2\delta _{i1,r}2\delta _{i2,r}\right]$$
(19)
Rewriting (19) in terms of identity (E), symmetric (S) and anti-symmetric (A) matrices

$`S^{j+1}\tau \underset{l,n,m,k}{\mathrm{max}}\left|g_{m,k}^{n,l}\right|\underset{i,m,k}{\mathrm{max}}\left(\left|\theta _{x,i}^{m,j}\right|\left|\theta _{x,i}^{k,j}\right|\right),`$

$`A^{j+1}\frac{\tau }{h}\underset{l,n,m,k}{\mathrm{max}}\left|g_{m,k}^{n,l}\right|\underset{i,m,k}{\mathrm{max}}\left(\left|\theta _i^{m,j}\right|\left|\theta _i^{k,j}\right|\right)+\frac{3}{h^3}\underset{𝑛}{\mathrm{max}}\left|d_n\right|,`$

where $`\theta _{x,i}^j=\frac{\theta _{i+1}^j\theta _{i1}^j}{2h},n,m,k=1,2,..N,l=1,2,\mathrm{}5.`$

$`T^{j+1}^2=\left(T^{j+1}\right)^{}T^{j+1}=\left(EA^{j+1}+S^{j+1}\right)\left(E+A^{j+1}+S^{j+1}\right)`$

$`1+2S^{j+1}+\left(A^{j+1}+S^{j+1}\right)^2`$

$`e^{a(\tau ,h)\tau },`$

$$a(\tau ,h)=2\underset{l,n,m,k}{\mathrm{max}}\left|g_{m,k}^{n,l}\right|\underset{i,m,k}{\mathrm{max}}\left(\right|\theta _{x,i}^{m,j}\left|\right|\theta _{x,i}^{k,j}\left|\right)+\tau [\underset{l,n,m,k}{\mathrm{max}}|g_{m,k}^{n,l}|\underset{i,m,k}{\mathrm{max}}\left(\right|\theta _{x,i}^{m,j}\left|\right|\theta _{x,i}^{k,j}\left|\right)$$

$$+\frac{1}{h}\underset{l,n,m,k}{\mathrm{max}}\left|g_{m,k}^{n,l}\right|\underset{i,m,k}{\mathrm{max}}\left(\right|\theta _i^{m,j}\left|\right|\theta _i^{k,j}\left|\right)+\frac{3}{h^3}\underset{𝑛}{\mathrm{max}}\left|d_n\right|]^2$$

We have here a conditional stability. That is we require that $`\tau 0`$ more faster than $`h0.`$Namely we need

$$\tau \left(\text{cons}\mathrm{tan}\text{t}\right).h^6$$

## 8 Appendix B <br>The scheme convergence

We prove the convergence by proving that the norm of the error (between the difference solution and the exact solution) vanishes as the mesh is refined. Let $`\theta _i^j`$the difference solution of (17), u$`{}_{}{}^{j}{}_{i}{}^{}`$the exact solution. The error v$`{}_{}{}^{j}{}_{i}{}^{}`$ is given by $`v_i^j=\theta _i^ju_i^j.`$
The scheme converges when the norm $`V^j0`$ as $`\left(\tau ,\text{ }h\text{ }0\right)`$ where the norm is defined as $`V^j=\left(\underset{𝑖}{}\underset{𝑛}{}\left(v_i^{n,j}\right)^2h\right)^{1/2}`$

substitute in (17) by $`\theta _i^j=v_i^j+u_i^j`$keeping in mind that for $`\theta _i^j`$equation (17) is $`O\left(\tau +h^2\right)`$and using the operator T defined by

$$v_i^{n,j}\tau [\underset{m,k}{}(g_{m,k}^{n,1}(u_i^{m,j}\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}+v_i^{m,j}\frac{u_{i+1}^{k,j}u_{i1}^{k,j}}{2h})+g_{m,k}^{n,2}((u_i^{m,j})^2\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}$$

$$+(v_i^{m,j})^2\frac{u_{i+1}^{k,j}u_{i1}^{k,j}}{2h})+g_{m,k}^{n,3}(\frac{u_{i+1}^{m,j}u_{i1}^{m,j}}{2h}\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}+\frac{v_{i+1}^{m,j}v_{i1}^{m,j}}{2h}\frac{u_{i+1}^{k,j}u_{i1}^{k,j}}{2h})$$

$$+g_{m,k}^{n,4}(u_i^{m,j}\frac{v_{i+1}^{k,j}2v_i^{k,j}+v_{i1}^{k,j}}{2h}+v_i^{m,j}\frac{u_{i+1}^{k,j}2u_i^{k,j}+u_{i1}^{k,j}}{2h})+2g_{m,k}^{n,5}(u_i^{m,j}v_i^{k,j}$$

$$\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}+v_i^{m,j}u_i^{k,j}\frac{u_{i+1}^{k,j}u_{i1}^{k,j}}{2h}))+d_n\frac{v_{i+2}^{n,j}2v_{i+1}^{n,j}+2v_{i1}^{n,j}v_{i2}^{n,j}}{2h^3}]$$

$$=\underset{𝑟}{}T_{ir}^{j+1}v_r^{n,j}$$

So we obtain

$$v_i^{n,j+1}=\underset{𝑟}{}T_{ir}^{j+1}v_r^{n,j}+\tau f_{m,k,i}^{n,j}$$
(20)
where $`f_{m,k,i}^{n,j}=\underset{m,k}{}g_{m,k}^{n,1}v_i^{m,j}\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}+g_{m,k}^{n,2}(v_i^{m,j})^2\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}+g_{m,k}^{n,3}\frac{v_{i+1}^{m,j}v_{i1}^{m,j}}{2h}\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}+g_{m,k}^{n,4}v_i^{m,j}\frac{v_{i+1}^{k,j}2v_i^{k,j}+v_{i1}^{k,j}}{2h}+2g_{m,k}^{n,5}v_i^{m,j}v_i^{k,j}\frac{v_{i+1}^{k,j}v_{i1}^{k,j}}{2h}+O\left(\tau +h^2\right)`$

$`f^j=\left(\underset{𝑖}{}\left(f_{m,k,i}^{n,j}\right)^2h\right)^{1/2}`$
$`\frac{\left|g_{m,k}^{n,1}\right|_{\mathrm{max}}}{h^{3/2}}V^j^2+\frac{\left|g_{m,k}^{n,2}\right|_{\mathrm{max}}}{h^2}V^j^3+\frac{\left|g_{m,k}^{n,3}\right|_{\mathrm{max}}}{h^{5/2}}V^j^2+\frac{\left|g_{m,k}^{n,4}\right|_{\mathrm{max}}}{h^{5/2}}V^j^2`$
$`+\frac{\left|g_{m,k}^{n,5}\right|_{\mathrm{max}}}{h^2}V^j^3+O\left(\tau +h^2\right)`$

$`\frac{\left|g_{m,k}^{n,l}\right|_{\mathrm{max}}}{h^2}V^j^3+O\left(\tau +h^2\right),\left|g_{m,k}^{n,l}\right|_{\mathrm{max}}=\underset{n,m,k,1}{\mathrm{max}}g_{m,k}^{n,l}`$

Using Schwartz inequality so (20) becomes

$`V^{j+1}T^{j+1}V^j+\tau f^j`$

$`T^{j+1}T^jV^{j1}+\tau \left(T^{j+1}f^{j1}+f^j\right)`$

$`e^{a\tau j}V^o+\tau \left(e^{a\tau \left(j1\right)}f^o+e^{a\tau \left(j2\right)}f^1+\mathrm{}+f^j\right)`$

$`e^{a\tau j}V^o+M\left|g_{m,k}^{n,l}\right|_{\mathrm{max}}V^{j+1}^3+MO\left(\tau +h^2\right),M=\tau \frac{e^{a\tau j}1}{e^{a\tau }1}`$

Using $`V^o=0,`$ the above inequality has the solution

$`V^{j+1}P\left(M\right)O\left(\tau +h^2\right),P(M)`$ is a polynomial in $`M`$.
Since $`M`$ is bounded then $`V^{j+1}0as\tau ,h0`$and the convergence proved.


format: text

text:
# Virial expansion for charged colloids and electrolytes

## 1 Introduction

The two-component hard-core plasma (TCPHC) has been used for a long time as an idealized model for electrolyte solutions. In this model, also known as “primitive model,” the ions are spherical particles that interact with each other via the Coulomb potential and a hard-core potential, which avoids the collapse of oppositely charged particles onto each other. In the general asymmetric case, the positive ions have a charge $`q_+e`$ (where $`e`$ is the elementary charge) and ionic diameter $`d_+`$, while the negative ions have a charge $`q_{}e`$ and diameter $`d_{}`$. The particles are immersed in a structureless solvent whose presence is felt only through the value of the dielectric constant of the medium, and the system is (globally) electroneutral.

The TCPHC in this formulation is a gross simplification of real systems. When two ions are at distances of the order of the size of the solvent molecules, the assumption of a continuum solvent breaks down, giving rise to so-called solvation forcesisraelachvili:book . Such effects also lead to non-local contributions to the dielectric constant. Since the dielectric constant of ions or colloids typically differs from the surrounding aqueous medium, one also expects dispersion forces to act between these particles. All these effects, which contribute to “ion-specific” effectsninham-yaminsky:97 , and are possibly more relevant than previously thoughtalfridsson-al:2000 , are not accounted for in the TCPHC, which treats the solvent—usually water—as a structureless medium within which the charged particles are embedded, neglecting the molecular arrangement that occurs around the ions.

Nevertheless, even with such simplifications, the TCPHC is far from being amenable to an exact treatment. A better understanding of this model is a necessary step if one wishes to develop more realistic approaches to charged systems. In this article, we turn our attention to the low-density, or virial, expansion of the TCPHC. Since we have exact results for the thermodynamic properties at low concentrations using field-theoretic methods, we can obtain useful information on dilute systems with otherwise arbitrary hard-core radii and charge valences (like colloidal suspensions). In particular, we can test various approximations to treat the TCPHC, like the mapping of a colloidal solution on an effective one-component plasma.

Due to the long-range character of the Coulomb potential, it is not easy to obtain the thermodynamic behavior of the TCPHC through the usual methods of statistical mechanics. For example, it can be shownmayer:50 ; friedman:book ; mcquarrie:book  that the straightforward application of the cluster expansion to the TCPHC leads to divergent virial coefficients. Mayermayer:50  proposed a solution to this problem through an infinite resummation of the cluster diagrams, carried out such that the divergent contributions to the virial expansion are canceled. With this, he was able to obtain explicitly the first term in the virial (or low-density) expansion that goes beyond the ideal gas, which turns out to be the well-known Debye-Hückel limiting lawdebye-hueckel:23 . Hagahaga:53  carried the expansion further and went up to order $`5/2`$ in the ionic density. More or less at the same time, Edwardsedwards:59  also obtained the virial expansion of the TCPHC by mixing cluster expansion and field theory. In both cases only equally-sized ions were considered.

The aforementioned methods typically depend on drawing, counting and recollecting the cluster diagrams which give finite contributions to the expansion up to the desired order in the density. This can be quite a formidable task, and unfortunately it is easy to “forget” diagrams that are relevant to the series (see for instance the comment on pp. 222–223 of Ref. friedman:book ). Also, the generalization to ions with different sizes is rather complicatedfriedman:book . Besides, the final results are typically not obtained in closed form, i.e., the final expressions depend on infinite sums that usually have to be evaluated numerically, which reflects the infinite diagrammatic resummation.

We generalize here a novel field-theoretic techniqueroland-orland:tcp , introduced for the symmetric TCPHC ($`q_+=q_{}`$ and $`d_+=d_{}`$). We obtain the exact low-density expansion of the *asymmetric* (both in size and charge) TCPHC. This method does not use the cluster expansion (and resummation) and yields analytic, closed-form results. We go up to order $`5/2`$ in the volume fraction of a system where the sizes and the charge valences of positive and negative ions are unconstrained, that is, the results we obtain can be applied, without modifications, to both electrolyte solutions (where anions and cations have approximately the same size and valence) and to colloidal suspensions (where the macro- and counterions have sizes and valences that can be different by orders of magnitude). In the colloidal limit, we demonstrate how to map the TCPHC on an effective one-component hard-core plasma (OCPHC) and obtain the corrections due to the exclusion of the background (i.e. the counterions) from the colloidal particles. For electrolytes, we obtain effective ionic sizes in solution, or thermodynamic diameters, from experimental data for the mean activity. The method developed here allows, in principle, to measure the diameters of cations and anions as *independent* variables, given that the range of experimental data extends to low enough densities such that the expansion used here is valid. It should be noted that the diameters obtained in this way are fundamentally different from the hydrodynamic diameters.

This article is organized as follows. In section 2 we describe in detail the steps that lead to the low-density expansion. Readers which are not interested in this derivation can go directly to section 3, where we apply the expressions obtained to two particular problems, viz., (i) colloids, where one of the charged species is much smaller and much less charged than the other one and (ii) the mean activity coefficient (which is related to the exponential of the chemical potential) of electrolyte solutions, which is available from experiments and can be compared to our results. Finally, section 4 contains some concluding remarks.

## 2 The method

We begin our calculation by assuming a system with $`N_+`$ positively charged particles with charge valence $`q_+`$ and diameter $`d_+`$, and $`N_{}`$ negatively charged particles with charge valence $`q_{}`$ and diameter $`d_{}`$. Global electroneutrality of the system will be imposed at a later stage of the calculation. The canonical partition function $`𝒵`$ is given by

$$𝒵=\frac{1}{N_+!N_{}!}\underset{i=1}{\overset{N_+}{}}\frac{\mathrm{d}𝐫_i^{(+)}}{\lambda _{T,+}^3}\underset{j=1}{\overset{N_{}}{}}\frac{\mathrm{d}𝐫_j^{()}}{\lambda _{T,}^3}\mathrm{exp}\left(\frac{}{k_BT}\right)$$
(1)
where $`\lambda _{T,+}`$ and $`\lambda _{T,}`$ are the thermal wavelengths, and $`𝐫_i^{(+)}`$ and $`𝐫_i^{()}`$ are the positions of positively and negatively charged particles. The Hamiltonian $``$ is given by

$$\begin{array}{c}\frac{}{k_BT}=E_{\mathrm{self}}+\frac{1}{2}\underset{\alpha ,\beta =+,}{}d𝐫d𝐫^{}\widehat{\rho }_\alpha (𝐫)\omega _{\alpha \beta }(𝐫𝐫^{})\widehat{\rho }_\beta (𝐫^{})\hfill \\ \hfill +\frac{1}{2}d𝐫d𝐫^{}[\widehat{\rho }_+(𝐫)\widehat{\rho }_{}(𝐫)]v_c(𝐫𝐫^{})[\widehat{\rho }_+(𝐫^{})\widehat{\rho }_{}(𝐫^{})]\end{array}$$
(2)
where the charge-density operators of the ions are defined as

$$\begin{array}{cc}\hfill \widehat{\rho }_+(𝐫)& =q_+\underset{i=1}{\overset{N_+}{}}\delta (𝐫𝐫_i^{(+)})\hfill \\ \hfill \widehat{\rho }_{}(𝐫)& =q_{}\underset{i=1}{\overset{N_{}}{}}\delta (𝐫𝐫_i^{(+)}),\hfill \end{array}$$
(3)
and $`\delta (𝐫𝐫^{})`$ is the Dirac delta function. The indices $`\alpha `$ and $`\beta `$ in Eq. (2) stand for $`+`$ and $``$, and the sum over $`\alpha `$ and $`\beta `$ in Eq. (2) runs over all possible permutations (viz. $`++`$, $``$, $`+`$ and $`+`$), i.e., we consider a different short-ranged potential $`\omega _{\alpha \beta }`$ for each combination. The Coulomb potential is given by $`v_c(𝐫)=\mathrm{}_B/r`$, where $`\mathrm{}_Be^2/(4\pi \epsilon k_BT)`$ is the Bjerrum length, defined as the distance at which the electrostatic energy between two elementary charges equals the thermal energy $`k_BT`$. Finally, $`E_{\mathrm{self}}`$ is the self-energy of the system

$$E_{\mathrm{self}}=\frac{N_+q_+^2}{2}[\omega _{++}(0)+v_c(0)]+\frac{N_{}q_{}^2}{2}[\omega _{}(0)+v_c(0)]$$
(4)
which cancels the diagonal terms in Eq. (2).

We proceed by applying the Hubbard-Stratonovich transformation, the essence of which is given by

$$\begin{array}{c}\mathrm{e}^{\frac{1}{2}{\scriptscriptstyle d𝐫d𝐫^{}\widehat{\rho }(𝐫)v(𝐫𝐫^{})\widehat{\rho }(𝐫^{})}}=\hfill \\ \hfill \frac{𝒟\varphi \mathrm{e}^{\frac{1}{2}{\scriptscriptstyle d𝐫d𝐫^{}\varphi (𝐫)v^1(𝐫𝐫^{})\varphi (𝐫^{})}ı{\scriptscriptstyle d𝐫\varphi (𝐫)\widehat{\rho }(𝐫)}}}{𝒟\varphi \mathrm{e}^{\frac{1}{2}{\scriptscriptstyle d𝐫d𝐫^{}\varphi (𝐫)v^1(𝐫𝐫^{})\varphi (𝐫^{})}}}\end{array}$$
(5)
where $`v(𝐫)`$ is some general potential and $`𝒟\varphi `$ denotes a path integral over the fluctuating field $`\varphi `$. While this transformation can be used without problems when $`v(𝐫)`$ is the Coulomb potential, for a short-ranged potential this can be more problematic: for instance, a hard-core potential does not even have a well-defined inverse function. We will anyway take this formal step for the short-ranged potential, and, as we will see later, the way we handle the resulting expressions leads to finite (and consistent) results, viz., the virial coefficientsnote1 .

Applying Eq. (5) to Eq. (1) we obtain the partition function in field-theoretic form

$$𝒵=\frac{𝒟\psi _+𝒟\psi _{}}{𝒵_\psi }\frac{𝒟\varphi }{𝒵_\varphi }\mathrm{e}^{\overline{}_0}W_+W_{}$$
(6)
with the action

$$\begin{array}{c}\overline{}_0=\frac{1}{2}\underset{\alpha ,\beta =+,}{}d𝐫d𝐫^{}\psi _\alpha (𝐫)\omega _{\alpha \beta }^1(𝐫𝐫^{})\psi _\beta (𝐫^{})\hfill \\ \hfill +\frac{1}{2}d𝐫d𝐫^{}\varphi (𝐫)v_c^1(𝐫𝐫^{})\varphi (𝐫^{}).\end{array}$$
(7)
The inverse potentials are formally defined as the solution of the equation

$$\underset{\beta =+,}{}d𝐫^{}\omega _{\alpha \beta }(𝐫𝐫^{})\omega _{\beta \gamma }^1(𝐫^{}𝐫^{\prime \prime })=\delta _{\alpha \gamma }\delta (𝐫𝐫^{\prime \prime })$$
(8)
($`\delta _{\alpha \gamma }`$ is the Kronecker delta) and

$$d𝐫^{}v_c(𝐫𝐫^{})v_c^1(𝐫^{}𝐫^{\prime \prime })=\delta (𝐫𝐫^{\prime \prime }).$$
(9)
For the Coulomb potential, $`v_c^1(𝐫)=^2\delta (𝐫)/4\pi \mathrm{}_B`$. We also define

$$W_\alpha =\frac{1}{N_\alpha !}\left[\mathrm{e}^{q_\alpha ^2\left[\omega _{\alpha \alpha }(0)+v_c(0)\right]/2}\frac{\mathrm{d}𝐫}{\lambda _{T,\alpha }^3}\mathrm{e}^{ıq_\alpha \left[\psi _\alpha (𝐫)+\alpha \varphi (𝐫)\right]}\right]^{N_\alpha }$$
(10)
and the normalization factors

$$𝒵_\psi =𝒟\psi _+𝒟\psi _{}\mathrm{e}^{\frac{1}{2}\underset{\alpha \beta }{}{\scriptscriptstyle d𝐫d𝐫^{}\psi _\alpha (𝐫)\omega _{\alpha \beta }^1(𝐫𝐫^{})\psi _\beta (𝐫^{})}}$$
(11)
and

$$𝒵_\varphi =𝒟\varphi \mathrm{e}^{\frac{1}{2}{\scriptscriptstyle d𝐫d𝐫^{}\varphi (𝐫)v_c^1(𝐫𝐫^{})\varphi (𝐫^{})}}.$$
(12)
In order to make the calculations simpler we use the grand-canonical ensemble. This is achieved through the transformation

$$𝒬=\underset{N_+,N_{}=0}{\overset{\mathrm{}}{}}\mathrm{\Lambda }_+^{N_+}\mathrm{\Lambda }_{}^N_{}𝒵,$$
(13)
where $`\mathrm{\Lambda }_+`$ and $`\mathrm{\Lambda }_{}`$ are, respectively, the fugacities (exponential of the chemical potential) of the positively and negatively charged particles. We perform the sum over $`N_+`$ and $`N_{}`$ without constraints, i.e., without imposing the electroneutrality condition $`q_+N_+=q_{}N_{}`$. Imposing this condition before going to the grand-canonical ensemble makes the calculations much more difficult. Later, electroneutrality will be imposed order-by-order in the low-density expansion in a consistent way, any infinities arising from the non-neutrality of the system will then be automatically canceled.

In its full form, the grand-canonical partition function $`𝒬`$ reads

$$\begin{array}{c}𝒬=\frac{𝒟\psi _+𝒟\psi _{}}{𝒵_\psi }\frac{D\varphi }{𝒵_\varphi }\mathrm{exp}(\frac{\mathrm{\Lambda }_+}{\lambda _{T,+}^3}\mathrm{d}𝐫h_+(𝐫)\mathrm{e}^{iq_+\varphi (𝐫)}\hfill \\ \hfill +\frac{\mathrm{\Lambda }_{}}{\lambda _{T,}^3}\mathrm{d}𝐫h_{}(𝐫)\mathrm{e}^{iq_{}\varphi (𝐫)}\overline{}_0),\end{array}$$
(14)
where $`\overline{}_0`$ is given in Eq. (7). We defined the local non-linear operator

$$h_\alpha (𝐫)\mathrm{exp}\left(\frac{q_\alpha ^2}{2}\left[\omega _{\alpha \alpha }(0)+v_c(0)\right]ıq_\alpha \psi _\alpha (𝐫)\right)$$
(15)
(as before, $`\alpha `$ stands for both $`+`$ and $``$). At this point we rescale the fugacities such that $`\lambda _+=\mathrm{\Lambda }_+/\lambda _{T,+}^3`$ and $`\lambda _{}=\mathrm{\Lambda }_{}/\lambda _{T,}^3`$, i.e., the fugacities have from now on dimensions of inverse volume.

Introducing the Debye-Hückel propagator,

$$v_{\mathrm{DH}}^1(𝐫𝐫^{})=v_c^1(𝐫𝐫^{})+I_2\delta (𝐫𝐫^{})$$
(16)
with the ionic strength

$$I_2=q_+^2\lambda _++q_{}^2\lambda _{}$$
(17)
and after some algebraic manipulations, we finally obtain the grand-canonical free energy density. It is defined through $`g\mathrm{ln}\left(𝒬\right)/V`$, and reads

$$\begin{array}{c}g=\lambda _+\lambda _{}\frac{1}{2}I_2v_c(0)\frac{1}{V}\mathrm{ln}\left(\frac{𝒵_{\mathrm{DH}}}{𝒵_\varphi }\right)\hfill \\ \hfill \frac{1}{V}\mathrm{ln}\mathrm{e}^{\lambda _+{\scriptscriptstyle d𝐫Q_+(𝐫)}+\lambda _{}{\scriptscriptstyle d𝐫Q_{}(𝐫)}},\end{array}$$
(18)
where $`V`$ is the volume of the system and the brackets $`\mathrm{}`$ denote averages over the fluctuating fields $`\varphi `$ and $`\psi _\alpha `$ with the propagators $`\omega _{\alpha \beta }^1`$ and $`v_{\mathrm{DH}}^1`$. $`𝒵_{\mathrm{DH}}`$ is defined as

$$𝒵_{\mathrm{DH}}=𝒟\varphi \mathrm{e}^{\frac{1}{2}{\scriptscriptstyle d𝐫d𝐫^{}\varphi (𝐫)v_{\mathrm{DH}}^1(𝐫𝐫^{})\varphi (𝐫^{})}}$$
(19)
and the functions $`Q_+(𝐫)`$ and $`Q_{}(𝐫)`$ are defined by

$$Q_\alpha (𝐫)=h_\alpha (𝐫)\mathrm{e}^{ı\alpha q_\alpha \varphi (𝐫)}1+\frac{1}{2}q_\alpha ^2\varphi ^2(𝐫)\frac{1}{2}q_\alpha ^2v_c(0).$$
(20)
In the preceding steps we obtained the exact expression for the grand-canonical free energy density $`g`$, still without imposing electroneutrality. In order to obtain the low-density expansion of the free energy, we now (i) expand $`g`$ in powers of $`\lambda _+`$ and $`\lambda _{}`$ (up to a order $`5/2`$), (ii) calculate the concentrations of positive and negative particles and impose electroneutrality *consistently*, order-by-order, and (iii) make a Legendre transformation back to the canonical ensemble.

### 2.1 Expansion in powers of the fugacities

We start the expansion of $`g`$ by noting that the Fourier transform of the Coulomb potential is $`\stackrel{~}{v}_c(𝐤)=4\pi \mathrm{}_B/k^2`$. Using this, one is able to express $`v_c(0)`$ as

$$v_c(0)=\frac{\mathrm{d}𝐤}{(2\pi )^3}\frac{4\pi \mathrm{}_B}{k^2}.$$
(21)
It is easy to show that

$$\frac{1}{V}\mathrm{ln}\left(\frac{𝒵_{\mathrm{DH}}}{𝒵_\varphi }\right)=\frac{1}{2}\frac{\mathrm{d}𝐤}{(2\pi )^3}\mathrm{ln}\left(1+\frac{4\pi \mathrm{}_BI_2}{k^2}\right).$$
(22)
Using Eqs. (21) and (22), one finds

$$\frac{1}{2}I_2v_c(0)+\frac{1}{V}\mathrm{ln}\left(\frac{𝒵_{\mathrm{DH}}}{𝒵_\varphi }\right)=\frac{1}{12\pi }\left[4\pi \mathrm{}_BI_2\right]^{3/2}.$$
(23)
Introducing the dimensionless quantity

$$\mathrm{\Delta }v_0=\sqrt{4\pi \mathrm{}_B^3I_2}$$
(24)
and expanding the last term in the rhs of Eq. (18) in cumulants of $`Q_+(𝐫)`$ and $`Q_{}(𝐫)`$, one obtains the low-fugacity expansion

$$\begin{array}{c}g=\lambda _+\lambda _{}\frac{\mathrm{\Delta }v_0^3}{12\pi \mathrm{}_B^3}\lambda _+Z_1^+\lambda _{}Z_1^{}\frac{\lambda _+^2}{2}Z_2^{++}\hfill \\ \hfill \frac{\lambda _{}^2}{2}Z_2^{}\lambda _+\lambda _{}Z_2^++O(\lambda ^3),\end{array}$$
(25)
where the coefficients are given by

$$Z_1^+=\frac{1}{V}d𝐫Q_+(𝐫)$$
(26)
and an analogous formula for $`Z_1^{}`$, and

$$Z_2^{++}=\frac{1}{V}d𝐫d𝐫^{}\left\{Q_+(𝐫)Q_+(𝐫^{})Q_+(𝐫)Q_+(𝐫^{})\right\}$$
(27)
and similar formulas for $`Z_2^{}`$ and $`Z_2^+`$. The symbol $`O(\lambda ^3)`$ in Eq. (25) means that any other contribution to the expansion will be of order 3 or higher, i.e., with terms like $`\lambda _+^3`$, $`\lambda _+^2\lambda _{}`$, etc. Clearly, the expectation values in Eqs. (26) and (27) contain additional dependencies on the fugacity $`\lambda _\alpha `$ via the DH propagator, Eq.(16), but, and this stands at the very core of our method, all expectation values can itself be expanded with respect to $`\lambda _\alpha `$ and have finite values as $`\lambda _\alpha 0`$

In order to do a full expansion of Eq. (25), one needs first to calculate the coefficients $`Z_1^+`$, etc., in Eq. (25), for which the averages given in the appendix are needed, cf. Eq. (7882). Since $`v_c(0)v_{\mathrm{DH}}(0)=\mathrm{\Delta }v_0`$, we obtain

$$Z_1^\alpha =\mathrm{e}^{q_\alpha ^2\mathrm{\Delta }v_0/2}1\frac{1}{2}q_\alpha ^2\mathrm{\Delta }v_0$$
(28)
and

$$\begin{array}{c}Z_2^{\alpha \beta }=\mathrm{d}𝐫\{\mathrm{e}^{[q_\alpha ^2+q_\beta ^2]\mathrm{\Delta }v_0/2}[\mathrm{e}^{q_\alpha q_\beta [\omega _{\alpha \beta }(𝐫)+\alpha \beta v_{\mathrm{DH}}(𝐫)]}1]\hfill \\ \hfill +\frac{1}{2}q_\alpha ^2q_\beta ^2v_{\mathrm{DH}}^2(𝐫)[1\mathrm{e}^{q_\alpha ^2\mathrm{\Delta }v_0/2}\mathrm{e}^{q_\beta ^2\mathrm{\Delta }v_0/2}]\}.\end{array}$$
(29)
Note that

$$v_{\mathrm{DH}}(𝐫)=\frac{\mathrm{}_B}{r}\mathrm{e}^{\mathrm{\Delta }v_0r/\mathrm{}_B}$$
(30)
was defined (through its inverse function) in Eq. (16). We now introduce the hard-core through the short-range potentials

$$\omega _{\alpha \beta }(𝐫)=\{\begin{array}{cc}+\mathrm{}\hfill & \text{if }r<\left(d_\alpha +d_\beta \right)/2,\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}$$
(31)
where the indices $`\alpha `$ and $`\beta `$ stand again for $`+`$ and $``$; $`d_+`$ and $`d_{}`$ are respectively the (effective) ionic diameters of the positive and negative particles.

Since the expressions for $`Z_1^+`$, etc., do depend on $`\lambda _+`$ and $`\lambda _{}`$, and in order to have a consistent expansion of $`g`$ in the fugacities, one has to expand Eqs. (2829) in powers of $`\lambda _+`$ and $`\lambda _{}`$, up to the appropriate order, before inserting them into $`g`$, Eq. (25). By doing this consistently up to order $`5/2`$ in the fugacities, one obtains the rescaled grand-canonical free energy density

$$\begin{array}{c}\stackrel{~}{g}d_+^3g=\hfill \\ \hfill \stackrel{~}{\lambda }_+\stackrel{~}{\lambda }_{}m_1\stackrel{~}{\lambda }_+^2m_2\stackrel{~}{\lambda }_{}^2m_3\stackrel{~}{\lambda }_+\stackrel{~}{\lambda }_{}[n_1\stackrel{~}{\lambda }_+^2+\\ \hfill n_2\stackrel{~}{\lambda }_{}^2+n_3\stackrel{~}{\lambda }_+\stackrel{~}{\lambda }_{}]\mathrm{\Delta }v_0[p_1\stackrel{~}{\lambda }_+^2+p_2\stackrel{~}{\lambda }_{}^2+p_3\stackrel{~}{\lambda }_+\stackrel{~}{\lambda }_{}]\mathrm{ln}\mathrm{\Delta }v_0\\ \hfill [r_1\stackrel{~}{\lambda }_+^2+r_2\stackrel{~}{\lambda }_{}^2+r_3\stackrel{~}{\lambda }_+\stackrel{~}{\lambda }_{}]\mathrm{\Delta }v_0\mathrm{ln}\mathrm{\Delta }v_0[s_1\stackrel{~}{\lambda }_+\\ \hfill +s_2\stackrel{~}{\lambda }_{}]\mathrm{\Delta }v_0^2[t_0+t_1\stackrel{~}{\lambda }_++t_2\stackrel{~}{\lambda }_{}]\mathrm{\Delta }v_0^3+\mathrm{\Omega }_0[q_+\stackrel{~}{\lambda }_+q_{}\stackrel{~}{\lambda }_{}]\\ \hfill \mathrm{\Delta }v_0\{\mathrm{\Omega }_1[q_+\stackrel{~}{\lambda }_+q_{}\stackrel{~}{\lambda }_{}]^2\mathrm{\Omega }_0[q_+^4\stackrel{~}{\lambda }_+^2+q_{}^4\stackrel{~}{\lambda }_{}^2\\ \hfill 2q_+q_{}\left[\frac{q_+^2+q_{}^2}{2}\right]\stackrel{~}{\lambda }_+\stackrel{~}{\lambda }_{}]\}+O(\stackrel{~}{\lambda }^3),\end{array}$$
(32)
where we used dimensionless fugacities

$$\stackrel{~}{\lambda }_+d_+^3\lambda _+$$
(33)
and

$$\stackrel{~}{\lambda }_{}d_+^3\lambda _{}.$$
(34)
In this expansion, we utilized that $`\mathrm{\Delta }v_0`$ scales like $`\stackrel{~}{\lambda }^{1/2}`$. The coefficients $`m_1`$, etc. are given explicitly in the appendix, Eqs. (8392).

The coefficients $`\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_1`$ are given by the divergent integrals

$$\begin{array}{cc}\hfill \mathrm{\Omega }_0& =2\pi \mathrm{}_B\underset{0}{\overset{\mathrm{}}{}}drr\hfill \\ \hfill \mathrm{\Omega }_1& =2\pi \underset{0}{\overset{\mathrm{}}{}}drr^2.\hfill \end{array}$$
(35)
These terms are present in Eq. (32) because global charge neutrality has not been yet demanded. By imposing this condition, these divergent terms cancel exactly, as is shown next.

### 2.2 Imposing electroneutrality

The electroneutrality condition ensures that the global charge of the system is zero, i.e., $`q_+N_+=q_{}N_{}`$. In the grand-canonical ensemble, $`N_+`$ and $`N_{}`$ are no longer fixed numbers but average values. This means that the electroneutrality condition in the grand-canonical ensemble is given by

$$q_+N_+=q_{}N_{}.$$
(36)
Defining the rescaled ion density

$$\stackrel{~}{c}_+=d_+^3N_+/V,$$
(37)
it is easy to show that

$$\stackrel{~}{c}_+=\stackrel{~}{\lambda }_+\frac{\stackrel{~}{g}}{\stackrel{~}{\lambda }_+}$$
(38)
with an analogous formula for $`\stackrel{~}{c}_{}`$; $`\stackrel{~}{c}_+`$ is the volume fractions of the positive ionsnote\_volfrac . As one imposes Eq. (36), the fugacities will depend on each other in a non-trivial way such that the system is, on average, neutral. If the system were totally symmetric (i.e., $`q_+=q_{}`$ and $`d_+=d_{}`$), this dependence would be given by the relation $`\stackrel{~}{\lambda }_+=\stackrel{~}{\lambda }_{}`$roland-orland:tcp . However, this is not the case here: one needs to find the relation between the fugacities order-by-order.

First, assume that $`\stackrel{~}{\lambda }_{}`$ can be expanded in terms of $`\stackrel{~}{\lambda }_+`$ such that

$$\begin{array}{c}\stackrel{~}{\lambda }_{}=a_0\stackrel{~}{\lambda }_++a_1\stackrel{~}{\lambda }_+^{3/2}+a_2\stackrel{~}{\lambda }_+^2+a_3\stackrel{~}{\lambda }_+^2\mathrm{ln}\stackrel{~}{\lambda }_+\hfill \\ \hfill +a_4\stackrel{~}{\lambda }_+^{5/2}+a_5\stackrel{~}{\lambda }_+^{5/2}\mathrm{ln}\stackrel{~}{\lambda }_++O(\stackrel{~}{\lambda }_+^3).\end{array}$$
(39)
This is inspired by the expanded form of the grand-canonical free energy $`\stackrel{~}{g}`$, Eq.(32).

After calculating $`\stackrel{~}{c}_+`$ and $`\stackrel{~}{c}_{}`$ from $`\stackrel{~}{g}`$, Eq. (32), according to Eq.(38), and insertion into the electroneutrality condition Eq. (36), we substitute the fugacity $`\stackrel{~}{\lambda }_{}`$ by its expanded form Eq. (39). This leads to the expanded form (up to order $`\stackrel{~}{\lambda }_+^{5/2}\mathrm{ln}\stackrel{~}{\lambda }_+`$) of the electroneutrality condition. Solving it consistently, order-by-order, yields the values of the coefficients in Eq.(39), which ensure electroneutrality order by order. For instance, at linear order in $`\stackrel{~}{\lambda }_+`$, the expanded form of Eq. (36) reads

$$q_+a_0q_{}=0,$$
(40)
which naturally gives $`a_0=q_+/q_{}`$. With the knowledge of $`a_0`$, one can solve the next-order term (in this case $`\stackrel{~}{\lambda }_+^{3/2}`$) and obtain the value of $`a_1`$, and so on. The resulting coefficients $`a_0`$ up to $`a_5`$ are given in the appendix—cf. Eqs. (9398).

As this order-by-order neutrality condition is imposed, one notices that the terms $`\mathrm{\Omega }_0`$ and $`\mathrm{\Omega }_1`$ in Eq. (32) are exactly canceled in a natural way, without any further assumptions. The resulting expression for $`\stackrel{~}{g}`$, now expanded only in one of the fugacities (in this case $`\stackrel{~}{\lambda }_+`$), is then a well behaved expansion (we omit this expression here since it is quite lengthy). With this, one can finally obtain the canonical free energy (as a density expansion) through a Legendre transform.

### 2.3 The canonical free energy

In order to obtain the free energy in the canonical ensemble, we use the back-Legendre-transformation

$$\stackrel{~}{f}=\stackrel{~}{g}+\stackrel{~}{c}_+\mathrm{ln}\left(\stackrel{~}{\lambda }_+\right)+\stackrel{~}{c}_{}\mathrm{ln}\left(\stackrel{~}{\lambda }_{}\right),$$
(41)
where the dimensionless canonical free energy density is

$$\stackrel{~}{f}d_+^3F/Vk_BT$$
(42)
(F is the canonical free energy). Note that at this point $`\stackrel{~}{\lambda }_{}`$ is a function of $`\stackrel{~}{\lambda }_+`$, according to Eq. (39).

The first step to obtain $`\stackrel{~}{f}`$ is to invert the expression given in Eq. (38) such that $`\stackrel{~}{\lambda }_+`$ is obtained as an expansion in $`\stackrel{~}{c}_+`$, neglecting any terms of order $`\stackrel{~}{c}_+^3`$ or higher. With this, we obtain $`\stackrel{~}{\lambda }_+=\stackrel{~}{\lambda }_+(\stackrel{~}{c}_+)`$: plugging this into Eq.(41), we finally obtain $`\stackrel{~}{f}`$, which reads

$$\begin{array}{c}\stackrel{~}{f}=\stackrel{~}{f}_{\mathrm{id}}+B_{\mathrm{DH}}\stackrel{~}{c}_+^{3/2}+B_2\stackrel{~}{c}_+^2+B_{2\mathrm{l}\mathrm{o}\mathrm{g}}\stackrel{~}{c}_+^2\mathrm{ln}\stackrel{~}{c}_+\hfill \\ \hfill +B_{5/2}\stackrel{~}{c}_+^{5/2}+B_{5/2\mathrm{l}\mathrm{o}\mathrm{g}}\stackrel{~}{c}_+^{5/2}\mathrm{ln}\stackrel{~}{c}_++O\left(\stackrel{~}{c}_+^3\right).\end{array}$$
(43)
Defining the valence ratio parameter

$$\eta q_{}/q_+,$$
(44)
the diameter ratio parameter

$$\xi d_{}/d_+,$$
(45)
and the coupling parameter

$$ϵ_+q_+^2\mathrm{}_B/d_+$$
(46)
(which is the ratio between the Coulomb energy at contact between two positive ions and the thermal energy $`k_BT`$), the coefficients in Eq. (43) can be explicitly written as

$$\stackrel{~}{f}_{\mathrm{id}}=\stackrel{~}{c}_+\mathrm{ln}\stackrel{~}{c}_++\frac{\stackrel{~}{c}_+}{\eta }\mathrm{ln}\left(\frac{\stackrel{~}{c}_+}{\eta }\right)\left[1+\frac{1}{\eta }\right]\stackrel{~}{c}_+$$
(47)
which is the ideal contribution to the free energy,

$$B_{\mathrm{DH}}=\frac{2}{3}\sqrt{\pi ϵ_+^3[1+\eta ]^3}$$
(48)
which is the coefficient of the Debye-Hückel limiting law term (order $`3/2`$ in $`\stackrel{~}{c}_+`$). It is useful to define the function

$$H(x)=\frac{11}{6}2\gamma +\frac{1}{x^3}\mathrm{e}^x\left[2x+x^2\right]\mathrm{\Gamma }(0,x)\mathrm{ln}x,$$
(49)
where $`\gamma `$ is the Euler’s constant and $`\mathrm{\Gamma }(a,b)`$ is the incomplete Gamma-functionabramowitz:book . Using $`H(x)`$, the higher order coefficients can be explicitly written as

$$\begin{array}{c}B_2=\frac{\pi }{3}ϵ_+^3\{H(ϵ_+)\mathrm{ln}ϵ_++2\eta ^2[H(\frac{2\eta ϵ_+}{1+\xi })\hfill \\ \hfill +\mathrm{ln}\left(\frac{2\eta ϵ_+}{1+\xi }\right)]\eta ^4[H\left(\frac{\eta ^2ϵ_+}{\xi }\right)+\mathrm{ln}\left(\frac{\eta ^2ϵ_+}{\xi }\right)]\\ \hfill 2\eta ^2[1\eta ^2]\mathrm{ln}\eta +\frac{1}{2}[1\eta ^2]^2\mathrm{ln}\left(36\pi ϵ_+^3[1+\eta ]\right)\}\end{array}$$
(50)
$$B_{2\mathrm{l}\mathrm{o}\mathrm{g}}=\frac{\pi }{6}ϵ_+^3\left[1\eta ^2\right]^2,$$
(51)
$$\begin{array}{c}B_{5/2}=\frac{2}{3}\left[\pi ϵ_+^3\right]^{3/2}[1+\eta ]^{1/2}\{\frac{5}{8}+H\left(ϵ_+\right)+\mathrm{ln}\left(ϵ_+\right)\hfill \\ \hfill +\eta ^6[\frac{5}{8}+H\left(\frac{\eta ^2ϵ_+}{\xi }\right)+\mathrm{ln}\left(\frac{\eta ^2ϵ_+}{\xi }\right)]+2\eta ^3[\frac{5}{8}+H(\frac{2\eta ϵ_+}{1+\xi })\\ \hfill +\mathrm{ln}\left(\frac{2\eta ϵ_+}{1+\xi }\right)]+\frac{1}{8}[1+\eta \left]^2\right[512\eta +17\eta ^212\eta ^3+5\eta ^4]\\ \hfill 2\eta ^3[1+\eta ^3]\mathrm{ln}\eta \frac{1}{2}[1+\eta ^3]^2\mathrm{ln}\left(64\pi ϵ_+^3[1+\eta ]\right)\}\end{array}$$
(52)
and

$$B_{5/2\mathrm{l}\mathrm{o}\mathrm{g}}=\frac{1}{3}\left[\pi ϵ_+^3\right]^{3/2}\left[1+\eta \right]^{1/2}\left[1+\eta ^3\right]^2.$$
(53)
The free energy Eq. (43) is the exact low density expansion of the asymmetric TCPHC and constitutes the main result of this paper. The only parameter that is demanded to be small is the ion density $`\stackrel{~}{c}_+`$, that means, this result is non-perturbative in the coupling $`ϵ_+`$, charge ratio $`\eta `$ and size ratio $`\xi `$.

We chose the positive ions as the “reference species” (i.e., the expansion is done with respect to $`\stackrel{~}{c}_+`$) without any loss of generality, since the relation between $`\stackrel{~}{c}_+`$ and $`\stackrel{~}{c}_{}`$ is fixed through the electroneutrality condition. As consistency checks, we notice that our expression for $`\stackrel{~}{f}`$ is symmetric, as expected, with respect to the simultaneous exchange of $`d_+`$ with $`d_{}`$ and $`q_+`$ with $`q_{}`$. Also, in the limit $`d_+=d_{}`$ and $`q_+=q_{}`$, we recover the same expression as previously calculated in Ref. roland-orland:tcp  for totally symmetric systems. Finally, as one turns off the charges in the system (or equivalently, as one takes the limit $`ϵ_+0`$), the pure hard-core fluid is recovered, i.e., $`\stackrel{~}{f}`$ becomes the usual virial expansion with $`B_{3/2}`$ and $`B_{5/2}`$ equal to zero and $`B_2`$ given by the second virial coefficient of a two-component hard-core gas, $`B_2^{\mathrm{HC}}=B_2(ϵ_+0)`$, which in our units reads

$$B_2^{\mathrm{HC}}=\frac{\pi }{3}\left[2+\frac{2\xi ^3}{\eta ^2}+\frac{[1+\xi ]^3}{2\eta }\right].$$
(54)
This limit can be also understood as the high-temperature regime: as the thermal energy largely exceeds the Coulomb energy at contact, the hard core interaction becomes the only relevant interaction between the particles.

## 3 Results

### 3.1 The virial coefficients

The behavior of the coefficients $`B_2`$ and $`B_{5/2}`$ as functions of the coupling parameter $`ϵ_+`$ are depicted in Figs. 1 and 2, the behavior of $`B_{2\mathrm{l}\mathrm{o}\mathrm{g}}`$ and $`B_{5/2\mathrm{l}\mathrm{o}\mathrm{g}}`$ is rather trivial and not shown graphically. In Fig. 1 the ionic diameters of positive and negative ions are kept equal, $`d_+=d_{}`$, and the ratio between the charge valences, $`\eta `$, is varied, while in Fig. 2 the charge valences are equal and the ratio between the ionic diameters is varied. These figures highlight the fact that both coefficients diverge as $`ϵ_+`$ goes to infinity. In this limit,

$$B_2\frac{\pi }{4}\frac{[1+\xi ]^4}{\eta ^2}\frac{1}{ϵ_+}\mathrm{exp}\left(\frac{2\eta ϵ_+}{1+\xi }\right)$$
(55)
and

$$B_{5/2}\frac{\pi ^{3/2}}{2}\frac{[1+\xi ]^4\sqrt{1+\eta }}{\eta }\sqrt{ϵ_+}\mathrm{exp}\left(\frac{2\eta ϵ_+}{1+\xi }\right).$$
(56)
Note that $`B_{5/2}`$ diverges faster (and with opposite sign) than $`B_2`$, reflecting that this low-density expansion is badly converging, as we will discuss further below. The exponentially divergent behavior of $`B_2`$ and $`B_{5/2}`$ when $`ϵ_+\mathrm{}`$ is due to the increasing importance of the interaction between oppositely charged particles (ionic pairing) as the coupling parameter increasesbjerrum:26 ; fisher-levin:93 ; yeh-zhou-stell:96 , corresponding for instance to lower temperatures. This is confirmed by noting that the argument in the exponential occurring in both asymptotic forms Eqs. (55) and (56), viz. $`2\eta ϵ_+/[1+\xi ]`$, can be re-expressed as $`2q_+q_{}\mathrm{}_B/[d_++d_{}]`$, which is the coupling between positive and negative ions (in this case, the ratio between the Coulomb contact energy between oppositely charged ions and the thermal energy $`k_BT`$).

To estimate roughly the ionic density $`\stackrel{~}{c}_+`$ up to which the expansion in Eq. (43) is expected to be valid, we use the following simple criterion: the terms proportional to $`B_{\mathrm{DH}}`$ and $`B_{5/2}`$ give the same contribution to the free energy when $`|B_{\mathrm{DH}}|\stackrel{~}{c}_+^{3/2}=|B_{5/2}|\stackrel{~}{c}_+^{5/2}`$ or, in other words, at a critical density $`\stackrel{~}{c}_+=|B_{\mathrm{DH}}/B_{5/2}|`$. Note that we set up this criterion separately for the integer and fractional coefficients, since the scaling behavior for the two different classes is very different and no meaningful results can be obtained by mixing them. The analogous criterion for the integer terms leads to a critical density $`\stackrel{~}{c}_+=(1+1/\eta )/|B_2|`$ where for the ideal contribution in Eq. 47 we replace the logarithmic term by a linear one (which is the ideal contribution to the osmotic pressure). We tacitly assume that the higher-order terms (which we have not calculated) show the same relative behavior, an assumption which seems plausible to us but which we cannot check. In Fig. 3 we show the estimate for the critical densities obtained from a) the ratio of the integer terms and from b) the ratio of the fractional terms. In a) the critical density settles at a finite value for vanishing coupling parameter $`\epsilon _+0`$, and decreases to zero for increasing coupling parameter (the divergence in the curves at finite value of $`\epsilon _+`$ is not significant since it is caused by a change in sign of $`B_2`$.). This means that for small coupling constants the integer coefficients are expected to show regular convergence behavior for not too large concentrations. For large values of $`\stackrel{~}{c}_+`$ on the other hand, the range of densities which can be described by the expansion quickly goes to zero. In b) the displayed behavior is more complex. The expansion of the coefficient $`B_{5/2}`$ for small values of the coupling parameter $`\epsilon _+`$ leads to

$`B_{5/2}`$ $`=`$ $`\pi ^{3/2}\sqrt{1+\eta }(\xi 1)^2(\xi +1)\epsilon _+^{3/2}`$
$`\pi ^{3/2}\sqrt{1+\eta }(2+2\eta ^2\xi ^2+\eta (1+\xi )^2)\epsilon _+^{5/2}`$
plus terms with scale as higher powers in $`\epsilon _+`$. The important point is that except for the singular case $`\xi =1`$, that is for particles of identical radius, the leading term of $`B_{5/2}`$ scales as $`\epsilon _+^{3/2}`$ and thus identical to the lower leading term $`B_{\mathrm{DH}}`$. In the case $`\xi =1`$ the leading term scales as $`\epsilon _+^{5/2}`$ and thus shows a different scaling behavior than $`B_{\mathrm{DH}}`$. In Fig. 3b) this difference is cleary seen: For $`\xi =1`$ the curves diverge as $`\epsilon _+0`$ and the convergence of the fractional density terms in the series is expected to be guaranteed. For $`\xi =0.1`$ the behavior is similar to the results in a), showing a saturation at a finite value as $`\epsilon _+0`$ (and an unsignificant divergence at finite $`\epsilon _+`$). All curves go to zero as the coupling parameter grows. In conclusion, the reliability of the low-density expansion becomes worse as the coupling parameter increases, but one can always find a window of small densities within which the expansion should work.

Finally, we mention that in a different calculation scheme, based on an integral-equation-procedure within the so-called mean-spherical-model (MSM) approximation, similar results for the free energy and other thermodynamic functions have been obtainedwaisman-lebowitz:72 ; Palmer ; Blum . Those results also reproduce the limiting laws, namely the hard-core behavior as the coupling parameter goes to zero, and the leading Debye-Hückel correction at low density, and give a quite satisfactory description even for much larger values of the density. It is important to note, however, that the next-leading (beyond Debye-Hückel) terms in the low-density expansion are not correctly reproduced by the MSM approximation.

### 3.2 The colloidal limit

In colloidal suspensions, flocculation or coagulation (driven by attractive van der Waals interaction between colloidal particles) can be prevented by the presence of repulsive electrostatic forces. These suspensions are typically quite dilute, with volume fractions of colloidal particles usually not higher than a few percent. The macro-particles normally have dimensionsottewill:review ; gonnet-al:94 ; kurth-al:2000  ranging from 10 to 1000 nm and charges of several thousands $`e`$ (elementary charge), with much smaller counterions that have a charge of a few $`e`$. In such systems the charge and size asymmetry between ions and counterions is immense; the TCPHC with unconstrained charge and size asymmetry is a suitable model for dilute colloidal suspensions, and the free energy Eq. (43) can be used to study the thermodynamic behavior of such systems.

With this in mind, let us take the following limit: assume the parameters describing the positive ions (representing the colloids) $`d_+`$, $`q_+`$ and $`\stackrel{~}{c}_+`$ fixed, and make both $`\eta q_{}/q_+`$ and $`\xi d_{}/d_+`$ vanishingly small (this has to be done with some care, since the limit $`\eta =0`$ is not well-defined). One can then rewrite the free energy Eq. (43) up to second order in $`\stackrel{~}{c}_+`$ as

$$\begin{array}{c}\stackrel{~}{f}^{\mathrm{cl}}=\stackrel{~}{f}_{\mathrm{id}}++B_{\mathrm{DH}}^{\mathrm{cl}}\stackrel{~}{c}_+^{3/2}+(B_2^{\mathrm{HC},\mathrm{cl}}+B_2^{\mathrm{cl}})\stackrel{~}{c}_+^2\hfill \\ \hfill +B_{2\mathrm{l}\mathrm{o}\mathrm{g}}^{\mathrm{cl}}\stackrel{~}{c}_+^2\mathrm{log}\stackrel{~}{c}_++O(\stackrel{~}{c}_+^{5/2})\end{array}$$
(58)
where $`\stackrel{~}{f}_{\mathrm{id}}`$ is the ideal term Eq. (47) and

$$B_2^{\mathrm{HC},\mathrm{cl}}=\frac{\pi }{3}\left[\frac{2\xi ^3}{\eta ^2}+\frac{[1+\xi ]^3}{2\eta }\right]$$
(59)
is the hard core contribution due to the counterion–counterion and macroion–counterion interactions only, without the macroion–macroion contributions, which explains the difference to the full hard-core virial coefficient in Eq.(54). It is necessary to treat this term separately, as it does not have a well-defined behavior in the double limit $`\eta 0`$ and $`\xi 0`$. The other coefficients in Eq.(58) follow from Eqs.(48), (50), (51) by performing the limits $`B_{\mathrm{DH}}^{\mathrm{cl}}=B_{\mathrm{DH}}(\eta 0,\xi 0)`$, $`B_2^{\mathrm{cl}}=B_2(\eta 0,\xi 0)B_2^{\mathrm{HC},\mathrm{cl}}`$, $`B_{2\mathrm{l}\mathrm{o}\mathrm{g}}^{\mathrm{cl}}=B_{2\mathrm{l}\mathrm{o}\mathrm{g}}(\eta 0,\xi 0)`$, and are given by

$$B_{\mathrm{DH}}^{\mathrm{cl}}=\frac{2}{3}\sqrt{\pi ϵ_+^3},$$
(60)
$$B_2^{\mathrm{cl}}=\frac{\pi }{3}ϵ_+^3\left\{H\left(ϵ_+\right)\frac{1}{2}\mathrm{ln}\left(36\pi ϵ_+\right)\right\}+\frac{\pi }{2}ϵ_+,$$
(61)
$$B_{2\mathrm{l}\mathrm{o}\mathrm{g}}^{\mathrm{cl}}=\frac{\pi }{6}ϵ_+^3.$$
(62)
It is to be noted that the contributions in $`\stackrel{~}{f}_{\mathrm{id}}`$ and $`B_2^{\mathrm{HC},\mathrm{cl}}`$ contain terms that scale with $`1/\eta `$. For dilute colloidal systems therefore the ideal contribution of the counterions dominates the free energy and cannot be neglected. Effects due to the electrostatic interaction between the ions are corrections to the ideal behavior.

At this point we briefly introduce yet another model that is also widely used to describe charged systems. It is the one-component plasma (OCP), which in its simplest form consists of a collection of $`N`$ equally charged point-like particles immersed in a neutralizing background that ensures the global charge neutrality of the system (in the TCPHC electroneutrality is ensured by oppositely charged particles). The OCP, or its quantum mechanical counter-part (“jellium”) has been used in different contexts in physics, as for instance to describe degenerate stellar matter (the interior of white dwarfs or the outer layer of neutron stars) and the interior of massive planets like Jupiter. Another example comes from condensed matter physics, where jellium is often used as a reference state when calculating the electronic structure of solids. For reviews see Refs. ishimaru:82 ; baus-hansen:80 ; abe:59 .

When the particles have a hard core, the OCP is called one-component hard-core plasma (OCPHC): what we will see next is that the electrostatic contribution to the free energy in dilute colloidal suspensions can be almost described through the OCPHC. If we compare the coefficients Eqs. (6051) with the ones previously obtainedroland-orland:tcp  for the low density expansion of the OCPHC

$$\stackrel{~}{f}_{\mathrm{OCPHC}}=\stackrel{~}{c}\mathrm{ln}\stackrel{~}{c}\stackrel{~}{c}+d_{3/2}\stackrel{~}{c}^{3/2}+d_2\stackrel{~}{c}^2+d_{\mathrm{ln}2}\stackrel{~}{c}^2\mathrm{ln}\stackrel{~}{c}+O(\stackrel{~}{c}^{5/2}),$$
(63)
where $`\stackrel{~}{c}`$ is the volume fraction of particles in the OCPHC, we find (cf. Eqs. (51), (56) and (57) of Ref. roland-orland:tcp )

$$B_{\mathrm{DH}}^{\mathrm{cl}}=d_{3/2},$$
(64)
$$B_2^{\mathrm{cl}}=d_2+\frac{\pi }{2}ϵ_+$$
(65)
and

$$B_{2\mathrm{l}\mathrm{o}\mathrm{g}}^{\mathrm{cl}}=d_{\mathrm{ln}2},$$
(66)
where the notation used in Ref. roland-orland:tcp  is kept in Eq. (63) and on the rhs of Eqs. (6466). The comparison between $`B_2^{\mathrm{cl}}`$ and $`d_2`$ is shown in Fig. 4, where we rescaled both coefficients by their value at vanishing coupling, $`d_2(ϵ_+=0)=B_2^{\mathrm{cl}}(ϵ_+=0)=2\pi /3`$. This shows that, including the order $`\stackrel{~}{c}_+^2\mathrm{ln}\stackrel{~}{c}_+`$ (i.e., for very dilute colloidal suspension), the OCPHC is almost recovered, except an additional term in Eq.(65). This additional term can be quite important at intermediate values of $`ϵ_+`$ since it changes the sign of the virial coefficient, as seen in Fig. 3.

This is in fact what one would intuitively expect: each macroion has around it a very large number of small neutralizing counterions which act like a background. Our results show that the background formed by counterions, which is a deformable background since the counterion distribution is not uniform, acts almost like a rigid homogeneous background, as it is assumed in OCP calculations. However, there is a small difference between the TCPHC in this limit and the OCPHC: in the latter, the background penetrates the particles, while in the TCPHC it cannot (for a thorough discussion of this difference see Hansen ). In our calculation, this is reflected in the $`\stackrel{~}{c}_+^2`$ term, where $`\pi ϵ_+/2`$ is the positive extra cost in the free energy that the OCPHC has to pay (at this order) to expel the background from the hard-core particles. The first consequence of the non-penetrating background is that the effective charge of the colloids increases by exactly the amount of background charge that is expelled from the colloidal interior. A simple calculation shows that the increased effective charge $`q_+^{\mathrm{eff}}`$ of colloids for the non-penetrating background turns out to be

$$q_+^{\mathrm{eff}}=\frac{q_+}{1\pi \stackrel{~}{c}_+/6}.$$
(67)
In the low-density expansion of the OCPHC, Eq.(63), one would have to reexpand all coefficients with respect to the density, using Eq.(67). However, since the leading term which depends on the charge $`q_+`$ goes like $`\stackrel{~}{c}_+^{3/2}`$, the effect would come in at order $`\stackrel{~}{c}_+^{5/2}`$ and therefore is not responsible for the additional term in Eq.(65). The reason for the extra term in Eq.(65) has to do with an increase of the self-free-energy of a colloidal particle, which can be understood in the following way: Assume that the OCPHC is very dilute, such that each colloidal particle and its neutralizing background form a neutral entity (in the spirit of the cell model, see for instance Ref. alexander-al:84 ) that can be regarded independent from the other particles. The free energy difference *per colloidal particle* between a system without penetrating background and with penetrating background (Fig. 5) is then given by two contributions: one coming from the entropy lost by the background (formed by $`N_{}`$ counterions) since it cannot penetrate the macroions, given by

$$\frac{\mathrm{\Delta }F_{en}}{N_+k_BT}=\frac{N_{}}{N_+}\mathrm{ln}(1\pi \stackrel{~}{c}_+/6)=\frac{N_{}}{N_+}\frac{\pi \stackrel{~}{c}_+}{6}+O(\stackrel{~}{c}_+^2).$$
(68)
This is one of the contributions to the second virial coefficient due to the hard-core interaction, and already included in Eq. (59). This term does not depend on the coupling $`ϵ_+`$, and therefore has nothing to do with the extra term in Eq.(65). The second contribution to the self-energy is electrostatic in origin. Defining the charge distributions for the situations in Fig. 4 as $`\rho _a(𝐫)=q_+\delta (𝐫)\rho _0`$ when the background can enter the colloids and $`\rho _b(𝐫)=q_+\delta (𝐫)\rho _0[1\theta (d_+/2r)]`$ when the background cannot enter the colloids (where $`\rho _0=c_+q_+`$ is the background charge density and $`\theta (x)=1`$ if $`x>0`$ and $`0`$ otherwise), the electrostatic self-energy difference between the two cases reads

$$\begin{array}{c}\frac{\mathrm{\Delta }F_{el}}{N_+k_BT}=\frac{1}{2}\mathrm{d}𝐫\mathrm{d}𝐫^{}\{\rho _b(𝐫)v_c(𝐫𝐫^{})\rho _b(𝐫^{})\hfill \\ \hfill \rho _a(𝐫)v_c(𝐫𝐫^{})\rho _a(𝐫^{})\}\\ \hfill =\rho _0q_+4\pi \mathrm{}_B_0^{d_+/2}drr+O(c_+^2)\\ \hfill =\frac{\pi }{2}\stackrel{~}{c}_+ϵ_++O(c_+^2).\end{array}$$
(69)
It follows that the free energy difference per volume is given by

$$\mathrm{\Delta }\stackrel{~}{f}=\frac{\stackrel{~}{c}_+\mathrm{\Delta }F_{el}}{N_+k_BT}=\frac{\pi \stackrel{~}{c}_+^2ϵ_+}{2}.$$
(70)
Note that this is a positive contribution to the $`\stackrel{~}{c}^2`$ term with the same coefficient as the extra term in Eq. (65). Therefore, we conclude that this extra term is due to the electrostatic energy associated with expelling the background from the colloidal volume. This is in accord with more general results on the difference between OCPHC models with penetrating and excluded neutralising backgroundsHansen . In summary, whenever using the OCPHC to describe highly asymmetric charged systems such as colloids, one has to to take into account the exclusion of the background from the macroions. As we demonstrated, if this is taken into account, then the TCPHC maps exactly (at least up to order $`\stackrel{~}{c}_+^2`$) onto the OCPHC.

### 3.3 Ionic activity and diameters

In electrochemistry, it is usual to define the mean activity $`\lambda _\pm `$ of an electrolyte as

$$\lambda _\pm \left[\lambda _+^q_{}\lambda _{}^{q_+}\right]^{1/[q_++q_{}]},$$
(71)
where $`\lambda _+`$ and $`\lambda _{}`$ are respectively the fugacities of the positive and the negative ions. The mean activity coefficient $`\mathrm{f}_\pm `$ is the ratio between the mean activity of the electrolyte and that of an ideal gas, in general given by

$$\mathrm{f}_\pm =\mathrm{exp}\left(\frac{q_{}}{q_++q_{}}\frac{\stackrel{~}{f}_{\mathrm{ex}}}{\stackrel{~}{c}_+}\right),$$
(72)
for a two-component system where the positively charged particles are used as reference species in the same way as in Eq. (43). $`\stackrel{~}{f}_{\mathrm{ex}}`$, the excess free energy, is the difference between the free energy of the interacting systems and the free energy of an ideal gas, i.e., $`\stackrel{~}{f}_{\mathrm{ex}}=\stackrel{~}{f}\stackrel{~}{f}_{\mathrm{id}}`$.

There are different ways of measuring $`\mathrm{f}_\pm `$, as for instance, through the change of the freezing point of the solvent (usually water) with the addition of saltrobinson-stokes:book , by measuring the change of the potential difference between the electrodes of a concentration cell as salt is addedbutler-roy:91 ; robinson-stokes:book  (potentiometry), or by direct measurement of the solvent activity through vapor exchange between a solution with known activity and the samplerard-platford:91 ; robinson-stokes:book  (isopiestic). Although dating from the early nineteen hundreds, these are still the most common techniques used today, especially the potentiometry, which is regarded as the most precise technique of all. The values of $`\mathrm{f}_\pm `$ are tabled as function of the salt concentration for many different electrolytesconway:book ; parsons:book ; robinson-stokes:book .

From the free energy Eq. (43) and the definition Eq. (72) we can obtain the low density expansion for the mean activity coefficient of a $`q_+:q_{}`$ salt. To compare with experimental results, it is useful to note that $`\stackrel{~}{c}_+=6.022\times 10^4q_{}d_+^3\varrho `$, where $`d_+`$ is in Ångströms and $`\varrho `$ is the salt concentration in moles/liter. After the appropriate expansion, $`\mathrm{f}_\pm `$ reads

$$\begin{array}{c}\mathrm{f}_\pm =1+\nu _{\mathrm{DH}}\varrho ^{1/2}+\nu _1\varrho +\nu _{1\mathrm{l}\mathrm{o}\mathrm{g}}\varrho \mathrm{ln}\varrho +\nu _{3/2}\varrho ^{3/2}\hfill \\ \hfill +\nu _{3/2\mathrm{l}\mathrm{o}\mathrm{g}}\varrho ^{3/2}\mathrm{ln}\varrho +O(\varrho ^2)\end{array}$$
(73)
(the order $`3/2`$ in the mean activity coefficient is the one consistent with a free energy given up to $`5/2`$). The coefficients $`\nu _{\mathrm{DH}}`$, $`\nu _1`$, etc. are given in the appendix, Eqs. (99103). Experimentally determined activity coefficients are typically measured at constant pressure, while the theoretical results are obtained for constant volume. Note that in the limit of vanishing density, the distinction between the activities calculated in the Lewis-Randall or in the McMillan-Mayer descriptions (constant pressure versus constant volume ensemble) becomes negligiblepailthorpe-al:82  in comparison to the first corrections to ideal behavior. The experimental data at low density can be directly compared with the theory without the need to convert between the two ensembles.

At infinite dilution, the mean activity coefficient Eq. (73) goes to $`1`$, which is the prediction for an ideal gas. The first correction to the ideal behavior is the term $`\nu _{\mathrm{DH}}\varrho ^{1/2}`$, which is the prediction one obtains from the Debye-Hückel limiting law (DHLL), and is independent of the ionic diameters, see Eq. (99). This means that there is always a range of concentrations where different salts (but with the same $`q_+`$ and $`q_{}`$) will deviate from ideality, but have the same activity. At higher concentrations, other terms have to be taken into account, and the ionic sizes begin to play an important role. In fact, as one fixes $`\mathrm{}_B`$ ($`7.1`$ Å in water at 25 <sup>o</sup>C) and the valences $`q_+`$ and $`q_{}`$, the only free parameters in the rhs of Eq. (73) are the ionic diameters $`d_+`$ and $`d_{}`$, which can be used in the theoretical predictions to fit the experimental values. This leads to effective equilibrium values of the ionic diameters (when in solution), which we call the “thermodynamic diameter”, in contrast to the bare diameterisraelachvili:book  (obtained through crystallographic methods) and the hydrodynamic diameter (obtained from mobility measurementswennerstroem:book ). By construction, it is this diameter which should be used in equilibrium situations if one wants to describe an electrolyte solution as a TCPHC, as for example computer simulations of electrolytes.

#### 3.3.1 Fitting assuming one mean diameter

We now show the fitting procedure assuming that $`d_+=d_{}=d`$, where $`d`$ is the mean diameter. We also restrict here this procedure to $`1:1`$ salts, where more experimental data is available at reasonably low densities (below $`0.05`$ mole/liter). For asymmetric salts, this method demands precise data at the range below $`0.01`$ mole/liter, where Debye-Hückel is often assumed to account for all effects and few experimental points are available.

The assumption of equal ionic sizes has been often used in the past to fit activities to modified Debye-Hückel theoriesrobinson-stokes:book  to account for the ionic sizes (as previously mentioned, the Debye-Hückel limiting law is insensitive to it). Since we force the two diameters to be equal, we only need the expansion for $`f_\pm `$ up to linear order since this will determine uniquely the single unknown parameter. From the experimental data $`\mathrm{f}_\pm ^{exp}`$sheldlovsky:50 ; conway:book  we subtract the diameter-independent term and obtain the difference

$$\mathrm{\Delta }\mathrm{f}_\pm \mathrm{f}_\pm ^{exp}\nu _{\mathrm{DH}}\varrho ^{1/2}=\nu _1^{exp}\varrho +O(\varrho ^{(3/2)}).$$
(74)
The fit to this function leads to the coefficient of the linear term $`\nu _1^{exp}`$, and that can be used to determine the diameter by solving the equation $`\nu _1(d)=\nu _1^{exp}`$ (with $`d`$ as unknown). In Fig. 6 we show the experimentalsheldlovsky:50 ; conway:book  $`f_\pm `$ for HCl, NaCl and KCl and the theoretical results (up to linear order in $`\varrho `$) using, respectively, the diameters $`4.0`$ Å, $`3.7`$ Å and $`3.6`$ Å. These values come from the fitting described above applied to the experimental data in the range $`\varrho =0.01`$ to $`0.05`$ mole/liter. The diameters obtained are very close to the ones obtained in Ref. sheldlovsky:50  with a similar fitting, but using the DH theory with an approximate way to incorporate the ionic sizes.

The fact that $`d_{\mathrm{HCl}}>d_{\mathrm{NaCl}}>d_{\mathrm{KCl}}`$, which is the opposite to the sequence of bare diameters, is a consequence of the existence of hydration shells around the ions which tend to be larger for smaller bare ion sizeisraelachvili:book . Notice however that the values obtained here lie between the bare diameters and the hydrated values available in the literature. This is not surprising since the effective hard-core size simultaneously reflects both the presence of the hydration shells and the “softness” of the outer water layer as two oppositely charged ions approach each other.

#### 3.3.2 Fitting assuming two diameters

If we now assume that both $`d_+`$ and $`d_{}`$ are unconstrained, we have to use one more term in the expansion of the activity coefficient and solve the coupled equations

$$\begin{array}{cc}\hfill \nu _1(d_+,d_{})& =\nu _1^{exp}\hfill \\ \hfill \nu _{3/2}(d_+,d_{})& =\nu _{3/2}^{exp}\hfill \end{array}$$
(75)
where the coefficients $`\nu _1^{exp}`$ and $`\nu _{3/2}^{exp}`$ are obtained by fitting the function

$`\mathrm{\Delta }\mathrm{f}_\pm `$ $`=`$ $`\mathrm{f}_\pm ^{exp}\nu _{\mathrm{DH}}\varrho ^{1/2}\nu _{1\mathrm{l}\mathrm{o}\mathrm{g}}\varrho \mathrm{ln}\varrho `$ (76)
$`=`$ $`\nu _1^{exp}\varrho +\nu _{3/2}^{exp}\varrho ^{3/2}+O(\varrho ^{3/2}\mathrm{ln}\varrho ).`$
For the actual fitting we divide by the density and obtain

$$\frac{\mathrm{\Delta }\mathrm{f}_\pm }{\varrho }=\nu _1^{exp}+\nu _{3/2}^{exp}\varrho ^{1/2}.$$
(77)
Note that we subtracted the term proportional to $`\varrho \mathrm{ln}\varrho `$, since like the DH term it is independent of the ionic diameters. What is interesting about this fitting procedure is that the two ionic diameters are independent parameters, that is, the effective sizes obtained through this fitting do not depend on the size of a “reference” ion, contrary to what happens when using crystallographic methodsschwister-al:book . The method we show here can in principle lead to useful results, as long as the experimental data is accurate enough at very low densities, as we now demonstrate.

Fig. 7 shows the mapping between the parameter space of ionic diameters and the two coefficients $`\nu _1`$ and $`\nu _{3/2}`$ (for a $`1:1`$ salt in water at room temperature). The mapping is restricted to the range $`3`$ Å $`<d_+<27`$ Å and $`3`$ Å $`<d_{}<27`$ Å. This allows the “inverse mapping” between the ($`\nu _1`$;$`\nu _{3/2}`$)–space and the ($`d_{}`$;$`d_+`$)–space: any system with ionic diameters within the values above should have the coefficients $`\nu _1`$ and $`\nu _{3/2}`$ inside the region in the ($`\nu _1`$;$`\nu _{3/2}`$)–space, as shown in Fig. 7. One problem becomes obvious: the coefficient $`\nu _{3/2}`$ is larger than the coefficient of the linear term, $`\nu _1`$, meaning that experimental data at very low concentrations are needed. As an example, Fig. 8 shows $`\mathrm{\Delta }\mathrm{f}_\pm /\varrho `$ for NaClconway:book  as a function of $`\varrho ^{1/2}`$, which should be asymptotically linear in the limit $`\varrho 0`$ and therefore give the coefficients $`\nu _1`$ and $`\nu _{3/2}`$ by a simple linear fit according to Eq.(77). The lines shown are two possible asymptotic linear fits to the experimental points (line $`a`$ is given by $`\mathrm{\Delta }f_\pm /\varrho =2.599.94\varrho ^{1/2}`$, and line $`b`$ by $`2.476.00\varrho ^{1/2}`$). The inset to Fig. 8 gives the corresponding positions in the ($`\nu _1`$;$`\nu _{3/2}`$)–space (cf. Fig.7) for the two lines. Clearly, only one of the fits (line $`b`$) leads to reasonable values for the diameters (between $`3`$ and $`9`$ Å as seen from the position of point $`b`$ in the inset, $`d_+=3.8`$ Å and $`d_{}=5.4`$ Å as obtained by solving Eq. (75)), while the other fit (line $`a`$) leads to unreasonable values for the diameters (outside the range $`3`$ to $`27`$ Å). In other words, the asymptotic extrapolation of the experimental data is very sensitive to small errors, and in order to obtain the diameters with this method one needs more accurate data at very low densities (which to the best of our knowledge is not available in the literature). Although NaCl was used as example, the situation is identical for other salts.

## 4 Conclusions

Using field theoretic methods we obtained the exact low density (“virial”) expansion of the TCPHC up to order $`5/2`$ in density. In its general form, the model can be applied to both electrolyte solutions and dilute colloidal suspensions (when the van der Waals forces are unimportant); the free energy derived here provides a unified way for handling both systems in the limit of low concentration. As the calculations show, the generalization to short-ranged potentials other than hard core is possible.

The behavior of the coefficients $`B_2`$ and $`B_{5/2}`$ suggests that the series is badly convergent, meaning that the inclusion of higher order terms does not necessarily extend the validity of the free energy to larger densities. We saw that the divergent behavior of these coefficients is related to the ionic pairingbjerrum:26 ; fisher-levin:93 ; yeh-zhou-stell:96  which is favored as the coupling increases (this is not present e.g. in the OCPHCroland-orland:tcp ). This also means that such a low-density expansion is not very useful to study the phase behaviorthree-component  of ionic systems: in the situation of phase-separation between a very dilute phase and a dense phase, typically the average density in the dense phase already falls outside the range of validity of the low density expansion.

In applying our results to colloidal systems we concluded that at low density the counterion entropic contribution dominates over the electrostatic contribution due to the macroions; the latter contribution can be described, in this limit and up to this order, by an OCPHC corrected by effects due to the exclusion of counterions from the colloidal particles.

Finally, we used the theoretical results for the mean activity coefficient to fit experimental data and extract effective ionic sizes. In the simplest fitting, where we assumed the two ionic diameters to be equal, we obtained sizes that are reasonable and close to what one would expect from the results obtained by other methods. For the more interesting case where the two ionic sizes are taken as free parameters and determined independently, we showed that one would need more experimental data for the mean activity coefficient at very low densities (which, to the best of our knowledge, is not available in the literature) to obtain the correct values for the diameters. With the proper experimental data it would then be a simple matter to obtain the effective thermodynamic ionic sizes, which could serve as useful input for computer simulations of two-component-plasmas with hard-cores to model electrolyte solutions.

###### Acknowledgements.

We thank N. Brilliantov and H. Löwen for useful discussions. AGM acknowledges the support from FCT through grant PRAXIS XXI/BD/13347/97 and DFG.

## Appendix

### Averages needed for $`Z_1`$ and $`Z_2`$

The following expressions were used to obtain Eqs. (2829). In order to have more compact formulas we use the Greek letters $`\alpha `$ and $`\beta `$ instead of $`+`$ or $``$. For instance, $`\alpha q_\alpha `$ means both $`+q_+`$ and $`q_{}`$.

$$h_\alpha (𝐫)=\mathrm{e}^{q_\alpha ^2v_c(0)/2},$$
(78)
$$\mathrm{e}^{ı\alpha q_\alpha \varphi (𝐫)}=\mathrm{e}^{q_\alpha ^2v_{\mathrm{DH}}(0)/2},$$
(79)
$$h_\alpha (𝐫)h_\beta (𝐫^{})=\mathrm{e}^{\omega _{\alpha \beta }(𝐫𝐫^{})+\left[q_\alpha ^2v_c(0)+q_\beta ^2v_c(0)\right]/2},$$
(80)
$$\begin{array}{cc}\hfill \mathrm{e}^{ı\alpha q_\alpha \varphi (𝐫)ı\beta q_\beta \varphi (𝐫^{})}& =\mathrm{e}^{\left[q_\alpha ^2v_{\mathrm{DH}}(0)+q_\beta ^2v_{\mathrm{DH}}(0)\right]/2}\hfill \\ & \times \mathrm{e}^{\alpha \beta q_\alpha q_\beta v_{\mathrm{DH}}(𝐫𝐫^{})}\hfill \end{array}$$
(81)
$$\varphi ^2(𝐫^{})\mathrm{e}^{ı\alpha q_\alpha \varphi (𝐫)}=\mathrm{e}^{q_\alpha ^2v_{\mathrm{DH}}(0)/2}\left[v_{\mathrm{DH}}(0)q_\alpha ^2v_{\mathrm{DH}}^2(𝐫𝐫^{})\right].$$
(82)
The brackets $`\mathrm{}`$ denote averages over the fields $`\varphi `$ and $`\psi _\alpha `$ where $`\omega ^1`$ and $`v_{\mathrm{DH}}^1`$ are the propagators.

### The coefficients in the grand-canonical free energy

We give here the explicit expressions for the coefficients $`\stackrel{~}{g}`$, Eq. (32). Using the function $`H(x)`$, defined in Eq.(49), the coefficients read

$$m_1=\frac{\pi }{3}\frac{q_+^6\mathrm{}_B^3}{d_+^3}\left\{H\left(\frac{q_+^2\mathrm{}_B}{d_+}\right)+\mathrm{ln}\left(\frac{3d_+}{\mathrm{}_B}\right)\right\},$$
(83)
$$m_2=\frac{\pi }{3}\frac{q_{}^6\mathrm{}_B^3}{d_+^3}\left\{H\left(\frac{q_{}^2\mathrm{}_B}{d_{}}\right)+\mathrm{ln}\left(\frac{3d_{}}{\mathrm{}_B}\right)\right\},$$
(84)
$$m_3=\frac{2\pi }{3}\frac{q_+^3q_{}^3\mathrm{}_B^3}{d_+^3}\left\{H\left(\frac{2q_+q_{}\mathrm{}_B}{d_++d_{}}\right)\mathrm{ln}\left(\frac{3[d_++d_{}]}{2\mathrm{}_B}\right)\right\},$$
(85)
$$n_1=\frac{\pi }{3}\frac{q_+^8\mathrm{}_B^3}{d_+^3}\left\{\frac{5}{8}2H\left(\frac{q_+^2\mathrm{}_B}{d_+}\right)+\mathrm{ln}\left(\frac{12d_+^2}{\mathrm{}_B^2}\right)\right\},$$
(86)
$$n_2=\frac{\pi }{3}\frac{q_{}^8\mathrm{}_B^3}{d_+^3}\left\{\frac{5}{8}2H\left(\frac{q_{}^2\mathrm{}_B}{d_{}}\right)+\mathrm{ln}\left(\frac{12d_{}^2}{\mathrm{}_B^2}\right)\right\},$$
(87)
$$\begin{array}{cc}\hfill n_3& =\frac{2\pi }{3}\frac{q_+^3q_{}^3\mathrm{}_B^3}{d_+^3}\{\frac{5}{8}q_+q_{}+\frac{[q_+q_{}]^2}{2}H(\frac{2q_+q_{}\mathrm{}_B}{d_++d_{}})\hfill \\ & +q_+q_{}\mathrm{ln}\left(\frac{2[d_++d_{}]}{\mathrm{}_B}\right)\frac{q_+^2+q_{}^2}{2}\mathrm{ln}\left(\frac{3[d_++d_{}]}{2\mathrm{}_B}\right)\},\hfill \end{array}$$
(88)
$$p_1=\frac{\pi }{3}\frac{q_+^6\mathrm{}_B^3}{d_+^3},p_2=\frac{\pi }{3}\frac{q_{}^6\mathrm{}_B^3}{d_+^3},p_3=\frac{2\pi }{3}\frac{q_+^3q_{}^3\mathrm{}_B^3}{d_+^3},$$
(89)
$$r_1=\frac{2\pi }{3}\frac{q_+^8\mathrm{}_B^3}{d_+^3},r_2=\frac{2\pi }{3}\frac{q_{}^8\mathrm{}_B^3}{d_+^3},$$
(90)
$$r_3=\frac{2\pi }{3}\frac{q_+^3q_{}^3\mathrm{}_B^3}{d_+^3}\left[q_+q_{}\frac{q_+^2+q_{}^2}{2}\right],$$
(91)
$$s_1=\frac{q_+^4}{8},s_2=\frac{q_{}^4}{8},t_0=\frac{d_+^3}{12\pi \mathrm{}_B^3},t_1=\frac{q_+^6}{48},t_2=\frac{q_{}^6}{48}.$$
(92)
### The coefficients in the fugacity

In Eq. (39) we define the fugacity of negative ions $`\stackrel{~}{\lambda }_{}`$ as an expansion in terms of $`\stackrel{~}{\lambda }_+`$. The coefficients $`a_0`$, $`a_1`$, etc. are determined by the condition of global electroneutrality. They are given by

$$a_0=\frac{q_+}{q_{}},$$
(93)
$$a_1=\frac{q_+^2}{q_{}}\left[q_+^2q_{}^2\right]\sqrt{\frac{\pi \mathrm{}_B^3}{d_+^3}\left[1+\frac{q_{}}{q_+}\right]},$$
(94)
$$\begin{array}{c}a_2=\frac{q_+}{q_{}^2}\left[2m_2q_++2m_1q_{}+m_3[q_1q_2]\right]+\frac{\pi }{6}\frac{q_+^2\mathrm{}_B^3}{q_{}d_+^3}\left[q_++q_{}\right]^2\hfill \\ \hfill \times \left[7q_+^39q_+^2q_{}+2q_{}^3\right]+\frac{\pi }{3}\frac{q_+^2\mathrm{}_B^3}{q_{}d_+^3}\left[q_+^5q_+^3q_{}^2+q_+^2q_{}^3q_{}^5\right]\\ \hfill \times \mathrm{ln}\left(\frac{4\pi \mathrm{}_B^3}{d_+^3}q_+\left[q_++q_{}\right]\right),\end{array}$$
(95)
$$a_3=\frac{\pi }{3}\frac{q_+^2\mathrm{}_B^3}{q_{}d_+^3}\left[q_+^5q_+^3q_{}^2+q_+^2q_{}^3q_{}^5\right],$$
(96)
$$\begin{array}{c}a_4=\frac{\mathrm{}_B^{3/2}}{24q_{}^3}\sqrt{\pi q_+\left(q_++q_{}\right)}\hfill \\ \hfill \times \{24n_1[5q_+q_{}^2q_{}^3]+24n_2[q_+^35q_+q_{}^2]\\ \hfill +72n_3[q_+^2q_{}q_+q_{}^2]+24m_1[q_+^2q_{}^33q_+q_{}^4]+24m_2[4q_+^4q_{}\\ \hfill q_+^3q_{}^2+7q_+^2q_{}^3]+12m_3[2q_+^4q_{}q_+^3q_{}^26q_+^2q_{}^3+5q_+q_{}^4]\\ \hfill +\frac{\pi \mathrm{}_B^3}{d_+^3}q_+^2q_{}^2[q_++q_{}]^2[26q_+^534q_+^4q_{}31q_+^3q_{}^2+67q_+^2q_{}^3\\ \hfill 45q_+q_{}^4+17q_{}^5]\}+\frac{\pi ^{3/2}}{6}\frac{q_+^{5/2}\mathrm{}_B^{9/2}}{q_{}d_+^{9/2}}[q_+q_{}\left]\right[q_++q_{}]^{5/2}\\ \hfill \times [10q_+^411q_+^3q_{}+13q_+^2q_{}^24q_+q_{}^3+3q_{}^4]\\ \hfill \times \mathrm{ln}\left(\frac{4\pi \mathrm{}_B^3}{d_+^3}q_+\left[q_++q_{}\right]\right),\end{array}$$
(97)
$$\begin{array}{c}a_5=\frac{\pi ^{3/2}}{6}\frac{q_+^{5/2}\mathrm{}_B^{9/2}}{q_{}d_+^{9/2}}[q_+q_{}][q_++q_{}]^{5/2}[10q_+^411q_+^3q_{}\hfill \\ \hfill +13q_+^2q_{}^24q_+q_{}^3+3q_{}^4]\mathrm{ln}\left(\frac{4\pi \mathrm{}_B^3}{d_+^3}q_+[q_++q_{}]\right).\end{array}$$
(98)
Notice that $`m_1`$, $`m_2`$, etc. were defined in Eqs. (8392).

### The coefficients in the mean activity coefficient

In Eq. 73 we obtained the low-density expansion of the mean activity coefficient of ionic solutions where the ions have valences $`q_+`$ and $`q_{}`$ and effective diameters $`d_+`$ and $`d_{}`$ (in Ångströms). Defining $`\omega 6.022\times 10^4\mathrm{}_B^3q_+^6q_{}`$ (with $`\eta q_{}/q_+`$ and $`\xi d_{}/d_+`$), the coefficients in Eq. 73 read

$$\nu _{\mathrm{DH}}=\eta \sqrt{\pi \omega [1+\eta ]},$$
(99)
$$\begin{array}{c}\nu _1=\frac{\pi \omega \eta }{6[1+\eta ]}\{1+3\eta +8\eta ^2+3\eta ^3\eta ^4+4[H\left(ϵ_+\right)+\mathrm{ln}\left(ϵ_+\right)]\hfill \\ \hfill +4\eta ^4\left[H\left(\frac{\eta ^2ϵ_+}{\xi }\right)+\mathrm{ln}\left(\frac{\eta ^2ϵ_+}{\xi }\right)\right]8\eta ^2\left[H\left(\frac{2\eta ϵ_+}{1+\xi }\right)+\mathrm{ln}\left(\frac{2\eta ϵ_+}{1+\xi }\right)\right]\\ \hfill +8\eta ^2[1\eta ^2]\mathrm{ln}(\eta )[24\eta ^2+2\eta ^4]\mathrm{ln}\left(36\pi \omega [1+\eta ]\right)\},\end{array}$$
(100)
$$\nu _{1\mathrm{l}\mathrm{o}\mathrm{g}}=\frac{\pi \omega \eta }{3}[1\eta ]^2[1+\eta ],$$
(101)
$$\begin{array}{c}\nu _{3/2}=\frac{\eta [\pi \omega ]^{3/2}}{24\sqrt{1+\eta }}\{426\eta 14\eta ^2+68\eta ^314\eta ^46\eta ^5+42\eta ^6\hfill \\ \hfill +40[1\frac{2\eta }{5}][H\left(ϵ_+\right)+\mathrm{ln}\left(ϵ_+\right)]+40\eta ^6[1\frac{2}{5\eta }][H\left(\frac{\eta ^2ϵ_+}{\xi }\right)\\ \hfill +\mathrm{ln}\left(\frac{\eta ^2ϵ_+}{\xi }\right)]+112\eta ^3[H(\frac{2\eta ϵ_+}{1+\xi })+\mathrm{ln}\left(\frac{2\eta ϵ_+}{1+\xi }\right)]\\ \hfill 16\eta ^3\left[72\eta ^2+5\eta ^3\right]\mathrm{ln}(\eta )+8\eta \left[1\eta ^2\right]^2\mathrm{ln}\left(36\pi \omega [1+\eta ]\right)\\ \hfill 20[1+\eta ^3]^2\mathrm{ln}\left(64\pi \omega [1+\eta ]\right)\},\end{array}$$
(102)
$$\begin{array}{c}\nu _{3/2\mathrm{l}\mathrm{o}\mathrm{g}}=\frac{\eta [\pi \omega ]^{3/2}}{6}\sqrt{1+\eta }\hfill \\ \hfill \times \left[57\eta +7\eta ^2+7\eta ^37\eta ^4+5\eta ^5\right].\end{array}$$
(103)