Thermodynamic decomposition of Debye–Hückel activity coefficients: resolving attractive, repulsive, and entropic components

Abstract

Ionic activity coefficients are central to the thermodynamics of electrolyte solutions, yet their thermodynamic structure at the level of individual ions is rarely made explicit. In this work, Debye–Hückel theory is formulated within classical partial molar thermodynamics, starting from the linearised Poisson–Boltzmann free energy functional. The resulting framework provides an explicit decomposition of the excess free energy into attractive and induced-repulsive energetic contributions, as well as an entropic contribution. These contributions are thermodynamically consistently split and defined at the level of excess partial molar quantities. This formulation reveals physical information that remains implicit in the conventional global free energy construction, while preserving the classical Debye–Hückel result. The derived expressions reproduce both the limiting law and extended Debye–Hückel equations, and establish a transparent thermodynamic connection between the mean-field description of ionic atmospheres and activity coefficients. More generally, the framework provides a systematic basis for analysing energetic and entropic contributions to electrolyte thermodynamics within mean-field theory.

---

1. Introduction

Activity coefficients quantify deviations from ideal behaviour that, in electrolyte solutions, arise primarily from interionic interactions. Since the pioneering work of Debye and Hückel,¹ the linearised Poisson–Boltzmann (PB) theory has provided a remarkably successful mean-field framework for describing these effects in the dilute limit. Despite its well-known limitations, the Debye–Hückel (DH) theory remains a cornerstone of electrolyte theory and continues to be a reference point for more elaborate models.^2,3

The DH activity coefficient can be derived from an excess free-energy functional within the linearised PB approximation. Although the classical DH theory is usually derived via mean-field electrostatics and Boltzmann statistics,^1,4 more formal treatments based on equilibrium statistical mechanics also exist. In particular, Kirkwood and Poirier provided an early statistical-mechanical foundation for DH theory that relates its approximations to systematic expansions of the potential of mean force in ionic systems.⁵ Consequently, much of the subsequent development of electrolyte theory has focused on extending the underlying PB framework to incorporate additional physical effects, such as asymmetry in ion sizes,⁶ ion association,⁷ relative permittivity decrement,^8,9 or electrostatic correlations.¹⁰

In standard derivations,¹¹ the activity coefficient of a given ion is obtained from the excess electrostatic free energy associated with the interaction between a central ion and its surrounding ionic atmosphere. Within this global construction, key aspects of the energetic and entropic contributions inherent in the mean-field description remain implicit when expressed in macroscopic thermodynamic form. The literature contains only occasional references to specific thermodynamic aspects of the DH free energy. For example, the entropic term, obtained as the temperature derivative of the limiting DH expression, has been adopted intermittently as a reference benchmark when analysing more complex theoretical models.^12,13

In the present work, we address a long-standing conceptual gap in the classical DH description of activity coefficients by analysing their thermodynamic origin at the level of individual ions. Although the DH result is well established, its internal thermodynamic structure has, to the best of our knowledge, not yet been examined in a systematic and explicitly decomposed form. The key conceptual step adopted here differs from conventional approaches in that it exploits the full spatial information contained in the radial charge density distribution around a central ion, rather than reducing the ionic atmosphere to an effective electrostatic potential acting at the ion position. This distinction is essential: while the latter representation captures the net ion–atmosphere interaction, it obscures the internal electrostatic structure of the ionic atmosphere itself. By retaining the complete charge-density description, both the direct interaction between the central ion and its atmosphere and the induced interactions within the ionic atmosphere can be treated on an equal footing.

The formulation is developed from the outset within the rigorous framework of excess partial molar quantities in classical thermodynamics, where the excess chemical potential of a species is defined as the derivative of the free energy with respect to particle number in the thermodynamic limit.¹⁴ This approach enables a direct mapping between the mean-field electrostatic description and thermodynamic quantities at the single-ion level. In particular, it provides a natural route for decomposing the excess chemical potential into distinct contributions with clear physical meaning.

A central feature of this construction is that it separates contributions that are only implicitly included in the conventional free energy functional. Specifically, the induced electrostatic interactions within the ionic atmosphere, arising from the redistribution of charge in response to the central ion, are identified explicitly alongside the direct ion–atmosphere interaction. When combined with the corresponding entropic contribution obtained from the temperature dependence of the free energy, this decomposition yields a transparent energy–entropy balance within the DH framework.

While the sum of all contributions recovers the classical DH expression for the activity coefficient, their explicit separation reveals the underlying thermodynamic structure and clarifies the relative roles of attractive, induced-repulsive, and entropic effects. Because the analysis is grounded in general principles of partial molar thermodynamics, it also provides a conceptual framework that could potentially be extended to more elaborate PB–based models, including those incorporating spatially varying dielectric response, excluded-volume effects, or alternative geometries.^10,15

The remainder of this article is organized as follows. Section 2 introduces the thermodynamic framework and notation. Section 3 presents the partial molar construction of excess electrostatic energy, Helmholtz free energy, and entropy within the DH model. Section 4 presents the results of the decomposition graphically and discusses differences from classical derivations. The conclusions, which include plans for future work and possible extensions, are given in Section 5.

2. Theoretical framework and definitions

2.1. Description of electrolyte solutions within Debye–Hückel theory

We adopt a partial molar formulation, expressing all quantities per particle. This choice enables a consistent and transparent decomposition of excess thermodynamic quantities, including Gibbs and Helmholtz free energies, electrostatic energy, and entropy, within a single conceptual framework, without introducing any additional thermodynamic assumptions. The excess partial molar quantities used here therefore provide a fully equivalent formulation of the chemical potential that is particularly well suited to the present analysis. For notational simplicity, all partial molar quantities are expressed per particle rather than per mole; they therefore differ from their conventional definitions only by a factor of the Avogadro constant N_A.

While classical thermodynamics determines the excess chemical potential uniquely, its decomposition into distinct energetic and entropic contributions can be explicitly determined within the linearised PB mean-field framework. This is possible once the excess free energy is expressed in terms of partial molar quantities with respect to the conventional ideal-solution reference state (a hypothetical 1 mol dm⁻³ ideal solution) and the standard charging process is applied. Activity coefficients are accordingly defined relative to this reference state. The excess (partial molar) quantities therefore represent the difference between the real solution and the hypothetical ideal solution at the same composition, temperature and pressure¹⁶ (or volume, when the Helmholtz free energy is considered). This formulation allows the interaction contributions to the chemical potential to be identified unambiguously while retaining the conventional definition of activities. In the calculations below, the evaluated excess chemical potential thus corresponds to the interaction contribution relative to infinite dilution, whereas the activity coefficient itself refers to the 1 mol dm⁻³ ideal reference state.

We consider a strong electrolyte solution described within the PB framework. The solvent is treated as a homogeneous dielectric continuum with concentration-independent permittivity ε₀ε, where ε₀ is the vacuum permittivity and ε is the relative permittivity of the solvent. Ions are modelled as spherical particles of finite radius, defining a distance of closest approach a, or, alternatively, as point-like charges. Any ion j, characterized by charge q_j = z_je₀ (z_j – charge number of ion j; e₀ – elementary charge) and radius r_j, may be chosen as the central ion, while all remaining ions form its ionic atmosphere. The ionic atmosphere is represented by a spherically symmetric continuous charge density, corresponding to a statistical average over microscopic ionic configurations.

Within the DH framework, electrostatic interactions between ions are described by the linearised PB equation. For sufficiently dilute solutions of strong electrolytes, where the electrostatic interaction energy q_jψ_IA (ψ_IA – electric potential of the ionic atmosphere interacting with ion j) is small compared to the thermal energy k_BT, the PB equation can be linearised. This approximation is valid primarily for dilute solutions of monovalent ions and progressively deteriorates with increasing concentration or ionic charge.

In the extended Debye–Hückel (EDH) formulation, a distance of closest approach a between the central ion and the ionic atmosphere is introduced. This accounts approximately for finite ion size effects and thus extends the concentration range where the equation can be effectively used.This distance a is commonly interpreted as the sum of the radii of oppositely charged ions. In the domain of very dilute solution, one may simplify the equation by neglecting finite-size of ions, thus entering in the domain of point-like ions (a = 0), governed by the Debye–Hückel limiting law (DHLL).

Solution of the linearised PB equation introduces the screening parameter κ, defined as¹

where I is the ionic strength,In the equation above, c_i and z_i denote the molar concentration and charge number of ion species i, respectively.

The inverse screening parameter κ⁻¹, often referred to as the Debye screening length, characterizes the spatial extent of the ionic atmosphere. The thermodynamic properties of dilute electrolyte solutions therefore depend primarily on the ionic strength.¹⁷

For distances r ≥ a, the total electrostatic potential around the central ion, ψ_tot(r), arising from all ions in the solution, is given by the classical DH expression¹ expressed for EDH formulation as

and asfor DHLL formulation.

The exponential decay reflects electrostatic screening by the ionic atmosphere and limits the validity of the DH theory to dilute solutions characterized by a single decay length.¹⁰

2.2. Separation of central-ion and ionic-atmosphere contributions

By virtue of the additivity (superposition) of the electrostatic potential, the total potential ψ_tot(r) can be decomposed into the contribution of the isolated central ion, ψ_j(r), ψ_j(r) = z_je₀/4πε₀εr, and that of the ionic atmosphere, ψ_IA(r),¹

ψ_tot(r) = ψ_j(r) + ψ_IA(r). (4)

Using eqn (3a) and (4), the electrostatic potential of the ionic atmosphere for r ≥ a can be written explicitly asIn the following, equations are written in the EDH form; the DH limiting law is recovered by setting a = 0.

Since the charge density of the ionic atmosphere vanishes for r < a, Gauss’ law implies that the electric field generated by the ionic atmosphere inside the exclusion sphere is zero and the potential ψ_IA is therefore constant for r ≤ a.

Evaluating ψ_IA at r = a gives

The central ion is located at r = 0 so eqn (6) defines the potential of the ion atmosphere that interacts with the charge of the central ion. The term κ/(1 + κa) = (κ⁻¹ + a) may be interpreted as an effective distance at which a point charge of opposite sign would reproduce the electrostatic potential of the ionic atmosphere at the position of the central ion. In dilute solutions, where κ⁻¹ ≫ a, this distance reduces approximately to the Debye screening length.

The electrostatic energy of the direct attractive interaction between the central ion and its ionic atmosphere is therefore given by

where ψ_IA(0) denotes ψ_IA(r = 0).

This interaction energy provides the energetic basis in standard derivations of the Debye–Hückel activity coefficient.^1,4 In the following sections, we depart from these conventional treatments and instead analyse the thermodynamic contributions underlying the activity coefficient in a more explicit manner.

3. Calculation of the activity coefficient from excess partial molar properties

3.1. Activity coefficient of a central ion

Activity coefficients quantify deviations from ideal solution behaviour. Within the DH framework, nonideality arises exclusively from ion–ion electrostatic interactions. The activity coefficient of a spescies j, γ_j, is related¹⁸ to its excess partial molar Gibbs free energy ,where G^ex stands for excess Gibbs free energy of the system, N_j and N_i denote the numbers of particles of species j and all other species, respectively.

While the activity coefficient is defined in terms of derivatives at constant pressure and temperature, the explicit calculation is carried out in the canonical (N, V, T) ensemble using excess partial molar Helmholtz free energies. As emphasized by Hill, partial molar quantities acquire a rigorous meaning only in the thermodynamic limit, where different thermodynamic potentials give equivalent excess partial molar quantities up to well-defined volume contributions.¹⁴ Neglecting electrostriction effects, the excess partial molar Gibbs free energy may therefore be identified with the excess partial molar Helmholtz free energy, .¹⁹

The latter can be obtained from the excess partial molar electrostatic energy via a standard charging process,¹⁹

where λ denotes a coupling parameter that uniformly scales the ionic charges.

3.2. Excess partial molar electrostatic energy of a central ion

To elucidate the origin of the excess partial molar electrostatic energy of a central ion, , we consider the formal insertion of a single ion j into an infinitely large volume of electrolyte solution characterized by a screening parameter κ (Fig. 1). This insertion is a thermodynamic construction employed to evaluate partial molar quantities in the thermodynamic limit and does not correspond to a physical particle exchange process. The insertion is a conceptual device used to isolate the contribution of a single ion. It does not represent a statistical-mechanical test-particle insertion in the sense of the Widom method,^20,21 although both constructions ultimately address a partial molar free energy contribution. This approach enables identification of induced-repulsive electrostatic interactions in the ionic atmosphere, existing in the model DH solution. These repulsive interactions are already implicitly included in the solution of PB equation but are practically never explicitly treated when DH equation is discussed. Yet, the explicit recognition of these induced interactions as well as of entropic changes, as illustrated schematically in Fig. 1, are needed when thermodynamic decomposition of DH equation is sought.

Fig. 1 (a) Time-averaged picture of a model DH electrolyte solution around a non-perturbing reference point. The instantaneous image (inset) shows one of possible configurations that can be expected in such a system. The time-averaged charge density throughout the system is 0. The time-averaged electric potential at the reference point is equal to the time-averaged electric potential at all other points in the system (here set to zero). The overall excess electrostatic energy of the system is E^ex₁, and the overall excess Helmholtz free energy is A^ex₁; (b) electrolyte solution around the positively charged central ion located at r = 0. A typical instantaneous image (inset) taken at a sufficiently large radius r again shows one of possible configurations that can be expected in a model DH electrolyte solution. The concentrically shaded area around the reference point at r = 0 indicates the time-averaged charge density of the negative charge from the ionic atmosphere. The higher the charge density, the darker the area. The overall excess electrostatic energy of the system is now E^ex₂, and the overall Helmholtz free energy is A^ex₂.

Prior to insertion, the charge is uniformly distributed all over the system and no local excess charge exists. Consequently, in the time-averaged sense the electrostatic potential throughout the solution is equal. Upon insertion of the ion, a spherically symmetric ionic atmosphere forms around the inserted ion, while the bulk properties of the solution remain unaffected due to the infinite system size.

The change in excess electrostatic energy associated with this formal insertion of ion j defines the excess partial molar electrostatic energy of the inserted (central) ion,

where E^ex₁ and E^ex₂ denote the excess electrostatic energies of the system before and after insertion, respectively, both defined relative to the corresponding ideal, non-interacting reference state. At this point we reiterate that the physical insertion of the chosen ion into the system is not required at all to obtain (and in continuation ). The same picture (or result) would be obtained simply by comparing the time-averaged view of the system from the position of the central ion with that from a non-perturbing reference point that does not disturb the rest of the system.

The difference between E^ex₂ and E^ex₁ may be decomposed into two distinct electrostatic contributions: (i) the direct attractive interaction between the central ion and its ionic atmosphere, u_j-IA, and (ii) the repulsive electrostatic interactions within the ionic atmosphere itself, u_IA–IA,

The first contribution is given by eqn (7). The second contribution is induced by the presence of the central ion. It arises from the electrostatic energy associated with the continuous charge distribution in the ionic atmosphere:

This integration is restricted to r ≥ a, i.e. outside the exclusion sphere of the central ion, which is defined by the contact distance a.

Using the linearised PB equation for r ≥ a, the charge density ρ(r) of the ionic atmosphere can be expressed as

Evaluation of eqn (12) using eqn (5) and (13) givesSummation of u_j-IA (eqn (7)) and u_IA–IA (eqn (14)) gives the excess partial molar electrostatic energy of the central ion:The excess partial molar electrostatic energy of the central ion, , is conceptually and numerically distinct from the effective central ion–ionic atmosphere interaction energy u_j-IA commonly used in standard derivations of the excess energy of DH solution. The central ion–ionic atmosphere energy represents the effective electrostatic interaction energy of the central ion with its surrounding ionic cloud. Half the sum of these energies for all ions in the system gives the electrostatic energy of the system. In contrast, represents the excess electrostatic energy of the last ion added to the system. In the following section, is used as the starting point for the calculation of the corresponding excess partial molar Helmholtz free energy.

3.3. Excess partial molar Helmholtz free energy and activity coefficient

Application of the charging process (eqn (9)), following the approach used by Marcus,¹⁹ gives the excess partial molar Helmholtz free energy of ion j, , via charging process on excess quantities. A similar route for calculating the excess free energy from the excess electrostatic energy has been employed by, for example, Varela et al.²² (eqn (17) in the cited work). In our case, can be decomposed into contributions arising from the interaction between the central ion and its ionic atmosphere, ,and from the self-interaction of the ionic atmosphere, ,Summing of eqn (16) and (17) gives the overall excess partial molar Helmholtz free energy:which coincides with the classical DH result.

The excess partial molar entropy can be now obtained thermodynamically as the temperature derivative of the excess partial molar Helmholtz free energy (eqn (18))

This result for the excess partial molar entropy satisfies also the thermodynamic relation connecting the Helmholtz free energy (eqn (18)), the electrostatic energy (eqn (15)) and the entropic contribution (eqn (19))This confirms thermodynamic consistency of our decomposition of the DH free energy within the framework of excess partial molar quantities.

The excess partial molar entropy (eqn (19)) reflects the loss of configurational entropy of the system. In the ideal reference state, ions are randomly distributed, while in the model electrolyte solution, the presence of the central ion imposes spatial correlations, reducing the number of accessible microstates. The entropic penalty is a direct measure of the ordering effect of the electrostatic field.

Finally, the individual activity coefficient of ion j, taking into account the considerations outlined in Section 3.1, follows directly from the excess partial molar Helmholtz free energy, :

It can be expressed as a sum of contributions associated with the direct central ion–atmosphere interaction, the self-interaction of the ionic atmosphere, and the effective entropic term,

ln γ_j = ln γ_j(j-IA) + ln γ_j(IA–IA) + ln γ_j(entr), (22)

where the individual terms are given by the following equationsThis decomposition clarifies the distinct energetic and entropic contributions to the nonideality of dilute electrolyte solutions. It retains as the fundamental thermodynamic quantity within the linearised PB framework, for both the DHLL and EDH equations. Interestingly, the sum of the repulsive energetic contribution (eqn (24)) and the entropic contribution (eqn (25)) to the activity coefficient is exactly one half (in absolute value) of the attractive energetic contribution (eqn (23)). This observation applies to both the DHLL and EDH formulations.

4. Results and discussion

4.1. Decomposition of the Debye–Hückel activity coefficient

The central result of the present analysis is an explicit decomposition of the Debye–Hückel activity coefficient into three distinct contributions: (i) a direct attractive interaction between the central ion and its ionic atmosphere, (ii) an induced-repulsive electrostatic interaction within the ionic atmosphere itself, and (iii) an entropic contribution arising from the formation of the ionic atmosphere around the central ion (Fig. 2). These contributions are denoted in Fig. 2 by subscripts in γ_j(x), where (x) = (i), (ii), or (iii), while (iv) refers to the total activity coefficient obtained as their sum. This decomposition provides direct physical insight into ionic non-ideality. It shows how competing energetic and entropic effects combine to produce the observed activity coefficients. In the following, the decomposition is first illustrated and then analysed as a function of the screening parameter κ and ion size.

Fig. 2 Decomposition of contributions to the individual activity coefficient γ_j of a central ion j at 25 °C in an aqueous 1 : 1 electrolyte solution as a function of the DH screening parameter κ. The notation (x) in γ_j(x) denotes the source of each contribution. To provide better insight into how individual contributions affect the activity coefficient, the figure shows their direct impact on γ_j rather than on ln γ_j (eqn (23)–(25)). The value κ = 0.2 nm⁻¹ corresponds approximately to a 0.004 mol dm⁻³ solution of 1 : 1 electrolyte, while κ = 1 nm⁻¹ corresponds approximately to a 0.1 mol dm⁻³ solution. (a) Point-ion limit (a = 0). The individual contributions are shown from (i) attractive interactions between the central ion and its ionic atmosphere, (ii) induced-repulsive electrostatic interactions within the ionic atmosphere, and (iii) the entropic component arising from the formation of the ionic atmosphere, together with (iv) their sum, which corresponds to the conventional DH activity coefficient. (b) Same as in (a), but for ions of finite size with a distance of closest approach a = 0.5 nm. In this case, the repulsive and entropic contributions no longer coincide.

Fig. 2 shows the decomposition of the excess partial molar Helmholtz free energy for an individual ion in terms of γ_j, derived from the solution of the linearised PB equation. Besides the well-known attractive interaction between the central ion and its ionic atmosphere, the partial molar framework reveals two additional contributions: an induced-repulsive electrostatic term from interactions within the ionic atmosphere, and an entropic term due to the redistribution of mobile charges around the central ion in response to its electric field.

In the point-ion limit (a = 0), the induced-repulsive electrostatic and entropic contributions are exactly equal across all ionic strengths (Fig. 2(a)). Both increase the excess free energy, partially counteracting the attractive central ion–atmosphere interaction, which lowers the free energy. This equality arises naturally from the linearised PB framework when excess free energy is expressed in partial molar terms, reflecting the quadratic structure of the mean-field free-energy functional. The attractive contribution dominates in absolute terms, yielding activity coefficients below unity that increasingly deviate from ideality as ionic strength rises.

When ions are treated as finite-sized particles, the absolute magnitudes of all contributions decrease, with the reduction becoming more pronounced as the distance of closest approach a increases. Fig. 2(b) illustrates this behaviour for a = 0.5 nm. In contrast to the point-ion case, the induced-repulsive and entropic contributions no longer coincide. The entropic term, which arises from the loss of configurational entropy due to ionic ordering (Section 3.3), exhibits the strongest response to increasing ion size. This sensitivity on a arises because the maximal charge density of the ionic atmosphere diminishes as the contact distance increases; the surplus charge becomes more evenly distributed, reducing the induced entropic change. Conversely, the induced-repulsive electrostatic interaction within the ionic atmosphere shows the weakest dependence on a.

The relationships between the contributions of attractive, induced-repulsive, and entropic terms – plotted in Fig. 3 and 4 as a function of the screening parameter κ – remain constant only in the point-ion limit (a = 0). These relationships vary when finite ion sizes are considered, with energetic and entropic contributions developing distinct dependences on κ and the distance of closest approach a.

Fig. 3 (a) Relative contributions of attractive (full line, eqn (7)) and repulsive (dashed line, eqn (14)) electrostatic interactions compared to the overall excess partial molar electrostatic energy (eqn (15)) as function of κ, for varying distances of closest approach a. (b) Similarly to (a), but now the ratio of the electrostatic energies for the attractive and repulsive interactions is shown.

Fig. 4 (a) Ratio of the entropic contribution to excess partial molar Helmholtz free energy as a function of κ and contact distance a. (i) Point-ion limit (a = 0); (ii) a = 0.1 nm; (iii) a = 0.2 nm; (iv) a = 0.5 nm. (b) Ratio of the entropic contribution to excess partial molar electrostatic energy under the same conditions.

Fig. 3(a) illustrates this variation by showing the relative magnitudes of the attractive and induced-repulsive electrostatic interaction energies, normalised by the overall excess partial molar electrostatic energy, for different values of a.

For point-like ions (a = 0), the relative energies of the attractive and induced-repulsive contributions are fixed at 4/3 and −1/3, respectively, independent of κ, reflecting the analytical structure of the DH limiting law. For finite ion sizes, increasing κ at fixed a or increasing a at fixed κ causes modest deviations from these limiting values, indicating a gradual redistribution of electrostatic energy between attractive and induced-repulsive components.

Fig. 3(b) presents the ratio of attractive to repulsive electrostatic interaction energies. In the point-ion limit, this ratio remains constant at −4, but decreases smoothly with increasing κ and a for finite-sized ions. At κ = 1 nm⁻¹ and a = 0.5 nm, corresponding to moderately concentrated (approximately 0.1 mol dm⁻³) aqueous solutions of 1 : 1 electrolytes, the ratio falls to about −3, demonstrating that finite ion size systematically, but moderately, alters the balance between attractive and repulsive interactions.

As shown in Fig. 4(a), the relative entropic contribution to the excess partial molar Helmholtz free energy decreases with increasing a and κ. This ratio is exactly 0.5 for the point-ion limit at all κ, but falls to approximately 1/3 for a = 0.5 nm at κ = 1 nm⁻¹. Since the excess partial molar electrostatic energy always exceeds the excess partial molar Helmholtz free energy in absolute terms, the ratio of the entropic contribution to the excess partial molar (Fig. 4(b)) is smaller yet follows the same pattern, having value of exactly 1/3 for point-like DH solutions at infinite dilution (κ = 0 nm⁻¹).

4.2. Comparison with other charging processes used to derive Debye–Hückel activity coefficients

Several charging procedures for DH activity coefficients have been developed within the linearised PB framework. Though equivalent under standard assumptions, they differ in their thermodynamic pathways to the excess Helmholtz free energy, leading to distinct physical interpretations (Table 1).

Table 1 Comparison of key elements in standard derivations of the DH activity coefficients with the approach presented here

Charging process	Charging target	Ionic atmosphere model during the charging process	κ (screening parameter)	Self-energy of the central ion	Thermodynamic level
Debye	All ions in the system	Infinitely thin spherical shell	Varies with λ	Subtracted after the charging process	Helmholtz free energy of all ions in the solution
Güntelberg	Single central ion (other ions are precharged)	Infinitely thin spherical shell	Constant	Subtracted after the charging process	Helmholtz free energy of the central ion
MacInnes^23 (shortened Güntelberg; found in most modern textbooks)	Interaction energy within ionic atmosphere (“excess potential of the ions”)	Infinitely thin spherical shell	Constant	Evaded through subtraction before the charging process	Helmholtz free energy of the central ion
present (excess partial molar)	(i) Interaction energy between the single central ion and induced ionic atmosphere; (ii) interaction energy within the induced ionic atmosphere	Full spatial distribution (linear PB solution)	Varies with is	Absent	Helmholtz free energy of the central ion

---

In the original DH construction, all ions in the solution are charged simultaneously. By contrast, the Güntelberg charging procedure charges only a selected central ion, while all other ions remain fully charged and maintain a constant κ. When using these two methods, the charge-integrated values include self-energy, which must then be subtracted at infinite dilution. The remaining interaction contribution is the quantity relevant for the calculation of the activity coefficient. When the solvent relative permittivity is concentration-independent and no volume changes occur, both charging routes (DH and Güntelberg) yield identical activity coefficients for equally sized ions with a common contact distance.¹¹

Our approach begins with the interaction energy of a single central ion, explicitly excluding self-energy contributions. This new approach eliminates the need for subsequent subtraction and provides a more transparent thermodynamic pathway.

The present procedure resembles the Güntelberg approach in that it focuses on a single central ion, but important differences exist. In standard Güntelberg-type derivations (including the shortened version used by, for example, MacInnes²³), the Debye screening parameter κ is held constant during charging at its final value. The interaction between the central ion and its ionic atmosphere is simplified to that with an infinitely thin spherical shell.

In contrast, our method retains the full interaction with the spatially distributed ionic atmosphere throughout the charging process. Here, κ varies with the charging parameter λ, reflecting the gradual development of ionic correlations. This allows the Helmholtz free energy to be decomposed into the three distinct contributions mentioned above.

Recent literature has introduced the so-called “ionic cloud charging process” when discussing DH and Güntelberg charging under conditions where the solvent relative permittivity depends on electrolyte concentration.²⁴ That approach incorporates changes in solvent properties during the DH charging process, in order to reconcile differences between the classical charging routes when the relative permittivity varies with κ.

The present work assumes a constant relative permittivity and an incompressible medium, thereby neglecting solvent–ion interactions. As a result, it is unrelated to the ionic cloud charging process. Instead, our charging procedure explicitly describes the energetic and entropic cost of forming the ionic atmosphere from an initially uniform charge distribution.

5. Conclusions

To the best of our knowledge, this work presents the first systematic dissection of Debye–Hückel theory, a prototypical mean-field description of electrolyte solutions, which has traditionally focused solely on net attractive electrostatic interactions between a central object and its surroundings, by explicitly separating these from the induced-repulsive interactions within the surrounding environment. While the approach presented here, as applied to Debye–Hückel theory, offers a new perspective on the thermodynamic origin of activity coefficients, it may also provide a basis for extension to other mean-field theories.

The extension of this decomposition to more accurate approaches within the restricted primitive model framework – such as Monte Carlo simulations, the hypernetted-chain (HNC) approximation, or the mean spherical approximation (MSA), which yield radial distribution functions – appears feasible.

In this regard, we have performed preliminary calculations using the HNC equation, from which the radial dependence of the electrostatic potential of the ionic atmosphere was obtained via Poisson integration of the charge density constructed from radial distribution functions.²⁵ The results demonstrate not only the feasibility of this approach, but also show that the excess partial molar electrostatic energies of the attractive (u_j-IA) and repulsive (u_IA-IA) interactions derived within the Debye–Hückel framework are in good agreement with those obtained from HNC calculations under conditions where the Debye–Hückel theory is expected to be valid. It should be noted that, when using HNC calculations (and other theoretical methods that provide radial distribution functions), an accurate evaluation of these functions over a sufficiently large radial range is essential, as they approach unity only asymptotically. In particular, careful treatment of the long-range tails is required²⁶ to maintain electroneutrality and to avoid systematic shifts in the separate attractive and induced-repulsive contributions.

With this in mind, we compared energetic contributions obtained using Debye–Hückel theory and the HNC equation for a 1 : 1 aqueous electrolyte at 25 °C with ion radii of 0.1 nm.²⁷ Within the Debye–Hückel framework, the values of u_j-IA and u_IA–IA were −180.1 and 45.9 J mol⁻¹ at 0.001 M, and −545.5 and 144.8 J mol⁻¹ at 0.01 M, respectively, while the corresponding values obtained using the HNC equation were −187.6, 46.1, −597.7, and 147.6 J mol⁻¹. The latter values yield ratios u_j-IA/u_IA–IA of −4.07 at 0.001 M and −4.05 at 0.01 M, which are consistent with those obtained within the Debye–Hückel framework, while suggesting that the variation of this ratio may be weaker in the HNC framework. A more detailed analysis of this behaviour is planned for a separate study.

One important motivation for the thermodynamic analysis presented here arises from studies of polyelectrolyte solutions, where polyion activities are often derived from the non-linear Poisson–Boltzmann equation within the cell model. In these systems, strong electrostatic interactions cause a significant fraction of counterions to remain associated with the polyion, even at very low concentrations.²⁸ Consequently, the reference state for polyelectrolyte solutions differs from that of simple electrolytes.²⁹ Due to this convention, the numerical values of activity coefficients – especially for polyions – often lack direct physical meaning, and discussions usually focus on their concentration dependence.^30,31 Reinterpreting polyion activity coefficients using the thermodynamic decomposition presented here could, in principle, provide further insight into the relative roles of energetic and entropic contributions in polyelectrolyte solutions. However, such an analysis would also need to ensure consistency between thermodynamically and operationally defined entropic contributions when using the non-linear Poisson–Boltzmann equation.³² Additionally, the role of electrostatic repulsive interactions between polyions may require consideration.³³ Given the complexity, this analysis will be addressed in a separate study, which may also examine differences between cylindrical³⁴ and spherical³⁵ systems. In this context, spherically symmetric polyvalent fullerenehexamalonates could be particularly interesting, as the non-linear Poisson–Boltzmann equation provides a reliable estimate of the electric potential around these polyions³⁶ and activity coefficients³⁵ determined from experimentally measured osmotic coefficients³⁷ are available. Alternatively, at least for alkali fullerenehexamalonate salts (1–12 electrolytes) such a thermodynamic decomposition could also be attempted using, for example, a modified HNC equation, in which the correlation functions are corrected to account for the pronounced asymmetry in both the charge and size of the ions.³⁸

While transferring the present excess-partial-molar-quantity approach to polyelectrolyte solutions by solving the non-linearised Poisson–Boltzmann equation within the standard cell model is conceptually straightforward, applying it to more advanced electrolyte theories beyond the standard extended Debye–Hückel formulation³⁹ may prove more challenging. Further extensions of Debye–Hückel theory – such as those explicitly accounting for finite ion size, relative permittivity decrement, or additional solvation effects – could also be analysed using this framework. However, the mathematical complexity of these models would require decomposing significantly more expressions into energetic and entropic contributions. A systematic analysis of these models is beyond the scope of this work but could be explored in future studies after the versatility of the present thermodynamic decomposition has been further established within restricted primitive models.

By bridging electrostatic theory and thermodynamic formalism, this work deepens understanding of electrolyte solutions and demonstrates the versatility of the approach through its successful application to Debye–Hückel theory. These findings suggest that similar decompositions may be applied to other mean-field theories of electrolyte solutions, providing a broader framework for analysing Coulomb interactions in complex charged systems. Exploring such extensions would further help to establish the generality and utility of the proposed methodology.

Conflicts of interest

The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Data availability

The data that support the findings of this study are contained within the article. The results are based on analytical expressions provided in full detail in the manuscript, from which all figures can be reproduced. Representative numerical values obtained from hypernetted-chain (HNC) calculations are included in the text. No additional datasets or code were generated.

Acknowledgements

We acknowledge the support from the Slovenian Research Agency (research core funding no. P1-0201). The author of this study would like to thank his colleagues, Prof. Jurij Reščič and Prof. Miha Lukšič, for fruitful discussions on HNC calculations and for providing radial distribution functions using their HNC simulation codes.

References

1. P. Debye and E. Hückel, Phys. Z., 1923, 24, 185–206.
2. H. Falkenhagen, Theorie der Elektrolyte, S. Hirzel Verlag, Stuttgart, 1971.
3. J. P. Simonin and O. Bernard, Fluid Phase Equilib., 2023, 571, 113805.
4. E. Güntelberg, Z. Phys. Chem., 1926, 123, 199–247.
5. J. G. Kirkwood and J. C. Poirier, J. Phys. Chem., 1954, 58, 591–596.
6. D. Zuckerman, M. Fisher and S. Bekiranov, Phys. Rev. E, 2001, 64, 011206.
7. H. Yokoyama and H. Yamatera, Bull. Chem. Soc. Jpn., 1975, 48, 1770–1776.
8. I. Y. Shilov and A. K. Lyashchenko, J. Phys. Chem. B, 2015, 119, 10087–10095.
9. Q. Lei, B. Peng, L. Sun, J. Luo, Y. Chen, G. M. Kontogeorgis and X. Liang, AIChE J., 2020, 66, e16651.
10. R. Kjellander, Phys. Chem. Chem. Phys., 2020, 22, 23952–23985.
11. G. M. Silva, X. Liang and G. M. Kontogeorgis, Mol. Phys., 2022, 120, e2064353.
12. B. B. Laird and A. Haymet, J. Chem. Phys., 1994, 100, 3775–3779.
13. E. S. Velázquez and L. Blum, J. Mol. Liq., 1997, 73-4, 75–89.
14. T. L. Hill, Thermodynamics of small systems, Dover Publications, Incorporated, Mineola, New York, 2013.
15. R. Netz and H. Orland, Eur. Phys. J. E, 2000, 1, 203–214.
16. H. L. Friedman, J. Chem. Phys., 1960, 32, 1351–1362.
17. G. N. Lewis and M. Randall, J. Am. Chem. Soc., 1921, 43, 1112–1154.
18. K. S. Pitzer, J. Phys. Chem., 1973, 77, 268–277.
19. R. A. Marcus, J. Chem. Phys., 1955, 23, 1057–1068.
20. B. Widom, J. Chem. Phys., 1963, 39, 2808–2812.
21. B. R. Svensson and C. E. Woodward, Mol. Phys., 1988, 64, 247–259.
22. L. Varela, M. García and V. Mosquera, Phys. A: Stat. Mech. Appl., 2003, 323, 75–87.
23. D. A. MacInnes, The principles of electrochemistry, Reinhold Publishing Corporation, New York, 1939.
24. G. M. Silva, X. D. Liang and G. M. Kontogeorgis, J. Chem. Eng. Data, 2023, 69, 291–306.
25. T. S. Sørensen, Mol. Simul., 1993, 11, 267–303.
26. J. Zhang, D. J. Searles and T. Duignan, J. Chem. Theory Comput., 2024, 20, 6957–6970.
27. E. Gutiérrez-Valladares, M. Lukšič, B. Millán-Malo, B. Hribar-Lee and V. Vlachy, Condens. Matter Phys., 2011, 14, 33003.
28. A. Zhang, X. Yang, F. Yang, C. Zhang, Q. Zhang, G. Duan and S. Jiang, Molecules, 2023, 28, 2042.
29. A. Katchalsky, S. Lifson and J. Mazur, J. Polym. Sci., 1953, 11, 409–423.
30. N. Ise and T. Okubo, J. Phys. Chem., 1965, 69, 4102–4109.
31. D. Dolar, Z. Phys. Chem. N. F, 1968, 58, 170–180.
32. D. A. McQuarrie, Statistical Mechanics, Harper & Row, New York, 1976.
33. M. Muthukumar, Macromolecules, 2017, 50, 9528–9560.
34. D. Dolar, in Polyelectrolytes, ed. E. Sélégny, M. Mandel and U. P. Strauss, D. Reidel, Dordrecht, 1974, vol. 1, pp. 97–113.
35. J. Škerjanc, Acta Chim. Slov., 2012, 59, 559–563.
36. J. Cerar and T. Urbic, J. Phys. Chem. B, 2008, 112, 12240–12248.
37. J. Cerar and J. Škerjanc, J. Phys. Chem. B, 2003, 107, 8255–8259.
38. J. Reščič, V. Vlachy and A. Haymet, J. Am. Chem. Soc., 1990, 112, 3398–3401.
39. G. M. Kontogeorgis, B. Maribo-Mogensen and K. Thomsen, Fluid Phase Equil., 2018, 462, 130–152.

---