Ion-ion association is lost by linearizing the Poisson-Boltzmann equation when deriving the Debye-Hückel equation

In this work, we demonstrate how the ion association constant can be attributed to the difference between the full Poisson-Boltzmann equation and its linearized version in very dilute solutions. We follow a pragmatic approach first by deriving an analytical approximated solution to the Poisson-Boltzmann equation, then calculating its respective Helmholtz free energy and activity coefficient, and then finally comparing it to the contribution from the mass action law principle. The final result is the Ebeling association constant. We conclude that electrostatic ion-ion interaction models miss the ion association contribution naturally introduced in higher-order electrostatic theories. We also demonstrate how the negative deviations from the Debye-Hückel limiting law can be physically attributed to the ion association phenomenon.