Statistical theory of individual activity coefficients of electrolytes including multiple ionic charges

In previous work we developed a new statistical method for calculating the individual activities of ions including the association of ions. Here we study multi-particle electrostatic interactions connected within higher cluster integrals and identify the ionization constants of the mass action law of associating ion clusters. In contrast to Bjerrum and Fuoss, our concept of association is not based on spatial criteria, but instead on the strength of interaction measured in powers of the Bjerrum parameter ($e^2 / D_0 k_{\text{B}} T a$ ; $a$ is contact) and defined by asymptotic properties of the cluster integrals. For ion pair formation our mass action constant is the classical counterpart of Planck's famous hydrogenic partition function. As a rule, the new association constants are smaller than traditional expressions, e.g., by Fuoss and Kraus, in the interesting regions of interaction parameters about fifty percent. Several examples including CaCl$_2$, MgCl$_2$, Na$_2$SO$_4$, K$_2$SO$_4$, LaCl$_3$ and a model of seawater are studied. For several associating electrolytes and seawater, reasonable agreement with experiments and Monte Carlo results is achieved.


Introduction
The problem of individual ionic activities, osmotic coefficients and ion association are key problems in the theory of electrolytes [1][2][3][4][5][6][7][8]. Usual approaches to association effects in electrolytes with multiple charged ions are based on the classical concepts of Bjerrum, Fuoss and Kraus [1], which define pairs and triples as special spatially defined configurations. Our concept of electrostatic association is completely different from that, we do not use spatial criteria, but the strength of interaction measured in powers of the parameter 0 = ( 2 / 0 B ) ( is contact distance of ions). We follow the concepts of Onsager, who stated in 1968 at a conference in Montpellier: "Bjerrums choice is good but we could vary it within reason. In a complete theory this would not matter; what we remove from one side of the ledger would be entered elsewhere with the same effect" [3]. In our concept for the definition of pairs we use this freedom and assume that pair-associates are formed by higher order (negative) contributions of binary charge interactions 0 with ⩾ 4 to the pressure and other thermodynamic functions. Triple and quadruple association is generated by (negative) contributions of three or four charges with opposite signs contributing higher orders in 0 to the pressure. Such a definition of association may seem less transparent in comparison with spatial definitions, although it is easier in the light of statistical thermodynamics. The background of the present work are the cluster expansions for systems with Coulombic interactions based on the fundamental papers of Joseph Mayer since 1949 which were developed in the 50th and 60th by E. Haga, E. Meeron, H. Falkenhagen, G. Kelbg, I. R. Yukhnovskii, H. Friedman and others. Note that the cluster expansions were developed in two versions, based on density and those based on fugacity expansion [9,10]. Both versions are connected with diagrammatic expansions with respect to the interaction strength 2 or in the case of hard charged spheres with the Bjerrum parameter 0 = 2 / 0 B . An important role played personal meetings in Rostock and Lviv between H. Falkenhagen, G. Kelbg, I. R. Yukhnovskii, H. Friedman and the present authors who witnessed these discussions.
The basic concepts which we use here were presented first in [4,5,[11][12][13][14][15] and were developed later in [6,[16][17][18]. These concepts are mainly based on a mathematical analysis of the cluster functions of the Mayer-theory of ionic solutions [4,5,7]. This basic theory of ionic solutions led to cluster expansions, which were rederived and generalized on the basis of collective variables by Yukhnovskii [8] and Kelbg [14]. The connections between cluster expansions and association theory were established in [4,11,12]. The main idea of this concept is that the contributions of ion pairs and ion triples are given by certain relevant parts of the second and third cluster integrals in the pressure expansion. The parts relevant to bound states are in the low temperature asymptotically large and negatively definite, this way strongly decreasing the pressure. These contributions are of higher order in the interaction parameter , e.g., ⩾ 8 for pairs. These concepts lead to quite transparent definitions of the ionization constants for pairs and triples as we showed in the foregoing work [4,5,18,19]. Alternative concepts of electrostatic association were considered in many works [4, 24, 29-33, 39, 40].
The influence of electrostatic pair and triple association of ions is of high relevance for many problems where electrolytes with higher charges such as MgSO 4 , MgCl 2 or Na 2 SO 4 play a role, e.g., in seawater research [1]. In recent works [18][19][20], the present authors summarized the existing theoretical knowledge on the activity coefficients and the individual conductivities of electrolytes in such systems. We responded this way to urgent needs for extending and deepening the knowledge on the physico-chemical properties of the components of seawater and other complex natural and technological electrolytes which include higher charges. For example in seawater in addition to univalent ions, such as Na + , K + , Cl − also the double-charged ions Mg 2+ , Ca 2+ and SO 2− 4 are of importance for seawater properties. In many technologies several multiple-charged ions like the vanadium ions are becoming of substantial interest, e.g., for modern battery development. Here, we develop a new look at the analytical theory of ions with higher charges, where the differences between the individual and the mean activities are large [1,2,6,18,26]. We devote special attention to the consequences of charge asymmetry and higher charges to the individual ionic activities. In previous work, triple association was mainly neglected, while here we give a more systematic approach. In our statistical approach we use the ion densities (in particle numbers per cm 3 ) or the molarities (in moles of ions per dm 3 = liter) as the basic primary quantities. For seawater we also use the chlorine molality ( Cl , in moles of chlorine per mass of pure water), and the salinity ( is the mass of dissolved sea salt per mass of seawater). Using the -scale, the activity coefficients are denoted by and the so-called practical (or molal) activity coefficients are defined using molalities, [6,7] where the molalities are measured in moles per mass of solvent [7]. Further we use the partial osmotic coefficients derived from the osmotic pressures of species . We use standard methods for estimating individual activities and osmotic coefficients as the virial expansion of the thermodynamic functions [4-7, 18, 27, 28, 43]. The approach given here is restricted to ion association problems, where less than about 1/4 of the ions are in associated states. This is due to our quasi-linear approach to the mass action law and to the first order approach to rational pressure expansions. The idea behind our approach is that a complete physical description of any higher interaction orders includes everything and in particular the chemical effects. The problem is, however, to identify the responsible higher order terms and work out these contributions relevant for association. We have described this here in some detail for pairing effects and in less detail for triple and quadruple association effects. A restriction of our method is that the nonlinearities contained already in a simple mass action law, are included in our procedure only in a quasi-linear way. As shown in figure 1 for a simple model, the results for the degree of ionization of individual ions using our simple method which is half way between physical and chemical description, is for weak association in reasonable agreement with the complete mass action law. The price to pay for the simplicity of our method is that we have to stay in the region of weak association effects, say that the degrees of ionization are higher than 75 percent.

Interaction potential and association functions
We follow here basically the concepts developed in [18]. First we define the potential of the mean force between the ions and in the solutions as . The average potentials and forces are split into a Coulombic and a short-range part The electrical part is given by Coulomb's law where ( , ) is the relative dielectric constant of pure water and ℓ the Coulomb length (also called Landau length or with the pre-factor 1/2 the Bjerrum length). Both are functions of temperature and pressure.
In what follows we perform all calculations for the temperature = 298.15 K (i.e., 25 degrees Celsius). and assume for the relative dielectric constant the value = 78.36. Then, we get for the Coulomb length ℓ = 715.4 pm. In our model this is the only parameter which is temperature and pressure dependent, so the transition to other values of , is reduced to changing ℓ. The short range forces are of hard-core type, where are the contact distances. These are the most important key data in our approach. Several values for the contact distance are given in table 1. Most of the values were given already in [18]. We are of opinion that the most important data about contact distances may be obtained from MC and MD simulations as presented, e.g., in [22,24]. We have added a few not so well studied ions as La and Cd. The crystallographic radius for Cd 2+ is with 95 pm just a few pm higher than that for Mg 2+ which is 86 pm. Therefore, we may assume that the contact distances in solution are also close, we took CdCl = 420 pm. For La 3+ , we know that the crystallographic radius is smaller than that for Mg 2+ . Following canonical MC simulations by Valisko and Boda we assume for those ions in water LaLa = 430 pm and LaCl = 270 pm corresponding to a quite large Bjerrum parameter pm = 7.95. The pair association constants given in table1 stem from an earlier theory, the triple association constants from new calculations.
Consequently, we expect the asymptotic This result corresponds to the estimates proposed by Friedman and one of these authors [15] by using the results of mathematical studies of cluster integrals [7,11]. Including the pre-factors, this estimate of the asymptotic reads [7,11,15] ≃ 8π 2 6 This estimate was obtained in Kelbg's early work [11] by using the assumption that the integrands have a cusp at contact, which provides the asymptotic. A qualitative procedure to treat triple association of the type (− + −) or (+ − +) is the method of effective charges as known from the treatment of helium formation as bound state of two electrons and one two times positively charged alpha-particle. The bound 23602-3 state energy of an electron-alpha pair is 2 where = 2 is the charge of the He-nucleus and the corresponding pair association constant is In one of his latest but quite important works, Max Planck showed that for Coulombic systems, here for hydrogen-like bound states, the exponential function should be replaced by a cropped exponential, where the first two terms of the exponential are cropped Plancks partition function begins with the order 8 similar to the Falkenhagen partition function [4]. For a pair of single charged ion and a Z-fold charged ion, the Falkenhagen theory provides with 0 = ℓ/ +− the expression [4] The arguments which lead to this series on the Bjerrum parameter are completely different from Plancks arguments. Therefore, we consider the relations between the quantum theory of Planck and the classical approaches by Bjerrum and Falkenhagen [4] more in detail. The fundamental work of Max Planck published in 1924 is notoriously difficult to understand, which is the reason why this work is not so often cited as the earlier papers by Planck. In order to apply Planck's ideas to classical systems, we come back to our own paper, which is an exercise about Planck's work with more detailed calculations [13]. The question to be studied in detail is, how the definition of the Falkenhagen partition function equation (8) is related to the Planck concept equation (7). In Planck's quantum statistics, the basic assumption is where +− ( ; ) is the bound state part of the quantum pair probability. This term is proportional to the sum of probabilities to find an opposite charge in a bound state around the central charge, which is expressed here by the diagonal part of the density matrix. In order to find the classical counterpart of this function first we find the probability density that a pair (+−) is in distance in states with negative relative energy Here, 0 = √︁ 2 2 / is the momentum where the pair energies change the sign from negative to positive values ( is relative mass). Carrying out the integral over all momenta < 0 , first we get an integral over the error functions and then explicitly the following series [13] The integral over the distances of this complicated function is divergent, so we have to omit according to Plancks arguments the first two terms which correspond to states close to = 0 and treat them together with the free states. Since ( ) ∼ ℓ, these first terms correspond to low orders in the interaction 2 . The remaining higher order terms ( ) +1/2 with > 2 provide us with the wanted bound state contribution

23602-4
The integral over 2 d is now convergent and gives a finite integral which is the direct classical counterpart of Plancks equation (7). Investigating this integral over the radial coordinate in detail, we find out that it has a special property [13] 4π We may conclude for the same integral with the lower limit +− : From the point of physics, this result is remarkable, since it says that Falkenhagen's mass action constant which originally was based on more mathematical arguments than simplicity, is indeed a complete analog of Planck's famous hydrogenic mass action constant. Further this result means that equation (9) is valid in the classical as well as in the quantum-statistical case, i.e., the approach is equivalent for the classical as well as for quantum Coulombic pairs. For the formation of He-like triples, the mass action constant in zeroth order, according to the idea of effective charge, can be approximated by a product where˜≃ 1.8. The decrease of the charge from 2 to 1.8 is due to the screening of an interacting charge pair by the third charge in the triple. Here we transfer this effective charge approach from the helium theory to the 2-1-association problem, assuming an effective charge of the double-charged ions of about˜≃ 1.8. The method of effective charge also works for electrolytes and provides a simple semi-quantitative approach to triple ionization, which provides a bit more flexibility and is in reasonable agreement with the other approaches.
We start here our investigation from the free energy in a physical description of the ions by the canonical ensemble [11,18]. General expressions from statistical thermodynamics for the cluster contributions ( ) to the negative free excess energy of hard charged spheres read [5,7,8,11,12] where ( ) is the so-called ring function in Debye-Hückel approximation and 0 the Onsager relative screening factors defined by The corresponding chemical potential is The sums are to be extended over the species of ions and all orders of clusters . The strong coupling contributions of ions read in the cluster order = 2, 3, 4, ... [5,7,8,11,12,40] (2) = 1 2 Here the strong coupling function which is of higher order ( 4 ) in the interactions is defined by We use here the so-called nonlinear Debye-Hückel approximation for charged hard spheres proposed by Onsager and Fuoss and developed by these authors [3,18,21]. The contributions of the cluster theory given here including a nonlinear DH-screening contain the full 2nd virial coefficient in good agreement with available HNC data. Including the third virial coefficient for associating electrolytes, e.g., for MgSO 4 with = 420 pm brings the theory even to a quite good agreement with MC data up to nearly 1 mol/liter [18]. An earlier comparison of our approximation with MC calculations [35] also gives a rather good agreement. The present approximation gives the values for the free energy which are close to MC results up to about ≃ 1 mol/liter [25]. This way we see that the contributions provided by the 2nd and 3rd virial coefficients are quite relevant for describing two-and three-particle association effects. The idea behind our physical theory of association is as follows: we connect the definition of the associates with the convergent, and in the asymptotic for strong interactions, dominant parts of the cluster integrals. Correspondingly, we define the degree of ions of kind bound in pairs or triples as the asymptotic convergent strong coupling parts of the cluster integrals which lead to big negative contributions to the free energy: , ( ); (2) , ( ) ∼ asy ∫ conv d2 (2) , , ( );

23602-6
Asymptotic and convergent means here that the given integrals should be treated as follows: the integrals should be extended over relative coordinates and performed after identification of the convergent and (at strong association) dominating positive definite parts of the cluster integrals. The results of these operations depend only on the temperatures and provide the corresponding association constants. We underline that this is the central point of our concept of association. The problem is that the asymptotic part is not uniquely defined. We have some asymptotic freedom. However, at this point we may use the freedom which gives us Onsagers statement, that in the ledger of the theory only the observable sums of mass action effects and long range interaction effects are fixed and controllable by experiments [40]. We note that a systematic approach to the derivation of mass action constants from cluster integrals was developed already in [41] on the basis of the grand canonical ensemble.
Here, we will explain our concept on examples. An essential part of the concept is that for Coulombic association, the classical exponential factors, appearing, e.g., in equation (6) should be replaced by some kind of cropped exponential functions which do not contain the first two terms. This will be demonstrated now. For the case of pair formation, we get the estimate of Falkenhagen given already above [4][5][6] (2) ( ) = 8πℓ 3 2 ( ) = 8π 3 Here, the pair association function ( ) is related to the ( ) functions and to the so-called Kirkwood function [4,5]. More difficult is the question how to define the association constant for triple associates as (+) (−−) (+) and (−) (++) (−). We follow mathematical approaches to the pair and triple cluster integrals [15,35]. The idea behind the so-called constant factor approach is as follows: one or more factors ( ) in the integrands of the relevant clusters belong on the case of associating clusters to repulsive interactions. For these factors , the exponent is positive expressing the repulsion of two equally charged ions. The essential point for the integration is that repulsive ions on average cannot come closer than certain minimal distance 2 +− for triples. In any case, these factors expressing repulsion are slowly changing. Fixing the values of the repulsing functions at the most probable values, we may take the factors out of the integral. For triple clusters, we find, e.g., assuming a double charged ion at 1 and and two single charged ions at 2 and 3, the estimate We consider the association constant of a triple consisting of two univalent ions , with a multiple charged ion imbedded in between. The repulsive factor in the cluster integral for 3 ions may be approximated by the first term in the Taylor expansion (ℓ −− /2 +− ) 2 which including a steric 2 provides just a constant factor. Further, in order to have more flexibility, we may introduce an effective chargeĩ n agreement with previous studies of the properties of the triple cluster integrals [11,12,35]. For 2-1 associates, we have˜2 ≃ 3.5;˜≃ 1.8. We see that our definition of mass action constants does not contain contributions in lower orders of 2 , the pair association constant starts with 8 . The contributions of lower orders are missing since otherwise, we would have conflicts between contributions from the mass action terms with the screening contributions, in particular with the extended Debye limiting law [17]. The degree of free ions of kind which is that part of the ions which are not associated in pairs, triples, quadruples etc., is defined as the relation between the number of free ions to the total number of those ions This is the typical mathematical structure we have in the semi-chemical approach suggested first by Justice et al. [31][32][33]. The (3) and (4) we find in our approach as the asymptotically big negative 23602-7 definite contributions from the 3rd and 4th virial coefficients providing relevant contributions to triple and quadruple association. According to our estimates in the previous work, the pair contribution is of order 4 , the triple contribution is of order 10 [19]. The corresponding thermodynamical functions as, e.g., the activities have, as a rule, a typical structure of rational functions which stems from the nonlinear Debye-Hückel and MSA theories. One way to arrive at these structures is an extension of the first corrections to the limiting law for the activities which reads [18] As already mentioned, there is a close relation to an approach developed by Justice et al. [31][32][33]. Note that instead of the relative screening factor 0 , used already by Onsager (however, denoted by the letter ), we worked in some earlier papers [20] with the half quantity = 0 /2. For a binary electrolyte, the calculation of the 1 , 2 is particularly simple We use here Onsagers notation, but changing his letter which we use for the chemical potential. For for some examples of binary electrolyte, the calculation of the 1 , 2 give, e.g., for MgCl 2 , CdCl 2 the pair (1/2, 2/3) and for LaCl 3 the pair (1/4, 3/4). We use here, in the simplest approximation for the zeroth order, the Debye-Hückel-approximation; more advanced possibilities are the Mean Spherical Approximation (MSA) and the related Henderson-Smith approximation (HSA). The nonlinear Debye-Hückel-approximations (DHA) take into account the first-and second-order terms 1 ( ), 2 ( ),˜1( ),˜2( ) [20]. We will not explain the extension of this theory to the MSA approximation in detail and refer to other works [43,45,[47][48][49]. Our favorite approach is based on a generalization of the Henderson-Smith formula for the pair distribution to arbitrary contact distances [18]. According to this approach, the transition from the DHA to the MSA consists formally in replacing the Debye-Hückel parameter by [18]: Note that the more advanced expression includes higher orders in and that this way programs written for DHA may be easily extended to the Henderson version of the MSA. The central point in our theory is taking into account the differences of the diameters of the ions, since the differences between the individual activities and conductivities depend strongly on the differences between the contact distances [18,20]. The degrees of ionization as well as the conductivities depend mainly on the contact distance of oppositely charged ions.

Association theory using rational pressure expansions
In the second section we identified the association contributions as the big positively definite terms in the virial expansions leading to negative contributions to the free energy. In an alternative approach based on virial expansions of the osmotic pressure, we construct now rational representations of the pressure, where the negatively definite terms in pressure expansions are moved to the denominator of rational expressions. We do not need here explicit definitions of mass action constants. For the partial osmotic pressure contributed by the ionic species we get, by introducing the distribution function into the virial formula, The electrical energy density can be expressed by Euler functions [18,20]. We define weakly and moderately coupling, separated from strongly coupling and hard core parts of the osmotic pressure The contribution of weak coupling 0 is given by a Debye-like osmotic function and 1 is a correction for asymmetric ionic systems [18].
The G-functions represent the different orders in and [18]. The first correction˜1 is relevant for asymmetric ionic systems and leads to logarithmic terms The strong coupling contributions to the pressure sc = sc2 + sc3 + sc4 + ... which are all strictly negative, stem from the second, third, forth, etc. virial coefficients obtained from the virial expansion as derivatives of the free energy cluster contributions in the virial formula given above [18,34]. The strong coupling contributions increase strongly with 2 since = ( 4 ). As a result, the pressure may eventually become negative which would be unphysical. Quite formally we may, in order to avoid negative values of , transform the expression by replacing 1 + sc which can get negative values by the strictly positive Padé-like expression 1/(1 − sc ) which would correspond to a geometric series and leads to The first term is now interpreted as the bound part of the ions of species i, which means that in the present approximation there appears, as a degree of ionization, the expressioñ The˜( ) are interpreted as degrees of association of -clusters. This leads finally to a rational virial representation of the osmotic pressure This way without changing the accuracy in the linear order in density ( ), we arrive at an mathematically equivalent but in the context of association theory quite useful rational expressions for the pressure. We note that we may construct more powerful Padé-like rational expressions by including the higher association coefficients in nominator and denominator [50]. Finally, we note that the terms in the denominator may be interpreted as the effective degrees of association, Here, we included all contributions of higher order (not only their asymptotic) into the expression for the degree of association, which makes some difference to the results of the previous section. Consequently, the expressions for and the new here˜as well as the old˜( ) and the new˜( ) are not identical, although they are in effect very close and converge to each other at small densities, the differences being within the "Onsager freedom in a ledger". We remember the statement that chemical species are not uniquely defined, but depend weakly on the context, i.e., on the imbedding into a theory [40]. The expressions for˜have a quasi-chemical meaning, since they may be interpreted as degrees of association. Formally, these new expressions may be obtained by a nonlinear extension of the canonical ensemble expansions assuming that these terms are part of a geometric series. The correctness of such an extension is justified by the studies within the grand-canonical ensemble [15,50].
The terms˜2,˜3 stemming from second, third virial terms are here included into the strong coupling sc since they describe the association effects and give in our new approach the "degree of association" of pair formatioñ Here, the circle is closing since in the lowest approximation restricting to 4 -term, we get a result close to the previous result for a weak binding: For small densities, we find by summing up the series in , an expression valid for both versions of (2) -definitions for the pair formation which provides the degree of low-density ionization: This means that we have at low density a full agreement with the results of the previous section. We remember that 2 ( ) was previously discussed and expressed by ( ) and the Kirkwood function. An application of the expressions (40) and (41) to a 6-component model of seawater is given in figure 2.
The pressure approach leads to alternative expressions describing a weak pair association. It does not include assumptions about symmetries of the charges and may be applied to any neutral mixture of charges. For symmetrical systems, the odd terms cancel out and we get back familiar expressions used in many papers [4,29]. The theory presented so far does not include triple and higher association. For MgSO 4 in standard seawater, the degree of association was measured based on the attenuation of sound by Fisher [51]. Based on these data, Fisher concluded that about 9.2 percent of the total Mg in seawater exists as MgSO 4 . Kester and Pytkowicz estimated the mass action constants and for found for Mg, Ca a degree of association of about 10 percent [52]. Our calculation gives for normal ocean salinity for Mg a degree of association of about 12 percent, for SO 4 about 9 percent and for Ca about 4 percent. Looking at the uncertainties of the experiment and the theory, the agreement seems to be satisfactory. Now, we propose a similar procedure for calculating the degree of association of triples˜( 3) . The easiest way to find an estimate for the 3-particle association function is the effective charge method. For the special case of a (+) (2−) (+) or (−) (2+)(−) triple association, we find the degree of association of triples just as a product two -particle terms modifying the charge. Further corrections arise by introducing some factors for the interaction of equal charges. For example, for triples formed by an ion is with some higher charge and the effective charge˜. In the lowest weak binding approximation, which includes only the terms of order 10 : We note that in the order 10 this new pressure approach is compatible with equation (27). According to this estimate, the triple association is as a rule, e.g., for seawater electrolytes, only a small correction, and the pair association dominates. A first application of our model to the triple association of ions in sea water was presented also in figure 2 (right-hand panel). We estimated the degree of triple association in seawater for the triples MgCl 2 and Na 2 SO 4 . We find, according to out estimates, that the formation of triples in seawater is rather seldom. The predicted degrees of association to MgCl 2 and Na 2 SO 4 obtained from our estimate equation (43) are in the range of 10 −4 to 10 −3 . Qualitatively, the shape seems to be reasonable. The maxima or minima, respectively, located as a rule near salinities between 10 and 15, demonstrate the typical screening effects. A stronger screening destroys Coulombic binding.
In the next section we are going from our extended physical approach as given by equations (40)(41)(42)(43) where association constants do not appear in explicit way, to an alternative semi-chemical picture where association constants ... -s are used explicitly.
Summarizing, we developed here an "extended physical description" of electrostatic association based on rational expressions including higher cluster integrals. This way we also include "binding effects" in some order and find agreement with the semi-chemical formulae found in the last section. Our recipe is in brief as follows: (i) Identify in the usual virial expansion of the pressure the "binding contributions" of plus-minus interaction which provide big negative and for big 0 asymptotically dominant effects. As we have shown such contributions are contributed by terms of order 4 +− for pairing and by 2 −− 8 +− or 2 ++ 8 +− for corresponding triple formation. (ii) Bring the "binding contributions" to the denominator anticipating that the corresponding geometric series stem from the grand-canonical ensemble [9,15,50].
We mention that the relation between the representations in the canonical and in the grand-canonical ensemble, which is not discussed here in detail, is one of the keys for a deeper understanding of binding and non-binding contributions in cluster expansions. The idea is that strong attracting contributions are better represented in the grand ensemble and repulsive contributions are better represented in the canonical ensemble. The term "better" means here, that the convergence of the series is improved [50].
Here, we were able to treat association without defining any mass action laws and mass action constants. Further, the present method predicts degrees of ionization and degrees of association for pairs˜ ( 2) and triples, etc., for each of the ions. This allows us to describe a variety of electrolytes with multi-charged ions including models of seawater. A comparison of our theory for the corresponding activities with MC calculations by Ulfsbo et al. [42] shows a reasonable agreement. As a comment, we mention that the meaning of the word "extended physical description" is here that we specify in a physical description the main higher order contributions which are responsible for binding effects. Our idea is that a complete physical description includes everything and in particular also chemical effects. However, the problem is in this view, to identify and work out the contributions relevant for association or chemical effects. We have described here this approach in some detail for pairing effects. An extension to triple formation cannot be given at the same level of rigor due to the lack of detailed studies of the 3rd virial coefficient. Our comparison with the existing data provided some hints to the existence of triples. We found that triples may be responsible for some corrections up to a few percent.
This new approach allows us to go further in improving a semi-chemical description and to formulate the mass action laws. We will not go here the full way up to a complete chemical description. Instead, we go the half way to a more simple semi-chemical description which works with simplified mass action laws and is comparable or even equivalent to our extended physical description. Our simplified semi-chemical approach is an approximation which works for the case that association effects are weak. The treatment of strong association in a fully chemical approach without simplification of the mass action law was discussed in [16,29].

Semi-chemical treatment of pair and triple association
Here, we show that the extended pressure approach suggested above is indeed half-way to a chemical description. We study again only the model of charged hard spheres and follow the concepts developed previously in [4,5,11,15,29,[31][32][33].
The main effect of pair and triple association is the decrease of effective ion numbers from to where is the degree of ionization. The following associates are of the main interest These estimates show that we should not expect the formation of triples in the case of univalent ions, since the formation of two separate pairs gives a lower energy than the formation of one triple. However, in the case of divalent ions, the formation of triples is of advantage. For the mass action constant of such triples, we expect the asymptotic These results correspond to early estimates proposed by Kelbg, Friedman and these authors [11,12,15] by using the results of mathematical studies of cluster integrals [7,11,34,35]. This estimate was derived by using the assumption that the integrands have a sharp cusp at ion contact, which provides the asymptotic. This is an approximate approach to triple ionization which gives at least the correct asymptotic.
We consider now the configurations of 4 ions like in MgCl 2 solutions. In the case of quadruple formation, we estimate an amount of 8 negative energy units stemming from 4 attractive Mg 2+ -Cl − interactions. Positive contributions come from the repulsion of the equally charged ions on opposite edges which are approximately at a distance √ 5 +− . This estimate shows that these 4 ions would get in a solution more energy namely 4 0 by forming two pairs instead of a quadruple due to the loss by positive repulsive contributions. A first estimate of the energy of a LaCl 3 -quadruple with La in the center and 3 Cl at the edges gives Evidently, the formation of this triple is of some advantage in comparison to forming a pair with energy 3 0 or a triple with energy 6 − (1/2) 0 = 5.5 0 . For a transition from a chemical to a semi-chemical approach, we first discuss the pair association. In the standard approach [4,5,18,29], the degree of ionization is given by the classical mass action law which we write in a form which is appropriate for iterative solutions The semi-chemical approach works with the first iteration following the zeroth approximation (0) = 1 The range of validity ends if more than 1/4 of the ions are associated (see figure1). This is a strict assumption which still leads to a great simplification of the mathematics and is justified for many interesting systems as, e.g., seawater. For the activity coefficients ± which appear in the mass action law, we use the standard expressions for the electrical parts We use here the so-called opposite-charge approximation (OPA), which is a specific property of Coulombic systems, based on the fact that in the region of stronger interactions (larger Bjerrum parameters), the encounter of opposite charges dominates [50]. The main effect of pair and triple association is the decrease of effective particle numbers which leads to a decrease of the osmotic coefficients and the conductivities. We solve the MAL given above by iteration beginning with = 1. Our first and linear approximation, which works only for the regions where is close to one, i.e., we are close to a full ionization, is equivalent to the expressions in our semi-physical approximation. This interesting result means that physical and chemical expressions meet here at half way. In case the association includes more than two ions, the MAL for the degrees of ionization looks as follows: Again, the first order solution may be obtained by iteration starting with = 1 in the denominator. This way, by extending our treatment of pair formation to higher association in first semi-chemical approximation, we get the general expression 23602-13 where the are the activities in electrical approximation: The definition of a mass action constant introduced above, applies to all charge-symmetrical ionic associates including, e.g., Mg 2+ -SO 2− 4 . We have shown above that the association constant of pair formation starts with the order 8 and for triple formation with 20 . The restriction to these lowest orders gives the so-called weak binding constants of one multiple charged ions > 1 with univalent ions (2) , ( ) = π 2 24 2 These results correspond to the findings within our pressure approach in the previous section. We see that both approaches agree in the lowest approximation. Introducing a new association function ( ) = ( )/ 3 , we may write down the results from the second section for triples in constant factor approximation (2) We study now some properties of the association functions useful for practical calculations [19].
The mass action function 2 ( ) which determines the connection between the interaction parameters and the association constants for for pair and triple formation increases first nearly linearly in +− then for +− < 11 like the given polynomial and for +− > 11 the mass action constant starts growing exponentially. We may define an 3 -function which depends on the ion parameters on 4 parameters in a complicated way A graphical representation of the triple association functions depending on the plus-minus Bjerrum parameter is shown in figure 4. Note that we need two functions for triple association in the symmetric  and in the asymmetric case for different combinations of the charges. For the relevant asymmetric charge combination with and without an approximation for the factor of repulsion, leads to the expression A graphical representation of the functions 2 ( ), 3 ( ) is given in figure 4. Since a symmetrical combination of 3 charges is energetically not favorable we concentrate here on the asymmetrical configurations (+) (2−) (+) and (−) (2+) (−). Figure 4 shows that for triple ionization, the constant factor approximation and the effective charge approximation approximately agree for < 6; this is the region where triple association is relevant for our applications, e.g., to K 2 SO 4 , Na 2 SO 4 and MgCl 2 , CaCl 2 . In what follows we combine both methods which leads to the curves in between.

Discussion of associating electrolytes with applications to K 2 SO 4 , CaCl 2 and LaCl 3
Here, we treat the pair and triple association in aqueous solutions based on the full nonlinear second and third virial coefficients. We may expect a corresponding contribution stemming from the strong coupling and negative definite part of the 4th virial coefficient. Relevant for pair association are large negative parts from the 2nd virial coefficient, and important for triple association are the big negative definite contributions from the 3rd virial coefficient. Following the general results from statistical thermodynamics [5,7,8,11,34,35] we have shown that the key quantities for association are the asymptotically dominant parts of the strong coupling terms in the cluster integrals. Several examples Table 1. Table of contact distances for several ion pairs including alkaline earth metal ions, sulfate ions and adapted "ideal" seawater ions according to [18] with a new value for CaCl = 500. In the last but one column we give our new estimates for the ionization constants of several electrostatic pairs and triples in water at 25℃, The values for 3 and 3 are corrected by a factor due to the effect of charge asymmetry.  of the resulting values of the association constants are given in table 1. Note that in comparison to earlier work [19] we introduced some corrections improving the dependence on the ℓ ++ , ℓ +− , ℓ −− parameters. We compared here two different ways of estimating the triple association constant, the constant factor approximation and the effective charge approximation. Both methods are in reasonable agreement for +− < 6, and then they start to disagree. In the constant factor approximation, the triple association constants for salts like K 2 SO 4 , CaCl 2 and ions like (LaCl 2 ) + are estimated again expanding the repulsive factor which leads to the expressions Several numerical values of mass action constants for pair and triple formation estimated this way are given in  [24]. Again, the agreement of our results with the data of other workers for CaCl 2 is sufficient, although not quantitative. For LaCl 3 , our theory is not capable of reproducing the pronounced minimum around √ ∼ 0.4(mole/liter) 0.5 predicted in [24]. The reason is possibly that we did not include so far the fourth virial coefficient. We underline that our approach provides also the individual activities.
For the fourth cluster integral and the corresponding association constants, a consequent statistical analysis is still missing. Therefore, we restrict ourselves here to an estimate following the lines valid for the third virial coefficient. In order to estimate the quadruple association constant for a salt like LaCl 3 , we fix the Cl-Cl-distances at some energetically favorite distance like √ 3 LaCl , which we guess from figure 2. Then, the factors for the repulsive terms may be taken out, which leads to factorization of the remaining terms in the integral and again we get approximating the repulsive term by the quadratic order We estimated several values of the quadruple association constant following equation (63) and found now ≃ 1/3456. These and other values for the association constant that we obtained this way are given in   [18] and [19]. We added a few not so well studied ions as Cd and La. The crystallographic radius for Cd 2+ is with 95 pm just a few pm higher than that for Mg 2+ which is 86 pm. Therefore, we may assume that the contact distances in a solution are also close, we took CdCl = 420 pm. For La 3+ , we know that the crystallographic radius is smaller than that for Mg 2+ . Following canonical MC simulations by Valisko and Boda (2017, 2018), we assume for those ions in water LaLa = 430 pm and LaCl = 270 pm corresponding to a quite large Bjerrum parameter ± = 7.95. Several applications to the binary electrolytes CaCl 2 and LaCl 3 are presented in figure 6. So far, the agreement with available data is not yet quantitative [24,53]. Note that we differ between association constants ,... in the density scale of statistical mechanics particle number/cm 3 , and ... in the usual chemical concentration scale mol/liter. In order to switch between the figures in both of the scales, we remember the relations [ −3 ] = 6.023 · 10 20 [mol/liter]; 2 = 6.022 · 10 −4 2 .
Using these factors of recalculation we derived the numbers given in the last column of table 1.

Conclusions
In the present survey we discuss ion association and activity coefficients in 1-1, 1-2, and 2-2 electrolytes. For 1-3-electrolytes, we restrict ourselves to some qualitative analysis. The association constants for triple and quadruple association are estimated in the effective charge and constant factor approximation. We repeat and correct several fully analytical results to calculate the degrees of weak association and activity coefficients [18,19]. In particular, we discuss extended physical methods based on rational extensions of pressure virial expansions, which do not use mass action laws in an explicit way. Further, we use semi-chemical methods which define the mass action constants for weak association but simplify the mass action laws. This way the restriction to weak association allows us to avoid the use of full (nonlinear) mass action laws. We estimate the ionization constants for pair and triple electrostatic association from the cluster integrals and calculate the degree of ionization for electrolytic mixtures including seawater in the regions of weak association, in general 20-30 percent smaller. Several applications to the electrolytes CaCl 2 and LaCl 3 are presented in figure 6.

23602-17
Summarizing our findings: based on the results of statistical physics, we recommend in addition to standard methods of calculating the individual activities, new statistical tools for the calculation of individual and mean activities and degrees of association of ions from lower concentrations up to moderate concentrations/salinity. The methods are based on the model of hard spheres with non-additive radii in combination with the nonlinear Debye-Hückel (or mean spherical) approximations for screening. Association effects are included by rational virial expressions for the osmotic pressure, which take into account higher order terms stemming from the grand canonical ensemble. These methods avoid the solution of nonlinear mass action laws and use instead rational polynomials which include terms from the fugacity expansions of the pressure. For the first time, pair and triple association is taken into account on same footing in a systematic way; an extension to quadruple association is proposed. A table of most relevant association constants for of electrolytes including seawater ions is given. We use for all calculations only hard-charged sphere models and, as parameters, the charges and non-additive contact distances. The proposed formulae are fully analytical and results can be obtained on a home computer. We stay within the traditional physical and chemical approaches and develop a rational virial approach to the pressure and a semi-chemical description of associates. In the purely physical approach to the osmotic pressure based on rational expressions avoiding explicit definitions of bound pairs or triples. In our semi-chemical approach, explicit association constants are defined and calculated; this is based on the concept that binding is due to higher powers in the interaction determining the asymptotic properties of the cluster coefficients. In both approaches, we concentrate on the relevant contributions of strong Coulombic interactions to association through the higher virial coefficients which determine association effects. Note that both approaches are in full agreement at low densities. The input parameters needed in our theory are, beside dielectric constants, the contact distances of all ion pairs.