Potential of mean force between oppositely charged nanoparticles: A comprehensive comparison between Poisson–Boltzmann theory and Monte Carlo simulations

Ion-mediated interactions between like-charged polyelectrolytes have been paid much attention, and the Poisson–Boltzmann (PB) theory has been shown to fail in qualitatively predicting multivalent ion-mediated like-charge attraction. However, inadequate attention has been paid to the ion-mediated interactions between oppositely charged polyelectrolytes. In this work, the potentials of mean force (PMF) between oppositely charged nanoparticles in 1:1 and 2:2 salt solutions were investigated by Monte Carlo simulations and the PB theory. Our calculations show that the PMFs between oppositely charged nanoparticles are generally attractive in 1:1 and 2:2 salt solutions and that such attractive PMFs become weaker at higher 1:1 or 2:2 salt concentrations. The comprehensive comparisons show that the PB theory can quantitatively predict the PMFs between oppositely charged nanoparticles in 1:1 salt solutions, except for the slight deviation at very high 1:1 salt concentration. However, for 2:2 salt solutions, the PB theory generally overestimates the attractive PMF between oppositely charged nanoparticles, and this overestimation becomes more pronounced for nanoparticles with higher charge density and for higher 2:2 salt concentration. Our microscopic analyses suggest that the overestimation of the PB theory on the attractive PMFs for 2:2 salt solutions is attributed to the underestimation of divalent ions bound to nanoparticles.

Importantly, the aggregation/assembly of oppositely charged polyelectrolytes is strongly dependent on the effective interaction between them. The existing experiments and theories including the PB theory and weak-/ strong-coupling theory have shown that interactions between oppositely charged polyelectrolytes can be modulated by monovalent and multivalent salt ions 34,[68][69][70][71][72][73][74][75][76] . In asymmetrical salt solution, repulsive interactions between oppositely charged polyelectrolytes were observed in experiments and simulations at very high salt concentrations [73][74][75][76] . Moreover, Borkovec et al. employed the PB theory with involving charge regulation to calculate the effective interactions between oppositely charged colloids in monovalent and multivalent ion solutions, and they compared the results with their experiments. They concluded that the PB theory can well describe the effective interactions between oppositely charged polyelectrolytes mediated by multivalent ions at large separation and low salt concentration [68][69][70] . The charge regulation involved in their PB calculations can be attributed to the change of effective charges on colloids due to ion binding or releasing in the inner layers of colloids [68][69][70] , and the charge regulation parameters were obtained through fitting their PB calculations to the experimental measurements for the corresponding symmetric colloid systems [68][69][70] . Therefore, it still remains unclear whether the PB theory in the original version can describe the effective interaction between oppositely charged polyelectrolytes, particularly in the presence of multivalent salt, given that the PB theory qualitatively fails to describe the multivalent-mediated attractions between like-charged polyelectrolytes 46,49,50 . In addition, there remains a lack of comprehensive understanding on the effective interaction between oppositely charged polyelectrolytes modulated by salt ions, especially multivalent salt ions.
In this work, we investigated the potential of mean force (PMF) between oppositely charged nanoparticles in symmetrical salt solutions by MC simulations and the PB theory. Our calculations were performed over wide ranges of ion conditions and charge densities on nanoparticles and would provide an overall illustration of ion-mediated PMFs between oppositely charged nanoparticles. In addition, the present work provides an extensive examination on when the PB theory works well for oppositely charged polyelectrolytes and microscopic analyses of the corresponding microscopic mechanism.

Results and Discussion
In this work, the PMFs ∆G(x) between oppositely charged nanoparticles immersed in 1:1 and 2:2 salt solutions were calculated by the PB theory and MC simulations, as displayed in Fig. 1. Our calculations were performed over wide ranges of 1:1 and 2:2 salt concentrations and charge densities on nanoparticles. We have comprehensively compared the PB theory and MC simulations, and analysed when the PB theory works well and the corresponding microscopic mechanism in detail.  Coulomb attraction, and higher ion concentration can cause stronger ionic screening for the oppositely charged nanoparticles owing to the lower entropy penalty for bound ions. Furthermore, with increasing charges on nanoparticles, the attractive PMFs become stronger. This is reasonable because the increasing charge density on nanoparticles can enhance the opposite-charge Coulomb attraction between nanoparticles and consequently enhance the ion-modulated attractive PMFs.
Figures 2(A-C) and S1A,B in the SM also show that the PMFs calculated from the PB theory are generally in quantitative accordance with those from the MC simulations for nanoparticles with different charges in 1:1 salt solutions, and the PB theory slightly overestimates the attractive PMFs between oppositely charged nanoparticles at very high 1:1 salt concentrations (e.g., ~0.3 M). This is not surprising even if the PB theory is generally considered to work well for polyelectrolytes in monovalent salt solutions. The PB theory assumes a continuous fluid-like ion distribution and consequently ignores the discrete properties of ions. For monovalent salt, although ion-ion Coulomb correlations are generally weak owing to the low valence, the exclusion volume correlation can become important at very high salt concentration 35,[42][43][44][45]50 , and this exclusion correlation is ignored in the PB theory. As a result, the deviation between the PB theory and MC simulations for nanoparticles with x = 22 Å is approximately zero at low and moderate 1:1 salt concentrations, and it becomes visible at very high 1:1 salt concentration (e.g., the relative deviation of PMF ~10% at 0.3 M 1:1 salt), as shown in Figs S1 and S2 in the SM.
To gain a deep understanding on the PMFs, the bound ions around oppositely charged nanoparticles were analysed through calculating the net ion charge fraction Q(r) within a distance r from the centre of each nanoparticle. In our cases, the electrostatic potential is zero at the plane in the middle of two oppositely charged nanoparticles in symmetric 1:1 (or 2:2) salt, and this middle plane splits the space into two parts, in which the bound ions are dominated by the respective charged nanoparticles. Therefore, we calculate the net ion charge fraction Q(r) in one side of the middle plane 49,50,77 where Z is the charge on a nanoparticle, and r is a radial distance from the centre of a nanoparticle. z i is the valence of ion species i, and c i (r) denotes its concentration at r. Here, the net ion charge fraction Q(r) represents the radial distribution of net bound ion charges and is generally calculated for the nanoparticles with separation  underestimates the net ion charge fraction, suggesting fewer bound counterions to charged nanoparticles from the PB theory than those from the MC simulations. The slightly fewer bound counterions contribute to the above described slight overestimation on the attractive PMFs of the PB theory, because counterions provide the major contribution to reduce the Coulomb attraction between oppositely charged nanoparticles. Overall, the deviation between the predictions from the MC simulations and PB theory is insignificant on both the PMFs and the ion neutralization fraction for 1:1 salt solutions; see Figs S1 and S2 in the SM.

2:2 salt solutions.
The PMFs between oppositely charged nanoparticles in 2:2 salt solutions were calculated over wide ranges of 2:2 salt concentrations and charge densities on nanoparticles. As shown in Figs 3A-C and S3A,B in the SM, in analogy to 1:1 salt, the PMFs between oppositely charged nanoparticles are always attractive in 2:2 salt solutions. This effective attraction becomes weakened with increasing 2:2 salt concentration, which is attributed to more bound ions at higher ion concentration.
Figures 3A-C and S3A,B in the SM show that compared with the PMFs from the MC simulations, the PB theory generally overestimates the effective attraction between oppositely charged nanoparticles in 2:2 salt solutions. This overestimation becomes more pronounced for nanoparticles with high charge densities and for high 2:2 salt concentrations. For example, for the nanoparticles with |Z| = 12e (e is the elementary charge) in 0.1 mM 2:2 salt solution, the relative deviation of PMFs at x = 22 Å between the PB theory and MC simulations is ~7% and can become as large as ~73% for nanoparticles with |Z| = 24e in 10 mM 2:2 salt solution. This suggests that for 2:2 salt, the PB theory generally qualitatively predicts the effective interactions between oppositely charged nanoparticles rather than quantitative predictions, and the deviation of the PB predictions can become very pronounced for nanoparticles with high charges in high 2:2 salt solutions.
To understand the severe deviations between the PMFs from the PB theory and MC simulations, we analysed the bound ions around oppositely charged nanoparticles by calculating the net ion charge fraction Q(r) within a distance r from the centre of each nanoparticle; see Eq. 1. As shown in Figs 3D-F and S3C,D in the SM, the net ion charge fraction Q(r) from the MC simulations is generally greater than that from the PB theory except for the nanoparticles with very low charge density. The underestimation of the PB theory on Q(r) becomes more pronounced for nanoparticles with higher charge density and for higher 2:2 salt concentration. For example, ∆Q(r) = (|Q MC (r) − Q PB (r)|) can become as large as ~0.25 at r = 14 Å for |Z| = 24e and 10 mM 2:2 salt concentration. The underestimation of the PB theory on Q(r) is in accordance with the above described overestimation on the attractive PMFs between oppositely charged nanoparticles in 2:2 salt solutions. The apparent underestimation of the PB theory on bound ions is attributed to the high charge of 2:2 salt ions, which can cause strong Coulomb correlations between divalent ions. As described above and previously 35  simulations explicitly account for the ion-ion correlations. The ion-ion correlations can drive ions to self-organize to low-energy micro-states and generally cause more bound ions to a charged polyelectrolyte 35,40,[42][43][44][45]50,78,79 . Therefore, the PB theory generally underestimates the bound divalent ions and consequently overestimates the attractive PMFs between oppositely charged nanoparticles in 2:2 salt solutions. With increasing charge density on nanoparticles or increasing 2:2 salt concentration, the ion-ion correlations near the surface of nanoparticles would become stronger, and consequently the deviation of the PB predictions on bound ions and the resultant PMFs would become more pronounced; see Figs S3 and S4 in the SM.
1:1 salt versus 2:2 salt. Physically, multivalent ions can interact more strongly with charged nanoparticles and consequently can bind more strongly to charged nanoparticles than monovalent ions. The larger Q(r) of bound multivalent ions relative to monovalent ions becomes more pronounced for nanoparticles with higher charges; see Figs 2D-F and 3D-F. As a result, compared with 1:1 salt, 2:2 salt is more efficient in screening the Coulomb attraction between oppositely charged nanoparticles, and this efficiency becomes more pronounced for nanoparticles with higher |Z|, as shown in Figs 2 and 3. For example, in 30 mM 1:1 salt and 1 mM 2:2 salt solutions, the attractive PMFs are almost identical for the nanoparticles with |Z| = 12e, whereas for the nanoparticles with |Z| = 18e, the PMF in 1 mM 2:2 salt solution becomes visibly less attractive than that in 30 mM 1:1 salt solution.
To quantitatively compare the efficiency of monovalent and divalent ions in screening the Coulomb attraction between oppositely charged nanoparticles, we calculated the equivalent 1:1 salt concentrations, which can cause (nearly) the same PMFs as those at 0.1 mM, 1 mM, and 10 mM 2:2 salt concentrations, respectively. Practically, the interpolation technique was used to obtain this equivalent 1:1 salt because the 1:1 concentrations cannot be continuously covered in our MC simulations, and we used ∆G(  Fig. 4A, divalent ions are much more efficient than monovalent ions, and this higher efficiency of divalent ions over monovalent ions becomes more pronounced for nanoparticles with higher charge density. For example, 0.1 mM and 10 mM 2:2 salts can cause (nearly) the same attractive PMFs as ~3 mM and ~80 mM 1:1 salts for nanoparticles with |Z| = 9e, respectively. However, for nanoparticles with |Z| = 24e, 0.1 mM and 10 mM 2:2 salts can cause the same PMFs as ~30 mM and ~300 mM 1:1 salts, respectively.
To understand the efficient role of divalent ions over monovalent ions in screening the Coulomb attraction between oppositely charged nanoparticles, we further calculated another kind of equivalent 1:1 salt concentration to 2:2 salt by comparing the net ion charge fraction Q(r); see Eq. 1. Practically, we also used the interpolation technique and Q(r = 22 Å; three layers of ions) with the separation x = 22 Å to compare 1:1 salts and 2:2 salts. Here, the large value of r = 22 Å was used to compare 2:2 salt with 1:1 salt of bound ions because monovalent ions bind loosely relative to divalent ions. As shown in Fig. 4B, the relationship between the equivalent 1:1 salt and 2:2 salts from Q(r = 22 Å) shares a similar trend to that from ∆G(x = 22 Å) in Fig. 4A; i.e., to achieve the same ionic neutralization, divalent ions are much more efficient than monovalent ions, and this apparently higher efficiency of divalent ions becomes more pronounced for nanoparticles with higher charge density. Therefore, the comparison between Fig. 4A and B clearly shows that the more efficient role of divalent ions than monovalent ions in

Deviation between the PB theory and MC simulations and an apparent parameter.
To obtain a direct understanding on the deviation between the PMFs from the PB theory and MC simulations, we further calculated the charge density distribution landscapes that include both cations and anions around the nanoparticles with |Z| = 18e and x = 22 Å. As shown in Fig. 5, the distributions of net ion charge density around nanoparticles from the PB theory and MC simulations share similar landscapes. For high 1:1 salt concentration, the deviation of ion net charge density appears slight, whereas such deviation becomes pronounced for high 2:2 salt concentration.
In the following, we quantitatively examine the relative deviation of the PMFs between the PB theory and the MC simulations, which is characterized by here, x 0 = 22 Å was used because the maximum deviation of PMFs between the PB theory and the MC simulations are generally located at x 0 = 22 Å; see Figs S2A-C and S4A-C in the SM. As shown in Fig. 6A,B, the relative difference ∆∆g between the PMFs from the PB theory and the MC simulations is generally rather small at low 1:1 salt concentrations and becomes visible at high 1:1 salt concentrations. For example, for |Z| = 24e, ∆∆g is ~3.5% in 0.03 M 1:1 salt and can reach ~15% in 0.3 M 1:1 salt. In contrast, for 2:2 salt, the relative deviation ∆∆g generally appears significant except for nanoparticles with very low |Z| and becomes more pronounced for nanoparticles with higher |Z| and for higher 2:2 salt concentrations. For example, in 0.1 mM 2:2 salt solution, ∆∆g is ≤~7% for the nanoparticles with |Z| ≤ 12e and becomes ~34% for those with |Z| = 24e. When the 2:2 salt concentration is increased to 10 mM, ∆∆g becomes as large as ~18% for the nanoparticles with |Z| = 12e and ~73% for those with |Z| = 24e. Therefore, for two oppositely charged nanoparticles in 1:1 salt solution, the PB theory generally predicts the PMFs between them reliably except for very high salt concentration (~0.3 M). However, for oppositely charged nanoparticles in 2:2 salt solution, the PB theory generally apparently overestimates the attractive PMFs except for the nanoparticles with very low |Z|. To further characterize the relative deviation ∆∆g of the PMFs between the PB theory and MC simulations, an apparent parameter |Z|∆Q * of the difference in total net ion charges was defined as |Z|∆Q(r) at r = 14 Å, where the maximum values of |Z|∆Q(r) are located at high 1:1 and 2:2 salts; see Figs S2D-F and S4D-F in the SM. Figure 6C showed the relative difference ∆∆g as a function of |Z|∆Q* for both 1:1 and 2:2 salts. It is interesting that ∆∆g appears tightly and positively coupled to |Z|∆Q * over wide ranges of salt conditions and charge densities on nanoparticles. This tight and positive coupling between ∆∆g and |Z|∆Q * is reasonable. Because bound ions mainly provide ionic neutralization to reduce the Coulomb attraction between oppositely charged nanoparticles, the large difference in bound ion charges should correspond to the large deviation in the PMFs between the PB theory and MC simulations. Therefore, |Z|∆Q * can be served as an apparent parameter to well characterize the deviation degree of the PMFs between the PB theory and MC simulations for oppositely charged nanoparticles.

Conclusions
In this work, both the PB theory and MC simulations were employed to investigate the PMFs between oppositely charged nanoparticles in 1:1 and 2:2 salt solutions. The work covers wide ranges of 1:1 and 2:2 salt concentrations and charge densities on nanoparticles. Through systematic calculations and detailed analyses, the following conclusions have been obtained: (1) For 1:1 salt solutions, the PMFs between oppositely charged nanoparticles are always attractive and become more attractive for nanoparticles with higher charge density and in lower 1:1 salt solution. The PMFs from the PB theory are generally in good agreement with those from the MC simulations except for the slight deviation at very high 1:1 salt concentrations, which suggests that the PB theory can generally describe the PMFs between oppositely charged nanoparticles in 1:1 salt solutions. (2) For 2:2 salt solutions, the PMFs between oppositely charged nanoparticles are always attractive, and this effective attraction becomes stronger for nanoparticles with higher charge density and lower 2:2 salt concentration. The PB theory generally overestimates the attractive PMFs between oppositely charged nanoparticles in 2:2 salt solutions, especially for highly charged nanoparticles and for high 2:2 salt concentrations. (3) Compared with 1:1 salt, 2:2 salt is more efficient in screening the Coulomb attraction between oppositely charged nanoparticles, and this efficient role of 2:2 salt solutions is more pronounced for more highly charged nanoparticles. (4) The overestimation of the PB theory on the attractive PMFs is attributed to the underestimation of bound ions around charged nanoparticles. The relative deviation in the PMFs is tightly and positively coupled to an apparent parameter |Z|∆Q* of the difference in the total net ion charges between the PB theory and MC simulation for oppositely charged nanoparticles.
Our conclusion that the PB theory generally overestimates the attractive PMFs between oppositely charged nanoparticles in 2:2 salt solutions is attributed to the higher ionic charge and the resultant strong ion correlations, as discussed above. This is not contradictory to the combination work of atomic force microscopy and the PB equation with charge regulation, which can be attributed to the change of effective charges on nanoparticles due to possible ion binding or releasing in the inner layer of nanoparticles [68][69][70] . The present work was mainly focused on the nanoparticles with constant charges, and our PB calculations based on the original version of the PB theory did not involve charge regulation in the inner layers of nanoparticles. Our conclusions do not exclude the validity of the PB equation with charge regulation on the PMFs between oppositely charged nanoparticles mediated by multivalent ions, especially for low salt concentrations and at large separations [68][69][70] , because the effective charges on nanoparticles can be apparently reduced by the ion binding in the inner layers of nanoparticles and consequently the electrostatic potential in the diffusive layers can be greatly reduced [68][69][70] . Moreover, in the present Here, x = 22 Å was used because the deviation between PMFs from the work, the PMFs between oppositely charged nanoparticles are calculated as the effective attractions in 1:1 and 2:2 salt solutions, which is also not contradictory to the effective repulsions observed in previous experimental and theoretical studies 42,43,[71][72][73][74][75][76] . The effective repulsions in previous studies might be attributed to the adopted very high asymmetrical salts and the resultant over-neutralization of bound multivalent ions to the nanoparticles 42,43,73,75,76 . However, the present work involves only symmetrical 1:1 and 2:2 salts and does not cover very high salt concentrations, thus only predicts effective attractions between oppositely charged nanoparticles.
The present work also involves some approximations and simplifications. First, the solvent molecules were modelled implicitly as a continuum medium with a high dielectric constant of 78, and for simplicity, the dielectric discontinuity at the interface between nanoparticles and solvent was ignored 35,50,80,81 . This assumption may slightly affect the bound ions at the nanoparticle surface and should be considered in future work. Second, the nanoparticles and ions were treated as hard-core spheres, which should not essentially affect our results and analyses. Finally, the present work did not cover the effective interactions at very large separation (e.g., ≫ nanoparticle size), because the assembly of nanoparticles generally involves interactions within a moderate separation range [68][69][70] . Nevertheless, the present work illustrates that the effective attraction between oppositely charged nanoparticles can be greatly modulated by salt ions, especially by multivalent ions. In addition, this work suggests that the PB theory can generally provide reliable predictions on the effective interactions between oppositely charged nanoparticles in 1:1 salt solutions, while it cannot quantitatively describe the PMFs between oppositely charged nanoparticles in multivalent salt solutions. The present work would be very helpful for understanding the ion-mediated assembly and complexation of oppositely charged polyelectrolytes [42][43][44][45][71][72][73][74][75][76] .

Model and Method
Model system for oppositely charged nanoparticles. In this work, for simplicity, we employed a model system of two oppositely charged nanoparticles immersed in 1:1 and 2:2 salt solutions to study the PMFs between oppositely charged nanoparticles [50][51][52] . In the model system, both charged nanoparticles and salt ions were represented by spheres with radii of 10 Å and 2 Å, and the solvent was modelled as a continuous medium with dielectric constant ε = 78; see Fig. 1. The charges on two oppositely charged nanoparticles were taken as +6e/−6e, +9e/−9e, +12e/−12e, +18e/−18e and +24e/−24e. In the model system, the interactions between nanoparticles and ions were simplified into Coulomb and hard-sphere repulsion interactions, and the interaction energy U ij between spheres i and j in the model system is given by Here, σ i and q i represent the radius and charge of particle i (nanoparticles or ions). r ij is the distance between the centres of two charged particles i and j (nanoparticles and ions). ε is the dielectric constant of solvent, and ε 0 is the permittivity of vacuum. Thus, our systems involved only the effects of Coulomb and exclusion volume interactions and did not involve possible chemical processes such as specific proton adsorption or surface ionization in some colloid systems 68-70 . Monte Carlo simulations for PMFs. For the model system, the MC simulations with the Metropolis algorithm have been employed to calculate the PMFs between oppositely charged nanoparticles in salt solutions. The simulation cell was generally a rectangular box, and to diminish the boundary effect, the box size was kept larger than the centre-to-centre separation x between oppositely charged nanoparticles by at least six Debye-Hückel lengths [35][36][37]50 .
Following previous studies 46,49,53 , we employed the pseudo-spring method to calculate PMFs 46,49,52 , and a spring with spring constant k = 9 nN/Å was added to link the centres of the two oppositely charged nanoparticles. In our MC simulations, one nanoparticle remains frozen, while the other can move along x-axis with a constraint of the added spring. The effective force F(x) between the two nanoparticles with a centre-to-centre separation x can be given by 46,49,53 where ∆x is the small deviation of spring length away from the original distance x at equilibrium. The convergence of the average separation between two nanoparticles is shown in Fig. S6 of the SM. Afterwards, the PMF ∆G(x) between the two nanoparticles can be calculated by the following integration 46,49,53 :  denotes the ion concentration of species i at r. It is noted that the electrostatic potential r ( ) ψ changes with separation x between two nanoparticles. The derivation of Eq. 7 is detailed in the SM; see Section 2 in the SM. Therefore, PB ref gives the PMF between two nanoparticles by the PB theory, where x ref is taken as 40 Å [35][36][37]46,49,50,53 . Thus, our PB calculations are based on the original version of the PB theory, without considering charge regulation due to ion binding or releasing in the inner layers of nanoparticles [68][69][70] .
In this work, the three-dimensional algorithm developed in the tightly bound ion theory was used to numerically solve the PB equation [35][36][37]50 . To account for the excluded volume layer of ions, a thin layer of ion radius was added to the nanoparticle surface, and the charge density ρ f of fixed charges on nanoparticles can be obtained approximately by partitioning the fixed charges to their eight nearest grid points when solving the PB equation on a cubic lattice with the finite-difference method [35][36][37]50 . The focusing process of three steps has been used to compute the detailed electrostatic potential near the nanoparticle surface [35][36][37]50 . In the first step, the PB equation was solved numerically on a cubic box with large grid size, which depends on salt concentration. Generally, we chose a box-size approximately six times larger than the Debye-Hückel length from the nanoparticle surface to include the salt effect in solution [35][36][37]50 . In the second step of focusing, the PB was solved on a smaller cubic box with smaller grid size, and the initial electrostatic potentials were obtained from the interpolation based on those in the first step. A similar procedure was repeated in the third focusing step to obtain the electrostatic potentials with high resolution near the nanoparticles. The resolution of the first step varies with the grid size to increase the efficiency of the iterative process. The grid sizes (L x , L y , L z ) were kept at (160 Å, 120 Å, 120 Å) and (100 Å, 60 Å, 60 Å) in the second and the third focusing steps, respectively, and the corresponding resolutions were 1.0 Å and 0.25 Å per grid size. As a result, the numbers of grid points were 161 × 121 × 121 in the second step and 401 × 241 × 241 in the third step. The iteration for each step was continued until the electrostatic potential change δψ(r) for an iteration was less than 10 −4 k B T/e [35][36][37]50 .