# Solving non-linear ODE for divalent solution at a 1-D surface boudary

I am trying to solve the following equation for a positively charged plane with charge density $\sigma$ at $z = 0$.

$$\phi''(z)=-\frac{e}{\epsilon \epsilon_0} \big(z_+n_{+} e^{-\beta z_+ e\phi}-z_{-} n_{-} e^{\beta z_- e\phi} \big)$$

So far I have determined that the electric field and potential must go to zero at infinity, and that, assuming there is no electric field at $z < 0$, I have the following boundary condition:

$$\phi'(0) = -\frac{\sigma}{\epsilon\epsilon_0}$$

I am not sure how to proceed with solving the equation given the two nonlinear terms.

• I'm a bit confused by your notation. On the right hand side of your first equation, I think you are using $z_-$ and $z_+$ to mean the valency/charge of the free ions, while $z$ is your coordinate. Is that correct? – Colin McFaul Apr 25 '14 at 14:32
• Yes, that is correct. – user34433 Apr 25 '14 at 19:35

## 1 Answer

The Poisson-Boltzmann equation is principally valid when dealing with 1:1 electrolytes i.e. $z_- = -1$ and $z_+ = +1$.Then, as soon as divalent ions enter the system things become a bit more tricky and people (physicists at least) might jump at your throat if you use Poisson-Boltzmann in this case. I will therefore assume in the following that you do have a 1:1 electrolyte.

The usual thing is to use dimensionless variables for the field. if you denote $\psi \equiv \beta e \phi$ then you get a rewriting of the PB equation as:

$\psi'' = \kappa^2 \sinh \psi \:\:\:\:\:\:\:\:\:\:\:\:\:\:\:\:(1)$

with the boundary condition $\psi'(0) = -4\pi l_B \overline{\sigma}$

where $\overline{\sigma}$ is the surface charge density expressed in electric charge $e$ per unit length squared, $l_B \equiv \frac{e^2}{4\pi \epsilon \epsilon_0 k_B T}$ is the Bjerrum length expressed in SI units and $\kappa^2 = 8 \pi l_B n_b$ is the square of the inverse Debye length where $n_b = n_+ = n_-$ is your bulk salt concentration.

Before trying to solve the non linear case, it is important to notice that $\psi$ tells you the typical electrostatic potential energy of an ion with respect to $k_B T$. Basically, if $\psi$ is small everywhere, then your charged plate doesn't perturb very much the homogeneous distribution of ions or at least it does so in a linear way.

As a matter of fact, the above rewritten PB equation becomes $\psi'' \simeq \kappa^2 \psi$ in this case which is called the Debye-Huckel equation or the linearized PB equation. Its solution is then a trivial decaying exponential.

Let us go back to the full non linear problem for one plate.

The idea is to solve it by quadrature. To do so we first notice that if we multiply (1) by $\psi'$ then we get:

$\left(\frac{1}{2}(\psi')^2 - \kappa^2 \cosh \psi \right)' = 0$

This implies that:

$\frac{1}{2}(\psi')^2 - \kappa^2 \cosh \psi = Cte = E$

I named this constant $E$ because if you are used to nonlinear mechanics then you call energy the function which, upon derivation with respect to the evolution variable, gives you back the equation of motion.

Let us assume there won't be any sign change in $\psi'$, we can then write (assuming the surface charge density is positive):

$\frac{d\psi}{dz} = -\sqrt{2E + 2\kappa^2 \cosh \psi}$

which leads to

$dz = -\frac{d\psi}{\sqrt{2E + 2\kappa^2 \cosh \psi}}$

Upon integration we get:

$z = -\int_{\psi=\psi(0)}^{\psi(z)} \:\frac{d\psi}{\sqrt{2E + 2\kappa^2 \cosh \psi}}$

In the general case, $E$ can be anything and the solution has to be expressed as an elliptic integral (which is what is written above) with $E$ as a parameter that has be found by solving an implicit non linear equation derived from the boundary conditions. And therefore, even the simple problem of an electrolyte confined in between two plates is still an active field of research in mathematical physics and colloidal science.

In our case though, we implicitly chose the gauge potential so that both $\psi$ and $\psi'$ vanish at infinity. Because $E$ has a constant value throughout space, it means that it has to be $-\kappa^2$ everywhere since it is $-\kappa^2$ at infinity. We therefore get the simpler integral:

$z = -\int_{\psi=\psi_0}^{\psi(z)} \:\frac{d\psi}{\sqrt{2}\kappa \sqrt{-1+ \cosh \psi}} = \frac{1}{\kappa}\ln \coth \frac{\psi(z)}{4} - \frac{1}{\kappa}\ln \coth \frac{\psi(0)}{4}$.

Let us denote $\kappa z_0 = \ln \coth \frac{\psi(0)}{4}$ we therefore have:

$\kappa(z+z_0) = \ln \coth \frac{\psi(z)}{4}$

we now use the fact that $\coth \frac{x}{2} = \frac{1+e^{-x}}{1-e^{-x}}$ to rewrite the above equation as:

$e^{\kappa(z+z_0)} = \frac{1+e^{-\psi(z)/2}}{1-e^{-\psi(z)/2}}$

Doing the algebra step by step yields:

Step 1- $e^{\kappa(z+z_0)}(1-e^{-\psi(z)/2})=1+e^{-\psi(z)/2}$

Step 2- $-e^{-\psi(z)/2}(1+e^{\kappa(z+z_0)})=1-e^{\kappa(z+z_0)}$

Step 3- $e^{-\psi(z)/2} = \frac{e^{\kappa(z+z_0)}-1}{1+e^{\kappa(z+z_0)}}$

Step 4- $e^{-\psi(z)/2} =\frac{e^{-\kappa(z+z_0)}}{e^{-\kappa(z+z_0)}} \frac{e^{\kappa(z+z_0)}-1}{1+e^{\kappa(z+z_0)}} = \frac{1-e^{-\kappa(z+z_0)}}{1+e^{-\kappa(z+z_0)}}$

Final step- $\psi(z) = 2 \ln \frac{1+e^{-\kappa(z+z_0)}}{1-e^{-\kappa(z+z_0)}}$

Which is the general solution to the problem. The unique solution corresponding to the problem is found by looking for the value of $z_0$ that satisfies the boundary condition at the beginning.

• Why does your expression have a minus sign in front $$\frac{d\psi}{dz} = -\sqrt{2E + 2\kappa^2 \cosh \psi},$$ when I solve for $\psi(z)$ I obtain $$\frac{d\psi}{dz} = \sqrt{2E + 2\kappa^2 \cosh \psi}$$How did you get $$z = -\int_{\psi=\psi_0}^{\psi(z)} \:\frac{d\psi}{\sqrt{2}\kappa \sqrt{-1+ \cosh \psi}} = -\frac{1}{\kappa}\ln \coth \frac{\psi(z)}{4} + \frac{1}{\kappa}\ln \coth \frac{\psi(0)}{4},$$ explicitly, how did you perform this integration? It seems your integration variable is $\psi$ which is also your limit of integration. – user34433 Apr 25 '14 at 16:14
• And also, is $\psi_0=\psi(0)$? Also, how do you then write $$\psi(z) = 4\ln \frac{1+e^{-\kappa z+ \kappa z_0}}{1-e^{-\kappa z+\kappa z_0}}?$$ What I obtain is $$z=-\frac{1}{\kappa} \ln \coth \frac{\psi(z)}{4}+z_0,\to \psi(z)=4\coth^{-1}\big(e^{-\kappa(z-z_0)} \big).$$ – user34433 Apr 25 '14 at 16:17
• $d\psi/dz$ is negative because of the boundary condition if I assume a positively charged plate. $\psi(0)=\psi_0$ yes and for the integral I was lazy to do it by hand, I therefore used wolfram alpha and it gave me a result a bit more complicated which included $\sqrt{\cosh x -1}/\sinh (x/2)$ which is equal to $\sqrt{2}$ for all positive values of $x$. For the last part, I may have omitted a factor 2 but the idea is that $\frac{1+e^{-x}}{1-e^{-x}}=\frac{e^{x/2}}{e^{x/2}} \frac{1+e^{-x}}{1-e^{-x}}=\frac{e^{x/2}+e^{-x/2}}{e^{x/2}-e^{-x/2}}=\coth \frac{x}{2}$ – gatsu Apr 25 '14 at 18:24
• I have one more question, about a possible sign error, thanks again. For the integral I obtain (using mathematica) $$z=-\frac{1}{\sqrt 2 \kappa}\int_{\psi_0}^{\psi(z)} \frac{dx}{\sqrt{\cosh x-1}}=\frac{1}{\kappa} \log \coth \frac{x}{4}\bigg|^{x=\psi(z)}_{x=\psi_0}=\frac{1}{\kappa}\left(\log \coth \frac{\psi(z)}{4}-\log \coth \frac{\psi_0}{4} \right).$$ However you have $$z=-\frac{1}{\kappa}\left(\log \coth \frac{\psi(z)}{4}-\log \coth \frac{\psi_0}{4} \right).$$ Thanks a lot! – user34433 Apr 25 '14 at 19:23
• you're right I have the wrong sign there. I will try to edit my answer during the weekend to patch it up of all these little mistakes towards the end. Do not hesitate to tell me if you have other concerns. – gatsu Apr 25 '14 at 20:37