Ionic structure around polarizable metal nanoparticles in aqueous electrolytes

Metal nanoparticles are receiving increased scientific attention owing to their unique physical and chemical properties that make them suitable for a wide range of applications in diverse fields, such as electrochemistry, biochemistry, and nanomedicine. Their high metallic polarizability is a crucial determinant that defines their electrostatic character in various electrolyte solutions. Here, we introduce a continuum-based model of a metal nanoparticle with explicit polarizability in the presence of different kinds of electrolytes. We employ several, variously sophisticated, theoretical approaches, corroborated by Monte Carlo simulations in order to elucidate the basic electrostatics principles of the model. We investigate how different kinds of asymmetries between the ions result in non-trivial phenomena, such as charge separation and a build-up of a so-called zero surface-charge double layer.


I. INTRODUCTION
The field of electrostatic interactions in classical softmatter and biological systems has a long and rich history, recognized by many intellectual challenges and ideas. [1][2][3] Maybe one of the most remarkable concepts is the approximation of implicit solvent, a continuum-level approach where a system composed of charged species and solvent molecules is simply treated as a gas of the charged species only, now with their interactions attenuated by the relative dielectric permittivity ε of the solvent. [4] Remarkably, this so-called primitive model or dielectric approximation works very well for simple ions in aqueous solutions down to only several layers of water molecules between the ions. [5,6] This is also one of the reasons for the efficiency of the Poisson-Boltzmann equation to describe monovalent ions in the water environment. Nonetheless, an implementation of this approximation can become technically involved if dielectric discontinuities are present in the system. This is unfortunately also one of the reasons why dielectric discontinuities are commonly neglected in (too) many studies. Neglecting them is not always justified, because in the presence of charges, a dielectric discontinuity leads to a polarization surface charge density at the boundary, which further influences the local electrostatic potential and interactions with surrounding charges. In the case of a planar dielectric discontinuity, the * matej.kanduc@helmholtz-berlin.de electrostatic potential can simply be expressed as the electrostatic potential arising from a fictive "image charge" residing on the other side of the discontinuity. Therefore, in this context, the dielectric discontinuity effects are sometimes referred to as image charges. The dielectric effects in double-layer problems of planar geometry have been elaborated by Torrie, Valleau, and Patey, [7] and by Bratko, Jonsson, and Wennerström [8] using computer simulations, and by Kjellander and Marčelja [9] and Outhwaite, Bhuiyan, and Levine [10] utilizing various theoretical frameworks. The image-charge concepts have been adapted to spherical symmetry by Linse. [1,11,12] He showed that approximating the exact mathematical expressions for the spherical geometry leads to a simplified picture in which the polarization is described by image charges as in planar cases. The image charges in the spherical geometry are of paramount importance, since a vast majority of the soft-matter electrostatics research in the recent decades has focused on colloidal and biological systems, where various macromolecular structures (e.g., colloids, proteins, polysaccharides, micelles) in water can be modeled as spherical entities with a lower dielectric interior ε (due to their predominantly hydrocarbon architectures) than the surrounding water environment (ε ε). [1,[13][14][15][16][17][18][19][20][21][22][23][24] The other side of the spectrum, containing spherical bodies of a much higher dielectric interior than water (i.e., ε ε), such as for example small metal particles in aqueous environments, has been much less explored. However, the interest in this field has boosted with recent advances in metal nanoparti-cle chemistry and physics, which have emerged as a broad new discipline in a subdomain of colloids and surfaces. [25,26] One of the most prominent discoveries was that gold nanoparticles (of a size 1-10 nm) are active catalysts for oxidation reactions. [27] This has triggered a tremendous research activity in nanocatalysis, which presently remains one of the fastest growing areas of nanoscience. [28][29][30] Furthermore, applications involving metal nanoparticles can for instance be found in electrochemistry for nanoelectrodes, [31] photovoltaic cells [32] electroosmosis, [33] or in biochemistry and nanomedicine for drug delivery, therapeutics, diagnostics, and bioimaging. [34][35][36][37][38] At the same time, experimental findings pointed out cytotoxic features of some metal nanoparticles. [39,40] Several studies suggested that metal nanoparticles interact with cell membranes in a complex way, [41][42][43][44] governed by electrochemical potentials and ion distributions around the membrane and a nanoparticle. These achievements emphasize the importance of a deeper theoretical understanding of the interface between a nanoparticle and the solvent, which acts as a determining factor for many properties of the nanoparticle and its complexes in aqueous environments. [45] On a simplified level of theoretical description, a basic elucidation involves an implicit-solvent treatment of the electrostatic double-layer problem, adopted from the well-established framework of colloidal science. Now of course, the high dielectric interior inverses the role of the image charges as compared to the case of low-dielectric colloidal particles, which thus become attractive and can trigger completely new physics. The attraction between the metal nanoparticle surface and ions can lead to their accumulation and adsorption and thus to a build-up of an electric double layer surrounding the particle, which crucially impacts the colloidal stability [46][47][48][49][50] and interactions with other molecules. This phenomenon is also of great importance in the catalysis by metal nanoparticles in liquid phase. [30,51] The reaction rates of surface-catalyzed bimolecular reactions depend on the concentration at the nanoparticle surface of both reactants, [52,53] which typically have asymmetric properties (charge, specific adsorption, [30] etc.). In this work, we employ theoretical approaches established in the colloidal electrostatics framework, and apply them to less investigated systems of neutral polarizable nanoparticles in different electrolytes. We corroborate the theoretical outcomes by Monte Carlo (MC) simulations, which enable us to assess their regimes of applicability. We show how different kinds of asymmetries between ions result 1. Schematic description of a neutral polarizable nanoparticle of radius a (yellow circle) immersed in an electrolyte solution of cations and anions (red and blue circles). The ions have a radius r0, which specifies the closest-approach distance to other ions and to the nanoparticle. The solvent is treated as a background continuum with the relative permittivity that is much smaller than the permitivity of the nanoparticle. In the MC simulations, the whole system is enclosed in a spherical simulation box of radius R with a reflecting boundary condition.
in non-trivial phenomena, such as charge separation and a build-up of net electrostatic potential and effective surface charge.

II. MODEL AND METHODS
We consider a metal nanoparticle as a neutral sphere with a radius a and a relative permittivity that is much larger than the permittivity of the surrounding electrolyte solution . The electrolyte comprises a mixture of cations with valency q + and anions with valency q − with bulk concentrations n (+) 0 and n (−) 0 , respectively, by which the electroneutrality condition, q + n (+) 0 +q − n (−) 0 = 0, has to be fulfilled. We express all distances by the Bjerrum length at ambient temperature, defined as B = e 2 0 /(4πεε 0 k B T ) (in water at room temperature, the value is B = 0.72 nm), where k B is the Boltzmann constant and T the absolute temperature. The ions are treated as spherical charges with the radius r 0 = 0.2λ B , which specifies the closest approach to other ions as well as to the nanoparticle surface (see Fig. 1). The presence of the dielectric inhomogeneity across the boundary of the sphere influences the electrostatic potential, which can be thus described by the Green's function connecting two points r and r out-side the sphere as u(r, r ) = u 0 (r, r ) + u im (r, r ) Here, u 0 is the direct standard Coulomb kernel in the absence of the dielectric inhomogeneities, and u im is the "image correction" term (sometimes referred to also as the "reaction field") due to the dielectric jump, just as in the case of a planar discontinuity. It turns out that for spherical dielectrics in the limit there is an elegant analytical solution for the electrostatic potential, which is obtained with the help of "image charges". [11,54] Namely, the electrostatic potential of a point charge q at distance r from the center of the sphere is the same as if there were two additional "image" charges instead of the sphere: one in the center of the sphere with the charge q(a/r ) and the other one with the charge −q(a/r ) dislocated by a 2 /r from the center on the line towards the real charge, With this exact Green's function at hand, we now turn to investigate the behavior of ions in the proximity of neutral metal nanoparticles in terms of analytical theories as well as MC simulations.

A. Theoretical approaches
The most common theoretical framework for treatment of electrostatics problems is the Poisson-Boltzmann (PB) equation, based on mean-field premises. [2,3] As such, it cannot describe any image-charge effects on its own. Therefore, for our setting, the PB equation yields a trivial result of non-perturbed ionic distributions around the nanoparticle. In order to account for the polarizability effects, we follow the original ideas of Onsager and Samaras [55]: We first calculate the self-energy (i.e., the potential of mean force) of an ion near the dielectric boundary and then combine it with the PB equation.

Onsager-Samaras self-energy
The simplest theoretical treatment to calculate the self-energy of the above introduced model would be to ignore interactions between ions and considering only the image-charge attraction of an ion with the nanoparticle as dictated by Eq. (3). The self-energy of the ion in this approach is given simply by the interaction potential of the ion with its own images, that is (1/2)e 2 0 u im (r, r). But due to the screening action of the "ionic atmosphere" caused by surrounding ions, the image force is considerable only within distances of the order of the Debye length from the surface, defined in terms of the screening coefficient κ (i.e., the inverse Debye length) as where the sum runs over all ion species. In order to heal the impairments stemming from the surrounding ions, Onsager and Samaras [55] proposed a "screening coefficient" to the image charge in the form of exp(−2κz), where z is the distance of the charge from the dielectric plane. The factor of 2 in the exponent arises because the total "actionreaction" screening distance from the ion to the surface and back to the ion is 2z. Besides, the distance 2z corresponds also to the distance between the ion and its virtual image charge in the planar geometry. Note that Onsager and Samaras originally proposed the correction for the planar geometry as a first approximation in order to simplify the laborious calculations by Wagner [56] who used a spatially varying screening length. Their primary aim was to compute the excess surface tension of electrolyte solutions by integrating the Gibbs adsorption equation. [55] In our first two approaches, we adopt the screening coefficient of Onsager and Samaras to derive an approximate image self-energy of a monovalent ion near a metal sphere. Yet, in the spherical geometry, we have at least two possibilities of adapting the screening distance 2z. In the first approach, we assume twice the distance between the ion and the sphere surface, 2(r − a), which gives the self-energy In the combination with the spherical non-screened image interaction u im given by Eq. (3), this yields a simple analytical expression (rescaled by the thermal energy β −1 = k B T ), which we will refer to as "Onsager-Samaras" (OS) image-charge interaction This interaction can be also seen as the adsorption potential of a monovalent ion to the metal nanoparticle.
In our second approach, we consider the screening distance as the separation between the ion and its images. In this case, each of the two induced image charges is screened by its own screening coefficient, which is exp(−κr) for the first and exp[−κr(r − a 2 /r)] for the second image term in Eq. (3). The resulting Onsager-Samaras* (OS*) expression of this approach (which we denote with an asterisk) is then (7) Note that both expressions, OS and OS*, have not been self-consistently derived but obtained by an adhoc "stitching" together the effects of dielectric discontinuity and ionic screening, and are therefore not exact. Consequently, it is also not a priori clear, which of the two approaches yields more accurate results.

Debye-Hückel self-energy
In our third approach, we base the image selfenergy on the exact Green's function u DH (r, r ) of the Debye-Hückel (DH) equation in the presence of a metal sphere, [50,57] viz.
The Green's function simultaneously accounts for dielectric and screening discontinuities at the surface of the metal sphere. The derivation details are provided in Appendix. The final result for the "DH image self-energy" of a monovalent ion reads Here, the primes denote derivatives of the spherical modified Bessel functions of the first and second kind, which are defined as where I l+1/2 (x) and K l+1/2 (x) are the conventional modified Bessel functions of the first and second kind, respectively.
In the limit of infinitely large radius a → ∞ (i.e., planar metal wall), both Eq. (6) and (9) simplify to e −2κz (12) where z is the distance of the ion from the wall. Exactly the same expression but with the opposite sign applies in the case of an ion near a planar wall with much lower dielectric interior than the electrolyte solution (ε ε). [58] 3

. Boltzmann distribution
In cases when cations and anions have symmetric properties, no electrostatic potential is generated, and their distribution around the nanoparticle is solely governed by the image self-energy. In a thermodynamic equilibrium, we therefore expect the ion densities to follow the Boltzmann distribution. Using the OS [Eq. (6)] and OS* [Eq. (7)] self-energies, this leads respectively to (13) and which we will term as "Boltzmann-Onsager-Samaras" approximations (B-OS and B-OS*, respectively). Similarly, using the DH form (9), gives us which we term the "Boltzmann-Debye-Hückel" (B-DH) approximation. Here, q i is the valency and n (i) 0 the bulk concentration of species i.

Modified Poisson-Boltzmann
If cations and anions redistribute dissimilarly around the nanoparticle, the resulting charge separation can lead to a net electrostatic potential, and thus Eqs. (13) and (15) become inaccurate. As already mentioned, the standard PB equation, which relates electrostatic potential and charge distributions, lacks the image self-energy term. A simple heuristic "remedy" to account for the image effects is to insert by hand the self-energy correction into the Boltzmann factor, thus leading to a modified Poisson-Boltzmann equation [59,60] in the form Here, the summation runs over all ion species i. The self-energy term is in principle given either by Eq. (6), (7), or (9). In our analyses, however, we will limit ourselves only to the DH-based form, Eq. (9). Once the potential φ(r) is known, the ion densities can be evaluated as which we term as the "Poisson-Boltzmann-Debye-Hückel" (PB-DH) approach in this paper. In the case of a symmetry between cations and anions, the electrostatic potential vanishes, φ = 0, and Eq. (17) reduces to Eq. (15). Note again, that the obtained expressions, Eqs. (13)(14)(15)(16)(17), cannot be considered as mean-field results, because they do not follow from the PB equation. Even though many studies [61][62][63][64][65] generalized the seminal work of Onsager and Samaras, it nevertheless remains widely misinterpreted what is the actual theoretical framework of their approach. In fact, these results extend beyond the mean-field level and can be deduced from the thermodynamic fluctuations of the instantaneous electric fields around the PB solution. [66] Alternatively, the PB-DH equation can be derived from a self-consistent variational analysis [67][68][69][70] by setting by hand the screening length κ to be location independent.

B. Monte Carlo simulations
In order to provide "exact" solutions to the introduced model, we perform MC simulations in the canonical NVT ensemble using the standard Metropolis algorithm. [71] The system with mobile ions is enclosed in a spherical simulation box with an outer radius R, containing N + cations and N − anions with valencies q + and q − , respectively, such that their amounts fulfill the electroneutrality condition N + q + + N − q − = 0, see Fig. 1. A reflecting boundary condition is applied to the external box boundary. As opposed to periodic boundary conditions, this treatment significantly simplifies the implementation and increases the performance of the simulations (as no Ewald summation is needed), whereas it distorts ionic distributions near the outer boundary.
In all simulations, the radius of the spherical box is set to R = 17 λ B , which is significantly larger than the largest Debye length of κ −1 ≈ 9 λ B in the study. This guarantees that the outer boundary does not impact the ionic behavior near the nanoparticle.

A. Symmetric case
We start our theoretical analysis by first considering symmetric electrolytes (q + = −q − ). In Fig. 2 we plot ion profiles for 1:1 (top) and 2:2 (bottom) cases at three different bulk concentrations (from left to right: 2.2, 22, and 220 mM) at a nanoparticle of size a = B as predicted by all three theoretical approaches Eqs. (13)(14)(15) and MC simulations. Since in this case the cations and anions have symmetric properties, their density distributions are equivalent. As seen from the plots, ions are considerably attracted to the metal nanoparticle surface due to attractive image charge interactions. The density peaks right at the surface vicinity (at r = a + r 0 ) and is therefore highly sensitive to the minimum approach distance of an ion to the surface, or to be more precise, to the dielectric boundary, in our model determined by the ion radius r 0 . From the OS equation (6), it can be easily appraised that the attractive adsorption energy at the surface (i.e., at r = a + r 0 ), scales as w 0 ∼ −1/r 0 , thus making it very sensitive to the choice of r 0 . Nevertheless, we will keep the value fixed at r 0 = 0.2 B for all further results in this study. For very low ion concentrations (2.2 mM), the screening length is considerably large (κ −1 = 9 B for the monovalent and 4 B for the divalent case), such that the interaction near the surface is predominantly governed by the unscreened part of the image charge interaction. Moreover, in the limit of vanishing salt concentration, all three theories become equivalent and exact. In the cases shown in the figure, all the theories agree very well even for salt concentrations up to 220 mM. The size of the nanoparticle is another important parameter that determines the strength of the image attraction. To demonstrate this effect, we plot in Fig. 3 the density profiles for a monovalent 1:1 electrolyte at 220 mM for different radii a of the nanoparticle. With an increasing size, the densities at the surface get higher. Larger metal nanoparticles have namely higher polarizability, thus attracting the ions more efficiently. From the plot it can  also be observed that all three theories are becoming equivalent as the particle size increases. Figure 3d shows a normalized ion density at the nanoparticle surface (i.e., at r = a + r 0 ) as a function of its radius a. In the limiting case of vanishing particle (a → 0), clearly, the polarizability vanishes and the density becomes bulk-like. The density increases with the radius and saturates at the limit of a planar wall (indicated by arrows), where the interaction is given by Eq. (12). In the limit of vanishing ionic strength (blue solid line), both theories become exact, since in that case the ion-ion interactions become rare and negligible. For higher concentrations (220 mM), the particle of size of a ∼ 2 B already nearly reaches the planar-wall limit. Here, it can also be noted that the theories (except B-OS* in the case of small particles) tend to slightly underestimate the densities near the surface compared with MC simulations. This can be of course attributed to several effects. One of them might be the absence of screening in the ion-free layer of the width r 0 around the surface, which has been for instance discussed by Levin and Mena. [72] Coming to the question of which of the three the-oretical approaches is the most accurate: It is of course expected that B-DH should predict more accurate results than either B-OS or B-OS*, because it properly takes into account the spherical geometry of the problem on the DH level. As can be seen from Figs. 2 and 3, the results of both approximate theories are very close to the results of B-DH. Interestingly, B-OS seems to consistently yield a bit lower densities than B-DH, meaning that it underestimates the overall attraction of the ion to the metal sphere. On the contrary, B-OS* predicts consistently slightly larger results than B-DH. It seems that for small spheres, the B-OS* performs slightly better than B-OS. However, this cannot be claimed for larger spheres, as shown in Fig. 3c, where both B-OS and B-OS* are approximately equally off, yet in opposite directions. However, the advantage of the approximate OS and OS* expressions is their much simpler mathematical form than B-DH.

B. Asymmetric case: specific adsorption
Continuum theoretical descriptions based on the dielectric approximation generally treat ions as equivalent point charges and neglect the nonelectrostatic interactions between ions and the particle surface, which occur in realistic systems. [73] The origin of these ion-specific interactions is still the subject of vivid debate, but in recent years, it has become well established that they are mainly influenced by three parameters: ion-surface, ion-water, as well as water-surface interactions, namely hydrophilicity and hydrophobicity. [74][75][76][77] Different ions are expected to bind to nanoparticle surfaces with different affinities, which typically follow the Hofmeister series. [77][78][79] Binding of ions to the nanoparticle significantly influences their surface charge and the surface potential, which are crucial for the stability of colloidal suspensions based on electrostatic repulsion. [48,49,79] On the continuum-level description, the specific effects can be phenomenologically incorporated via various approximate approaches. In the simplest approximation, the specifically adsorbed ions in the Stern layer close to the surface can be, for instance, treated as a fixed pre-determined surface charge, which is a concept adopted in many theoretical approaches. The main shortcoming of this approximation is that it neglects the dependence of the adsorbed amount of ions on the bulk concentration. Furthermore, it also neglects the influence of surface polarizability and ion correlations. Another approach, which we will adopt here, is to assume an additional attractive potential U s (r) between the ions and the nanoparticle. For simplicity, we use a square-well potential of depth ∆U = −2k B T and the range of r s = 0.3 B from the effective nanoparticle surface, as presented in Fig. 4. In order to introduce an asymmetry in our system, we apply this potential only to cations, while we assume no specific interactions for anions. We plug the potential U s (r) into the Boltzmann factor of Eqs. (15) and (17). This break of the symmetry, assumed in the previous section, has far-reaching consequences as we will see in the following. The ion distributions, shown in Fig. 5, now exhibit a distinct accumulation of cations due to the specific adsorption potential. Notably, the simple Boltzmann-based approach B-DH [Eq. (15)] already captures the densities sufficiently well at low concentrations, since the generated electrostatic potential has negligible influence on ions. But as we increase the concentration, the relative cation density n (+) (r)/n (+) 0 near the surface starts to decrease and anion density slightly to increase. Namely, the potential generated by the adsorbed cations is hindering further accumulation of cations. This behavior is well captured by PB-DH [Eq. (17)], whereas the simple B-DH starts breaking down. A crucial difference between PB-DH and B-DH shows up when zooming in to the far-field region [panel (d)] that extends beyond the specific adsorption potential. There, the ion distributions are considerably influenced by the generated potential. As can be noted, the anion density is higher than the cationic, since anions have to compensate the accumulated positive charge at the surface. An interesting comparison can be made when considering only a PB equation with the specific adsorption but without the image-charge self- energy, The resulting densities are shown in Fig. 5d by dashdotted lines and are in the vicinity of the surface expectedly flatter and lower than the other results due to the missing image-charge attraction. Nevertheless, the PB provides reasonable agreement for the ion densities for distances beyond the specific potential. By integrating the density profiles, we obtain the cumulative charge Z(r) contained within a sphere of radius r around the nanoparticle, which can be then used to evaluate the electrostatic potential (e.g., in MC simulations) as The electrostatic potential generated due to the specific adsorption is shown in Fig. 6a for 22 and 220 mM of a 1:1 electrolyte. The results of PB-DH compare excellent to the MC results. On the other hand, PB that neglects the polarization [Eq. (18)] yields a bit smaller potentials. It is interesting to examine, how does the surface potential, defined as the electrostatic potential at the nanoparticle surface, φ 0 = φ(a), evolve with the ionic strength. As shown in Fig. 6b, φ 0 first linearly rises with concentration, but the rise is becoming gradually weaker for higher concentrations. This slow-down can be attributed to higher repulsion due to accumulated ions and to more effective screening (larger κ) at higher concentrations. The potential stemming from the B-DH approximation (15) can be estimated via the cumulative charge integration, Eq. (20), as is also done for MC data. Since the B-DH approach neglects the potential, which is especially important for far-field behavior, its predictions are severely off compared with other approaches. The B-DH is therefore in this case only a useful predictor for local ion densities, but it fails for far-field. As is well established in colloid science, we expect the generated potential φ(r) to follow a well-known DH law in the far-field, where Z eff is the effective charge (normalized by the unit charge e 0 ) of the nanoparticle. Indeed, by fitting Eq. (21) with Z eff and κ as fitting parameters to the electrostatic potentials at large distances (shown in Fig. 6c for MC data), we obtain very good agreement. The effective charge hence arises as a result of charge separation around an otherwise neutral particle. Figure 6d further demonstrates that the effective charge rises almost linearly with the salt concentration. Both PB-based approaches PB-DH and PB predict very good results (comparable to MC) at low salinities, but tend to underestimate (PB slightly more) the values at higher concentrations. On the other hand, estimating Z eff from the B-DH approach is not possible, since the accumulated charge is effectively not screened by the electrolyte and the evaluated potential in this theory does not follow the DH form of Eq. (21). Instead, B-DH predicts a saturation of the cumulative charge Z(r) to a non-zero value, whereas it is realistically expected to vanish at r → ∞, as is the case for the other two theories and MC simulations. Even though, this is due to a deficiency of the B-DH approach, the value reflects the charge accumulation right at the surface, where B-DH performs reasonably well. Therefore, it is interesting to compare the total cumulative charge Z tot = Z(r → ∞) from B-DH to Z eff from other approaches. Indeed, in Fig. 6d we see that Z tot compares remarkably to Z eff .

C. Asymmetric case: valency
From specific adsorption we now turn our attention to a different kind of asymmetry, the asymmetry that stems from different ion valencies in an electrolyte. According to Eqs. (13)(14)(15), the self-image attraction of an ion to a metal sphere exhibits a square dependence on its charge, ∼ q 2 . Consequently, in cases of asymmetric electrolytes, such as 2:1 or 3:1, this valency dependence engenders strong differences in adsorption between both ion species.
As opposed to specific adsorption, where the particle polarizablity is only an accompanying effect to asymmetric adsorption and related phenomena, it is the main agent for similar phenomena in the case of asymmteric ion valencies. As before, we first look into the ion distributions, which are shown in Fig. 7 for asymmetric 2:1 and 3:1 electrolytes, and compare the theoretical approaches B-DH and PB-DH with MC simulations. As in the case of the ion-specific adsorption, the theories yield better results at low salt concentrations. At higher concentrations, they perform worse due to delicate ion-ion interactions, in particular for higher asymmetry (i.e., 3:1). This theoretical break-down is not unexpected, since multivalent ions are known for significant correlation effects, not accounted for on a mean-field level, a feature that is well established in the double layer literature. [80][81][82][83][84][85][86] As such, Fig. 7 demonstrates a dramatic influence of the valency on the local densities of ions. The relative ionic density at the surface, n(a + r 0 )/n 0 , scales namely as ∼ exp(const. × q 2 ), which for low ionic strengths leads "only" to around 2-fold enrichment of monovalent ions in our system (Fig. 2), 16-fold (∼ 2 4 ) of divalent, and an enormous 512-fold (∼ 2 9 ) enrichment of trivalent ions compared to bulk. This implies high ability of metal particles to take-up multivalent ions from a solution. Cases of highly asymmetric electrolytes are very relevant also in catalytic science, where one of well-studied benchmark "model reactions" involves the reduction of trivalent hexacyanoferrate (III) ions by monovalent borohydride ions catalyzed by metal nanoparticles. [30,87] The local density of the reactant at the surface is one of the governing factors that determines the reaction rate.
[52] In Fig. 8a and b we plot the electrostatic potentials generated by asymmetric electrolytes. While PB-DH gives satisfactory agreement at 22 mM of 2:1 salt, it becomes poorer at 220 mM, where the deviation reaches a factor of 2. The situation significantly worsens for 3:1 case. The surface potential φ 0 as a function of concentration is plotted in panel (c). The theoretical prediction, which is now only qualitative, predicts non-monotonic behavior. The surface potential first rapidly rises with concentration due to increased adsorption of ions. At larger concentrations, the rise of the adsorption slows down with increasing concentration due to electrostatic repulsion of already adsorbed ions. Additionally, increasing the salt concentration increases also the screening of the electrolyte, which eventually leads to a drop in the surface potential at high concentrations. Whereas PB-DH yields satisfactory agreement for the 2:1 case (deviating by factor of 2 from MC at large concentrations), it fails considerably for the 3:1 case. As predicted by the MC simulations, a 3:1 electrolyte creates approximately 20 mV of surface potential in the range of 20-220 mM. This is comparable to the specific-adsorption model discussed in the previous section. We now fit the DH theory [Eq. (21)] to the longdistance potential, which gives us the effective charge Z eff , shown in (d). Contrary to the specificadsorption model in the previous section, the effective charge in this case is notably a non-linear function of concentration. It first shows a rapid increase with the concentration that turns into a more gradual trend at higher concentrations. Consistently with the results for φ 0 in (c), the PB-DH theory underestimates the effective charge. Similarly as in the previous section, the total cumulative charge Z tot from the B-DH approach is very similar to Z eff from MC and PB-DH, with an exception for high concentrations of the 3:1 electrolyte. The last plot is revealing an immense influence of the valency asymmetry on the effective charge. According to the MC result, a neutral nanoparticle of a radius 1 B gains an effective charge of around 1 e 0 at 220 mM of 2:1 electrolyte, and an impressive 6 e 0 in a 3:1 electrolyte of the same concentration. Here we note that the expected effective charge scales with an increasing nanoparticle size faster than its surface, since, as we have seen in Fig. 3, larger particles adsorb ions more effectively due to their higher polarizability. In the limit of large nanoparticle sizes, we then expect Z eff ∼ a 2 . That means that in the case of a polydisperse solution of particle sizes, larger ones gain significantly larger charges than smaller ones. The presented model points to a practical relevance in the physical chemistry, namely the build-up of an electric double layer even in the absence of surface charge, solely because of the difference in cation and anion concentrations in the surface vicinity. The so-called "zero surface-charge double layer", a concept introduced by theoretical models a few decades ago, [88,89] helped to interpret several experimental facts, such as electrokinetic effects of uncharged colloids. [33,90,91] A charged nanoparticle surface enhances its chemical reactivity and consequently has a strong impact on its growth. [92] In real-ity, metal nanoparticles can also possess an intrinsic charge. Partially because nanoparticles can be contaminated with various compounds from electrolytes and oxidized material. [92,93] On the other hand, some syntheses techniques of gold nanoparticles (e.g., pulsed laser ablation) lead to partial oxidation (3.3-6.6% [94]) of surface atoms, forming a pH-dependent equilibrium of Au-OH/AuO − terminal groups, which thus contribute to the overall negative charge of gold nanoparticles.

IV. CONCLUSIONS
In this study, we revisited a continuum electrostatics problem of image-charge interactions and applied it to a model of a metal nanoparticle, featuring a high dielectric interior and hence high polarizability. We compared the predictions of various theoretical approaches, differing in their mathematical complexity and applicability regimes, with Monte Carlo simulations. Focusing first on the case of symmetric electrolyte, we found very good agreement of the theoretical approaches and MC simulations. Here, the polarizability effects lead to sizable ion accumulation near the nanoparticle surface, which further depends on the ionic strength as well as on the nanoparticle size. In addition, we investigated how an asymmetry in the adsorption affinities for cations and anions influences their distributions. We separately considered two different kinds of asymmetries, in one case stemming from an additional specific adsorption potential to one ionic species, and in the other case stemming from an asymmetric electrolyte (i.e., 2:1 and 3:1). The asymmetries, which give rise to asymmetric distributions of ionic profiles, engender a net electrostatic potential and an effective charge of the nanoparticle. Here, even the most simple approaches that neglect the generated potentials can nevertheless very satisfactorily predict local ion densities (i.e., in the surface vicinity). Of course, at larger distances, where ions tend to neutralize the accumulated charge, it is necessary to invoke a Poisson-Boltzmann description with implemented image-charge corrections. For very high charge asymmetries, such as in a 3:1 electrolyte, the theories face difficulties when compared with the "exact" solutions of MC simulations. The difficulties may be associated with correlation effects between multivalent ions, which are not captured within our theoretical framework. Still, the theories are able to capture the qualitative behavior considerably well and thus help to elucidate basic principles of elec-trostatics of metal nanoparticles in electrolyte solutions. Finally, we need to be aware of various conceptual challenges that occur in such systems containing metal-like particles in aqueous solutions. Due to high ionic adsorption affinities, the surface details become very important. This is in stark contrast to low-dielectric macromolecules, where ions are typically repelled from the surfaces, and therefore their molecular structure becomes less relevant. One of such details is for instance the exact geometry of the nanoparticles, which typically possess a welldefined atomic arrangement (e.g., resembling the face-centered cubic structure, [95]) and seems to be critical for the nanoparticle's activity. [96] A deeper understanding of fine details of metal nanoparticles calls for approaches beyond the idealized continuum model. In this context, in particular atomistic models that take the granularity of the nanoparticle surface and solvent into account are nowadays becoming the focus of sophisticated simulation approaches. [45,97,98]

APPENDIX: GREEN'S FUNCTION NEAR A METAL SPHERE
We decompose the Debye-Hückel Green's function u DH (r, r ) into the direct and indirect part [similarly as the Coulomb Green's function given by Eq. (1)], u DH (r, r ) = u DH 0 (r, r ) + u DH im (r, r ) We now express it in spherical coordinates with the center of the sphere located in the origin of the coordinate system, that is, r(r, θ, ϕ) and r (r , θ , ϕ ). A multipole expansion of the direct part yields a form, [99] u DH 0 (r, r ) = 8κ 4πεε 0 lm i l (κr < )k l (κr > )Y lm (θ , ϕ )Y * lm (θ, ϕ) (23) Here, r < and r > correspond respectively to the smaller and the larger radial value among r and r . The functions Y lm are spherical harmonics, the asterisk denotes the complex conjugate value, and the spherical modified Bessel functions are defined by Eq. (11). The summation in Eq. (23) runs over integer values l = 0, 1, 2 . . . and m = −l, . . . , +l. A general Ansatz for the second term in Eq. (22) is [50,57] u DH im (r, r ) = lm A lm k l (κr)Y * lm (θ, ϕ).
In order to determine coefficients A lm , we apply two boundary conditions. The first condition requires that the potential on the surface of the metal sphere with a radius a is constant, that is, independent of the solid angle Ω ∂u DH (r, r ) ∂Ω | r=a = 0 (25) which leads to k l (κr )Y lm (θ , ϕ ) for l, m = 0.
(26) Note that for l = m = 0, Y 00 = 1/ √ 4π, and consequently ∂Y 00 /∂Ω = 0, thus trivially satisfying Eq. (25). For that reason, A 00 has to be determined by an additional boundary condition. Since the sphere is electrically isolated, it is overall charge neutral. We apply the Gaussian law, E · dS = 0, where we integrate the electric field E = −e 0 ∇u DH over the sphere surface. This provides the second boundary condition ∂u DH ∂r r=a sin θdθdϕ = 0 (27) where we integrate over the entire solid angle. Using the identity Y lm (θ, ϕ) sin θdθdϕ = √ 4π δ l0 δ m0 , which eliminates all the terms but l = m = 0 when performing the integration in Eq. (27), provides the remaining coefficient The self-energy of a monovalent ion then follows as w DH 0 = (1/2)e 0 u DH (r, r), which is Eq. (9) in the main text.

CONFLICTS OF INTEREST
There are no conflicts of interest to declare.