{'fraction_non_alphanumeric': 0.04975896531452087, 'fraction_numerical': 0.02393885949441505, 'mean_word_length': 4.573525557011796, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Highly asymmetric electrolytes in the primitive model: Hypernetted chain solution in arbitrary spatial dimensions

The pair-correlation functions for fluid ionic mixtures in arbitrary spatial dimensions are computed in hypernetted chain (HNC) approximation. In the primitive model, all ions are approximated as non-overlapping hyperspheres with Coulomb interactions. Our spectral HNC solver is based on a Fourier-Bessel transform introduced by Talman [J. Comput. Phys., 29, 35 (1978)], with logarithmically spaced computational grids. Numeric efficiency for arbitrary spatial dimensions is a commonly exploited virtue of this transform method. Here, we highlight another advantage of logarithmic grids, consisting in efficient sampling of pair-correlation functions for highly asymmetric ionic mixtures. For three-dimensional fluids, ion size- and charge-ratios larger than one thousand can be treated, corresponding to hitherto computationally not accessed micrometer-sized colloidal spheres in 1-1 electrolyte. Effective colloidal charge numbers are extracted from our primitive model results. For moderately large ion size- and charge-asymmetries, we present Molecular Dynamics simulation results that agree well with the approximate HNC pair correlations.

The pair-correlation functions for fluid ionic mixtures in arbitrary spatial dimensions are computed in hypernetted chain (HNC) approximation. In the primitive model, all ions are approximated as non-overlapping hyperspheres with Coulomb interactions. Our spectral HNC solver is based on a Fourier-Bessel transform introduced by Talman [J. Comput. Phys., 29, 35 (1978)], with logarithmically spaced computational grids. Numeric efficiency for arbitrary spatial dimensions is a commonly exploited virtue of this transform method. Here, we highlight another advantage of logarithmic grids, consisting in efficient sampling of pair-correlation functions for highly asymmetric ionic mixtures. For three-dimensional fluids, ion size-and charge-ratios larger than one thousand can be treated, corresponding to hitherto computationally not accessed micrometer-sized colloidal spheres in 1-1 electrolyte. Effective colloidal charge numbers are extracted from our primitive model results. For moderately large ion size-and charge-asymmetries, we present Molecular Dynamics simulation results that agree well with the approximate HNC pair correlations.

I. INTRODUCTION
Most of the essential features of polyelectrolyte solutions can be efficiently modeled by a combination of charges and excluded volume, see [1][2][3][4] for reviews. In the so-called primitive model (PM), all specific properties of the solvent are neglected except for its dielectric constant. Research over the past decades has shown that some of the basic properties of polyelectrolytes (like screening and Coulomb association) are contained in the asymmetric PM of electrolytes which enables wide applications to charged colloidal suspensions, micelles and globular proteins. However, for high asymmetry in charge and size between the microions and macroions, as occurring in suspensions of charged colloidal particles, the PM is not easy to solve numerically in general. For example, structural correlations in the PM were obtained by numerically expensive computer simulations only up to charge and diameter asymmetries of about 1:100 [5][6][7][8], corresponding to the micellar rather than the colloidal regime.
In the present paper, liquid structure is computed by solving integral equations based on the Ornstein-Zernike (OZ) equation [9]. This approach requires an approximative closure for an explicit solution. A rather simple but successful closure is the hypernetted chain (HNC) scheme [10], which has been proven to be a realistic approximation for mixtures of charged particles. The numerical solution methods presented here can be straightforwardly generalized to more sophisticated, thermodynamically partially self-consistent OZ closure relations [11][12][13] like, e.g., the one proposed by Zerah and Hansen [14]. For the sake of simplicity, and since it was shown * mheinen@thphy.uni-duesseldorf.de that enforcing thermodynamic self-consistency leads to a weak accuracy improvement only [15], in the present work we restrict ourselves to the HNC approximation. A variety of liquid integral equation [15][16][17][18][19][20][21] or density functional studies [22][23][24] of the PM have been reported. The HNC equations have been solved for the asymmetric PM by Léger and Levesque [25], in case of non-zero macroion number densities for size asymmetries between micro-and macroions as high as 1:80 and charge asymmetries ranging up to 1:450.
The HNC scheme can be formulated in any spatial dimension d > 0. While d = 3 is the standard threedimensional situation, it is important to note that also two-dimensional and one-dimensional fluids occur in experiments with strong confinement between glass plates or at interfaces, or in one-dimensional channels [26]. Dimensions higher than d = 3 have no immediate realization. However, they play an important role in constructing and testing theories and are also helpful to find suitable mean-field-like approximations in lower dimensions [27][28][29][30]. Hence, there is a need to study charged systems also in d > 3. We formulate the solution method in the present paper for arbitrary d > 0, with explicit data presented for spatial dimensions d = 1, 2, 3, 4, 5, and 6.
Here we pick up a strategy of solving the PM-HNC equations that, for d = 3, has first been employed by Rossky and Friedman [31]. The key idea lies in using computational grids with logarithmic spacing in coordinate-and wavenumber space. For even number of dimensions d, in particular for d = 2, the use of logarithmic grids in a spectral OZ equation solver emerges as a necessary consequence of mapping the occurring Fourier-Bessel transforms to numerically efficient fast Fourier transform (FFT) methods [32][33][34][35][36]. Hence, logarithmic grids appear quite naturally in d = 2 liquid integral equation studies like, e.g., Refs. [37,38]. For odd number of dimensions, where FFT methods can be applied directly on uniformly spaced grids, using logarithmic grids is a less obvious approach which has nevertheless been employed in some d = 3 studies [39][40][41]. In none of these studies, however, extensive use has been made of an important virtue of the logarithmic grids, which is highlighted in our present paper: The simultaneous dense distribution of grid points at very different length scales renders logarithmic grids ideal for the discretization of pair-correlation functions in a highly asymmetric PM, at moderate numerical expense. Here we are presenting results for ion size-and charge-asymmetries as high as 1:1000 (both asymmetries simultaneously reached, for non-dilute suspensions), which represents well the colloidal regime. To our knowledge, no liquid integral equation or computer simulation studies have been published so far, where asymmetries of this magnitude have been reached.
We investigate the accuracy of the HNC paircorrelation functions by comparing to results of numerically expensive Molecular Dynamics (MD) simulations, for ion charge-and diameter asymmetries up to 1:500 and 1:250, respectively. The HNC results are found to be in overall good agreement with the MD results, except for a somewhat underestimated principal peak in the macroion radial distribution function at high macroion charge numbers.
Our results can be used to extract colloidal effective interaction potentials, including non-saturated effective charge numbers, at high numerical efficiency. As opposed to other approximate theories for colloidal effective charge numbers like the cell model [42,43] or the renormalized jellium model [44,45], pair correlations among all ionic species are treated on equal footing in our method, with the HNC entering as the only approximation. In a future study, the method described here could be augmented to include colloidal surface chemistry, described by a mass action balance that takes into account the local variations in the pH value near the colloidal surfaces. This would allow for parameter free ab initio calculation of pair correlations in colloidal suspensions with reactive electrolyte like, e.g., suspensions of silica spheres in N aOH with a reentrant fluid-crystal-fluid phase diagram, reported in Refs. [46,47].

II. METHODS
A. The d-dimensional primitive model Let n i , with 1 ≤ i ≤ m, denote the number density of ions of species i in an m-component fluid mixture of hyperspherical particles in an arbitrary positive number d of spatial dimensions. Then, n = i n i is the total particle number density and χ i = n i /n is the mole fraction of species i. Let Z i e denote the electric charge of a particle of species i, where e is the elementary charge. In order to prevent singular attractions between oppositely charged particles, ions of species i possess a hard-core diameter σ i . Hence, the particles of species i occupy a fraction (1) of the system hypervolume, where V (d) = π d/2 /Γ(d/2 + 1) is the d-dimensional unit hypersphere volume, Γ(x) denoting the Gamma function. All ions are assumed neutrally buoyant in an infinite structureless solvent, which is fully described by the solvent dielectric constant ǫ in the PM applied here. We express all real-space functions in units of the dimensionless particle-center to particlecenter distance x = rn 1/d . The functions u ij (x) = V ij (x)/(k B T ) are the pairpotentials of direct interaction, V ij (x), between ions of species i and j, divided by the Boltzmann constant k B and the absolute temperature T . The dimensionless pairpotentials can be decomposed as with pairwise additive hard-core diameters σ ij = 1/2(σ i + σ j ), and long-ranged (Coulomb) parts with coupling constants Γ ij ∝ Z i Z j . In d = 3 dimensions, connection to experimentally accessible systems can be made by choosing the coupling constant Γ ij in Eq. (4) as Γ ij = L B n 1/3 Z i Z j , involving the solvent-specific Bjerrum length L B = e 2 /(ǫk B T ) in Gaussian units. In the present paper, we investigate only such systems that obey the Berthelot mixing rule Γ 2 ij = Γ ii Γ jj [48]. One-component systems, m = 1, with pair-potential according to Eqs. (2)-(4), are commonly referred to as (d-dimensional) one component plasmas (OCPs). Global electroneutrality of an OCP implies the presence of a homogeneous background charge density that does not couple to the distribution of the correlated ions (like, e.g., an electron plasma at sufficiently high temperature). For systems with m > 1 components, global electroneutrality without a neutralizing charge background is enforced in all cases studied here, by requiring that i n i Z i = 0.

B. Hypernetted chain scheme
We compute the ion pair-correlations in the PM described in the previous section, by numerically solving the OZ equations [9] in combination with the approximate HNC closure [10][11][12]. In an isotropic, homogeneous fluid mixture, the coupled OZ equations may be written as (5) In Eqs. (5) and the rest of this paper we adhere to the Einstein summation convention. Equations (5) can be regarded as the definitions of the partial direct correlation functions c ij (x) in terms of the continuous partial indirect correlation functions . The latter identity comprises the total correlation functions h ij (x) and the partial radial distribution functions (rdf's), g ij (x), which are the conditional probabilities of finding a particle of species j at a dimensionless center-to-center distance x from a particle of species i [9].
An isotropic function f can be Fourier-transformed in d dimensions as with J n (x) denoting the Bessel function of the first kind and order n. Employing the convolution theorem, the OZ equations are Fourier-transformed into the space of dimensionless wavenumbers y, where they read γ ij (y) = χ kcik (y)c kj (y) + χ kcik (y)γ kj (y).
The OZ equations need to be supplemented by an appropriate closure relation. The HNC closure, which is known to be a good approximation for the PM [15], reads Numerical solution of the set of Eqs. (5) and (9) in combination with the long-ranged potentials in Eqs. (2)-(4) requires splitting the analytically known long-ranged asymptotic parts ∓u (l) ij (x) off the direct and indirect correlation functions as [9,25,49] and The so-defined functions c  ij (x) are considerably shorter in range than c ij (x) and γ ij (x). In terms of the short-ranged correlation functions and long-ranged potential parts, the OZ equations in wavenumber space can be written as the set of coupled algebraic equations with Kronecker-delta δ ij . The HNC closure in terms of the short-ranged correlation functions and the hard-core diameters is (13) with unit step function Θ(x).
The Fourier transform of the Coulombic part u (l) ij (x) of the potential is [50] u (l) where A(d) = 2π d/2 /Γ(d/2) denotes the surface of the d-dimensional unit hypersphere. We solve the closed set of Eqs. (12)- (14) by the numeric methods described in the following two sections. Our results are presented in form of the partial rdf's g ij (x), and partial static structure factors, S ij (y) = δ ij + √ χ i χ jhij (y) [9,12].

C. Numerical algorithm
In order to solve the set of Eqs. (12)-(14), we employ a generalized version of the numerically robust, quickly convergent solution method introduced by Ng in the appendix of Ref. [49]. The method shares great similarities with the direct inversion of iterative subspace (DIIS) method developed by Pulay [51,52], which is commonly used in the solution of quantum mechanical (Hartree-Fock) density functional problems (see Ref. [53] for a detailed analysis of the DIIS method). The DIIS method has been applied also to density functional theory of hard spheres [54], and reference interaction site model HNC equations of liquid water [55].
In conformity with Ng's notation, we formulate a fixed point problem ij (x), and the nonlinear fourth-rank operator A, that depends on all pair potentials u ij (x).
We solve Eq. (15) numerically by executing two nested instruction loops with iteration in- , provided that the iteration with respect to index n 2 converges. To avoid confusion with other kinds of indices or with particle number densities in the following, n 1 and n 2 are exclusively used to label entire operators, matrices or vectors (identified through bold font), and are never used in indexing scalars. Both n 1 and n 2 are always enclosed in curly brackets when used as an index, whereas other lower indices, like the species indices i, j, . . . stand always without brackets. For instance, c (s) {n1,n2} (x) is the m × m array of intermediate solutions for the short-ranged parts of partial direct correlation functions at iteration stage characterized by n 1 and n 2 , and the element of that array with particle-species indices i and j is the function (c (1)-(4). We employ the rules for potential ramp-up, where a near-optimal convergence rate, combined with good numerical stability of the outer loop iteration is achieved by a convergence-adaptive scaling parameter 0 < ε ≤ 1, which increases monotonically as a function of n 1 . After each outer iteration, the growth rate of ε(n 1 ) is increased if the previous inner loop took less than a certain threshold of iterations to converge, and decreased in the opposite case. Our experience shows that numerical stability of the algorithm benefits considerably from the superlinear form of Eq. (17), characterized by the (empirically chosen) exponent 1/10. This can be rationalized by considering the contact value, lim xցσij n 1/d u ij (x), of the pairpotential in Eqs. (1)-(4): The contact value decreases for increasing values of φ i and φ j , and increases for increasing Γ ij . The potential at particle contact has great influence on the strength of the undulations in the pairstructure functions, which, in turn, influence numerical stability. Therefore, it is favorable to choose potential ramp-up rules like Eqs. (16) and (17), where the φ i are increased quicker than the Γ ij .
That a potential ramp-up is necessary at all, is owed to critical dependence of the inner loop convergence on the quality of the inner loop seed c We proceed now with the discussion of the inner iteration loop, where n 1 is kept fixed. Operator A {n1} is defined by as the operator that transforms the input functions arrays, c {n1,n2} (x), for fixed indices n 1 and n 2 , into the corresponding output functions arrays, η {n1,n2} (x), the latter being defined further down the text of this subsection. Hence, A {n max 1 } = A, and equation (15) is equivalent to with function arrays d {n1,n2} (x) defined by In our implementation, iteration in the inner loop is stopped at a finite, n 1 -dependent value, n max 2 (n 1 ) > 1, of the index n 2 , once the convergence criterion with a small tolerance TOL(n 1 ), as specified below, has been fulfilled. In Eq. (21), the norm, f(x) , of an m× m function array f (x), is defined as and the bracket (f , g) denotes the inner product of two m×m function arrays f (x) and g(x) with elements f ij (x) and g ij (x), double dots indicating the contraction with respect to both particle species indices i and j. The interval [x 1 , x 2 ] should be chosen to contain the major structural features of the partial direct correlation functions. To obtain the results presented in this paper, we have used x 1 = min{0.5σ i n 1/d , i = 1 . . . m}, x 2 = 30, and values of TOL that decrease as a function of n 1 , with TOL(0) < 10 −4 and TOL(n max 1 ) < 10 −12 . In our algorithm, {n1,n2} (x) at fixed values of n 1 and n 2 by applying the following four steps: Step 1: {n1,n2} (y) by the fast transform method described in the following section, requiring logarithmically spaced grids in x-and y-space.
Step 2: The coupled OZ Eqs. (12) are solved to obtaiñ γ Step 3: The fast inverse transform method on logarithmic grids is applied to compute the inverse Fourier transform, γ Step 4: The elements of the function arrays η {n1,n2} (x) are calculated as the left-hand-sides of Eqs. (13) (the HNC closure), where the elements of the function arrays γ (s) {n1,n2} (x) from step 3 are entered to the right-handsides of Eqs. (13).
A straightforward way of selecting the input functions, c (s) {n1,n2+1} (x), for the next step of the inner loop, is the Picard-iteration scheme This simple scheme, however, converges only for weak pair-potentials (small values of |Γ ij | and φ i , in our case). Numerical stability of the Picard-iteration scheme can be somewhat improved, at the cost of increasing computational effort, if a fixed mixing parameter 0 < α < 1 is introduced in Eq. (24), which gives the alternative iteration rule Fixed point iterations on basis of Eq. (24) or (25) have been applied in a number of integral equation studies [21,25,31,37,38,56,57], where different strategies have been applied in computing the Fourier transforms, and various closure relations for the OZ-equations have been used, including the HNC closure. The mixing parameter α in Eq. (25) has been empirically determined in most cases.
Despite being numerically more robust than the Picard iteration in Eq. (24), the iteration scheme according to Eq. (25) still fails to converge for large values of |Γ ij | or φ i , especially if the number of components, m, is larger than one. We therefore use a generalized version of the fixed point iteration scheme proposed by Ng [49], which has proven to be numerically much more stable and efficient. Ng's iteration scheme, generalized to mcomponent mixtures and arbitrary number, M ≥ 0, of mixing coefficients, reads For M = 0, Eq. (26) reduces to the Picard iteration in Eq. (24). At every step of the iteration, the Mdimensional mixing coefficient vector α {n1,n2} is determined as the solution of the set of linear equations where the elements, and Solving for α {n1,n2} in Eqs. (27)-(30) is equivalent to solving the minimization problem with respect to α {n1,n2} , and the minimization in expression (31) can be motivated by approximating A {n1} as a locally linear operator [49,58]. It may be counterintuitive, but is worthwhile to note that one should not select M = n 2 , i.e., the maximum possible order at each inner loop iteration. Instead, numerical stability is increased if one chooses M to rise slower than possible, in our case as where the f ij (x) are the Mayer functions [9,57]. This results in a seed c (s) for n 1 = n 2 = 0, where u {n1} (x) denotes an array of pair-potentials (u {n1} ) ij (x) between particles of species i and j, obtained from inserting the potential parameters (Γ {n1} ) ij and (φ {n1} ) i into Eqs.
(1)-(4). To access HNC solutions at higher |Γ ij | and φ i , we construct the seeds for n 1 > 0 from the converged solutions of inner iterations corresponding to smaller values of n 1 . For n 1 = 1, we choose which is a Picard-iteration step in the outer loop. For n 1 > 1, we use the seed with coefficients and that extrapolate linearly on basis of the previous two converged inner iteration solutions. While the linear extrapolation in Eq. (35) is crucial for numerical stability at strong particle correlations, generalizing to quadratic and higher orders of polynomial extrapolation for n 1 > 2 seems to have a weakening effect on numerical stability. Variations of Ng's fixed point iteration method have been successfully applied in various liquid integral equation studies [17,20,49,55,[59][60][61][62][63][64]. General performance figures of the algorithm in Eqs. (26)-(37) are difficult to formulate, as its efficiency depends on the number of particle species and spatial dimensions, and, most importantly, on the pair-potential parameters. All individual HNC solutions presented in the present paper took few minutes or less to be computed on an inexpensive personal computer.
Note that the fixed point iteration scheme in Eqs. (26)-(37) is merely one among a wide variety of solution methods that have been developed. Part of the alternative solution methods have been reported to show superior numerical efficiency, at the cost of a more complicated implementation. Along with Ng's algorithm, Newton-Raphson-like fixed-point iteration schemes, first introduced by Gillan [65], Labík et al. [66], and Zerah [67], are routinely used in integral equation studies of liquids with strong pair-correlations [15,18,19,[68][69][70][71][72][73][74][75][76]. For an elaborate comparison of Ng's and Zerah's methods, including a formulation of the latter method for liquid mixtures, we refer to appendix A of Ref. [58]. Furthermore, highly elaborate Newton-GMRES (Krylov subspace) algorithms have been applied [77], and have been combined with multigrid techniques [78]. Yet another alternative approach is the vector extrapolation method [79]. For the sake of simplicity, in the present study we do not employ the methods laid out in Refs. [65][66][67][77][78][79]. In future studies however, use of such elaborate fixed point solution methods, in combination with the Fourier transform method in Sec. II D, might give access to liquid integral equation solutions of the PM for even larger ion-size and charge asymmetries than accomplished in the present work.

D. Fourier transform on logarithmic grids
In this subsection, we present our numerical method of choice to approximate the forward-and backward Fourier transforms in Eqs. (6) and (7), for functions f (x) and f (y) that are sampled on finite computational grids. The method employed here has been devised in essence by Talman [34], and constitutes a sophisticated version of the so-called quasi-fast Hankel transform method of Siegman [33]. It is based on the use of logarithmic variables, i.e., computational grids of the form with grid index n in the range −⌊N/2⌋ ≤ n ≤ ⌊N/2⌋, and L, N > 0. The use of logarithmic grids has been motivated by the work of Gardner et al. [32]. As demonstrated in Ref. [33], sampling on logarithmic grids allows to rewrite Eqs. (6) and (7) as discrete circular correlations, each of which can be treated by applying two subsequent fast Fourier transforms (FFTs). For a complete, and particularly clear-cut documentation, we refer to the work of Hamilton [35,36], where the method has been named FFTLog. Since we have closely followed Refs. [35,36] in our implementation of the FFTLog transform, we refrain from repeating all details here, and list only the essential expressions instead.
Adapting to the notation of Refs. [35,36], we define the primed sum symbol ′ through with weights w n = 1 for all n, except for w −N/2 = w N/2 = 1/2 if N is even. In a preprocess, preceding the many FFTLog transforms occurring in the fixed point iteration described in the previous section, lookup tables of the grid-specific coefficients are computed. As pointed out in Refs. [34,37], numerical evaluation of Eq. (40) for large arguments of the Gamma functions is considerably simplified by noticing that |Γ(x + iy)/Γ(x − iy)| = 1 for x, y ∈ R. The remaining problem of determining the complex phases of the Gamma functions in Eq. (40) is conveniently solved in resorting to Ref. [80]. An alternative way to determine the complex phases has been devised in Ref. [37], where the problem was tackled using recurrence relations of the Gamma function. The forward FFTLog transform from the space of dimensionless distances x to the space of dimensionless wavenumbers y can be evaluated as [35,36] f (y n ) = 2π and the inverse transform (from y-to x-space) reads [35,36] f with the star denoting the complex conjugate. Being nothing else but discrete, one-dimensional Fourier transforms, all four Eqs. (41a)-(42b) can be solved by numerically efficient FFT algorithms, which are available in a great variety of implementations. Here, we employ a mixed-radix FFT routine [80], imposing no constraints on the number of grid points 2⌊N/2⌋ + 1. For validity of Eqs. (41b) and (42b), it is necessary to fulfill the constraint u −⌊N/2⌋ = u ⌊N/2⌋ in choosing the parameters x 0 , y 0 , N , and L of the grids in Eq. (38). Fulfilling this constraint is equivalent to with an arbitrary integer number z, and Arg[c] denoting the phase of complex number c. In Refs. [35,36], the criterion in Eq. (43) has been named the low-ringing condition. In all calculations with results presented here, we have chosen x 0 = 1, and x 0 y 0 ≈ 1, while fulfilling Eq. (43).
An important virtue of the FFTLog transform in Eqs. (40)-(42b) is its computational efficiency for arbitrary dimensions d. Choosing the number of grid points as an integer power of 2 results in optimal performance of the transforms in Eqs. (41a)-(42b), each requiring O(N log 2 N ) arithmetic operations in that case. We note here that an alternative, numerically efficient method of calculating d-dimensional Fourier-Bessel transforms of the kind of Eqs. (6) and (7) has been used in Refs. [81][82][83][84]. In this alternative approach, which does not require logarithmic grids, the Fourier-Bessel transforms are replaced by a sequence of so-called hat transforms and FFTs.
For even d, in particular for d = 2, numerically less efficient methods for computing the transforms in Eqs. (6) and (7) have been reported [85][86][87], each of which requires O(N 2 ) arithmetic operations. Such numerically sub-optimal O(N 2 ) transforms have been applied in various liquid integral equation studies [21,57,88].
For d = 3 (odd number of dimensions, in general), using the transform in Eqs. (39)- (43) is not obvious since the standard FFT can be directly applied to functions sampled on grids with uniform spacing. Uniformly spaced grids, however, are not ideally suited for sampling the correlation functions of highly asymmetric PM fluids. Ion size-and charge-ratios of the order of 1:1000 in typical colloidal suspensions render it necessary to simultaneously resolve length scales that differ by a factor of more than one thousand, which requires huge numbers of grid points in uniformly spaced grids. For instance, in Ref. [25], 2 18 = 262144 points had to be used to sample correlation functions for ion diameter asymmetry of 1:80 and charge asymmetries up to 1:450 (at nonzero macroion number density), resembling only rather small macroions. Logarithmic grids, on the other hand, are ideally suited to capture the different length scales in (asymmetric) charged sphere systems, as it has been first pointed out by Rossky and Friedman [31]. Logarithmic grids have later been used in d = 3 liquid integral equation studies [39][40][41] but, to our knowledge, no liquid integral equation study has been conducted so far, where the advantages of the the transform in Eqs. (39)- (43) have been exploited in solving pair correlations of an extremely asymmetric PM. Our results presented here, for ion charge-and size-ratios both simultaneously as high as 1:1000, have been obtained using no more than 2 13 = 8192 grid points.

E. Molecular Dynamics Simulations
In this paper, we present MD simulation data for fluids in d = 3 spatial dimensions only. We have simulated globally electroneutral systems with three or four different ion species in a cubic box of edge length B with periodic boundary conditions in all three Cartesian directions. The MD simulation method used here is the same as in Refs. [8,89]. In order to handle the long-ranged Coulomb interactions, the Lekner summation method [90][91][92] is employed.
All parameters for the simulations presented here are listed in Table I, with N i denoting the total number of particles of species i in the simulation box, such that n i = N i /B 3 . All simulations are for a Bjerrum length of L B = 0.701 nm, corresponding to water at room temperature. A representative snapshot of particle positions for the four-component system (rightmost column of Table I) is shown in Fig. 1.
On average, one week of execution time on a 64-bit computer cluster is enough to get fairly good statistics for the macroion-macroion and the macroion-microion correlation functions, for the systems consisting of up to 32,208 charged particles. Achieving comparably good  microion-microion statistics would require even much longer execution times. Therefore, we restrict ourselves here to comparisons of macroion-macroion and macroionmicroion pair correlations obtained from MD simulations and the HNC scheme.

A. Thermodynamical properties and pair-correlations in one to six dimensions
As a first result, the HNC solutions computed in our numerical algorithm are in agreement with the ddimensional local electroneutrality (LEN) conditions (44) where τ i (x), defined as is the isotropic charge density around a test particle of species i. In Eqs. (44) and (45), sgn(x) is the sign function. The LEN condition states that the total charge of ions around a test particle cancels out with the test particle's charge. If a computational grid is chosen that extends to a very large outer radius x ⌊N/2⌋ , we find the LEN condition in Eq. (44)  ij (x) assume values too small to be resolved at machine precision. However, if taking the x → ∞ limit in Eq. (44) is replaced by insertion of an intermediately large value of x ≈ 100, where the oscillations in all g ij (x) have essentially died out, we find Eqs. (44) fulfilled to within good accuracy.
In the special case of m = 1, our algorithm allows to compute pair-correlations of the OCP in arbitrary dimensions, as illustrated in Fig. 2. Note here the good quality of the HNC solutions at small wavenumbers, where the vanishing compressibility of the OCP, with lim y→0 S 11 (y) = 0, is well described. The magnitude of the undulations in the S 11 (y) and g 11 (x) plotted in Fig. 2 is a non-monotonic function of the dimension d, with maximal undulations occurring for d = 3.
In Figs. 3 and 4, we compare the HNC rdf's g ij (r) for three-and four-component primitive models in d = 3 to the results of our MD simulations. Overall good agreement is observed between the HNC and MD results, the most prominent discrepancy being an underestimation (of up to about 20%) of the principal peak height in the HNC macroion-macroion rdf's, occurring at strong macroion correlations. Underestimation of the principal peak heights in the macroion-macroion pair-correlation functions is a known shortcoming of the HNC [9,93,94], which can be tackled by choosing an alternative, thermodynamically partially self-consistent integral equation scheme [11][12][13][14]. An alternative method to improve the accuracy of the HNC consists in using a tailored Ansatz [49,88] for the bridge function [9] at high coupling, which is neglected altogether in the HNC.
In Fig. 5 we display all HNC partial rdf's and static structure factors of three-component, globally electroneutral primitive models in all integer dimensions d One component plasma (m = 1) static structure factors, S11(y), and rdf's, g11(x), for coupling parameter Γ11 = 100. In each case, the hypervolume fraction φ1 has been chosen small enough to ensure a practically vanishing rdf contact value, g11(x = σ11n 1/d ) ≈ 0. Correlation functions for systems in all integer dimensions from d = 1 to d = 6 are plotted. Principal peak positions of S11(y) and g11(x) shift from left to right as d increases.
from 2 to 6. The dimensionless pair-potential parameters are the same for all systems in the figure, with coupling constants Γ 22 = Γ 33 = −Γ 23 = 0.01, Γ 13 = −Γ 12 = 25, and Γ 11 = 625, corresponding, for d = 3, to Z 1 = 250, Z 2 = −1, Z 3 = 1, and L B n 1/3 = 0.1. The hypervolume fraction of species 1 is φ 1 = 0.01 for all systems in Fig. 5, and n 3 = n 1 √ Γ 11 , which means that there is one salt coion per macroion-surface released counterion. A rather small ratio of macroion-to microion diameters, σ 1 /σ 2 = σ 1 /σ 3 = 5, has been chosen for the systems in Fig. 5 since, for larger size asymmetry, the solution in higher dimensions such as d = 6 becomes numerically very slowly convergent or divergent. Note that Fig. 5 features panels with double logarithmic as well as logarithmic-linear axes, and that the axes ranges vary from panel to panel, to exhibit simultaneously the details of the various plotted functions.
In Fig. 5, the g ij (x) and S ij (y) with the most pronounced oscillations around the asymptotic value one are for d = 2. For rising dimension and fixed dimensionless potential parameters, particle packing becomes less efficient and the decay of the Coulomb potentials becomes steeper. Therefore, the undulations in the paircorrelation functions get reduced for rising d until, for d = 6, the undulations have almost completely died out.
B. Application to the colloidal domain Figure 6 features all HNC g ij (x) and S ij (y) for a d = 3, ternary PM that resembles a realistic suspension of colloidal particles in aqueous (L B = 0.70 nm) electrolyte with a low concentration of dissociated N aCl. The diameter of the macroions is taken to be σ 1 = 1µm, corresponding to rather large colloidal particles, while the diameters σ 2 = 0.756 nm and σ 3 = 0.922 nm correspond to hydrated N a+ and Cl− ions, respectively. Here, we have chosen the concentration of Cl− coions (species 3) as n 3 = n 1 |Z 1 |, so that the suspension contains one coion per colloid-surface released counterion, and overall about twice as many counterions as coions, n 2 ≈ 2n 3 . With the assumed macroion charge Z 1 = −750, which is a realistic bare charge for a micron-sized colloidal sphere, this gives a coion concentration of n 3 = 0.71µM . This corresponds to an almost deionized aqueous solvent with a little amount of dissociated salt only (c.f., the number concentration, n = 0.1µM , of the water self-dissociation products H 3 O+ and OH− at neutral pH-value, which is a lower bound for the coion concentration).

C. Effective colloidal interactions
As a dimensionless effective pair potential between particles of the same species a in an m-component fluid mixture, one can define [69] u eff where c eff aa (x) is an effective direct correlation function between particles of species a. The Fourier transform of the latter isc eff aa (y) =h with the total correlation functionh aa (y) =γ aa (y) + c aa (y) =γ  (14). Equations (46) and (47) constitute an inversion of the HNC for species a only, meaning that a one-component fluid of particles with pair-potential βu eff aa (x), solved within the HNC approximation, shows exactly the same paircorrelation functions as component a of the m-component mixture.
Consider now a three-dimensional, ternary ionic liquid mixture of macroionic spheres (species a) with diameter σ a and charge number Z a , monovalent counterions, and monovalent coions. In this case, the repulsive part of the dimensionless Derjaguin-Landau-Verwey-Overbeek (DLVO) effective pair potential between two macroions at a non-overlap center-to-center distance x > σ a n 1/3 can be written as [95] βu DLVO aa (x) = Γ aa e kσ a n 1/3 1 + kσ a n 1/3 2 with Γ aa = L B n 1/3 Z 2 a , as defined further up this text. Equation (48) involves the dimensionless DLVO screening parameter k, which is given by with χ coion denoting the mole fraction of coions. The DLVO potential in Eq. (48) is valid for two macroions in a bath of microions whose distribution can be treated in the Debye-Hückel approximation. It is thus valid only at low macroion concentration, and for L B Z a /σ a 1, i.e., for ionic pair interactions that do not considerably exceed the thermal energy. Under conditions where L B Z a /σ a > 1, the potential in Eq. (48) is nevertheless a good approximation to the effective macroion pairpotential at sufficiently large particle separation, provided that the charge number Z a , entering via Γ aa and Eq. (49), is replaced by an effective charge number Z eff a < Z a [42,44,45,[96][97][98][99][100][101][102][103][104]. The effective charge Z eff a e has to be regarded as the net charge of a colloid dressed with closely associated counterions.
While theories for the saturation value, Z eff a (Z a → ∞), of the effective charge number are available [99,100,105], calculation of the non-saturated effective charge in a PM, in its full dependence on the concentrations and charges of all ionic species, remains a challenging task [42][43][44][45].
We determine the effective charge by replacing Z a with Z eff a in Eqs. (48) and (49), and tuning Z eff a until the re-sulting pair-potential optimally fits the HNC inversion effective potential in Eq. (46) at large x. Examples for systems of colloidal particles with σ a = σ 1 = 250 nm, suspended at φ a = φ 1 = 10 −4 in an aqueous 1-1 electrolyte with a fixed concentration, n 3 = 1µM , of salt coions, are shown in Fig. 7. Results for three different colloidal charge numbers, Z 1 = 100, 400, and 1000 are shown. Note that χ 3 = χ coion ≫ Z a χ a = Z 1 χ 1 for all three systems shown in Fig. 7, such that the DLVO screening length in Eq. (49) is not considerably altered by replacing Z a with Z eff a . Fitting the potential in Eq. (48) to the one obtained from Eq. (46) therefore corresponds to vertical translation of βu DLVO aa (x) in the linear-logarithmic plot of Fig. 7. For the lowest considered colloid charge number, Z 1 = 100, the DLVOpotential calculated according to Eqs. (48) and (49) (lowermost blue dotted curve) is in nearly perfect agreement with the HNC-inversion potential (lowermost red solid curve). Fitting the DLVO potential with effective charge number input to the HNC-inversion potential, results in Z eff 1 = 99.85, which is only slightly lower than the bare colloidal charge number. The same procedure, carried out for Z 1 = 400 and Z 1 = 1000, results in effective charge numbers of Z eff 1 = 390 and 825, respectively. This demonstrates the capability of the employed numerical methods to access HNC solutions of the PM for realistic suspensions of charged colloids, where charge renormalization plays an important role. In Fig. 8, we plot the effective charge of macroions (species 1) as a function of the macroion bare charge, for various concentrations of salt coions, n 3 = 1, 10, and 100µM . In agreement with a recent small angle x-ray scattering study for charged silica spheres in aqueous electrolyte [106], we find that association of counterions to the macroion surfaces is most efficient at low salinity. This is reflected in Fig. 8 by a steepening decay of Z eff 1 (Z 1 ), for increasing salinity n 3 .

IV. CONCLUSIONS
We have shown that a combination of a numerically robust fixed-point iteration scheme and logarithmically spaced computational grids allows efficient computation of HNC solutions of the d-dimensional primitive model, with explicit results shown for d ≤ 6. Logarithmic grids are ideally suited for the discretization of pair-correlation functions of ionic mixtures with large particle diameter and charge asymmetries. This has allowed us to access HNC solutions for primitive model parameters corresponding to realistic suspensions of micrometer-sized colloidal spheres with charge numbers as high as |Z| ≈ 1000, in an aqueous 1-1 electrolyte.
Numerical stability might be further improved in future studies, if another elaborate fixed-point iteration scheme is used [65][66][67][77][78][79]. We expect that this would give access to HNC solutions of the primitive model at charge and diameter asymmetries exceeding the ones reported here.
Future projects, based on the methods presented here, might include ab initio modeling of colloidal suspensions with reactive electrolytes such as N aOH [46,47], including chemical association-dissociation reactions influenced by the locally varying pH-Value near the colloidal particle's surfaces.
As opposed to molecular dynamics simulations of the asymmetric primitive model, with very long program execution times even for moderate ion size-and chargeasymmetries, the solution of the HNC equations with the method described here takes few minutes or less on an inexpensive, standard computer. In addition, the HNC is a good approximation for mixtures of charged particles with long-ranged pair-potentials, predicting paircorrelation functions in good agreement with the numerically expensive computer simulations. Despite continuing rapid progress in computer simulations, liquid integral equations therefore remain an indispensable approach in studying highly asymmetric electrolytes.

A. Acknowledgements
We gratefully acknowledge Martin Oettel and Thomas Palberg for helpful discussions. This work was funded by