Solution (source code)

= Solution

Let $\varepsilon=\varepsilon_0\varepsilon_r$ be the <permittivity>, $\beta=(k_BT)^{-1}$, $e>0$ the <elementary charge>, and $c$ the bulk salt formula-unit <number density>. The <Boltzmann distribution> gives
$$
n_{\mathrm{Mg}}(r)=c\,e^{-2\beta e\phi(r)},\qquad
n_{\mathrm{Cl}}(r)=2c\,e^{\beta e\phi(r)}.
$$
Consequently the exterior <Poisson-Boltzmann equation> is
$$
\frac1{r^2}\frac d{dr}\left(r^2\frac{d\phi}{dr}\right)
=-\frac{2ec}{\varepsilon}\left(e^{-2\beta e\phi}-e^{\beta e\phi}\right),\qquad r>a.
$$
For small $Q$, meaning $2\beta e|\phi(a)|\ll1$, linearize the exponentials. The <Debye–Hückel approximation> gives
$$
\rho_{\mathrm{ion}}\simeq-6\beta e^2c\phi,\qquad
\frac1{r^2}(r^2\phi')'=\kappa^2\phi,\qquad
\boxed{\kappa^2=\frac{6e^2c}{\varepsilon k_BT},\quad \lambda_D=\kappa^{-1}.}
$$
This is also $\kappa^2=\beta e^2\sum_i z_i^2 n_{i,\infty}/\varepsilon$: the divalent cation contributes $4c$ and the two monovalent anions contribute $2c$. In terms of molar salt concentration $C$ in $\mathrm{mol\,m^{-3}}$, replace $c$ by $N_A C$; if $C$ is in $\mathrm{mol\,L^{-1}}$, use $c=10^3N_A C$. The molar <ionic strength> is $3C$ in the same concentration units.

The boundary conditions are $\phi(\infty)=0$ and, by <Gauss's law>,
$$
-\varepsilon\phi'(a^+)=\frac{Q}{4\pi a^2}.
$$
Here the surface charge is uniform and the spherically symmetric interior field is zero. Set $u=r\phi$; the radial equation reduces to $u''=\kappa^2u$. Decay excludes the growing exponential, so write $\phi=B e^{-\kappa(r-a)}/r$. Then $\phi'(a)=-B(1+\kappa a)/a^2$, giving
$$
\boxed{\phi(r)=\frac{Q}{4\pi\varepsilon\,r(1+\kappa a)}e^{-\kappa(r-a)}.}
$$
It solves the \i[linearized] <Poisson-Boltzmann equation>, not the full nonlinear equation at finite $Q$. The radial ion profiles, consistently to first order, are
$$
\boxed{n_{\mathrm{Mg}}(r)\simeq c[1-2\beta e\phi(r)],\qquad
n_{\mathrm{Cl}}(r)\simeq2c[1+\beta e\phi(r)].}
$$
The exponential <Boltzmann distribution> formulas above express the full mean-field relation, but substituting a linearized potential into them is accurate only to the retained order. A useful <electroneutrality> check is
$$
4\pi\int_a^\infty\rho_{\mathrm{ion}}r^2dr
=-4\pi\varepsilon[r^2\phi']_a^\infty=-Q.
$$
The result is the <Debye-Hückel potential of a charged sphere>.