# Structural correlations and phase separation in binary mixtures of charged
and uncharged colloids

Elshad Allahyarov1,2,3, Hartmut Löwen2 1 Theoretical Department, Joint
Institute for High Temperatures, Russian Academy of Sciences (IVTAN), 13/19
Izhorskaya street, Moscow 125412, Russia
2 Institut für Theoretische Physik II: Weiche Materie, Heinrich-Heine
Universität Düsseldorf, Universitätstrasse 1, 40225 Düsseldorf, Germany
3 Department of Physics, Case Western Reserve University, Cleveland, Ohio
44106-7202, United States<EMAIL_ADDRESS>
###### Abstract

Structural correlations between colloids in a binary mixture of charged and
uncharged spheres are calculated using computer simulations of the primitive
model with explicit microions. For aqueous suspensions in a solvent of large
dielectric constant, the traditional Derjaguin-Landau-Vervey-Overbeek (DLVO)
theory of linear screening, supplemented with hard core interactions,
reproduces the structural correlations obtained in the full primitive model
quantitatively. However for lower dielectric contrast, the increasing Coulomb
coupling between the micro- and macroions results in strong deviations. We
find a fluid-fluid phase separation into two regions either rich in charged or
rich in uncharged particles which is not reproduced by DLVO theory. Our
results are verifiable in scattering or real-space experiments on charged-
uncharged mixtures of colloids or nanoparticles.

Keywords: Primitive Model simulations, DLVO theory, Binary colloidal mixture,
Entropic forces, Pair correlations, Structure factor, Charged macroion

##  1 Introduction

There are, in principle, two mechanisms to stabilize colloidal suspensions
against irreversible flocculation, namely charge-stabilization and steric
stabilization [1, 2, 3]. In the former case, the colloidal particles are
highly charged releasing counterions into the solvent such that they repel
each other by electrostatics which is traditionally described by a screened
Coulomb interaction. In the latter situation of steric stabilization, the
colloidal particles are typically coated with polymer brushes causing
repulsive entropic interaction forces which stabilize against flocculation.
While the effective interactions and structural correlations in strongly
interacting colloidal fluids are by now well-understood in one-component or
even polydisperse systems of either charged [4, 5, 6, 7, 8, 9, 10] or
uncharged colloidal spheres [11, 12, 13], much less is known about binary
mixtures of charged and uncharged particles. Such binary neutral-charged
systems occur frequently in mixtures of colloids with nanoparticles
culminating in the charged nanoparticle-halo effect around neutral colloids
which provides colloidal stabilization [15, 16, 17, 22, 18, 19, 20, 21, 14,
23, 24, 25, 26, 27, 28]. Moreover, in ordinary mixtures of charged colloids
the particle charge can be tuned by the pH of the solution [29, 30, 31] such
that one component stays charged and the other can become neutral close to the
isoelectric point realizing a charged-uncharged colloidal system. Further
examples are uncharged spherical vesicles exposed to charged colloids [32] and
mixtures of charged nanoparticles and neutral spherical bacteria [33].

About 30 years ago, first theoretical calculations were performed for a binary
mixture of charged and uncharged spheres [34]. The results were based on a
Yukawa model for the interaction between charged particles and a simple
excluded hard sphere interaction between the charged-neutral and neutral-
neutral spheres as predicted by standard DLVO-theory of linear screening
applied to such a mixture. Liquid integral equation closures were then used to
compute the partial pair correlation functions [34, 35, 36, 23, 31]. These
studies with effective pairwise interactions give some first insight into the
structural correlations but neglect nonlinear effects [37, 38, 39] beyond
linear screening. In general the latter cause effective many-body interactions
between the colloids [40, 41]. While earlier studies have employed approximate
Poisson-Boltzmann theory (see, e.g. Ref.[38]) and local classical density
functional theory for the inhomogeneous counterions in the field created by
the charged colloids [37, 39, 42, 43, 44], it has become possible by now to
calculate effective interactions and pair correlations with explicit
counterions based on the primitive model (PM) approach of electrolytes [45,
46, 47, 48, 49, 50, 51, 52, 53] which includes full nonlinear screening and
Coulomb correlations. Though there are simulations for charged mixtures [54,
55, 56, 57] and even for oppositely charged colloids [55, 58], to the best of
our knowledge a mixture of charged and neutral colloids has not yet been
simulated using the primitive model approach with explicit microions on a
large scale. The binary colloidal mixture considered in this work can be
assumed as a particular case of the charge regulated macroion system
investigated in Ref.[59] using mean-field formulation.

In this paper we close this gap and present computer simulations data for
charged and uncharged colloidal mixtures and extract their partial pair
correlation functions. Complementary to earlier studies of the colloidal halo
effect [15, 16, 17, 22, 18, 19, 20], we focus on the case of comparable hard
core radii of the two colloidal spheres. Here, we find that the traditional
Yukawa-hard sphere model is sufficient to describe the pair correlations in
aqueous suspensions where the dielectric constant (or relative permittivity)
$\epsilon$ of the solvent is pretty high (about 80 for water at room
temperature). However, in less polar solvents when $\epsilon$ is reduced by an
order of magnitude, the Coulomb coupling without screening between the charge
species is getting much stronger resulting in nonlinear screening effects. In
this study we show that for $\epsilon=8$, the standard Yukawa-hardcore
interaction model as proposed in [34] cannot be applied any longer to a
charged-uncharged mixture. There are significant deviations in the pair
correlations. Moreover we predict the existence of fluid-fluid phase
separation as documented by divergence in the partial static structure factors
at small wave vector. There are two regions either rich with charged or with
uncharged colloidal particles. This phase separation is absent within the
traditional DLVO-model when combined with excluded volume interactions.

The paper is organized as follows. In section 2 we describe the details of our
primitive model simulations for the binary colloidal system. The description
of the parameters used in the different runs is presented in section 3. The
results obtained for the partial pair correlation functions, as well as the
structure factors, are discussed in section 4 with supplementary snapshots
from the simulation boxes documenting phase separation for high Coulomb
couplings without screening. In section 5 we explore the role of the core size
of neutral colloids on the correlations and we conclude in section 6.

##  2  Details of the Primitive Model

We consider a three-component binary colloidal suspension consisting of
$N_{Z}$ macroions of charge $q^{(Z)}=Ze$ and size $\sigma_{Z}$ at positions
$\vec{r}_{i}^{\,(Z)}$ ($i$=1,…,$N_{Z}$), $N_{z}$ neutral colloids of zero
charge $q^{(z)}=0$ and size $\sigma_{z}$ at positions $\vec{r}_{j}^{\,(z)}$
($i$=1,…,$N_{z}$), and $N_{c}=ZN_{Z}$ monovalent counterions of charge
$q^{(c)}=-e$ and size $\sigma_{c}=\sigma_{Z}/600$ at positions
$\vec{r}_{\ell}^{\,(c)}$ ($\ell$=1,…,$N_{c}$). Here $e$ is the absolute value
of the electron charge. We fix the size ratio $\sigma_{Z}/\sigma_{z}$=1 to
unity such that $\sigma_{Z}=\sigma_{z}=\sigma$, and the number ratio
$N_{Z}/N_{z}$=1 to reduce parameter space. As a reference state we also
consider a case $q^{(z)}=q^{(Z)}$ which corresponds to a two component system
containing only one species of charged macroions and compensating counterions.

The pair interaction potentials between the species $\alpha$ and $\beta$ with
$\alpha,\beta\in\\{Z,z,c\\}$ are given as a combination of excluded volume and
Coulomb interactions (in SI units),

$V^{(\alpha\beta)}(r_{ij})=\left\\{\begin{array}[]{ l l
}\infty\,\,,&\quad\textrm{for}\,\,\,\,r_{ij}\leq\sigma_{\alpha\beta}\,\,,\\\
q^{(\alpha)}q^{(\beta)}/\Big{(}4\pi\varepsilon_{0}\varepsilon
r_{ij}\Big{)}\,\,,&\quad\textrm{for}\,\,\,\,r_{ij}>\sigma_{\alpha\beta}\,\,,\end{array}\right.$
(1)

where $\vec{r}_{ij}=\vec{r}_{j}^{\,(\alpha)}-\vec{r}_{i}^{\,(\beta)}$ with
$i\in 1,...,N_{\alpha}$ ($\alpha=Z,z,c$) and $j\in 1,...,N_{\beta}$
($\beta=Z,z,c$) is the distance between the two particles,
$\sigma_{\alpha\beta}=(\sigma_{\alpha}+\sigma_{\beta})/2$ is their additive
hard core diameter, $\varepsilon_{0}$ is the vacuum permittivity, and
$\varepsilon$ is the relative permittivity of the suspension. For simplicity,
we assume that $\varepsilon$ is the same throughout the system in order to
avoid image charge and dielectric boundary effects.

The following parameters characterize the intensity of interparticle
interactions and counterion screening effects in binary colloidal systems:
\- the total packing fraction $\eta$ of the macroions,
$\eta=\eta_{Z}+\eta_{z}$, where $\eta_{Z}=\pi N_{Z}\sigma^{3}/(6L^{3})$,
$\eta_{z}=\pi N_{z}\sigma^{3}/(6L^{3})$, and $L$ is the edge size of the cubic
simulation box.
\- the Debye-Hückel inverse screening length
$\kappa=\sqrt{n_{c}e^{2}/(\varepsilon_{0}\varepsilon k_{B}T)}$ of counterions,
where $k_{B}$ is the Boltzmann constant, $T$ is the temperature in the system,
and $n_{c}=N_{c}/L^{3}$ is the microion number density. High/low $\kappa$
values correspond to strong/weak counterion screening conditions.
\- the Bjerrum length, $\lambda_{B}=e^{2}/(4\pi\varepsilon_{0}\varepsilon
k_{B}T)$,
\- the average macroion-macroion separation distance in the charged and
neutral colloidal subsystems, $\bar{a}_{ZZ}=L\Big{(}6/(\pi
N_{Z})\Big{)}^{1/3}$ and $\bar{a}_{zz}=L\Big{(}6/(\pi N_{z})\Big{)}^{1/3}$. In
the current study, $\bar{a}_{zz}=\bar{a}_{ZZ}$ because we assume throughout
this paper $N_{z}=N_{Z}$.
\- the Coulomb coupling parameter in the charged subsystem with a screening
length $\kappa$, $\Gamma_{ZZ}$=$Z^{2}\exp(-\kappa\bar{a}_{ZZ})\lambda_{B}/a$,
and without screening, $\Gamma_{ZZ}^{*}$=$Z^{2}\lambda_{B}/a$.

##  3  Primitive Model Simulation Parameters

We have simulated globally electroneutral binary colloidal mixture in a cubic
box of an edge length $L$ with periodic boundary conditions in all three
Cartesian directions. The MD simulation method used here is the same as in
Refs. [60, 61, 56, 57]. In order to handle the long-ranged Coulomb
interactions, we use the Lekner summation method [62, 63, 64] which takes the
real-space particle coordinates as its only input.

We produced 4 different primitive model simulation run series $A_{i}$,
$B_{i}$, $C_{i}$, and $D_{i}$, $i$=1,..,3. Each run series consists of three
runs: a reference state run with $i$=1 where all macroions are charged, a
binary run with $i$=2 where one half of macroions are charged and the other
half is neutral, and another reference run with $i$=3 which is similar to the
run with $i$=2 but without neutral colloids. To distinguish the two reference
states, the reference state with $i$=1 will be referred as a “ground state”.

The run series differ in the values of $\eta$ and $\kappa$. In total, the run
series $A_{i}$ had low $\eta$ and low $\kappa$, the run series $B_{i}$ had
high $\eta$ and low $\kappa$, the run series $C_{i}$ had low $\eta$ and high
$\kappa$, and the run series $D_{i}$ had high $\eta$ and high $\kappa$.
Simulation parameters for all runs are collected in Table 1. To achieve high
values for $\kappa$ in the run series $C_{i}$ and $D_{i}$, i.e. low Debye
screening lengths $r_{D}=1/\kappa$, the relative permittivity $\varepsilon$ of
the solvent was decreased from 80 to 8. All simulations were carried out at
room temperature $T$=293 K, and the macroion diameter was fixed to
$\sigma$=100 nm. For the run series $A_{i}$ and $B_{i}$, the Bjerrum length
was $\lambda_{B}$=0.0071$\sigma$, whereas for the run series $C_{i}$ and
$D_{i}$, $\lambda_{B}$=0.071$\sigma$.

Table 1: Primitive model simulation parameters for the run series $A_{i}$, $B_{i}$, $C_{i}$, and $D_{i}$, $i$=1,2,3. The quantities listed in the first row are explained in the text. The binary composition is 1:1, thus each component has a packing fraction $\eta_{Z}=\eta_{z}=\eta$/2.  Runs | $Z$ | $z$ | $N_{Z}$ | $N_{z}$ | $\eta$ | $N_{c}$ | $\bar{a}_{ZZ}$ | $\kappa\sigma$ | $\Gamma_{ZZ}$ | $\Gamma_{ZZ}^{*}$ | $\varepsilon$
---|---|---|---|---|---|---|---|---|---|---|---
$A_{1}$ | 100 | 100 | 500 | 500 | 0.1 | 100000 | 2.15 | 1.31 | 1.98 | 33 | 80
$A_{2}$ | 100 | 0 | 500 | 500 | 0.1 | 50000 | 2.71 | 0.92 | 2.17 | 26 | 80
$A_{3}$ | 100 | - | 500 | 0 | 0.05 | 50000 | 2.71 | 0.92 | 2.17 | 26 | 80
$B_{1}$ | 100 | 100 | 500 | 500 | 0.2 | 100000 | 1.70 | 1.85 | 1.81 | 42 | 80
$B_{2}$ | 100 | 0 | 500 | 500 | 0.2 | 50000 | 2.15 | 1.32 | 1.94 | 33 | 80
$B_{3}$ | 100 | - | 500 | 0 | 0.1 | 50000 | 2.15 | 1.32 | 1.94 | 33 | 80
$C_{1}$ | 100 | 100 | 500 | 500 | 0.1 | 100000 | 2.15 | 4.14 | 0.05 | 330 | 8
$C_{2}$ | 100 | 0 | 500 | 500 | 0.1 | 50000 | 2.71 | 2.93 | 0.09 | 262 | 8
$C_{3}$ | 100 | - | 500 | 0 | 0.05 | 50000 | 2.71 | 2.93 | 0.09 | 262 | 8
$D_{1}$ | 100 | 100 | 500 | 500 | 0.2 | 100000 | 1.70 | 5.87 | 0.02 | 418 | 8
$D_{2}$ | 100 | 0 | 500 | 500 | 0.2 | 50000 | 2.15 | 4.15 | 0.04 | 330 | 8
$D_{3}$ | 100 | - | 500 | 0 | 0.1 | 50000 | 2.15 | 4.15 | 0.04 | 330 | 8
Table 2: Yukawa-DLVO model simulation parameters corresponding to the PM run series $A_{i}$, $B_{i}$, $C_{i}$, and $D_{i}$, $i$=1,2,3. The runs are labeled with an additional “Y”. Runs | $Z$ | $z$ | $N_{Z}$ | $N_{z}$ | $\eta$ | $\kappa\sigma$ | $\varepsilon$
---|---|---|---|---|---|---|---
$A_{1}Y$ | 100 | 100 | 500 | 500 | 0.1 | 1.31 | 80
$A_{2}Y$ | 100 | 0 | 500 | 500 | 0.1 | 0.92 | 80
$A_{3}Y$ | 100 | - | 500 | 0 | 0.05 | 0.92 | 80
$B_{1}Y$ | 100 | 100 | 500 | 500 | 0.2 | 1.85 | 80
$B_{2}Y$ | 100 | 0 | 500 | 500 | 0.2 | 1.32 | 80
$B_{3}Y$ | 100 | - | 500 | 0 | 0.1 | 1.32 | 80
$C_{1}Y$ | 100 | 100 | 500 | 500 | 0.1 | 4.14 | 8
$C_{2}Y$ | 100 | 0 | 500 | 500 | 0.1 | 2.93 | 8
$C_{3}Y$ | 100 | - | 500 | 0 | 0.05 | 2.93 | 8
$D_{1}Y$ | 100 | 100 | 500 | 500 | 0.2 | 5.87 | 8
$D_{2}Y$ | 100 | 0 | 500 | 500 | 0.2 | 4.15 | 8
$D_{3}Y$ | 100 | - | 500 | 0 | 0.1 | 4.15 | 8

For each PM run listed in Table 1, we also carried out corresponding binary
Yukawa-DLVO system simulations [34] with no explicit counterions. The pair
interaction potentials are given within DLVO-like theory as,

$V^{(\alpha\beta)}(r_{ij})=\left\\{\begin{array}[]{ l l
}\infty\,\,,&\quad\textrm{for}\,\,\,\,r_{ij}\leq\sigma\,\,,\\\
\frac{q^{(\alpha)}\exp{(\kappa\sigma/2)}}{1+\kappa\sigma/2}\,\frac{q^{(\beta)}\exp{(\kappa\sigma/2)}}{1+\kappa\sigma/2}\,\frac{\exp{(-\kappa
r_{ij})}}{4\pi\varepsilon_{0}\varepsilon
r_{ij}},&\quad\textrm{for}\,\,\,\,r_{ij}>\sigma\,\,,\end{array}\right.$ (2)

where the values for the bare colloidal charges $q^{(\alpha)}$ and
$q^{(\beta)}$, $\alpha,\beta=Z,z$, and for the inverse screening length
$\kappa$ are given in Table 2.

##  4 PM Simulation Results

###  4.1 Run series $A_{i}$ and $B_{i}$ for weak counterion screening

$A_{1}$ $A_{2}$ $B_{1}$ $B_{2}$

$A_{1}Y$ $A_{2}Y$ $B_{1}Y$ $B_{2}Y$

Figure 1:  Simulation snapshots from the PM runs $A_{1}$, $A_{2}$, $B_{1}$,
$B_{2}$ (upper row), and their counterpart Yukawa-DLVO runs $A_{1}Y$,
$A_{2}Y$, $B_{1}Y$, $B_{2}Y$ (bottom row). A color gradient from blue to
orange is used for the charged macroions according to their altitude in the
cell. Neutral colloids are shown as pink spheres. The PM counterions are shown
as scattered red dots.  Figure 2:  Simulation snapshot from the PM run
$B_{2}$.

The PM and Yukawa-DLVO simulation snapshots for the run series $A_{i}$ and
$B_{i}$ with low $\kappa\sigma\leq 1.85$, i.e. weak counterion screening, are
shown in Figure 1. The PM snapshot for the run $B_{1}$ is separately shown in
Figure 2 to make the counterions visible, which otherwise are hardly visible
in Figure 1. A visual inspection of the PM and the corresponding Yukawa-DLVO
runs reveals that in both cases the charged and neutral macroions are randomly
distributed across the system boundaries. This is an indirect indication of
the mutual repulsion between the charged macroions.

A detailed picture of the macroion distribution around a target macroion can
be achieved by calculating the partial pair correlation functions $g_{ij}(r)$,
$i,j$=Z,z, for the charged macroions and neutral colloids. These functions for
the run series $A_{i}$, $B_{i}$, $A_{i}Y$, and $B_{i}Y$ are shown in Figure 3.
Interestingly, both PM and the corresponding Yukawa-DLVO runs show similar
trends for the charged macroion $g_{ZZ}(r)$, such as the position of the first
maximum $r_{max}$ in $g_{ZZ}(r)$ for the ground state system (see the green
line) shifts to the right for the binary system run (see the black line). The
reference runs (see the pink line) have their $r_{max}$ in $g_{ZZ}(r)$ nearly
at the same distance $r$ for the binary runs. Such shifting of $r_{max}$ to
the right from the ground state run to the binary and reference runs clearly
signals about a stronger repulsion between the charged colloids in the latter
compared to the former run. The stronger repulsion originates from the
decreasing of the inverse screening length $\kappa$, i.e. from the weaker
screening of the macroions in the binary and reference runs.

$A_{1},A_{2},A_{3}$ $A_{1}Y,A_{2}Y$,$A_{3}Y$ $B_{1},B_{2},B_{3}$
$B_{1}Y,B_{2}Y$,$B_{3}Y$

Figure 3:  (Color in online)  Partial pair correlations $g_{ij}(r)$ for the PM
run series $A_{i}$ (a) and $B_{i}$ (c), and Yukawa-DLVO run series $A_{i}Y$
(b) and $B_{i}Y$ (d). Green lines- $g_{ZZ}(r)$ for the ground state runs
$A_{1}$, $A_{1}Y$, $B_{1}$, $B_{1}Y$. Black lines- $g_{ZZ}(r)$, red lines-
$g_{zz}(r)$, and blue lines- $g_{Zz}(r)$ for the binary charged-neutral runs
$A_{2}$, $A_{2}Y$, $B_{2}$, $B_{2}Y$. Pink lines- $g_{ZZ}(r)$ for the
reference runs $A_{3}$, $A_{3}Y$, $B_{3}$, $B_{3}Y$.

The PM and Yukawa-DLVO simulated $g_{Zz}(r)$ and $g_{zz}(r)$ for the charged-
neutral and neutral-neutral colloids, see the blue line and red lines Figure
3, correspondingly, look practically the same, except having major differences
at smaller separations $r$. At $r<1.1\sigma$, the PM simulated $g_{zz}(r)$ and
$g_{Zz}(r)$ show an upturn and reach higher contact values at $r=\sigma$,
whereas the related contact values in the Yukawa-DLVO results are around one.
The higher contact value in the PM simulated $g_{zz}(r)$ is an indication of a
pair-like clustering between neutral colloids. We assume that it is the
entropic force $F_{ent}(r)$ of the counterions acting on the neighboring
neutrals, which pushes them together. The definition of the counterion
entropic force acting on the macroions is briefly introduced in Appendix A.
Note that this force is absent in the Yukawa-DLVO model.

There is an upturn in $g_{Zz}(r)$ at small separations in the PM binary
simulations, which does not appear in the Yukawa-DLVO binary simulations. To
understand the origin of such short-range association between charged and
neutral colloids, we analyze the radial distribution of counterions,
$\rho_{c}^{(i)}(r)$, around the macroions in Appendix B. As can be seen from
Figure 9, both the charged and neutral colloids are wrapped by the screening
counterions. This gives rise to effective attraction between charged and
neutral particles via counterion depletion in the binary PM runs $A_{2}$ and
$B_{2}$ in Figure 3.

In total, because the upturn in $g_{zz}(r)$ and $g_{Zz}(r)$ appears at very
small separations, such pairing (or association) is not easy to detect
visually in the simulation snapshots in Figure 1.

$A_{1},A_{2}$ $A_{1}Y,A_{2}Y$ $B_{1},B_{2}$ $B_{1}Y,B_{2}Y$

Figure 4:  (Color in online)  Structure factors $S_{ij}(k)$ for the PM runs
$A_{1}$ and $A_{2}$ (a), $B_{1}$ and $B_{2}$ (c), and Yukawa-DLVO runs
$A_{1}Y$ and $A_{2}Y$ (b), $B_{1}Y$ and $B_{2}Y$ (d). Green lines- $g_{ZZ}(r)$
for the ground state runs $A_{1}$, $B_{1}$, $A_{1}Y$, $B_{1}Y$. Black lines-
$g_{ZZ}(r)$, and red lines- $g_{zz}(r)$ for the binary runs $A_{2}$, $A_{2}Y$,
$B_{2}$, $B_{2}Y$.

In Figure 4 we analyze the partial macroion-macroion structure factors
$S_{ij}(k)$. The colors used for the lines are the same as in Figure 3. There
are only from light to negligible differences between the PM and Yukawa-DLVO
data, which indicates that the structural correlations between the charged
macroions and between the neutral colloids, are nearly the same in both
simulation protocols. The calculated $S_{ij}(k)$ have non-diverging values at
small-$k$, which is pertinent to the random (i.e. no long-range clustering)
distribution of colloids in the system and indicates a finite compressibility.

###  4.2 Run series $C_{i}$ and $D_{i}$ for stromg counterion screening

$C_{1}$ $C_{2}$ $D_{1}$ $D_{2}$

$C_{1}Y$ $C_{2}Y$ $D_{1}Y$ $D_{2}Y$

Figure 5:  Simulation snapshots from the PM runs $C_{1}$, $C_{2}$, $D_{1}$,
$D_{2}$ (upper row), and the corresponding Yukawa-DLVO runs $C_{1}Y$,
$C_{2}Y$, $D_{1}Y$, $D_{2}Y$ (bottom row). Neutral macroions are shown as pink
spheres. The PM counterions are shown as scattered red dots.

The PM and Yukawa-DLVO simulation snapshots for the run series $C_{i}$ and
$D_{i}$ with high inverse screening $\kappa\sigma\geq 2.93$, i.e. strong
counterion screening, are shown in Figure 5. The PM and Yukawa-DLVO snapshots
look similar for the ground state runs $C_{1}$ and $C_{1}Y$, and for $D_{1}$
and $D_{1}Y$, however they strongly differ for the binary runs $C_{2}$ and
$C_{2}Y$, and for $D_{2}$ and $D_{2}Y$. The PM snapshots exhibit demixing and
clustering in the binary system, whereas such clustering is absent in the
binary Yukawa-DLVO runs.

Partial pair correlation functions $g_{ij}(r)$, presented in Figure 6, also
confirm the clustering in the binary PM runs. First, the height of the first
maximum in the PM ground state $g_{ZZ}(r)$ is smaller than its value in the
corresponding Yukawa-DLVO ground state, compare the green lines in Figure 6(a)
and (b), and in (c) and (d). We believe that this happens because of much
stronger macroion charge screening in the PM runs. As a result, strongly
screened macroions can easily approach each-other up to small separations in
the PM ground state runs. The same scenario is valid for the $g_{ZZ}(r)$ in
the binary PM runs, compare the black lines in Figure 6(a) and (b), and in (c)
and (d). The position of the first maximum $r_{max}$ in $g_{ZZ}(r)$ in the PM
binary run is shifted to the left compared to the corresponding Yukawa-DLVO
binary run. Moreover, the PM simulated black lines are above the line $g=1$
for the entire range of the separation distance $r$ for the run $D_{2}$, and
for $r>r_{max}$ for the run $C_{2}$, which is a direct indication of the
clustering of binary PM runs.

Second, as seen from the blue lines for $g_{Zz}(r)$ in Figure 6(a) and (c),
there is a repulsion between the charged and neutral colloids resulting in
$g_{Zz}(r)<1$ for $r\leq 4\sigma$. This can be viewed as another confirmation
of the demixing in the PM binary runs. In opposite to it, the Yukawa-DLVO
binary runs exhibit the macroion-neutral colloid association at low separation
distances where $g_{Zz}(r)>1$, see Figure 6(b) and (d). A similar association
was observed for the run series $A_{i}$ and $B_{i}$ in previous section.

Third, the binary PM results for the neutral-neutral pair correlation function
$g_{zz}(r)$, see the red curves in Figure 6(a) and (c), completely differ from
the corresponding Yukawa-DLVO results in Figure 6(b) and (d). Whereas the
latter shows no clustering at all, the former indicates a strong coagulation
effect in the system of neutral colloids with high contact values
$g_{zz}(\sigma)\approx 6$.

To analyze the origin of the observed clustering in binary PM runs $C_{2}$ and
$D_{2}$, additional two-macroion simulations were carried out for the runs
$D_{1}$ and $D_{2}$, see Appendix C for details. The calculated pair
interaction potentials show a repulsive interaction $U_{ZZ}(r)$ between the
charged macroions for the PM runs $D_{1}$ and $D_{2}$, see the the upper row
in Figure 10(a), and even much stronger repulsion $U^{Y}_{ZZ}(r)$ in the
Yukawa-DLVO runs $D_{1}Y$ and $D_{2}Y$, see the bottom row in Figure 10(a).

As seen from the upper row in Figure 10(b), the calculated cross interaction
potential $U_{Zz}(r)$ between the charged macroion and a neutral colloid is
always positive. It has a shallow minimum at small separations and a repulsive
barrier at larger distances. This repulsive barrier is the main reason why
neutral colloids avoid the neighborhood of the charged macroions in the PM
binary runs $C_{2}$ and $D_{2}$, as seen from the blue lines for $g_{Zz}(r)$
in Figure 6(a) and (c). Opposite to it, the Yukawa-DLVO cross term interaction
potential $U^{Y}_{Zz}(r)$ is always attractive, as seen from Figure 10(b),
bottom row. This explains the charged macroion-neutral colloid association for
the binary Yukawa-DLVO runs in Figure 6(b) and (d).

The neutral-neutral interaction potential $U_{Zz}(r)$ is attractive (by about
3 $k_{B}T$) in the PM run $D_{2}$, whereas $U^{Y}_{Zz}(r)$ is very weak in the
Yukawa-DLVO run $D_{2}Y$, see Figure 10(c) upper and bottom rows. This result
explains the observed clustering of the neutrals in the snapshots in Figure 5
for the binary PM runs $C_{2}$ and $D_{2}$.

$C_{1},C_{2},C_{3}$ $C_{1}Y,C_{2}Y,C_{3}Y$ $D_{1},D_{2},D_{3}$
$D_{1}Y,D_{2}Y,D_{3}Y$

Figure 6:  (Color in online)  Partial pair correlations $g_{ij}(r)$ for the PM
run series $C_{i}$ (a) and $D_{i}$ (c), and for the Yukawa-DLVO run series
$C_{i}Y$ (b ) and $D_{i}Y$ (d), $i$=1,2,3. Green lines- $g_{ZZ}(r)$ for the
ground state runs $C_{1}$, $C_{1}Y$, $D_{1}$, $D_{1}Y$. Black lines-
$g_{ZZ}(r)$, red lines- $g_{zz}(r)$, and blue lines- $g_{Zz}(r)$ for the
binary charged-neutral runs $C_{2}$, $C_{2}Y$, $D_{2}$, $D_{2}Y$. Pink lines-
$g_{ZZ}(r)$ for the reference runs $C_{3}$, $C_{3}Y$, $D_{3}$, $D_{3}Y$.

$C_{1},C_{2}$ $C_{1}Y,C_{2}Y$ $D_{1},D_{2}$ $D_{1}Y,D_{2}Y$

Figure 7:  (Color in online)  Structure factors $S_{ij}(k)$ for the PM run
series $C_{i}$ and $D_{i}$ ((a) and (c)), and for the Yukawa-DLVO run series
$C_{i}Y$ and $D_{i}Y$ ((b) and (d)). Green lines- $g_{ZZ}(r)$ for the ground
state runs $C_{1}$, $C_{1}Y$, $D_{1}$, $D_{1}Y$. Black lines- $g_{ZZ}(r)$, and
red lines- $g_{zz}(r)$ for the binary charged-neutral runs $C_{2}$, $C_{2}Y$,
$D_{2}$, $D_{2}Y$.

The PM and Yukawa-DLVO simulated structure factors $S_{ij}(k)$ are plotted in
Figure 7, The PM simulated $S_{zz}(k)$ for thew neutral component has strong
up-turn in the low-$k$ region for the binary runs $C_{2}$ and $D_{2}$. Such
divergence, though still finite in value, in practice is an indication of a
phase separation in the neutral component. Since all partial structure factor
are coupled, we expect that this indicates a global phase separation in a
charged macroion-rich and charged macroion-poor fluid. At the same time, no
such clustering (nor any sign of phase separation) is detected in the binary
Yukawa-DLVO runs $C_{2}Y$ and $D_{2}Y$ in Figure 7.

Radial distribution functions $\rho_{c}^{(Z)}(r)$ and $\rho_{c}^{(z)}(r)$ of
the counterions around the colloids for the runs series $C_{i}$ and $D_{i}$
are presented in Figure 9. Due to the strong Coulomb coupling of counterions
to the charged macroions, there are less counterions near neutral colloids
than in the bulk, which is another indication of the demixing in the charged-
neutral binary system.

##  5  Neutral colloid size effect on the clustering effect

Figure 8:  (Color in online)  PM simulated pair correlations functions
$g_{zz}(r)$ for neutral colloids for the binary runs $A_{2}$ (red line),
$E_{1}$ (pink line), $E_{2}$ (blue line), and $E_{3}$ (black line). The run
parameters are given in Table 3.

The entropic force of the screening counterions acting on neutral colloids,
see Eq.(3) in Appendix A, explicitly depends on the core radius $\sigma_{z}$
of the latter. Therefore, the bigger the neutral colloid, the stronger the
entropic force they endure. When a pair or several neutral colloids are
surrounded by the charged macroions, it is natural to expect that the entropic
force of the counterions will force the neutrals to stick to each other. For
smaller $\sigma_{z}$ the entropic force will be weak, and, as a result,
neutral colloids might evade the entropic pressure of the counterions and
attain random distribution across the system. In order to verify this
suggestion, the PM binary run $A_{2}$, for which a pair-like clustering in
$g_{zz}(r)$ was detected at low separations, see the red line in Figure 3(a),
is chosen as a reference system with $\sigma/\sigma_{z}$=1, and three
additional runs $E_{1}$-$E_{3}$ with the size ratio $\sigma/\sigma_{z}$=1.5,
2, and 4, see Table 3, were carried. While the size of neutral colloids is
decreased, their number $N_{z}$ was increased in order to keep their total
packing fraction fixed to $\eta/2$=0.05, provided that the charged macroions
keep their packing fraction fixed to $\eta/2$=0.05.

Table 3:  PM binary run parameters for $A_{2}$, $E_{1}$, $E_{2}$, and $E_{3}$.  Runs | $Z$ | $z$ | $N_{Z}$ | $N_{z}$ | $\eta$ | $N_{c}$ | $\bar{a}_{ZZ}$ | $\kappa\sigma$ | $\Gamma_{ZZ}$ | $\sigma/\sigma_{n}$ | $\varepsilon$
---|---|---|---|---|---|---|---|---|---|---|---
$A_{2}$ | 100 | 0 | 500 | 500 | 0.1 | 50000 | 2.71 | 0.92 | 2.17 | 1.0 | 80
$E_{1}$ | 100 | 0 | 500 | 1700 | 0.1 | 50000 | 2.71 | 0.92 | 2.17 | 1.5 | 80
$E_{2}$ | 100 | 0 | 500 | 4000 | 0.1 | 50000 | 2.71 | 0.92 | 2.17 | 2.0 | 80
$E_{3}$ | 100 | 0 | 500 | 32200 | 0.1 | 50000 | 2.71 | 0.92 | 2.17 | 4.0 | 80

Simulated pair correlation functions $g_{zz}(r)$, presented in Figure 8,
confirm that the smaller the neutral colloid, the weaker is the short-range
clustering in the neutral component. When the neutral colloid is 4 times
smaller than the charged macroion, the short-range clustering completely
disappears, see the black line in Figure 8.

##  6 Conclusions

To summarize, we have calculated the pair correlations in a binary mixture of
charged and neutral colloids using the primitive model with explicit
counterions and compared the results to the Yukawa-DLVO model. We found that
the traditional DLVO based approach is sufficient to describe the pair
correlations in aqueous suspensions with high dielectric constant. However, in
less polar solvents with reduced permittivity, the Coulomb coupling without
screening between the charge species is getting much stronger which results in
nonlinear screening effects. In this case the Yukawa-DLVO approach completely
fails to describe the PM results for the pair correlations. We predict, for
the first time, indications for a fluid-fluid phase separation in such
strongly nonlinear systems. Our simulations show that there are two regions
either rich with charged or with neutral colloidal particles.

Additional two macroion simulations were staged to calculate interaction
potentials between colloids in the binary runs that show fluid-fluid phase
separation. We found that the interaction between the charged colloids is
strongly repulsive, whereas it is moderately repulsive between the charged and
neutral colloids. More interestingly, the interaction between neutral colloids
is strongly attractive. We believe that the latter is the main reason for the
fluid-fluid phase separation in the simulated binary system with strong
Coulomb coupling.

Acknowledgments

The work of H.L. was supported by the DFG within project LO418/23-1. E.A.
thanks the financial support from the Ministry of Science and Higher Education
of the Russian Federation (State Assignment No. 075-01056-22-00).

##  Appendix A  Electrostatic and entropic forces in binary mixtures

In the PM simulations, the excluded volume of the macroions initiates entropic
forces arising from the contact density of counterions and neutral colloids on
the macroion surface. The entropic force acting on a $i$-th macroion of
species $\alpha$ at the position $\vec{r}_{i}^{\,(\alpha)}$ with $i\in
1,...,N_{\alpha}$ ($\alpha=Z,z$) is defined as [60, 61, 65, 66, 67, 68].

$\vec{F}_{ent}^{\,(\alpha)}(\vec{r}_{i}^{\,(\alpha)})=-k_{B}T\int_{S_{i}^{(\alpha)}}d\vec{f}\,\rho_{\beta}(\vec{r})$
(3)

where $\vec{f}$ is a surface normal vector pointing outwards the macroion’s
core, $S_{i}^{(\alpha)}$ is the surface of the hard core of the $i$-th
macroion centered around $\vec{r}_{i}^{\,(\alpha)}$ with diameter
$(\sigma_{\alpha}+\sigma_{\beta})/2$, $\beta$=$z$ for neutral colloids, and
$\beta$=$c$ for the counterions, and $\rho_{\beta}(\vec{r})$ is the density of
particles of sort $\beta$ around the $i$-th macroion. The entropic force,
usually neglected in weakly charged macroion systems, strongly modifies the
macroion interactions in highly charged and dense colloidal systems.

Likewise, the canonically averaged electrostatic force acting on the $i$-th
macroion of species $\alpha$ is defined as,

${\vec{F}_{elec}^{(\alpha)}(\vec{r}_{i}^{\,(\alpha)})}=\left<\sum_{\beta=Z,z,c}\,\,\sum_{j=1}^{N_{\beta}}\vec{F}^{(\alpha\beta)}(\vec{r}_{i}^{\,(\alpha)}-\vec{r}_{i}^{\,(\beta)})\,\Big{(}1-\delta_{\alpha\,\beta}\,\delta_{i\,j}\Big{)}\right>_{c}$
(4)

where $\alpha=Z,z$ and $\beta=z,c$. The Kronecker delta functions in this
expression nullify the self-interaction of macroions. Clearly, in Eq.(4) the
electrostatic pair interaction forces $\vec{F}^{(\alpha\beta)}$ are defined
as,

$\vec{F}^{(\alpha\beta)}(\vec{r}_{ij})=-\vec{\nabla}_{\vec{r}_{ij}}V^{(\alpha\beta)}(r_{ij})=\frac{1}{4\pi\varepsilon_{0}}\frac{q^{(\alpha)}q^{(\beta)}}{\varepsilon
r_{ij}^{2}}\frac{\vec{r}_{ij}}{r_{ij}}\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{ for
\,}r_{ij}>\sigma_{\alpha\beta}$ (5)

where $\vec{r}_{ij}=\vec{r}_{i}^{\,(\alpha)}-\vec{r}_{i}^{\,(\beta)}$.

##  Appendix B  Radial distribution of counterion around macroions

Radial distribution of counterions, $\rho_{c}^{(i)}(r)$, around the charged,
$i$=$Z$, and neutral, $i$=$z$, macroions is defined as,

$\rho_{c}^{(\alpha)}(r)=\Big{<}\sum_{i=1}^{N_{\alpha}}\rho_{c}\Big{(}\vec{r}-\vec{r}_{i}^{\,(\alpha)}\Big{)}\Big{>}_{\alpha}$
(6)

where the averaged counterion density field,

$\rho_{c}(\vec{r})=\Big{<}\sum_{\ell=1}^{N_{c}}\delta\Big{(}\vec{r}-\vec{r}_{\ell}^{\,(c)}\Big{)}\Big{>}_{c}$
(7)

parametrically depends on the fixed macroion positions
$\rho_{c}(\vec{r})\equiv\rho_{c}(\vec{r},\\{\vec{r}_{i}^{\,(Z)},\vec{r}_{j}^{\,(z)}\\})$,
$\\{\vec{r}_{i}^{\,(Z)},\vec{r}_{j}^{\,(z)},i=1,...,N_{Z};j=1,...N_{z}\\}$,
where the counterion averaging $\big{<}\cdot\cdot\cdot\big{>}_{c}$ and the
macroion averaging $\big{<}\cdot\cdot\cdot\big{>}_{\alpha}$ ($\alpha=Z,z$ )
are done for fixed macroion positions.

As seen from Figure 9, where the run series $A_{i}$ and $B_{i}$ are analyzed,
the counterion density at the surface of neutral colloids in the PM binary
runs is larger than in the bulk. This effect appears to be more solid for
$\eta$=0.2 simulations. Because no counterions are originally associated with
neutral colloids, it is the counterions stemming from the charged macroions
that exist near neutral colloids. This conclusion is supported by the
appearance of the minimum in the counterion distribution at around
$r=1.2\sigma$, see the red line, below which the macroion-neutral association
develops in the PM binary runs in Figure 3, see the blue line there.

$A_{1},A_{2}$ $B_{1},B_{2}$ $C_{1},C_{2}$ $D_{1},D_{2}$

Figure 9:  (Color in online)  Normalized radial distribution function
$\rho_{c}^{(q)}(r)$ of counterions around the macroions for the run series
$X_{i}$, $i=$1,2, $X$=$A$, $B$, $C$, $D$. Green lines- $q^{(Z)}$=100$e$ for
the ground state runs $X_{1}$, black lines- $q^{(Z)}$=100$e$, and red line-
$q^{(z)}$=0 for the binary runs $X_{2}$.

For the run series $C_{i}$ and $D_{i}$, there are less counterions near
neutral colloids than in the bulk, see the monotonic increase of
$\rho_{c}^{(z)}(r)$ from the distance $r=0.5\sigma$ to $r=3\sigma$. In other
words, the counterions are mostly localized to the regions of charged
macroions. This is another indication of the demixing in the charged-neutral
system. Another interesting factor is a much stronger screening of the
macroion charge in the binary PM run $D_{2}$ compared to the ground state run
$D_{1}$, see the gap between these two lines in Figure 9. The reason for such
stronger screening in $D_{2}$ is the stronger macroion clustering effect in
the $D_{2}$ run. Larger clusters of charged macroions practically trap all
their compensating counterions in and around them for effective screening.
This suggestion is also supported by the differences in the black line
($g_{ZZ}$ for the binary PM run $D_{2}$) and green line ($g_{ZZ}$ for the
ground state PM run $D_{1}$) in Figure 6.

##  Appendix C  Macroion pair interaction potentials of the averaged force

We calculated macroion-macroion interaction potentials explicitly in the PM
simulation runs $D_{1}$, and $D_{2}$, and in the Yukawa-DLVO runs $D_{1}Y$,
and $D_{2}Y$. For this purpose, two macroions of charge $q^{(Z)}$=100$e$ were
placed along the diagonal of the simulation box of size $L$ at a separation
distance $r$. First, the canonically averaged total macroion-macroion
interaction force $F_{ZZ}((r)$ was calculated for a set of $r$ varied between
1$\sigma$ and 4$\sigma$. This force includes the direct interaction between
the fixed macroions, the electrostatic interaction with other macroions and
counterions, and the entropic force arising from the collisions with other
macroions and counterions. As a remark, a full average has been performed here
over all remaining particles. The electrostatic and entropic forces are
defined in Appendix A. In the Yukawa-DLVO runs with no counterions, the total
macroion-macroion interaction force $F^{Y}_{ZZ}$ includes only the direct
interaction between the macroions, and the electrostatic and entropic
interaction with other macroions. Second, the obtained total force is
integrated over the separation distance $r$ to get the interaction potential
of the averaged force $U_{ZZ}(r)$ in the PM runs, and $U^{Y}_{ZZ}(r)$ in the
Yukawa-DLVO runs. Third, we repeat this procedure for another pair of fixed
charges $q^{(Z)}$=100$e$ and $q^{(z)}$=0, as well as for $q^{(z)}$=0 and
$q^{(z)}$=0, to calculate the interactions $U_{Zz}(r)$ and $U_{zz}(r)$ for the
PM runs, and $U^{Y}_{Zz}(r)$ and $U^{Y}_{zz}(r)$ for the Yukawa-DLVO runs. The
resulting macroion-macroion interaction potentials for the averaged force are
plotted in Figure 10.

Figure 10:  (Color in online)  Upper row- pair interaction potentials
$U_{ij}(r)/(k_{B}T)$ between two fixed macroions in the PM runs. (a)
Interaction potentials between the fixed macroions with $q^{(Z)}$=100$e$ in
the run $D_{1}$ (red line with circles), and in the run $D_{2}$ (blue line
with squares). (b) The interaction potential $U_{Zz}(r)$ between the fixed
macroions with $q^{(Z)}$=100$e$ and $q^{(z)}$=0 in the run $D_{2}$. (c) The
interaction potential $U_{Zz}(r)$ between the fixed macroions with $q^{(z)}$=0
and $q^{(z)}$=0 in the run $D_{2}$. Bottom row- the same as in the upper row,
but now for the Yukawa-DLVO runs.

## Literature

## References

  * [1] P. N. Pusey, Colloidal suspensions in Liquids, Freezing and Glass Transition, Les Houches Session 51, edited by J.-P. Hansen, D. Levesque, and J. Zinn-Justin, (North-Holland, Amsterdam) 2, 763-931 (1991).
  * [2] A. Ivlev, H. Löwen, G. Morfill, and C. P. Royall, Complex plasmas and colloidal dispersions: particle-resolved studies of classical liquids and solids, Series in Soft Condensed Matter 5, ISBN: 978-981-4350-06-8, (2012), https://doi.org/10.1142/8139.
  * [3] L. F. Rojas-Ochoa, R. Castaneda-Priego, V. Lobaskin, A. Stradner, F. Scheffold, and P. Schurtenberger, Density dependent interactions and structure of charged colloidal dispersions in the weak screening regime, Phys. Rev. Lett. 100, 178304 (1-4) (2008), https://doi.org/10.1103/PhysRevLett.100.178304.
  * [4] J. M. Mendez-Alcaraz, B. Daguanno, and R. Klein, The structure of binary-mixtures of charged colloidal suspensions, Physica A - Statistical Mechanics and its applications 178, 421-443 (1991), https://doi.org/10.1016/0378-4371(91)90031-7.
  * [5] B. D’Aguanno, R. Krause, J. M. Mendez-Alcaraz, and R. Klein, Structure factors of charged bidispersed colloidal suspensions, J. Phys. Condens. Matter 4, 3077-3086 (1992), https://doi.org/10.1088/0953-8984/4/12/007.
  * [6] P. Wette, H. J. Schoepe, and T. Palberg, Enhanced crystal stability in a binary mixture of charged colloidal spheres, Phys. Rev. E 80, 021407 (1-12) (2009), https://doi.org/10.1103/PhysRevE.80.021407.
  * [7] S. D. Finlayson and P. Bartlett, Non-additivity of pair interactions in charged colloids, J. Chem. Phys. 145, 034905 (1-10) (2016), https://doi.org/10.1063/1.4959122.
  * [8] L. Javidpour, B. A. Losdorfer, R. Podgornik, and A. Naji, Role of metallic core for the stability of virus-like particles in strongly coupled electrostatics, Scientific Reports 9, 3884 (1-13) (2019), https://doi.org/10.1038/S41598-019-39930-8.
  * [9] F. Smallenburg, N. Boon, M. Kater, M. Dijkstra, and R. van Roij, Phase diagrams of colloidal spheres with a constant zeta-potential, J. Chem. Phys. 134, 074505 (1-8) (2011), https://doi.org/10.1063/1.3555627.
  * [10] J. Dzubiella, J. Chakrabarti, and H. Löwen, Tuning colloidal interactions in subcritical solvents by solvophobicity: Explicit versus implicit modeling, J. Chem. Phys. 131, 044513 (1-5) (2009), https://doi.org/10.1063/1.3193557.
  * [11] M. Dijkstra, R. van Roij, and R. Evans, Phase diagram of highly asymmetric binary hard-sphere mixtures, Phys. Rev. E 59, 5744-5771 (1999), https://doi.org/10.1103/physreve.59.5744.
  * [12] R. Roth, Fundamental measure theory for hard-sphere mixtures: a review, J. Phys. Condens. Matter 22, 063102 (1-18) (2010), https://doi.org/10.1103/PhysRevE.91.052121.
  * [13] N. Heptner and J. Dzubiella, Equilibrium structure and fluctuations of suspensions of colloidal dumbbells, Molecular Physics 113, 2523-2530 (2015), http://dx.doi.org/10.1080/00268976.2015.102260.
  * [14] S. A. Barr and E. Luijten, Effective interactions in mixtures of silica microspheres and polystyrene nanoparticles, Langmuir 22, 7152-7155 (2006), https://doi.org/10.1021/la061291d.
  * [15] A. T. Chan and J. A. Lewis, Electrostatically tuned interactions in silica microsphere- polystyrene nanoparticle mixtures, Langmuir 21, 8576-8579 (2005), https://doi.org/10.1021/la0510073.
  * [16] V. Tohver, J. E. Smay, A. Braem, P. V. Braun, and J. A. Lewis, Nanoparticle halos: A new colloid stabilization mechanism, Proc. Natl. Acad. Sci. USA 98, 8950-8954 (2001), https://doi.org/10.1073/pnas.151063098.
  * [17] V. Tohver, A. Chan, O. Sakurada, and J. A. Lewis, Nanoparticle engineering of complex fluid behavior, Langmuir 17, 8414-8421 (2001), https://doi.org/10.1021/la011252w.
  * [18] J. Liu and E. Luijten, Rejection-free geometric cluster algorithm for complex fluids, Phys. Rev. Lett. 92, 035504 (1-4) (2004), https://doi.org/10.1103/PhysRevLett.92.035504.
  * [19] J. Liu and E. Luijten, Generalized geometric cluster algorithm for fluid simulation, Phys. Rev. E 71, 066701 (1-12) (2005), https://doi.org/10.1103/PhysRevE.71.066701.
  * [20] J. Liu and E. Luijten, Stabilization of colloidal suspensions by means of highly charged nanoparticles, Phys. Rev. Lett. 93, 247802 (1-4) (2004), https://doi.org/10.1103/PhysRevLett.93.247802.
  * [21] S. Karanikas and A. A. Louis, Dynamic colloidal stabilization by nanoparticle halos, Phys. Rev. Lett. 93, 248303 (1-4) (2004), https://doi.org/10.1103/PhysRevLett.93.248303.
  * [22] A. A. Louis, E. Allahyarov, H. Löwen, and R. Roth, Effective forces in colloidal mixtures: from depletion attraction to accumulation repulsion, Phys. Rev. E 65, 061407 (1-11) (2002), https://doi.org/10.1103/PhysRevE.65.061407.
  * [23] E. N. Scheer and K. S. Schweizer, Haloing, flocculation, and bridging in colloid-nanoparticle suspensions, J. Chem. Phys. 128, 164905 (1-15) (2008), https://doi.org/10.1063/1.2907721.
  * [24] D. Herman and J. Y. Walz, Stabilization of Weakly Charged Microparticles Using Highly Charged Nanoparticles, Langmuir 29, 5982-5994 (2013), https://doi.org/10.1021/la400699g.
  * [25] H. Huang and E. Ruckenstein, Repulsive force between two microparticles decorated with highly charged nanoparticles, Colloids and Surfaces A: Physicochemical and Engineering Aspects 436, 862-867 (2013), https://doi.org/10.1016/j.colsurfa.2013.08.024.
  * [26] N. M. Zubir, A. Badarudin, S. N. Kazi, M. Misran, A. Amiri, R. Sadri, and S. Khalid, Experimental investigation on the use of highly charged nanoparticles to improve the stability of weakly charged colloidal system, J. Colloid and Interface Science 454, 245-255 (2015), https://doi.org/10.1016/j.jcis.2015.05.019.
  * [27] B. M. Weight and A. R. Denton, Structure and stability of charged colloid-nanoparticle mixtures, J. Chem. Phys. 148, 114904 (1-11) (2018), https://doi.org/10.1063/1.5004443.
  * [28] A. R. Denton, Effective electrostatic interactions in colloid-nanoparticle mixtures, Phys. Rev. E 96, 062610 (1-14) (2017), https://doi.org/10.1103/PhysRevE.96.062610.
  * [29] H. J. M. Hanley, G. C. Straty, and P. Lindner, Order in a simple colloidal mixture suspension, Physica A 174, 60-73 (1991), https://doi.org/10.1016/0378-4371(91)90417-B.
  * [30] M. E. Leunissen, C. G. Christova, A. P. Hynninen, C. P. Royall, A. I. Campbell, A. Imhof, M. Dijkstra, R. van Roij, and A. van Blaaderen, Ionic colloidal crystals of oppositely charged particles, Nature 437, 235-240 (2005), https://doi.org/10.1038/nature03946.
  * [31] G. J. Ojeda-Mendoza, A. Moncho-Jorda, P. Gonzalez-Mozuelos, C. Haro-Perez, and L. F. Rojas-Ochoa, Evidence of electrostatic-enhanced depletion attraction in the structural properties and phase behavior of binary charged colloidal suspensions, Soft Matter 14, 1355-1364 (2018), https://doi.org/10.1039/c7sm02220d.
  * [32] S. Savarala, S. Ahmed, M. A. Ilies, and S. L. Wunder, Stabilization of Soft Lipid Colloids: Competing Effects of Nanoparticle Decoration and Supported Lipid Bilayer Formation, ACS Nano 5, 2619-2628 (2011), https://doi.org/10.1021/nn1025884.
  * [33] B. Q. Yang, F. Gao, Z. Y. Li, M. Y. Li, L. Chen, Y. Guan, G. Liu, and L. H. Yang, Selective Entropy Gain-Driven Adsorption of Nanospheres onto Spherical Bacteria Endows Photodynamic Treatment with Narrow-Spectrum Activity, J. Phys. Chem. Lett. 11, 2788-2796 (2020), https://doi.org/10.1021/acs.jpclett.0c00287.
  * [34] J. M. Mendez-Alcaraz, B. D’Aguano, and R. Klein, Structure of binary colloidal mixtures of charged and uncharged spherical particles, Langmuir 8, 2913-2920 (1992), https://doi.org/10.1021/la00048a012.
  * [35] C. Caccamo, M. Varisto, M. A. Floriano, E. Caponetti, R. Triolo, and G. Lucido, Phase stability of dense multicomponent charged and uncharged hard-sphere fluid mixtures, J. Chem. Phys. 98, 1579-1587 (1993), https://doi.org/10.1063/1.464274.
  * [36] R. Garibay-Alonso, J. M. Mendez-Alcaraz, and R. Klein, Phase separation of binary liquid mixtures of hard spheres and Yukawa particles, Physica A 235, 159-169 (1997), https://doi.org/10.1016/S0378-4371(96)00337-8.
  * [37] H. Löwen, P. A. Madden, and J. P. Hansen, An ab initio description of counterion screening in colloidal suspensions, Phys. Rev. Lett. 68, 1081-1084 (1992), https://doi.org/10.1103/PhysRevLett.68.1081.
  * [38] M. Fushiki, Molecular-dynamics simulations for charged colloidal dispersions, J. Chem. Phys. 97, 6700-6713 (1992), https://doi.org/10.1063/1.463676.
  * [39] H. Löwen, J. P. Hansen, and P. A. Madden, Nonlinear counterion screening in colloidal suspensions, J. Chem. Phys. 98, 3275-3289 (1993), https://doi.org/10.1063/1.464099.
  * [40] H. Löwen and E. Allahyarov, Role of effective triplet interactions in charged colloidal suspensions, J. Phys. Condens. Matter 10, 4147-4160 (1998), https://doi.org/10.1088/0953-8984/10/19/003.
  * [41] C. Russ, H. H. von Grünberg, M. Dijkstra, and R. van Roij, Three-body forces between charged colloidal particles, Phys. Rev. E 66, 011402 (1-12) (2002), https://doi.org/10.1103/PhysRevE.66.011402.
  * [42] H. Löwen and G. Kramposthuber, Optimal effective pair potential for charged colloids, Europhys. Lett. 23, 637-678 (1993), https://doi.org/10.1209/0295-5075/23/9/009.
  * [43] M. C. Barbosa, M. Deserno, and C. Holm, A stable local density functional approach to ion-ion correlations, Europhys. Lett. 52, 80-86 (2000), https://doi.org/10.1209/epl/i2000-00407-y.
  * [44] M. C. Barbosa, M. Deserno, C. Holm, and R. Messina, Screening of spherical colloids beyond mean field: A local density functional approach, Phys. Rev. E 69, 051401 (1-9) (2004), https://doi.org/10.1103/PhysRevE.69.051401.
  * [45] J.-P. Hansen and H. Löwen, Effective interactions between electric double-layers, Annual Reviews of Physical Chemistry 51, 209-242 (2000), https://doi.org/10.1146/annurev.physchem.51.1.209.
  * [46] L. Belloni, Colloidal interactions, J. Phys. Conf. Matter 12, R549-R587 (2000), https://doi.org/10.1088/0953-8984/12/46/201.
  * [47] P. Linse, Simulation of charged colloids in solution, Advanced Computer Simulation Approaches Soft Matter Sciences II, Book Series: Advances in Polymer Science 185, 111-162 (2005), https://doi.org/10.1007/b136795.
  * [48] A. Torres, A. Cuetos, M. Dijkstra, R. and van Roij, Sedimentation of charged colloids: The primitive model and the effective one-component approach, Phys. Rev. E 75, 041405 (1-8) (2007), https://doi.org/10.1103/PhysRevE.75.041405.
  * [49] M. Dijkstra, Computer simulations of charge and steric stabilized colloidal suspensions, Current Opinion in Colloid & Interface Science 6, 372-382 (2001), https://doi.org/10.1016/S1359-0294(01)00106-6.
  * [50] M. Quesada-Perez, E. Gonzales-Tovar, A. Martin-Molina, M. Lozada-Cassou, and R. Hidalgo-Alvarez, Overcharging in Colloids: Beyond the Poisson-Boltzmann Approach, Chem. Phys. Chem. 4, 234-248 (2003), https://doi.org/10.1002/cphc.200390040.
  * [51] R. Messina, Electrostatics in soft matter, J. Phys. Condens. Matter 21, 113102 (1-18) (2009), https://doi.org/10.1088/0953-8984/21/11/113102.
  * [52] E. Allahyarov and H. Löwen, Influence of solvent granularity on the effective interaction between charged colloidal suspensions, Phys. Rev. E 63, 041403 (1-13) (2001), https://doi.org/10.1103/PhysRevE.63.041403.
  * [53] R. Castaneda-Priego, L. F. Rojas-Ochoa, V. Lobaskin, and J. C. Mixteco-Sanchez, Macroion correlation effects in electrostatic screening and thermodynamics of highly charged colloids, Phys. Rev. E 74, 051408 (1-6) (2006), https://doi.org/10.1103/PhysRevE.74.051408.
  * [54] E. Allahyarov, H. Löwen, and S. Trigger, Effective forces between macroions: the cases of asymmetric macroions and added salt, Phys. Rev. E 57, 5818-5824 (1998), https://doi.org/10.1103/PhysRevE.57.5818.
  * [55] M. Dijkstra, J. W. Zwanikken, and R. van Roij, Sedimentation of binary mixtures of like- and oppositely charged colloids: the primitive model or effective pair potentials?, J. Phys. Condens. Matter 18, 825-836 (2006), https://doi.org/10.1088/0953-8984/18/3/005.
  * [56] E. Allahyarov and H. Löwen, Nonadditivity in the effective interactions of binary charged colloidal suspensions, J. Phys. Condens. Matter 21, 424117 (1-7) (2009), https://doi.org/10.1088/0953-8984/21/42/424117.
  * [57] E. Allahyarov, H. Löwen, and A. R. Denton, Structural correlations in highly asymmetric binary charged colloidal mixtures, Phys. Chem. Chem. Phys. 24, 15439-15451 (2022), https://doi.org/10.1039/D2CP01343F.
  * [58] A. F. Demirors, J. C. P. Stiefelhagen, T. Vissers, F. Smallenburg, M. Dijkstra, A. Imhof, and A. van Blaaderen, Long-ranged oppositely charged interactions for designing new types of colloidal clusters, Phys. Rev. X 5, 021012 (1-12) (2015), https://doi.org/10.1103/PhysRevX.5.021012.
  * [59] Y. Avni, R. Podgornik, and D. Andelman, Critical behavior of charge-regulated macro-ions, J. Chem. Phys. 153, 024901 (1-11) (2020), https://doi.org/10.1063/5.0011623.
  * [60] E. Allahyarov, I. D’Amico, and H. Löwen, Attraction between like-charged macroions by Coulomb depletion, Phys. Rev. Lett. 81, 1334-1337 (1998), https://doi.org/10.1103/PhysRevLett.81.1334.
  * [61] E. Allahyarov, G. Gompper, and H. Löwen, Attraction between DNA molecules mediated by multivalent ions, Phys. Rev. E 69, 041904 (1-13) (2004), https://doi.org/10.1103/PhysRevE.69.041904.
  * [62] J. Lekner, Summation of dipolar fields in simulated liquid-vapour interfaces, Physica A 157, 826-838 (1989), https://doi.org/10.1016/0378-4371(89)90068-X.
  * [63] J. Lekner, Summation of Coulomb fields in computer-simulated disordered systems, Physica A 176, 485-498 (1991), https://doi.org/10.1016/0378-4371(91)90226-3.
  * [64] M. Mazars, Lekner summations, J. Chem. Phys. 115, 2955-2965 (2001), https://doi.org/10.1063/1.1386904.
  * [65] J. Z. Wu, D. Bratko, H. W. Blanch, and J. M. Prausnitz, Monte Carlo simulation for the potential of mean force between ionic colloids in solutions of asymmetric salts, J. Chem. Phys. 111, 7084-7094 (1999), https://doi.org/10.1063/1.480000.
  * [66] J. Wu, D. Bratko, and J. M. Prausnitz, Interaction between like-charged colloidal spheres in electrolyte solutions, Proc. Natl. Acad. Sci. USA 95, 15169-15172 (1998), https://doi.org/10.1073/pnas.95.26.1516.
  * [67] Y.-R. Chen and C.-L. Lee, A cell-model study on counterion fluctuations in macroionic systems: effect of non-extensiveness in entropy, Phys. Chem. Chem. Phys. 16, 297-303 (2014), https://doi.org/10.1039/C3CP53318B.
  * [68] A. Gonzalez-Calderon, E. Gonzalez-Tovar, and M. Lozada-Cassou, Very long-range attractive and repulsive forces in model colloidal dispersions, Euro. Phys. 227, 2375-2390 (2019), https://doi.org/10.1140/epjst/e2018-800089-0.