Effect of Ionic Correlations on the Surface Forces in Thin Liquid Films: Influence of Multivalent Coions and Extended Theory

Experimental data for the disjoining pressure of foam films stabilized by anionic surfactant in the presence of 1:1, 1:2, 1:3, and 2:2 electrolytes: NaCl, Na2SO4, Na3Citrate, and MgSO4 are reported. The disjoining pressure predicted by the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory coincides with the experimental data in the case of a 1:1 electrolyte, but it is considerably greater than the measured pressure in all other cases. The theory is extended to account for the effects of ionic correlations and finite ionic radii. Original analytical expressions are derived for the local activity coefficient, electrostatic disjoining pressure, and asymptotic screening parameter. With the same parameter of counterion binding as for a 1:1 electrolyte, the curves predicted by the extended theory are in perfect agreement with the experimental data for 1:2 and 1:3 electrolytes. In comparison with the DLVO theory, the effect of ionic correlations leads to more effective screening of electrostatic interactions, and lower electric potential and counterion concentrations in the film’s midplane, resulting in lower disjoining pressure, as experimentally observed. The developed theory is applicable to both multivalent coions and multivalent counterions. Its application could remove some discrepancies between theory and experiment observed in studies with liquid films from electrolyte solutions.


Introduction
The Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [1,2] of surface forces in thin liquid films provides a quantitative description of the behavior of colloidal dispersions. Later, additional effects called non-DLVO forces have been taken into account [3]. In fact, one of the basic equations for the electrostatic component of disjoining pressure of a symmetric plane-parallel liquid film, Π el , was proposed earlier by Langmuir [4]: ppxq " kT rc 1 pxq`c 2 pxqs´k T 24π t e 2 ε 0 ε kT rz 2 1 c 1 pxq`z 2 2 c 2 pxqsu 3{2 (1.2) where ε 0 is the dielectric permittivity of vacuum; ε is the relative dielectric constant; e is the charge of positron; and z 1 and z 2 are the valences of the respective ions. The last term in Equation (1.2), which takes into account the effect of ion-ion correlations, gives an essential contribution to the value of p even for moderately high electrolyte concentrations, but this effect is not included in the DLVO theory [7]. Indeed, only the contribution of the first term in the right-hand side of Equation (1.2) has been taken into account in Equation (1.1). Nevertheless, the ionic correlations are expected to influence the electrostatic interactions in the thin liquid films and colloidal dispersions, especially in the presence of multivalent ions.
In the case of single electric double layer, Outhwaite and Bhuiyan [8] proposed an approach based on a modified Poisson-Boltzmann equation, which takes into account the contribution of ionic correlations to the double-layer potential and the coion and counterion profiles. In the case of two overlapping double layers (thin liquid film), Attard et al. [9] proposed an analytical approach, which extends the DLVO theory based on the Poisson-Boltzmann equation by including the effect of images and ionic correlations in the case of symmetric electrolytes. Ninham et al. [10,11] derived an analytical expression for the decay length of the asymptotic form of the double-layer force in mixed and asymmetric electrolytes, which is markedly different from the classical Debye length. The latter expression was verified experimentally by Pashley et al. [12] by using the method of colloidal-probe atomic force microscopy (CP-AFM). These authors concluded that the various proposed corrections to the electrostatic screening lengths in solutions of multivalent electrolytes are either invalid or beyond the limits of measurement of the used technique. It was found that surface force profiles measured by CP-AFM in solutions with monovalent counterions and multivalent coions were accurately described by the DLVO theory [13]. However, differences between the predictions of the Monte Carlo simulations and the DLVO theory have been reported in several studies [14][15][16][17] for solutions of asymmetric electrolytes.
In the case of 2:1 electrolytes, the divalent counterions adsorb strongly on the film surfaces and decrease the double-layer repulsion [18]. In the case of multivalent counterions, even charge reversal was observed [19][20][21]. The charge reversal is related to the effect of counterion adsorption, termed also counterion binding or condensation. Besteman et al. [22,23] showed that, in the presence of sufficiently high concentrations of multivalent salt, two originally oppositely charged surfaces could repel each other. Later, Monte Carlo simulations showed that this repulsion appears to be a consequence of ion-ion correlations in the film [24]; see also [25]. By theoretical analysis, strong evidence was found that ion-ion correlations can be the origin of charge inversion in the electric double layer near the interface between mercury and an aqueous MgSO 4 solution [26].
At present, the results about the importance of the effect of ionic correlations for the electrostatic double-layer interactions in thin liquid films appear to be controversial. In part, this is due to the fact that in many cases the correlation effect could be incorporated in an effective surface charge (or potential) [9], which can be determined as an adjustable parameter, thus obtaining an apparent agreement with the DLVO theory. The binding energies of counterions of various valences, which are frequently unknown, lead to the appearance of additional adjustable parameters. To overcome this problem, in the present study we measured the disjoining pressure of free thin liquid films stabilized by the anionic surfactant sodium dodecyl sulfate (SDS). Disjoining pressure isotherms were measured in the presence of 1:1; 1:2, and 1:3 electrolytes, viz. NaCl, Na 2 SO 4 , and Na 3 Citrate. In other words, the counterion, Na + , is the same, whereas the valences of the coions are different. It is very important that the binding energy of the Na + ion to the sulfate headgroup of the SDS molecule, characterized by the Stern constant K St,2 , is known from surface-tension studies [27][28][29][30] and has not been used as an adjustable parameter in the present study. In such a case, the effect of the ionic correlations on the surface force (if any) should be detected as a deviation from the predictions of the DLVO theory in regime of charge regulation, characterized by the same counterion binding constant, K St,2 , for the 1:1, 1:2, and 1:3 electrolytes. Disjoining pressure isotherm has been measured also in the presence of MgSO 4 to compare the ionic-correlation effects in the presence of 1:1 and 2:2 electrolytes. Section 2 describes the used materials and methods. In Section 3, the DLVO theory is extended to take into account the effect of ionic correlations, henceforth referred to as DLVO-IC for brevity. In Section 4, the experimental data are compared with the extended theory; the effect of ionic correlations is identified and interpreted in terms of ion concentrations and electrostatic potential inside the film. In Section 5, an analytical expression is derived for the asymptotic decay length of double-layer interactions affected by the ionic correlations. In Appendix A, a general expression for the local activity coefficient in electric double layers is derived, which is used in the DLVO-IC theory. We should mention in advance that a considerable effect of ionic correlations was detected in the cases of 1:2, 1:3 and 2:2 electrolytes, and that DLVO-IC accurately describes the experimental disjoining pressure isotherms.

Materials, Methods and Experimental Results
The chemicals used were: sodium dodecyl sulfate (SDS) from Across Organics (Geel, Belgium); sodium chloride and sodium sulfate (NaCl and Na 2 SO 4 ) from Merck (Kenilworth, NJ, USA); and sodium citrate and magnesium sulfate (Na 3 Citrate and MgSO 4 ) from Sigma-Aldrich (Taufkirchen, Germany). All solutions were prepared with deionized water of specific resistivity 18.2 MΩ¨cm (Elix purification system, Millipore, Molsheim, France). The natural pH of the Na 3 Citrate solutions was 7.0˘0.1, which guarantees that the citrate anion is present in its trivalent form. In all experiments, the concentration of SDS was 1 mM, whereas the ionic strength due to the added electrolyte was 30 mM. The last value corresponds to 30 mM NaCl, 10 mM Na 2 SO 4 , and 5 mM Na 3 Citrate. All experiments were carried out at a room temperature of 25˘1˝C.
The dependence of disjoining pressure, Π, versus the film thickness, h w , was measured by means of a Mysels-Jones (MJ) experimental cell [31], which was further modified and termed thin film pressure balance [32,33]; see Figure 1. The film formed in the center of a cylindrical hole drilled in a porous-glass plate. The two film surfaces are forced against each other at pressures of up to 9000 Pa in our experiments. The film thickness was determined by an interferometric method from the intensity of the monochromatic light (λ = 546 nm) reflected from the film [32,34,35]. The value of film thickness was calculated assuming that the film is a uniform plane-parallel layer with refractive index of water, n w = 1.33, and for this reason the determined optical thickness was called equivalent water thickness, h w . The reflected light was supplied to a photomultiplier and computer, and h w was recorded continuously during the experiment. The applied pressure was measured electronically by a pressure transducer. A detailed description of the used version of the MJ cell can be found elsewhere [36].
Materials 2016, 9,145 3 of 23 has been measured also in the presence of MgSO4 to compare the ionic-correlation effects in the presence of 1:1 and 2:2 electrolytes. Section 2 describes the used materials and methods. In Section 3, the DLVO theory is extended to take into account the effect of ionic correlations, henceforth referred to as DLVO-IC for brevity. In Section 4, the experimental data are compared with the extended theory; the effect of ionic correlations is identified and interpreted in terms of ion concentrations and electrostatic potential inside the film. In Section 5, an analytical expression is derived for the asymptotic decay length of double-layer interactions affected by the ionic correlations. In Appendix A, a general expression for the local activity coefficient in electric double layers is derived, which is used in the DLVO-IC theory. We should mention in advance that a considerable effect of ionic correlations was detected in the cases of 1:2, 1:3 and 2:2 electrolytes, and that DLVO-IC accurately describes the experimental disjoining pressure isotherms.

Materials, Methods and Experimental Results
The chemicals used were: sodium dodecyl sulfate (SDS) from Across Organics (Geel, Belgium); sodium chloride and sodium sulfate (NaCl and Na2SO4) from Merck (Kenilworth, NJ, USA); and sodium citrate and magnesium sulfate (Na3Citrate and MgSO4) from Sigma-Aldrich (Taufkirchen, Germany). All solutions were prepared with deionized water of specific resistivity 18.2 MΩ•cm (Elix purification system, Millipore, Molsheim, France). The natural pH of the Na3Citrate solutions was 7.0 ± 0.1, which guarantees that the citrate anion is present in its trivalent form. In all experiments, the concentration of SDS was 1 mM, whereas the ionic strength due to the added electrolyte was 30 mM. The last value corresponds to 30 mM NaCl, 10 mM Na2SO4, and 5 mM Na3Citrate. All experiments were carried out at a room temperature of 25 ± 1 °C.
The dependence of disjoining pressure, Π, versus the film thickness, hw, was measured by means of a Mysels-Jones (MJ) experimental cell [31], which was further modified and termed thin film pressure balance [32,33]; see Figure 1. The film formed in the center of a cylindrical hole drilled in a porous-glass plate. The two film surfaces are forced against each other at pressures of up to 9000 Pa in our experiments. The film thickness was determined by an interferometric method from the intensity of the monochromatic light (λ = 546 nm) reflected from the film [32,34,35]. The value of film thickness was calculated assuming that the film is a uniform plane-parallel layer with refractive index of water, nw = 1.33, and for this reason the determined optical thickness was called equivalent water thickness, hw. The reflected light was supplied to a photomultiplier and computer, and hw was recorded continuously during the experiment. The applied pressure was measured electronically by a pressure transducer. A detailed description of the used version of the MJ cell can be found elsewhere [36].  [31], where a foam film is formed in the center of a hole in a porous plate, which is connected to a capillary filled with the surfactant solution and to a pressure transducer.
During the experiment, the applied pressure was increased in small steps, of 40 Pa for Π < 3000 Pa, and of 100 Pa for Π > 3000 Pa. After each pressure increase, the computer's indications for the film thickness were observed (from 5 to 30 min) until an equilibrium thickness was established, which was recorded together with the respective value of pressure. After that, the  [31], where a foam film is formed in the center of a hole in a porous plate, which is connected to a capillary filled with the surfactant solution and to a pressure transducer.
During the experiment, the applied pressure was increased in small steps, of 40 Pa for Π < 3000 Pa, and of 100 Pa for Π > 3000 Pa. After each pressure increase, the computer's indications for the film thickness were observed (from 5 to 30 min) until an equilibrium thickness was established, which was recorded together with the respective value of pressure. After that, the pressure was increased again to obtain the next experimental point. For each of the three electrolytes, NaCl, Na 2 SO 4 and Na 3 Citrate, the experimental Π(h w ) dependence was measured in four independent runs. The experimental points shown in Figure 2 are average values from these four experiments. The error of the mean value of h w is˘0.4 nm for Π < 1 kPa;˘0.3 nm for 1 ď Π ď 4 (kPa); and˘0.15 nm for Π > 4 kPa.
Materials 2016, 9,145 4 of 23 pressure was increased again to obtain the next experimental point. For each of the three electrolytes, NaCl, Na2SO4 and Na3Citrate, the experimental Π(hw) dependence was measured in four independent runs. The experimental points shown in Figure 2 are average values from these four experiments. The error of the mean value of hw is ±0.4 nm for Π < 1 kPa; ±0.3 nm for 1 ≤ Π ≤ 4 (kPa); and ±0.15 nm for Π > 4 kPa.  was the same; I = 31 mM: 1 mM from the SDS and 30 mM from the added electrolyte; zi are the valences of the ions; and ci∞ are their bulk concentrations. Despite the fixed I, the Π(hw) isotherms measured in the presence of 1:1, 1:2, and 1:3 electrolytes differ essentially from each other: Π increases with the rise of coion valence. The predictions of the conventional DLVO theory coincide with the experimental Π(hw) curve only in the case of NaCl (computational details can be found below). In contrast, for the 1:2 and 1:3 electrolytes, the Π calculated by DLVO theory is considerably greater than the experimental Π. As demonstrated below, this discrepancy can be overcome if the DLVO theory is extended to account for the effect of ionic correlations. In the presence of asymmetric electrolytes, the ionic correlation effect leads to a stronger screening of the repulsive electrostatic interactions and to lower values of disjoining pressure, as experimentally observed ( Figure 2). The effect is stronger for a 1:3 electrolyte. For a 1:2 electrolyte the effect is weaker, but noticeable, similarly to the difference between the Monte Carlo data and the DLVO curve in [14]. The greater reduction of the electrostatic repulsion with the rise of coion valence found in the present study is in qualitative agreement with the inverse Schulze-Hardy rule established in [37].  As demonstrated below, this discrepancy can be overcome if the DLVO theory is extended to account for the effect of ionic correlations. In the presence of asymmetric electrolytes, the ionic correlation effect leads to a stronger screening of the repulsive electrostatic interactions and to lower values of disjoining pressure, as experimentally observed ( Figure 2). The effect is stronger for a 1:3 electrolyte. For a 1:2 electrolyte the effect is weaker, but noticeable, similarly to the difference between the Monte Carlo data and the DLVO curve in [14]. The greater reduction of the electrostatic repulsion with the rise of coion valence found in the present study is in qualitative agreement with the inverse Schulze-Hardy rule established in [37].

Expression for Disjoining Pressure Extended to Account for Ionic Correlations
Let us consider a free liquid film formed from a solution of surfactant and electrolyte; see Figure 3a. The film surfaces are electrically charged and the electric potential in the film, ϕ, obeys the Poisson equation: here, x is Cartesian coordinate normal to the film surfaces ( Figure 3b); ρ b is the local charge density in the electric double layer (EDL); c i = c i (x) are the local concentrations of the ions; and n is the number of ionic components. Under static conditions, the Navier-Stokes equation reduces to dp dx where p is the local osmotic pressure; and ρ b has been substituted from the Poisson Equation (3.1). Equation (3.2), which expresses the local force balance in the fluid inside the film, implies: Materials 2016, 9, 145 5 of 23

Expression for Disjoining Pressure Extended to Account for Ionic Correlations
Let us consider a free liquid film formed from a solution of surfactant and electrolyte; see Figure 3a. The film surfaces are electrically charged and the electric potential in the film, ϕ, obeys the Poisson equation:  Because the electrochemical potential, µ i , of the i-th ionic species is constant throughout the EDL, we obtain: where γ i is the activity coefficient. Upon summation in Equation (3.4) over all ionic species, its last term becomes ρ b dϕ, which in view of Equation (3.3) yields: The effect of ionic correlations on the local pressure, p, is incorporated in the activity coefficient, In Appendix A, using the Debye-Hückel theory extended to a nonuniform electrolyte solution (like an electric double layer), we have derived the following expression for γ i : here, r i is the radius of the respective ion; κ is the local Debye screening parameter; and L B is the Bjerrum length: The integration of Equation (3.5), along with Equations (3.6) and (3.7), yields: By differentiation, one could verify that the expression for p in Equation (3.8) satisfies Equation (3.5). In the limit κr i Ñ 0, Equation (3.8) reduces to Equation (1.2), generalized for the case of n ionic components. Moreover, the form of the last term in Equation (3.8) indicates that the correlation effect, related to the screening parameter, κ, is coupled with the effect of finite ionic radius, r i . The value of p in the bulk of solution can be calculated from Equations (3.7) and (3.8) with c i = c i8 : where c i8 (i = 1, . . . , n) are the respective bulk concentrations. The constancy of the electrochemical potential in the electric double layer, leads to Boltzmann relations for the activities of the ions: here, i = 1, . . . , n; γ i is the activity coefficient in the EDL given by Equation (3.6); γ i8 is the activity coefficient in the bulk of the solution, which is to be calculated from Equation (3.6) with c i = c i8 ; ϕ and Φ are the dimensional and dimensionless electric potential in the EDL; for convenience, Φ is defined as a nonnegative quantity. where Φ m and p(Φ m ) are the values of Φ and p(Φ) in the film's midplane ( Figure 3b). Equation (3.11), along with Equations (3.6)-(3.10), represents an expression for the electrostatic disjoining pressure extended to account for the effect of ionic correlations. It has to be noted that Equation (3.8) should also contain an additive constant of integration (background pressure), which is canceled in Equation (3.11) and does not affect Π el .

Calculation of the Electrostatic Disjoining Pressure with Ionic Correlations
By integration of Equation (3.2), one derives: A subsequent integration from the film midplane (x = 0) to the film surface (x = h/2) yields: where h is the film thickness; and Φ s is the dimensionless potential at the film surface. The boundary condition for the electric field at the film surface reads: here, ρ s is the surface charge density and Γ i is the adsorption of the respective ionic species, i.e., the number of adsorbed ions per unit surface area. (We have Γ i = 0 if component i is not present at the film surface.) The elimination of the derivative of electric potential between Equations (3.12) and (3.14) yields: As an example, our experimental system (see Section 2) contains three kinds of ions: component 1-surfactant ion (dodecyl sulfate anion); component 2-counterion (Na + ), and component 3-coions (Cl´, SO 4 2´o r Citrate 3´) . In this case, Γ 1 denotes the surfactant adsorption, Γ 2 denotes the number of counterions bound per unit surface area, and Γ 3 = 0, i.e., the negatively charged coions are not expected to bind to the like charged interface. For our computations, it is convenient to present the surfactant adsorption isotherm [27] in the form: here, K 1 is the surfactant adsorption constant; K St,2 is the Stern constant of Na + binding to the surfactant headgroups in the adsorption layer; Γ 8´1 is the excluded area per surfactant molecule in the adsorption layer; and β is a parameter that characterizes the attraction between the surfactant tails in the adsorption layer. All these parameters have been determined from the fits of surface tension isotherms at various salt concentrations [28]. Equation (3.16) corresponds to the van der Waals adsorption model, which turns out to be adequate for describing surfactant adsorption at a fluid interface [38]. In addition, the counterion binding can be described by the Stern isotherm [27]: The use of Equation (3.17) means that we are working in the regime of charge regulation, i.e., in the regime of constant electrochemical potential, which is the real experimental situation. Indeed, Equation (3.17) is derived by setting the electrochemical potential of the counterions in the bulk equal to that of the adsorbed counterions (bound in the Stern layer), where the Langmuir adsorption isotherm has been used. Note that Equations (3.16) and (3.17) are thermodynamically compatible, i.e., they satisfy the Euler condition for exact differential [27]: For a given value of Φ s , from Equations (3.6) and (3.10) we determine the respective values of the concentrations and activity coefficients, c i = c is and γ i = γ is , by iterations. The procedure is analogous to that in point 1 above.
The integral in Equation (3.13) is solved numerically. For this goal, at each given Φ from Equations (3.6) and (3.10) we determine the respective values of the concentrations and activity coefficients, c i and γ i , by iterations. The procedure is analogous to that in point 1 above. Next, from Equations (3.7) and (3.8) we calculate p(Φ s ). As a result of integration, we determine h(Φ m ). 7.
The value of Φ m is varied to obtain the dependences Π el = Π el (Φ m ) and h = h(Φ m ).

Conventional DLVO Theory of Electrostatic Disjoining Pressure
For the conventional DLVO theory, one can find the dependences Π el = Π el (Φ m ) and h = h(Φ m ) by using the same procedure (Section 3.2), but instead of Equations (3.8) and (3.9) the following simpler expressions have to be used to calculate the pressures: In addition, the activity coefficients have to be expressed from the Debye-Hückel limiting law [5]: The only equations, where the activity coefficients γ 18 and γ 28 appear, are the adsorption isotherms, Equations (3.16) and (3.17). Upon substituting γ i = γ i8 in Equation (3.10), the activity coefficients disappear from this equation. Then, Equation (3.10) represents an explicit expression for the concentrations c i and iterations are not necessary at steps 1, 3, and 6 of the procedure in Section 3.2.

Calculation of the Total Disjoining Pressure
For the considered free foam films, the total disjoining pressure, Π, is the sum of Π el and the van der Waals component of disjoining pressure, Π vw : As before, h is the thickness of the aqueous core; and d accounts for the thickness of the two surfactant adsorption layers. In our calculations, we are using the Hamaker parameter A H for water and h w = h + d is the experimental "water thickness" of the film. We recall that the experiment yields h w insofar as the intensity of light reflected from the film is determined, assuming that the film has a uniform refractive index equal to that of water; see Section 2.
The Hamaker parameter, A H , is calculated from the formula [39]: where the retardation effect is taken into account; h P = 6.63ˆ10´3 4 J¨s is the Planck constant; ν e « 3.0ˆ10 15 Hz is the main electronic absorption frequency; n i and n j are the refractive indexes, respectively, of the outer and inner phase; ε i and ε j are the respective relative dielectric constants; the dimensionless thickness r h is defined by the expression r h " 2 πν e h w n j pn 2 i`n 2 j q 1{2 {c 0 ; c 0 = 3.0ˆ10 8 m/s is the speed of light in vacuum; and ξ is an integration variable.

Numerical Results and Discussion
4.1. Results for 1:1, 1:2, and 1:3 Electrolytes For our system, n i = 1 and ε i = 1 (air); n j = 1.33 and ε j = 78.25 (water). A H calculated from Equation (3.22) decreases from 3.0ˆ10´2 0 J at h w = 5 nm to 1.7ˆ10´2 0 J at h w = 20 nm because of the electromagnetic retardation effect. In the experimental range of film thicknesses, Π wv is relatively small, but it is not negligible. For example, for the leftmost experimental point of the curve for Na 2 SO 4 in Figure 2 we have h w = 13.2 nm, Π = 7.10 kPa, and Π el = 7.59 kPa, whereas |Π wv | = 0.49 kPa, i.e. 6.5% of Π el .
For the hydrated ionic radii, r i , we used the values 0.36 nm for Na + and 0. 33  . Among all ions, the sodium counterions are present at a higher concentration in the EDL and, consequently, Π el is the most sensitive to the radius of the Na + ions. Conversely, the concentrations of the coins (repelled by the like charged film surfaces) are lower in the film, and Π el is much less sensitive to the values of their radii.
The solid lines in Figure 2 are drawn by using the DLVO theory with ionic correlations (DLVO-IC), as described in Section 3.2. A single adjustable parameter was used, the effective water thickness of the surfactant adsorption layers, d in Equation (3.21); see also Figure 3. The three experimental isotherms in Figure 2 were simultaneously fitted with the DLVO-IC theory and d = 1.7 nm was obtained; the agreement between theory and experiment is excellent.
The dashed lines in Figure 2 represent the predictions of the conventional DLVO theory (see Section 3.3) with the same d = 1.7 nm in Equation (3.21). All constant geometrical and physical parameters used in the calculations with DLVO theory are the same as with DLVO-IC. In the case of a 1:1 electrolyte, NaCl, the prediction of DLVO theory coincides with that of DLVO-IC and with the experimental curve. However, in the cases of Na 2 SO 4 and Na 3 Citrate the conventional DLVO theory predicts considerably greater values of Π in comparison with the experimental values and the predictions of DLVO-IC. This result confirms that the ionic correlation effects are really significant in the cases of divalent and trivalent coions.
It should be also noted that agreement between the conventional DLVO theory and experiment in the case of NaCl can be achieved only if we keep the activity coefficients γ 18  If an approximation for an ideal solution (γ i8 = 1) is used in the adsorption isotherms, the conventional DLVO theory predicts a Π(h w ) isotherm, which is markedly above the experimental isotherm with NaCl in Figure 2.
Because all parameters in theory are known, we can calculate various properties of the liquid film in the framework of the DLVO-IC theory. Thus, in Figure 4a the dependences of the magnitudes of surface electric potential, |ϕ s |, and surface charge density, |ρ s /e| = Γ 1´Γ2 , are plotted vs. the bulk concentration of counterions, c 28 , in the asymptotic case of large film thickness, h Ñ 8, i.e., for a single interface. Three theoretical curves have been calculated, corresponding to the cases of added NaCl, Na 2 SO 4 , and Na 3 Citrate, i.e., electrolytes with the same counterion, Na + , but with coions of different valence. It is remarkable that the three theoretical curves for |ϕ s | vs. c 28 , calculated for monovalent, divalent, and trivalent coions, practically coincide. This result, which is numerically very close to the prediction of DLVO, illustrates the known paradigm that the surface potential is determined by the counterion, and is insensitive to the valence of coions, which, in turn, has been used to predict the critical coagulation concentration (CCC) and to explain the Schulze-Hardy rule by the DLVO theory [1][2][3]. The decrease of |ϕ s | with the rise of c 28 can be explained with the enhanced Debye screening. The increase of |ρ s /e| = Γ 1´Γ2 with the rise of c 28 could be attributed to the rise of surfactant (SDS) adsorption at the air/water interface, Γ 1 , upon the addition of electrolyte, owing to the suppression of the electrostatic repulsion between the surfactant ions at the interface. As seen in Figure 4a, a relatively weak effect of the coion valence on ρ s is present at the higher values of c 28 . It should be also noted that agreement between the conventional DLVO theory and experiment in the case of NaCl can be achieved only if we keep the activity coefficients γ1∞ and γ2∞ in the adsorption isotherms, Equations (3.16) and (3.17). For I = 31 mM, Equation (3.20) yields γ1∞ = γ2∞ = 0.845. If an approximation for an ideal solution (γi∞ = 1) is used in the adsorption isotherms, the conventional DLVO theory predicts a Π(hw) isotherm, which is markedly above the experimental isotherm with NaCl in Figure 2.
Because all parameters in theory are known, we can calculate various properties of the liquid film in the framework of the DLVO-IC theory. Thus, in Figure 4a the dependences of the magnitudes of surface electric potential, |ϕs|, and surface charge density, |ρs/e| = Γ1 − Γ2, are plotted vs. the bulk concentration of counterions, c2∞, in the asymptotic case of large film thickness, h → ∞, i.e., for a single interface. Three theoretical curves have been calculated, corresponding to the cases of added NaCl, Na2SO4, and Na3Citrate, i.e., electrolytes with the same counterion, Na + , but with coions of different valence. It is remarkable that the three theoretical curves for |ϕs| vs. c2∞, calculated for monovalent, divalent, and trivalent coions, practically coincide. This result, which is numerically very close to the prediction of DLVO, illustrates the known paradigm that the surface potential is determined by the counterion, and is insensitive to the valence of coions, which, in turn, has been used to predict the critical coagulation concentration (CCC) and to explain the Schulze-Hardy rule by the DLVO theory [1][2][3]. The decrease of |ϕs| with the rise of c2∞ can be explained with the enhanced Debye screening. The increase of |ρs/e| = Γ1 − Γ2 with the rise of c2∞ could be attributed to the rise of surfactant (SDS) adsorption at the air/water interface, Γ1, upon the addition of electrolyte, owing to the suppression of the electrostatic repulsion between the surfactant ions at the interface. As seen in Figure  4a, a relatively weak effect of the coion valence on ρs is present at the higher values of c2∞.     Figure 4b shows the dependencies of surface potential and charge on the solution's ionic strength, I, predicted by the DLVO-IC theory for the same 1:1, 1:2, and 1:3 electrolytes. In this case, the values of ϕ s and ρ s are different for the coions of different valence at the same ionic strength, I. This difference is almost completely due to the fact that at the same I, the counterion concentration, c 28 , is different in the 1:1, 1:2, and 1:3 electrolytes. For example, at 1 mM SDS and ionic strength I = 31 mM, we have c 28 = 31, 21, and 16 mM in the solutions of NaCl, Na 2 SO 4 , and Na 3 Citrate, respectively. In view of Figure 4a different values of c 28 correspond to different ϕ s values.
In Figure 5, we compare the predictions of the DLVO and DLVO-IC theories with respect to the dependences of surface charge density and surface potential, ρ s and ϕ s , on the film thickness, h. In the case of DLVO theory, Figure 5a,b, the results indicate the following: (i) ρ s and ϕ s are weakly sensitive to the valence of the coion. (ii) The dependencies of both surface charge and potential on h are rather weak: ρ s varies with about 5%, whereas ϕ s varies with only 1.4%. Hence, in this case the regime of charge regulation is closer to the regime of fixed surface potential. (iii) The magnitude of surface charge density, |ρ s |, decreases with the decrease of h (Figure 5a), which can be explained by the enhanced counterion binding in the thinner films; (iv) Conversely, the magnitude of the surface potential, |ϕ s |, increases with the decrease of h (Figure 5b), which is counterintuitive in view of the decrease in surface charge density |ρ s |. This difference is almost completely due to the fact that at the same I, the counterion concentration, c2∞, is different in the 1:1, 1:2, and 1:3 electrolytes. For example, at 1 mM SDS and ionic strength I = 31 mM, we have c2∞ = 31, 21, and 16 mM in the solutions of NaCl, Na2SO4, and Na3Citrate, respectively. In view of Figure 4a different values of c2∞ correspond to different ϕs values.
In Figure 5, we compare the predictions of the DLVO and DLVO-IC theories with respect to the dependences of surface charge density and surface potential, ρs and ϕs, on the film thickness, h. In the case of DLVO theory, Figure 5a,b, the results indicate the following: (i) ρs and ϕs are weakly sensitive to the valence of the coion. (ii) The dependencies of both surface charge and potential on h are rather weak: ρs varies with about 5%, whereas ϕs varies with only 1.4%. Hence, in this case the regime of charge regulation is closer to the regime of fixed surface potential. (iii) The magnitude of surface charge density, |ρs|, decreases with the decrease of h (Figure 5a), which can be explained by the enhanced counterion binding in the thinner films; (iv) Conversely, the magnitude of the surface potential, |ϕs|, increases with the decrease of h (Figure 5b), which is counterintuitive in view of the decrease in surface charge density |ρs|. As seen in Figure 5c,d, the predictions of the DLVO-IC theory essentially differ from those of the conventional DLVO theory: (i) The surface charge and potential are sensitive to the kind of the coion; (ii) the dependencies of both surface charge and potential on h are significant: |ρs| increases about three times, and |ϕs| increases from 90 mV up to 120-135 mV, depending on the coion; (iii) both surface charge, |ρs|, and surface potential, |ϕs|, decreases with the decrease of h, which can be explained with enhanced counterion binding in the thinner films. With respect to the behavior of |ϕs| vs. h, the DLVO and DLVO-IC theories predict opposite tendencies; compare Figure 5b,d. To illustrate the differences between DLVO and DLVO-IC with respect to the ionic concentrations inside the film, in Figure 6 we have plotted the distributions of surfactant ions c 1 (x), counterions c 2 (x), and coions c 3 (x), vs. the distance to the film surface, h´x, predicted by the two theories for a film of thickness h = 14 nm formed from a solution of 1 mM SDS and 5 mM Na 3 Citrate. More significant relative differences between the two theories are seen for the surfactant and citrate coions, c 1 (x) and c 3 (x). Note, however, that the concentrations of coions in the film are much lower than those of the counterions, c 1 , c 3 << c 2 , so that Π el is dominated by the effect of counterion concentration, c 2 .
To illustrate the differences between DLVO and DLVO-IC with respect to the ionic concentrations inside the film, in Figure 6 we have plotted the distributions of surfactant ions c1(x), counterions c2(x), and coions c3(x), vs. the distance to the film surface, h−x, predicted by the two theories for a film of thickness h = 14 nm formed from a solution of 1 mM SDS and 5 mM Na3Citrate. More significant relative differences between the two theories are seen for the surfactant and citrate coions, c1(x) and c3(x). Note, however, that the concentrations of coions in the film are much lower than those of the counterions, c1, c3 << c2, so that Πel is dominated by the effect of counterion concentration, c2. Our computations show that the difference between the Π(h) curves predicted by the DLVO and DLVO-IC theories for 1:2 and 1:3 electrolytes (Figure 2) is due mostly to the difference between the counterion concentrations, c2m, in the film midplane. To illustrate this fact, the two upper curves for c2 in Figure 6a (which seem to be almost coinciding) are plotted in magnified scale in Figure 6b. In this figure, c2m is the value of c2 at h − x = 7 nm. The difference between the c2m values predicted by the two models is Δc2m = 20.5 − 19.2 = 1.3 mM. In view of the first term in the right-hand side of Equation  Our computations show that the difference between the Π(h) curves predicted by the DLVO and DLVO-IC theories for 1:2 and 1:3 electrolytes (Figure 2) is due mostly to the difference between the counterion concentrations, c 2m , in the film midplane. To illustrate this fact, the two upper curves for c 2 in Figure 6a (which seem to be almost coinciding) are plotted in magnified scale in Figure 6b. In this figure, c 2m is the value of c 2 at h´x = 7 nm. The difference between the c 2m values predicted by the two models is ∆c 2m = 20.5´19.2 = 1.3 mM. In view of the first term in the right-hand side of Equation In other words, ∆p gives the main contribution in ∆Π. The relatively small difference between ∆p = 3225 Pa (estimated from the first term in Equation (3.8)) and the experimental ∆Π = 3154 Pa, is due to the last term in Equation (3.8). The values of p 8 predicted by DLVO and DLVO-IC are very close; see Equations (3.9) and (3.11). In such a case, ∆Π is practically equal to the difference, ∆p m , between the midplane pressures, p(Φ m ), which is plotted vs. h in Figure 6c. As already mentioned, ∆p m is dominated by the contribution of ∆c 2m that, in turns, is related to the difference between the midplane potentials ϕ m , which is illustrated in Figure 6d.
In Figure 7a, we have plotted the experimental data and DLVO-IC theoretical curves from Figure 2 as Π el = Π´Π vw vs. h in semi-logarithmic scale. Π vw was calculated as explained in Section 3.4 with the same d = 1.7 nm, as above. We recall that the three experimental curves correspond to the same ionic strength, I = 31 mM, but to different coion valences, z 3 =´1,´2 and´3. As seen in Figure 7a, at a large h value the experimental curves follow a linear dependence, which corresponds to an exponential asymptotic law: Π el 9 expp´κ a hq for large h (4.1) where κ a is an asymptotic screening parameter. The values of κ a determined from the computed curves at h Ñ 8 are given in Table 1 In other words, Δp gives the main contribution in ΔΠ. The relatively small difference between Δp = 3225 Pa (estimated from the first term in Equation (3.8)) and the experimental ΔΠ = 3154 Pa, is due to the last term in Equation (3.8). The values of p∞ predicted by DLVO and DLVO-IC are very close; see Equations (3.9) and (3.11). In such a case, ΔΠ is practically equal to the difference, Δpm, between the midplane pressures, p(Φm), which is plotted vs. h in Figure 6c. As already mentioned, Δpm is dominated by the contribution of Δc2m that, in turns, is related to the difference between the midplane potentials ϕm, which is illustrated in Figure 6d.
In Figure 7a, we have plotted the experimental data and DLVO-IC theoretical curves from Figure 2 as Πel = Π − Πvw vs. h in semi-logarithmic scale. Πvw was calculated as explained in Section 3.4 with the same d = 1.7 nm, as above. We recall that the three experimental curves correspond to the same ionic strength, I = 31 mM, but to different coion valences, z3 = −1, −2 and −3. As seen in Figure 7a, at a large h value the experimental curves follow a linear dependence, which corresponds to an exponential asymptotic law: where κa is an asymptotic screening parameter. The values of κa determined from the computed curves at h → ∞ are given in Table 1  It is important to note that the experimental points obtained by the MJ cell do not belong to the asymptotic region of large h, where the slope of the curves in Figure 7a is characterized by κa in accordance with Equation (4.1). In fact, the experimental points correspond to smaller h values, where the curves for 1:1, 1:2, and 1:3 electrolytes are almost parallel (Figure 7a). To illustrate that, the points belonging to an experimental curve in Figure 7a (except the last one or two points with   Figure 7a and the theoretical asymptotic screening parameter, κ a , calculated from the asymptotic slope of the theoretical curves or from Equation (5.12). It is important to note that the experimental points obtained by the MJ cell do not belong to the asymptotic region of large h, where the slope of the curves in Figure 7a is characterized by κ a in accordance with Equation (4.1). In fact, the experimental points correspond to smaller h values, where the curves for 1:1, 1:2, and 1:3 electrolytes are almost parallel (Figure 7a). To illustrate that, the points belonging to an experimental curve in Figure 7a (except the last one or two points with smaller Π, which are measured with lower accuracy) were fitted with linear regression. The fits are good (regression coefficients between 0.9964 and 0.9994); the slope of each regression yields an apparent screening parameter, κ exp , which is also given in Table 1. The values of κ exp are close to κ 8 . Hence, the slope of the semi-logarithmic plot of the experimental data cannot serve as an indicator for the significance of the effect of ionic correlations on the electrostatic surface force. The close values of κ exp and κ 8 could be the reason for some authors to conclude that the force profiles in the presence of multivalent ions can be accurately described by the conventional DLVO (Poisson-Boltzmann) theory [12,13]. As demonstrated in Figure 2, the main effect of ionic correlations is the vertical shift of the experimental curves toward smaller Π values as compared with the DLVO prediction.

System
As demonstrated in Section 5, the asymptotic region characterized by slope κ a in semi-logarithmic scale, see Equation (4.1), corresponds to the region with small midplane potential, Φ m << 1. From this viewpoint, our experimental data (Figure 7a) correspond to a range of h, where the condition Φ m << 1 is not satisfied (Figure 6d), and for this reason the local slope, κ exp , is different from κ a ; see the data for Na 3 Citrate in Table 1. Figure 7b compares the predictions of the DLVO and DLVO-IC for the same three systems, but in wider ranges of variation of h and Π el . The curves in Figure 7b indicate that the effect of ionic correlations affects Π el at both short and long distances between the film surfaces. In the special case of a symmetric 1:1 electrolyte (NaCl), the differences between the predictions of the two theories are essential only at the smallest values of h, but for the asymmetric electrolytes these differences are essential for all h values. It is also remarkable that at the smallest h the effect of the different coion valences vanishes and the respective curves calculated by DLVO-IC almost coincide.

Results for a 2:2 Electrolyte
In the case of a symmetric 1:1 electrolyte (I = 31 mM), no differences between DLVO and DLVO-IC have been observed; see Figure 2. It is interesting to verify whether there are differences between DLVO and DLVO-IC in the case of a symmetric 2:2 electrolyte. For this goal, we measured experimentally Π vs. h w for foam films by using the MJ cell at the same ionic strength, I = 31 mM, which is due to 7.5 mM MgSO 4 , and 1 mM SDS; see Figure 8a. Because the radius of the hydrated Mg 2+ ion and the Stern constant characterizing the binding of the Mg 2+ ion to the surfactant headgroup are known (see below), the theoretical curves for a 2:2 electrolyte in Figure 8a predicted by the DLVO and DLVO-IC theories have been drawn without using any adjustable parameters. Details are following.
smaller Π, which are measured with lower accuracy) were fitted with linear regression. The fits are good (regression coefficients between 0.9964 and 0.9994); the slope of each regression yields an apparent screening parameter, κexp, which is also given in Table 1. The values of κexp are close to κ∞. Hence, the slope of the semi-logarithmic plot of the experimental data cannot serve as an indicator for the significance of the effect of ionic correlations on the electrostatic surface force. The close values of κexp and κ∞ could be the reason for some authors to conclude that the force profiles in the presence of multivalent ions can be accurately described by the conventional DLVO (Poisson-Boltzmann) theory [12,13]. As demonstrated in Figure 2, the main effect of ionic correlations is the vertical shift of the experimental curves toward smaller Π values as compared with the DLVO prediction.
As demonstrated in Section 5, the asymptotic region characterized by slope κa in semi-logarithmic scale, see Equation (4.1), corresponds to the region with small midplane potential, Φm << 1. From this viewpoint, our experimental data (Figure 7a) correspond to a range of h, where the condition Φm << 1 is not satisfied (Figure 6d), and for this reason the local slope, κexp, is different from κa; see the data for Na3Citrate in Table 1. Figure 7b compares the predictions of the DLVO and DLVO-IC for the same three systems, but in wider ranges of variation of h and Πel. The curves in Figure 7b indicate that the effect of ionic correlations affects Πel at both short and long distances between the film surfaces. In the special case of a symmetric 1:1 electrolyte (NaCl), the differences between the predictions of the two theories are essential only at the smallest values of h, but for the asymmetric electrolytes these differences are essential for all h values. It is also remarkable that at the smallest h the effect of the different coion valences vanishes and the respective curves calculated by DLVO-IC almost coincide.

Results for a 2:2 Electrolyte
In the case of a symmetric 1:1 electrolyte (I = 31 mM), no differences between DLVO and DLVO-IC have been observed; see Figure 2. It is interesting to verify whether there are differences between DLVO and DLVO-IC in the case of a symmetric 2:2 electrolyte. For this goal, we measured experimentally Π vs. hw for foam films by using the MJ cell at the same ionic strength, I = 31 mM, which is due to 7.5 mM MgSO4, and 1 mM SDS; see Figure 8a. Because the radius of the hydrated Mg 2+ ion and the Stern constant characterizing the binding of the Mg 2+ ion to the surfactant headgroup are known (see below), the theoretical curves for a 2:2 electrolyte in Figure 8a predicted by the DLVO and DLVO-IC theories have been drawn without using any adjustable parameters. Details are following.
(a)  As seen from the experimental values of hw (Figure 8a), in the presence of MgSO4 the equilibrium foam films are considerably thinner than in the presence of NaCl at the same ionic strength. Moreover, the maximal values of Π that can be achieved experimentally are not so high because a transition from primary (common black) film to secondary (Newton black) film takes place; see e.g., Figure 7.17 in reference [42].
Theoretically, in this case we are dealing with four ionic components, which are numbered as follows: 1-surfactant ion (dodecyl sulfate); 2-Na + ion (from SDS); 3-SO4 2− anion and 4-Mg 2+ cation (from MgSO4). All equations in Section 3 remain the same, except the surfactant adsorption isotherm, Equation (3.16) and the counterion binding isotherm, Equation (3.17). Equation (3.16) has to be replaced with the following equation [27]: where c4∞ and γ4∞ are the bulk concentration and activity coefficient of the Mg 2+ ions and KSt,4 is their Stern constant. Because both Na + and Mg 2+ counterions can adsorb in the Stern layer, Equation (3.17) has to be replaced with the following two Stern isotherms [27]: where Γ4 is the adsorption of Mg 2+ ions. The Stern constant for Mg 2+ ions, KSt,4 = 1.92 × 10 −4 (mM) −1 , has been determined in reference [43] from data for the surface tension of SDS solutions in the presence of MgSO4. The radius of the Mg 2+ ion is also known, r4 = 0.43 nm [3]. The values of all other physical parameters are the same as in Section 4.1.
The most important conclusions from Figure 8a are as follows. (i) It is remarkable that the DLVO-IC theoretical curve, drawn without using any adjustable parameters, is in excellent agreement with the experimental data for 2:2 electrolytes. (ii) The values of Π predicted by DLVO theory in the case of a 2:2 electrolyte are considerably greater than the experimental ones. (iii) The Π(hw) isotherm for a 2:2 electrolyte, as compared to that for a 1:1 electrolyte (at the same ionic As seen from the experimental values of h w (Figure 8a), in the presence of MgSO 4 the equilibrium foam films are considerably thinner than in the presence of NaCl at the same ionic strength. Moreover, the maximal values of Π that can be achieved experimentally are not so high because a transition from primary (common black) film to secondary (Newton black) film takes place; see e.g., Figure 7.17 in reference [42].
Theoretically, in this case we are dealing with four ionic components, which are numbered as follows: 1-surfactant ion (dodecyl sulfate); 2-Na + ion (from SDS); 3-SO 4 2´a nion and 4-Mg 2+ cation (from MgSO 4 ). All equations in Section 3 remain the same, except the surfactant adsorption isotherm, Equation (3.16) and the counterion binding isotherm, Equation (3.17). Equation (3.16) has to be replaced with the following equation [27]: where c 48 and γ 48 are the bulk concentration and activity coefficient of the Mg 2+ ions and K St,4 is their Stern constant. Because both Na + and Mg 2+ counterions can adsorb in the Stern layer, Equation (3.17) has to be replaced with the following two Stern isotherms [27]:  [43] from data for the surface tension of SDS solutions in the presence of MgSO 4 . The radius of the Mg 2+ ion is also known, r 4 = 0.43 nm [3]. The values of all other physical parameters are the same as in Section 4.1.
The most important conclusions from Figure 8a are as follows. (i) It is remarkable that the DLVO-IC theoretical curve, drawn without using any adjustable parameters, is in excellent agreement with the experimental data for 2:2 electrolytes. (ii) The values of Π predicted by DLVO theory in the case of a 2:2 electrolyte are considerably greater than the experimental ones. (iii) The Π(h w ) isotherm for a 2:2 electrolyte, as compared to that for a 1:1 electrolyte (at the same ionic strength), is shifted considerably to the left, which indicates stronger suppression of the electrostatic repulsion in the thin liquid films.
To understand the reasons for the last effect, in Figure 8b,c we have plotted the surface charge and surface potential as predicted by DLVO-IC. The surface charge is greater in the case of a 2:2 electrolyte (Figure 8a), because the binding constant of the Mg 2+ ions, K St,4 , is 3.4 times smaller than that of the Na + ions, K St,2 (see above). This could be explained by the greater degree of hydration of the Mg 2+ ions. In contrast, the surface potential, |ϕ s |, is about 2 times lower in the case of a 2:2 electrolyte (Figure 8b). This is due to the higher valence of the Mg 2+ ions, which leads to a significant rise of their concentration in the subsurface layer, in accordance with the Boltzmann law, Equation (3.10). In other words, the more effective screening of the electrostatic interactions in the case of a 2:2 electrolyte is due to the increased concentration of divalent counterions in the diffuse part of the electric double layer, which also leads to enhancement of the ion-correlation effects, as evidenced by the difference between the DLVO and DLVO-IC disjoining pressure isotherms in Figure 8a.

Analytical Expression for the Asymptotic Screening Parameter
Here, we derive an analytical expression for the asymptotic screening parameter, κ a , which allows calculation of this parameter from the valences and radii of the ions present in the solution.
In terms of the dimensionless electric potential, Φ, the Poisson equation can be expressed in the form: where we assume that the bulk solution is electroneutral, i.e., Σ i z i c i8 = 0. At sufficiently large values of the film thickness, h, the potential near the film's midplane is low and we could write: where ε ci and ε γi are small quantities. Substituting the quantities in Equation (5.2) into Equation (3.10) and keeping the leading terms in the power expansion, we obtain: To simplify the further derivations, let us assume that all counterions can be characterized with a mean ionic radius, r m . For r i = r m (i = 1, . . . , n), without any approximations Equation (3.6) acquires a much simpler form: Because of the general character of Equation (5.4), it holds also in the bulk of solution: where κ 8 is the conventional Debye parameter; see Equation (3.20). Equation (5.5) represents a generalized form of the Debye-Hückel limiting law and has already been derived in a different way [44]. Substituting γ i « γ i8 p1`ε γi q in Equation (5.4) and expanding in series for small (κ´κ 8 ), in view of Equation (5.5) we derive: Further, we substitute Equation (5.6) into Equation (5.3) and multiply by c i8 : Substituting c i « c i8 p1`ε ci q in the definition of κ, Equation (3.7), and expanding in series for small values of ε ci and (κ´κ 8 ), we obtain: To find an expression for the last term in Equation (5.8), we multiply Equation (5.7) with z i 2 and sum up: The substitution of Equation (5.9) in the right-hand side of Equation (5.8) yields: Finally, the substitution of (κ´κ 8 ) from Equation (5.10) into Equation (5.7), followed by substitution of the obtained expression for (c i´ci8 ) in Equation (5.1), yields: where the screening parameter κ a is defined as follows: In the case of EDL at a single interface, the solution of Equation (5.11) reads: where x is the distance to the interface; and Φ a is a constant pre-exponential factor of the asymptotics of electric potential. Furthermore, setting r i = r m (i = 1, . . . , n), as above, and applying analogous series expansion for small Φ in Equation (3.8), one can derive: where A is a constant coefficient in this asymptotic expression for Π el . Finally, following the approach by Verwey and Overbeek [2], we apply the superposition approximation for low potentials in the film midplane, viz. Φ m = 2Φ 1 (h/2), which in view of Equations (5.13) and (5.14) yields: Π el " 4AΦ 2 a expp´κ a hq (5.15) Equations (5.13) and (5.15) imply that if the ionic correlations are taken into account, the asymptotic decay length of both the potential Φ 1 and disjoining pressure Π el is κ a´1 , where κ a is given by Equation (5.12). In other words, κ a is the generalized screening parameter that takes into account the effects of ionic correlations and finite ionic radius, r m . The form of Equation (5.12) leads to the following conclusions: (i) The screening parameter κ a depends not only on the sum ř i z 2 i c i8 that enters the expression for the conventional Debye parameter κ 8 , but also on the sums ř i z 3 i c i8 and ř i z 4 i c i8 . This leads to a considerable effect of ions of higher valence, including coions, on the Debye screening. (ii) Equation (5.12) shows that κ a ě κ 8 , which leads to a stronger screening and weaker electrostatic repulsion as a result of the ionic correlation effect: compare the DLVO-IC with the DLVO curves in Figure 7b. (iii) The effect of ionic correlations is a long-range effect-it influences the long-range asymptotics of Π el and can be significant when the film thickness is tens of nanometers, not only in films a few nanometers thick. (iv) If all dissolved electrolytes are symmetric, then ř i z 3 i c i8 = 0 and Equation (5.12) yields κ a = κ 8 , i.e., the ionic correlation effect on the asymptotic decay length disappears for symmetric electrolytes, supposedly the electric potential Φ m in the film's midplane is low. This fact is related to the coincidence of the predictions of DLVO and DLVO-IC theories for NaCl in Figure 2.
(v) The ionic radius, r m , essentially affects the value of the screening parameter κ a . If we formally set r m = 0 in Equation (5.12), the calculated κ a will be markedly different from the asymptotic slope of the exact numerical solution in Figure 7b.
To calculate the values of κ a in Table 1 from Equation (5.12), we have set the mean ionic radius equal to that of the hydrated Na + counterions, which are the most abundant in the film, i.e., r m = 0.36 nm.

Conclusions
In this article, we report experimental data for the disjoining pressure of foam films stabilized by anionic surfactant in the presence of 1:1, 1:2, and 1:3 electrolytes: NaCl, Na 2 SO 4 , and Na 3 Citrate. The disjoining pressure predicted by the DLVO theory coincides with the experimental data in the case of a 1:1 electrolyte, but it is considerably greater than the measured pressure in the cases of divalent and trivalent coions. To achieve agreement with the experiment, the DLVO theory is extended to account for the effects of ionic correlations and finite ionic radii. Original analytical expressions are derived for the local activity coefficient, electrostatic component of disjoining pressure, and asymptotic screening parameter; see Equations (3.6), (3.8), (3.11) and (5.12). With the same parameters of the counterion-binding and van-der-Waals-disjoining-pressure isotherms, as for a 1:1 electrolyte, the curves predicted by the extended theory are in perfect agreement with the experimental data for 1:2 and 1:3 electrolytes. The obtained formula for the screening parameter, κ a , exactly predicts the asymptotic decay length of disjoining pressure. In the investigated range of ionic strengths and surface charge densities, the effect of ionic correlations turns out to be essential if 1:2, 1:3 and 2:2 electrolytes are present in the solution.
In comparison with the DLVO theory, the effect of ionic correlations leads to more effective screening of the electrostatic interactions in the electric double layer, as evidenced by the lower values of electric potential in the middle of the film (Figure 6d). This results in a lower concentration of counterions in the film's midplane (Figure 6b) and lower disjoining pressure, as experimentally observed (Figure 2). The effect of ionic correlations on disjoining pressure is a result of a fine balance of electrostatic and osmotic effects, plus the effects of counterion binding (adsorption), which can be accurately described only by using the respective consistent system of equations; see Sections 3.1 and 3.2.
In the considered case, it is possible to formally fit the experimental data using the DLVO theory by adjusting an effective (empirical) surface potential for each separate experimental curve or by introducing different adjustable Stern constants for the different experimental curves. However, physically the Stern constant K St,2 must be the same for the three curves in Figure 2 insofar as the counterions are the same (Na + ) and the interface is the same (dense adsorption layer of SDS); moreover, K St,2 is known from independent surface tension measurements [27][28][29][30]. In such a case, the experimental data cannot be quantitatively interpreted by the DLVO theory without taking into account the effect of ionic correlations.
For a symmetric 1:1 electrolyte, no effect of ionic correlations was observed. In contrast, for a symmetric 2:2 electrolyte a considerable ionic-correlation effect is present (Figure 8a). The fact that in the latter case the calculated disjoining-pressure isotherm coincides with the experimental curve without using any adjustable parameters is an additional argument in favor of the correctness of the DLVO-IC theory.
The extended theory shows that the ionic-correlation effects are essential in the whole range of film thicknesses, from nanometers to tens of nanometers, i.e., it affects both the short-range and long-range surface forces. The effect of ionic radius, r i , is coupled with the effect of ionic correlations and it affects even the long-range asymptotics of disjoining pressure; see Equation (5.12). Although the pressure in the thin liquid film can be presented as a sum of terms related to the ionic concentrations and activity coefficients, see Equation (3.8), the disjoining pressure cannot be presented as a sum of electrostatic and ionic-correlation terms (e.g., Π el + Π cor ), because the ionic correlations affect not only the activity coefficients, but also the ionic concentrations.
The extended DLVO theory with ionic correlations (DLVO-IC) is applicable to both multivalent coions and counterions. One possible extension of the present study is its application to systems with other multivalent counterions. We hope that the data analysis by DLVO-IC could remove some discrepancies between theory and experiment observed in studies with liquid films from electrolyte solutions.